# Properties in Stage-Structured Population Models with Deterministic and Stochastic Resource Growth

Authors: Tin Nwe Aye; Linus Carlsson
Journal: Journal of Applied Mathematics
 (2022)

---

## Abstract

Modelling population dynamics in ecological systems reveals properties that are difficult to find by empirical means, such as the probability that a population will go extinct when it is exposed to harvesting. To study these properties, we use an aquatic ecological system containing one fish species and an underlying resource as our models. In particular, we study a class of stage-structured population systems with and without starvation. In these models, we study the resilience, the recovery potential, and the probability of extinction and show how these properties are affected by different harvesting rates, both in a deterministic and stochastic setting. In the stochastic setting, we develop methods for deriving estimates of these properties. We estimate the expected outcome of emergent population properties in our models, as well as measures of dispersion. In particular, two different approaches for estimating the probability of extinction are developed. We also construct a method to determine the recovery potential of a species that is introduced in a virgin environment.


---

## Body

## 1. Introduction

The dynamics in general ecological systems are very complicated, and it is quite often difficult to draw appropriate precise conclusions. In spite of this difficulty, some properties are well-defined and commonly studied, such as biodiversity, stability, and food webs; see, e.g., Bardgett and van der Putten [1], Brännström et al. [2], Flores et al. [3], Loreau [4], and Loreau and de Mazancourt [5]. In this article, we have used a small aquatic ecological system as the base model, but the results herein can be applied to other types of ecological systems. More precisely, we assume that our ecosystem consists of one fish species and an underlying food resource.The field of stage-structured population models has attracted much attention; see, e.g., Lundström et al. [6], Meng et al. [7], Ackleh and Jang [8], Aiello et al. [9], Liz and Pilarczyk [10], Soudijn and de Roos [11], and de Roos et al. [12]. The stage-structured population model can be derived from the general physiologically structured population model; see, e.g., de Roos et al. [13]. In this paper, the models we use are stage-structured models with and without starvation, in which we assume that the fish population is divided into two stages: the juveniles and the adults, and the food resource is treated as an unstructured entity.Emergent properties of the stage-structured population model are normally studied in deterministic models. The main purpose of this project is to develop methods to understand the emergent properties in stochastic models. One commonly uses a system of deterministic partial differential equations to model physiologically structured populations and from these equations derives stage-structured population models. In this paper, we introduce randomness in the stage-structured models by adding randomness in the growth rate models for the resource. There is of course a wide variety of possibilities to include randomness in the models, but for the sake of clarity, we restrict ourselves to this extension. Randomness in the growth rate of the food resource is a natural adaptation to capture environmental variations such as the daily changes in temperature and sun exposure. Stochastic stage-structured models have been studied for some time; see, e.g., Burgman and Gerard [14], Castañera et al. [15], and Scranton et al. [16]. A major advantage in the stochastic setting is that, in addition to the expected values, we also get the dispersion of the different emergent properties of the population and resource. We have used Monte Carlo methods to evaluate the resilience, the recovery potential, and the probability of extinction. In Appendix C, we have also studied estimates for the juvenile/adult/resource biomass, yield, impact on biomass, and impact on size structure.In deterministic population models, the population either goes extinct or it stays positive for all future periods of time. However, in stochastic population models, depending on the choice of parameters, the population goes extinct within a finite period of time, with probability0<p<1. We call this the probability of extinction, which is closely connected to the minimum viable population (MVP) size. This property (MVP) is also investigated in, e.g., Wang et al. [17], Flather et al. [18], and Shaffer [19]. The MVP formulation of the probability of extinction is based on Monte Carlo simulations. We have also constructed a more natural approach to find the probability of extinction by utilizing the fact that when the recovery potential is larger than one, there is no extinction. This formulation, called RP formulation, uses statistical methods on the recovery potential (Section 3.5). The RP formulation corroborates well with the MVP formulation (see Section 4.3). In addition, we compare the emergent properties of the stage-structured population model, when starvation mortality is included or excluded.As proved by Abrams [20], Abrams and Matsuda [21], in some structured population models, an increase of biomass of a certain species can be obtained by increasing the mortality rate, which is called the hydra effect. We find that in some of our models, a hydra effect is present (see Conclusions and Discussion). The hydra effect was first noted by Ricker [22] and has later been incorporated in different models; see, e.g., Abrams [20], Adhikary et al. [23], and Ghosh et al. [24].This article is structured as follows: in Section2, we present the stage-structured population models investigated in this paper, both with the deterministic and stochastic resource dynamics. In Section 3, we define and discuss the emergent population properties: yield, impact on biomass and size structure, resilience, recovery potential, and the probability of extinction. In Section 4, the simulation results of these emergent properties are presented; the yield and impact on biomass and size structure are presented in Appendix C. In Section 5, we compare the simulation results with and without starvation mortality rates for all emergent properties. In Section 6, we end the paper with a discussion and conclusions of our findings.


## 2. The Stage Structured Population Models

In this paper, we consider an ecological model that consists of a single species and a food resource. The model is a two-staged structured fish population model with an unstructured resource. First, we describe the model in a deterministic setting; then, we change the growth rate of the resource from deterministic to stochastic. In our models, the proportional harvesting rates for the juvenile and adult stages are assumed to be equal. We investigate the impact on the solutions to our model when starvation mortality on the juvenile and the adult stages is included or not. We also describe and investigate emergent properties such as yield, impact on biomass, impact on size structure, resilience, and the recovery potential. When considering the stochastic stage-structured models, the probability of extinction is investigated by using two different methods.


### 2.1. Population Stage Dynamics

The individuals of the species are modeled by a stage-structured biomass model which is derived by formulating a size-structured population model on the individuals’ life history. More precisely, Individuals are divided into two stages: juveniles and adults, based only on their size. For many relevant properties of the species, a two-stage model is often enough; see, e.g., Ackleh and Jang [8], de Roos et al. [13], Lundström et al. [6], and Meng et al. [7]. Individuals are assumed to have the same size sbirth at birth, and the maximum size of individuals is denoted by smax. Both stages forage for the shared resource, R=Rt, and the metabolic requirement, T, is assumed to be constant. The growth rate depends on the available resource ([25]) and varies with population density ([26, 27]). Juvenile and adult biomasses are denoted by J=Jt and A=At, respectively.The derivation of the stage-structured biomass model is investigated by formulating a size-structured population model. We give a brief description of the model and refer the reader to de Roos et al. [13] for further details and motivations. The population is composed of juveniles, J, and adults, A, which forages on a shared resource R. Juvenile individuals use all the consumed energy for growth, development, and maintenance, whereas adult individuals use all their energy for maintenance and reproduction (see, e.g., [6], [7], and [13]). The biomass production rate of juveniles and adults depend on the resource abundance. Both foraging rate and metabolic requirements increase with body size.The dynamics of our stage-structured population models consist of three ODEs describing the rates of change of juvenile biomass (J), adult biomass (A), and food resource (R) which are
(1)dJdt=wJR−vwJR−MJR−FJ+wARA,(2)dAdt=vwJRJ−MAR+FA,(3)dRdt=ΦJ,A,R.Here,wJR and wAR are the net biomass production per unit of body mass of juveniles and adults, respectively. Furthermore, MJR and MAR denote the mortality rates of juveniles and adults, respectively. The mortality rates are studied with the inclusion and exclusion of the additional starvation. The interaction function Φ described below depends on what resource dynamic model we choose.The stage-dependent harvesting rates (in this study, to demonstrate the results in an elegant way and for the simplicity of implementing the model, we will assumeuniform harvesting rates, i.e., equal proportional harvesting rates F=FJ=FA on both stages. Compare this with [6], in which they conclude that uniform harvesting rates are a good strategy) of juveniles and adults are both proportional to the corresponding biomasses; the probability constant is F. The maturation rate, vwJR, is the resource-dependent rate at which juveniles mature and become adults. The net biomass production rates for juveniles and adults are given by
(4)wJR=maxw~JR,0,wAR=maxw~AR,0,where
(5)w~JR=ϵImaxRH+R−T,w~AR=ϵqImaxRH+R−T,where ϵ is the efficiency coefficient for the assimilated ingested resource. Here, H is the half-saturation constant of consumers and Imax is the maximum ingestion rate for juveniles. The mass-specific metabolic rate is denoted by T. The ratio between the adult and juvenile feeding rates is denoted by q. The factor q phenomenologically captures stage-specific differences in resource availability and resource use between juveniles and adults; see de Roos and Persson [28] for a detailed explanation. Here, we use a default parameter q=0.85 given in [12]. We set z=sbirth/smax to reduce the number of parameters. In the stage-structured models, the juvenile maturation rate is given by
(6)vwJR=γwJR1−zγ,wJR≠MJ+FJ,−MJ+FJlnz,otherwise,where γ=γR=1−MJ+FJ/wJR. In equations (1)–(6), the constant H is measured in biomass per unit of volume, while the constants T, r, and Imax as well as the mortality rates MA and MJ are expressed per unit of time. All parameters as well as the biomass densities J, A, and R can be considered nondimensional after rescaling.Two models for mortalities are considered. In one model, mortality is not due to starvation which we call natural mortality. In the other model, mortality depends on natural cases as well as starvation. In the model with starvation, the mortality rate functions are given by (we use the same model as [13])
(7)MJR=M0,R>HϵImax/T−1,M0−w~JR,otherwise,(8)MAR=M0,R>HϵqImax/T−1,M0−w~AR,otherwise,where M0 is the natural mortality rate for juvenile and adult individuals. The model we use for morality without starvation (given in [6]) assumes constant mortality rates MJR=MAR=M0. The ecological-mathematical correlations and derivations of expressions (5)–(8) can be found in de Roos et al. [13] where these equations are explained in a clear and instructive manner.


### 2.2. Resource Dynamics

In most aquatic stage-structured models, the resource growth rate is assumed to be of eithersemichemostat growth or logistic growth [29–31]. However, logistic growth might not be suitable for ecological purposes [32], and in this paper, we have found that the simulated solutions become asymptotically periodic solutions under low harvesting rates and asymptotically stable solutions otherwise. The reason for this periodicity is that the logistic growth rate curve has a horizontal tangent line at the origin (see Figure 1). To handle the periodicity in the logistic case, we use the mean value of the juvenile biomass and the adult biomass over a timespan of one solution cycle. The same timespan is used to find the average of the resource density.Figure 1
Graphs for the semichemostat, logistic, threshold-logistic, and semilogistic curves, used in the growth rate models for the resource.Moreover, the obstacle of cyclic solutions that arise in the logistic growth model can be overcome in several ways. To overcome this obstacle, we present and investigate two variations of the logistic model: (1) we assume that the fish are not able to find the resource if the resource density is less than a certain threshold, which we call athreshold-logistic growth model. (2) We assume that the resource is available in all densities with a small influx from the surrounding ecosystem. We call this variant as the semilogistic growth model.As we see in Figure1, all the curves, except the logistic growth curve, have a positive slope at the origin, which will stabilize the solutions in such a way that they approach a steady-state solution in finite time. In the remainder of this paper, our discussion will focus on the four basic growth rate models that are semichemostat, logistic, threshold-logistic, and semilogistic.In alignment with known stochastic growth models, we have converted the aforementioned deterministic growth rate models to stochastic growth rate models by adding a white noise. The stochastic semichemostat growth rate model is the same as the deterministic stochastic growth rate model with the addition of a white noise in terms of a Brownian motion multiplied with a constant volatility. As for the deterministic logistic and threshold-logistic growth rate models, we include stochasticity with the same white noise added, but with volatility proportional to the resource density. The stochastic model for semilogistic growth rate is created by adding a linear combination of the stochastic semichemostat and stochastic logistic growth rate models.


#### 2.2.1. Deterministic Setting

In the stage-structured population models, the growth rate models we consider are semichemostat, logistic, threshold-logistic, and semilogistic growth. The semichemostat growth rate is a natural model when plants or phytoplankton are considered as the food resource, but it may also be appropriate to use for zooplankton and insects when such food resources migrate into the relevant ecosystem. The logistic growth rate may be considered when the resource consists of zooplankton and insects in a closed ecosystem, but in general, the fish population will not be able to reach all of the resource, due to coverage or other types of inaccessibility (see, e.g., [32], for a deeper discussion of these phenomena in ecosystems).In our models, the population stages utilize one unstructured resource. We incorporate the resource models as derived by de Roos et al. [13], which are given by
(9)Semichemostat:dRdt=rRmax−R−ImaxRH+RJ+qA,Logistic:dRdt=rRRmax−R−ImaxRH+RJ+qA,where r is the resource turnover rate and Rmax is the carrying capacity of resource density. In addition, when we simulate the ecosystem with logistic growth, we get asymptotically periodic or stable solutions for the region of small harvesting rates. To handle this periodicity, the mean values of juvenile, adult biomasses, and resource density are used instead of the stable solutions for the deterministic case. To compensate for influx and inaccessibility of the resource, we have also defined the logistic growth rate with threshold and a new growth rate model, which we call logistic-threshold and semilogistic growth, respectively. The semilogistic growth model is a linear combination of the other two models, in which only a small portion of semichemostat growth is used.To stabilize the cyclic solutions that appear in the logistic growth model, we modify the dynamics (see Resource dynamics or Conclusions and Discussion, for a biological interpretation of these modifications) in the following ways:(10)Threshold–logistic:dRdt=rRRmax−R−ImaxR−R0H+R−R0J+qA,(11)Semilogistic:dRdt=prRmax−R+1−prRRmax−R−ImaxRH+RJ+qA.Note that in the threshold-logistic growth model, the functionR−R0 means that the resource available for consumers is the resource density above the threshold R0. When using the semilogistic growth rate model, we use the parameter 0<p<<1 to denote the proportion of the semichemostat growth and 1−p for the proportion of logistic growth.


#### 2.2.2. Stochastic Setting

In this section, the resource dynamics described in Section2.2.1 is modified by including environmental randomness; i.e., we include a stochastic term in the growth rate of the food resource. In the natural environment, it is reasonable to include the randomness on the weather, disease outbreaks, water stress, deforestation, overgrazing, and overcultivation. Basing on these factors, the dynamics of the stochastic food resource are introduced in this paper. We have presented the randomness by Brownian motions (see Section 6 for a discussion about this choice).We extend the abovementioned deterministic resource models by first rewriting the dynamics in differential forms, and secondly, we add appropriate noise terms. In the absence of consumers, four different types of stochastic resource dynamics are as follows:Semichemostat: this type is used in, e.g., Singh [33] (equation (10)). The semichemostat growth rate dynamics is given by
(12)dRt=rRmax−Rdt+ρdWt,where Wt is a Brownian motion and ρ is the standard deviation parameter (in finance, this model is often referred to as the Vasicek model).Logistic and threshold-logistic: these types are used in, e.g., Shah [34] and Roughgarden [35]. The growth rate dynamics for these models are given by
(13)dRt=rRRmax−Rdt+ρRdWt.Note that the same resource growth rate is used in both these models, but when consumers are introduced, not all of the resource is available in the threshold-logistic model (see in Table1 how the consumption of the resource differs in these models). The threshold idea is in line with the argumentations made by Persson et al. [32].Table 1
Summary of resource dynamic models in stochastic setting including consumption from population. The deterministic resource dynamics are given by analogous equations, where terms involvingdW are excluded; that is, we set the volatility ρ=0.
Semichemostat:dRt=rRmax−Rdt+ρdWt−ImaxR/H+Rj+qAdtLogistic:dRt=rRRmax−Rdt+ρRdWt−ImaxR/H+Rj+qAdtThreshold-logistic:dRt=rRRmax−Rdt+ρRdWt−ImaxR−R0/H+R−R0j+qAdtSemilogistic:dRt=r1−pR+pRmax−Rdt+pρdWt1+1−pρRdWt2−ImaxR/H+Rj+qAdtSemilogistic: the resource dynamics in this model is a linear combination of the stochastic semichemostat and the logistic growth rates. Given a constant 0<p≪1, we define the semilogistic dynamics as
(14)dRt=prRmax−Rdt+ρdWt1+1−prRRmax−Rdt+ρRdWt2=r1−pR+pRmax−Rdt+pρdWt1+1−pρRdWt2,where Wt1 and Wt2 are independent Brownian motions.When consumers are introduced to these models, the rate of change of the available resource biomasses,R, is given by the stochastic differential equations in Table 1.A bit of caution might be needed here. The stochastic models have quite a complicated feedback interaction between the juveniles, the adults, and the resource, and the expectation of the solutions in the stochastic models will thus not always be the solution of the corresponding deterministic models. If the models had a simpler setting, we would be able to compensate for this difference using techniques similar to the theory developed in, e.g., Giet et al. [36].


## 2.1. Population Stage Dynamics

The individuals of the species are modeled by a stage-structured biomass model which is derived by formulating a size-structured population model on the individuals’ life history. More precisely, Individuals are divided into two stages: juveniles and adults, based only on their size. For many relevant properties of the species, a two-stage model is often enough; see, e.g., Ackleh and Jang [8], de Roos et al. [13], Lundström et al. [6], and Meng et al. [7]. Individuals are assumed to have the same size sbirth at birth, and the maximum size of individuals is denoted by smax. Both stages forage for the shared resource, R=Rt, and the metabolic requirement, T, is assumed to be constant. The growth rate depends on the available resource ([25]) and varies with population density ([26, 27]). Juvenile and adult biomasses are denoted by J=Jt and A=At, respectively.The derivation of the stage-structured biomass model is investigated by formulating a size-structured population model. We give a brief description of the model and refer the reader to de Roos et al. [13] for further details and motivations. The population is composed of juveniles, J, and adults, A, which forages on a shared resource R. Juvenile individuals use all the consumed energy for growth, development, and maintenance, whereas adult individuals use all their energy for maintenance and reproduction (see, e.g., [6], [7], and [13]). The biomass production rate of juveniles and adults depend on the resource abundance. Both foraging rate and metabolic requirements increase with body size.The dynamics of our stage-structured population models consist of three ODEs describing the rates of change of juvenile biomass (J), adult biomass (A), and food resource (R) which are
(1)dJdt=wJR−vwJR−MJR−FJ+wARA,(2)dAdt=vwJRJ−MAR+FA,(3)dRdt=ΦJ,A,R.Here,wJR and wAR are the net biomass production per unit of body mass of juveniles and adults, respectively. Furthermore, MJR and MAR denote the mortality rates of juveniles and adults, respectively. The mortality rates are studied with the inclusion and exclusion of the additional starvation. The interaction function Φ described below depends on what resource dynamic model we choose.The stage-dependent harvesting rates (in this study, to demonstrate the results in an elegant way and for the simplicity of implementing the model, we will assumeuniform harvesting rates, i.e., equal proportional harvesting rates F=FJ=FA on both stages. Compare this with [6], in which they conclude that uniform harvesting rates are a good strategy) of juveniles and adults are both proportional to the corresponding biomasses; the probability constant is F. The maturation rate, vwJR, is the resource-dependent rate at which juveniles mature and become adults. The net biomass production rates for juveniles and adults are given by
(4)wJR=maxw~JR,0,wAR=maxw~AR,0,where
(5)w~JR=ϵImaxRH+R−T,w~AR=ϵqImaxRH+R−T,where ϵ is the efficiency coefficient for the assimilated ingested resource. Here, H is the half-saturation constant of consumers and Imax is the maximum ingestion rate for juveniles. The mass-specific metabolic rate is denoted by T. The ratio between the adult and juvenile feeding rates is denoted by q. The factor q phenomenologically captures stage-specific differences in resource availability and resource use between juveniles and adults; see de Roos and Persson [28] for a detailed explanation. Here, we use a default parameter q=0.85 given in [12]. We set z=sbirth/smax to reduce the number of parameters. In the stage-structured models, the juvenile maturation rate is given by
(6)vwJR=γwJR1−zγ,wJR≠MJ+FJ,−MJ+FJlnz,otherwise,where γ=γR=1−MJ+FJ/wJR. In equations (1)–(6), the constant H is measured in biomass per unit of volume, while the constants T, r, and Imax as well as the mortality rates MA and MJ are expressed per unit of time. All parameters as well as the biomass densities J, A, and R can be considered nondimensional after rescaling.Two models for mortalities are considered. In one model, mortality is not due to starvation which we call natural mortality. In the other model, mortality depends on natural cases as well as starvation. In the model with starvation, the mortality rate functions are given by (we use the same model as [13])
(7)MJR=M0,R>HϵImax/T−1,M0−w~JR,otherwise,(8)MAR=M0,R>HϵqImax/T−1,M0−w~AR,otherwise,where M0 is the natural mortality rate for juvenile and adult individuals. The model we use for morality without starvation (given in [6]) assumes constant mortality rates MJR=MAR=M0. The ecological-mathematical correlations and derivations of expressions (5)–(8) can be found in de Roos et al. [13] where these equations are explained in a clear and instructive manner.


## 2.2. Resource Dynamics

In most aquatic stage-structured models, the resource growth rate is assumed to be of eithersemichemostat growth or logistic growth [29–31]. However, logistic growth might not be suitable for ecological purposes [32], and in this paper, we have found that the simulated solutions become asymptotically periodic solutions under low harvesting rates and asymptotically stable solutions otherwise. The reason for this periodicity is that the logistic growth rate curve has a horizontal tangent line at the origin (see Figure 1). To handle the periodicity in the logistic case, we use the mean value of the juvenile biomass and the adult biomass over a timespan of one solution cycle. The same timespan is used to find the average of the resource density.Figure 1
Graphs for the semichemostat, logistic, threshold-logistic, and semilogistic curves, used in the growth rate models for the resource.Moreover, the obstacle of cyclic solutions that arise in the logistic growth model can be overcome in several ways. To overcome this obstacle, we present and investigate two variations of the logistic model: (1) we assume that the fish are not able to find the resource if the resource density is less than a certain threshold, which we call athreshold-logistic growth model. (2) We assume that the resource is available in all densities with a small influx from the surrounding ecosystem. We call this variant as the semilogistic growth model.As we see in Figure1, all the curves, except the logistic growth curve, have a positive slope at the origin, which will stabilize the solutions in such a way that they approach a steady-state solution in finite time. In the remainder of this paper, our discussion will focus on the four basic growth rate models that are semichemostat, logistic, threshold-logistic, and semilogistic.In alignment with known stochastic growth models, we have converted the aforementioned deterministic growth rate models to stochastic growth rate models by adding a white noise. The stochastic semichemostat growth rate model is the same as the deterministic stochastic growth rate model with the addition of a white noise in terms of a Brownian motion multiplied with a constant volatility. As for the deterministic logistic and threshold-logistic growth rate models, we include stochasticity with the same white noise added, but with volatility proportional to the resource density. The stochastic model for semilogistic growth rate is created by adding a linear combination of the stochastic semichemostat and stochastic logistic growth rate models.


### 2.2.1. Deterministic Setting

In the stage-structured population models, the growth rate models we consider are semichemostat, logistic, threshold-logistic, and semilogistic growth. The semichemostat growth rate is a natural model when plants or phytoplankton are considered as the food resource, but it may also be appropriate to use for zooplankton and insects when such food resources migrate into the relevant ecosystem. The logistic growth rate may be considered when the resource consists of zooplankton and insects in a closed ecosystem, but in general, the fish population will not be able to reach all of the resource, due to coverage or other types of inaccessibility (see, e.g., [32], for a deeper discussion of these phenomena in ecosystems).In our models, the population stages utilize one unstructured resource. We incorporate the resource models as derived by de Roos et al. [13], which are given by
(9)Semichemostat:dRdt=rRmax−R−ImaxRH+RJ+qA,Logistic:dRdt=rRRmax−R−ImaxRH+RJ+qA,where r is the resource turnover rate and Rmax is the carrying capacity of resource density. In addition, when we simulate the ecosystem with logistic growth, we get asymptotically periodic or stable solutions for the region of small harvesting rates. To handle this periodicity, the mean values of juvenile, adult biomasses, and resource density are used instead of the stable solutions for the deterministic case. To compensate for influx and inaccessibility of the resource, we have also defined the logistic growth rate with threshold and a new growth rate model, which we call logistic-threshold and semilogistic growth, respectively. The semilogistic growth model is a linear combination of the other two models, in which only a small portion of semichemostat growth is used.To stabilize the cyclic solutions that appear in the logistic growth model, we modify the dynamics (see Resource dynamics or Conclusions and Discussion, for a biological interpretation of these modifications) in the following ways:(10)Threshold–logistic:dRdt=rRRmax−R−ImaxR−R0H+R−R0J+qA,(11)Semilogistic:dRdt=prRmax−R+1−prRRmax−R−ImaxRH+RJ+qA.Note that in the threshold-logistic growth model, the functionR−R0 means that the resource available for consumers is the resource density above the threshold R0. When using the semilogistic growth rate model, we use the parameter 0<p<<1 to denote the proportion of the semichemostat growth and 1−p for the proportion of logistic growth.


### 2.2.2. Stochastic Setting

In this section, the resource dynamics described in Section2.2.1 is modified by including environmental randomness; i.e., we include a stochastic term in the growth rate of the food resource. In the natural environment, it is reasonable to include the randomness on the weather, disease outbreaks, water stress, deforestation, overgrazing, and overcultivation. Basing on these factors, the dynamics of the stochastic food resource are introduced in this paper. We have presented the randomness by Brownian motions (see Section 6 for a discussion about this choice).We extend the abovementioned deterministic resource models by first rewriting the dynamics in differential forms, and secondly, we add appropriate noise terms. In the absence of consumers, four different types of stochastic resource dynamics are as follows:Semichemostat: this type is used in, e.g., Singh [33] (equation (10)). The semichemostat growth rate dynamics is given by
(12)dRt=rRmax−Rdt+ρdWt,where Wt is a Brownian motion and ρ is the standard deviation parameter (in finance, this model is often referred to as the Vasicek model).Logistic and threshold-logistic: these types are used in, e.g., Shah [34] and Roughgarden [35]. The growth rate dynamics for these models are given by
(13)dRt=rRRmax−Rdt+ρRdWt.Note that the same resource growth rate is used in both these models, but when consumers are introduced, not all of the resource is available in the threshold-logistic model (see in Table1 how the consumption of the resource differs in these models). The threshold idea is in line with the argumentations made by Persson et al. [32].Table 1
Summary of resource dynamic models in stochastic setting including consumption from population. The deterministic resource dynamics are given by analogous equations, where terms involvingdW are excluded; that is, we set the volatility ρ=0.
Semichemostat:dRt=rRmax−Rdt+ρdWt−ImaxR/H+Rj+qAdtLogistic:dRt=rRRmax−Rdt+ρRdWt−ImaxR/H+Rj+qAdtThreshold-logistic:dRt=rRRmax−Rdt+ρRdWt−ImaxR−R0/H+R−R0j+qAdtSemilogistic:dRt=r1−pR+pRmax−Rdt+pρdWt1+1−pρRdWt2−ImaxR/H+Rj+qAdtSemilogistic: the resource dynamics in this model is a linear combination of the stochastic semichemostat and the logistic growth rates. Given a constant 0<p≪1, we define the semilogistic dynamics as
(14)dRt=prRmax−Rdt+ρdWt1+1−prRRmax−Rdt+ρRdWt2=r1−pR+pRmax−Rdt+pρdWt1+1−pρRdWt2,where Wt1 and Wt2 are independent Brownian motions.When consumers are introduced to these models, the rate of change of the available resource biomasses,R, is given by the stochastic differential equations in Table 1.A bit of caution might be needed here. The stochastic models have quite a complicated feedback interaction between the juveniles, the adults, and the resource, and the expectation of the solutions in the stochastic models will thus not always be the solution of the corresponding deterministic models. If the models had a simpler setting, we would be able to compensate for this difference using techniques similar to the theory developed in, e.g., Giet et al. [36].


## 2.2.1. Deterministic Setting

In the stage-structured population models, the growth rate models we consider are semichemostat, logistic, threshold-logistic, and semilogistic growth. The semichemostat growth rate is a natural model when plants or phytoplankton are considered as the food resource, but it may also be appropriate to use for zooplankton and insects when such food resources migrate into the relevant ecosystem. The logistic growth rate may be considered when the resource consists of zooplankton and insects in a closed ecosystem, but in general, the fish population will not be able to reach all of the resource, due to coverage or other types of inaccessibility (see, e.g., [32], for a deeper discussion of these phenomena in ecosystems).In our models, the population stages utilize one unstructured resource. We incorporate the resource models as derived by de Roos et al. [13], which are given by
(9)Semichemostat:dRdt=rRmax−R−ImaxRH+RJ+qA,Logistic:dRdt=rRRmax−R−ImaxRH+RJ+qA,where r is the resource turnover rate and Rmax is the carrying capacity of resource density. In addition, when we simulate the ecosystem with logistic growth, we get asymptotically periodic or stable solutions for the region of small harvesting rates. To handle this periodicity, the mean values of juvenile, adult biomasses, and resource density are used instead of the stable solutions for the deterministic case. To compensate for influx and inaccessibility of the resource, we have also defined the logistic growth rate with threshold and a new growth rate model, which we call logistic-threshold and semilogistic growth, respectively. The semilogistic growth model is a linear combination of the other two models, in which only a small portion of semichemostat growth is used.To stabilize the cyclic solutions that appear in the logistic growth model, we modify the dynamics (see Resource dynamics or Conclusions and Discussion, for a biological interpretation of these modifications) in the following ways:(10)Threshold–logistic:dRdt=rRRmax−R−ImaxR−R0H+R−R0J+qA,(11)Semilogistic:dRdt=prRmax−R+1−prRRmax−R−ImaxRH+RJ+qA.Note that in the threshold-logistic growth model, the functionR−R0 means that the resource available for consumers is the resource density above the threshold R0. When using the semilogistic growth rate model, we use the parameter 0<p<<1 to denote the proportion of the semichemostat growth and 1−p for the proportion of logistic growth.


## 2.2.2. Stochastic Setting

In this section, the resource dynamics described in Section2.2.1 is modified by including environmental randomness; i.e., we include a stochastic term in the growth rate of the food resource. In the natural environment, it is reasonable to include the randomness on the weather, disease outbreaks, water stress, deforestation, overgrazing, and overcultivation. Basing on these factors, the dynamics of the stochastic food resource are introduced in this paper. We have presented the randomness by Brownian motions (see Section 6 for a discussion about this choice).We extend the abovementioned deterministic resource models by first rewriting the dynamics in differential forms, and secondly, we add appropriate noise terms. In the absence of consumers, four different types of stochastic resource dynamics are as follows:Semichemostat: this type is used in, e.g., Singh [33] (equation (10)). The semichemostat growth rate dynamics is given by
(12)dRt=rRmax−Rdt+ρdWt,where Wt is a Brownian motion and ρ is the standard deviation parameter (in finance, this model is often referred to as the Vasicek model).Logistic and threshold-logistic: these types are used in, e.g., Shah [34] and Roughgarden [35]. The growth rate dynamics for these models are given by
(13)dRt=rRRmax−Rdt+ρRdWt.Note that the same resource growth rate is used in both these models, but when consumers are introduced, not all of the resource is available in the threshold-logistic model (see in Table1 how the consumption of the resource differs in these models). The threshold idea is in line with the argumentations made by Persson et al. [32].Table 1
Summary of resource dynamic models in stochastic setting including consumption from population. The deterministic resource dynamics are given by analogous equations, where terms involvingdW are excluded; that is, we set the volatility ρ=0.
Semichemostat:dRt=rRmax−Rdt+ρdWt−ImaxR/H+Rj+qAdtLogistic:dRt=rRRmax−Rdt+ρRdWt−ImaxR/H+Rj+qAdtThreshold-logistic:dRt=rRRmax−Rdt+ρRdWt−ImaxR−R0/H+R−R0j+qAdtSemilogistic:dRt=r1−pR+pRmax−Rdt+pρdWt1+1−pρRdWt2−ImaxR/H+Rj+qAdtSemilogistic: the resource dynamics in this model is a linear combination of the stochastic semichemostat and the logistic growth rates. Given a constant 0<p≪1, we define the semilogistic dynamics as
(14)dRt=prRmax−Rdt+ρdWt1+1−prRRmax−Rdt+ρRdWt2=r1−pR+pRmax−Rdt+pρdWt1+1−pρRdWt2,where Wt1 and Wt2 are independent Brownian motions.When consumers are introduced to these models, the rate of change of the available resource biomasses,R, is given by the stochastic differential equations in Table 1.A bit of caution might be needed here. The stochastic models have quite a complicated feedback interaction between the juveniles, the adults, and the resource, and the expectation of the solutions in the stochastic models will thus not always be the solution of the corresponding deterministic models. If the models had a simpler setting, we would be able to compensate for this difference using techniques similar to the theory developed in, e.g., Giet et al. [36].


## 3. Population Properties of the Stage-Structured Model

### 3.1. Yield

The stage-structured models are examined by introducing a wide range of equal harvesting rates of juveniles and adults with different consequences for the yield. The yield is defined as the continuous outtake of biomass by harvesting. The juvenile and adult biomasses at equilibrium are denoted byJ∗=J∗FJ,FA and A∗=A∗FJ,FA, respectively, where FJ and FA are the proportional harvesting rates. The yield objective function is defined by
(15)Yield=FJJ∗+FAA∗.In the absence of harvesting,Ju∗=J∗0,0 and Au∗=A∗0,0 denote the juvenile and adult biomasses, respectively, at equilibrium.In the stochastic setting, the yield is calculated, using equation (15), where J∗ and A∗ are replaced by the mean biomasses of juveniles and adults.


### 3.2. Impact on Biomass and Size Structure

For the models abovementioned, we investigate the impact of harvesting on the consumer population. In addition to harvesting, the biomasses of juveniles and adults also depend on the amount of resource which in turn decreases due to the consorted foraging of all consumers on it.The impact on biomass is a measure of how much biomass of the population has decreased with respect to the harvesting rate. It is defined as(16)Impactonbiomass=1−B∗Bu∗,where B∗=J∗+A∗ and Bu∗=Ju∗+Au∗. The value of the impact on biomass is zero when no harvesting is done and increases up to one when the harvesting rate depletes the population. In some of our models, the impact on biomass temporary decreases as harvesting rate increases, which implies an increase in the population biomass; that is, these models exhibit the hydra effect. Similarly, the impact on size structure measures the relative change of biomasses of the adult population versus the juvenile population with respect to different harvesting rates. It is evaluated by the expression
(17)Impactonsizestructure=J∗J∗+A∗Ju∗Ju∗+Au∗−1−1.If the impact on size structure is positive at a certain harvesting rate, then the fraction of juveniles in the population has increased compared to the nonharvesting scenario; cf. Lundström et al. [6].Moreover, the impact on biomass of harvesting in the stochastic case is measured by the expression(18)Impactonbiomass=1−EB∗EBu∗,where EB∗=EB∗F is the steady state of the expectation of the total biomass under a harvesting rate F and EBu∗=EBu∗0.The impact on size structure for stochastic case is explored in a similar way using equation (17).


### 3.3. Resilience

Resilience is one of the important components of stability of an ecosystem. It is a measure of how fast the ecosystem recovers after population perturbation. It is strongly influenced by the types of environmental fluctuations commonly encountered by an ecosystem [37].We consider the resilience of the population, for both deterministic and stochastic cases by measuring the reciprocal of the time needed for the population to recover to the positive equilibrium given a random perturbation [6]. We denote the equilibrium of biomasses as J∗,A∗,R∗. Let κ>0; this constant scales the maximum displacement of the population from the equilibrium. A trajectory of the stage model is started from a random point uniformly distributed in the cube 0,κJ∗×0,κA∗×0,κR∗.We then find the return time as the time needed for this trajectory to be close enough to the equilibrium in the sense that(19)J∗−JtJ∗2+A∗−AtA∗2+R∗−RtR∗21/2≤εd,for some small error bound εd>0 in the deterministic case. To find the resilience, we use the definition introduced in [6]. That is, after repeating this procedure N times, the resilience is defined by
(20)Resilience=1averagevalueofthereturntimes.The higher the resilience, the smaller the risk of extinction due to random drift [6]. To estimate the resilience in the stochastic case, we find the mean value of the biomasses of the resource, the juveniles, and the adults. This is done by performing a number of simulations using the stochastic model, in which the initial values are randomly picked in the cube defined above. The value of the error bound for the stochastic case, εs, has to be larger than the one for the deterministic case according to the variance of the mean values. We then find the return time to be the time it takes for the solution to come close enough to the equilibrium in the sense of equation (19) by using the mean values and εs in place of εd. This procedure is repeated to find the average value of return times. The average return time is used to find the resilience by again using equation (20). Finally, many samples of the resilience are drawn and used to find its mean and standard deviation.


### 3.4. Recovery Potential

We consider the basic reproduction ratio which represents the average number of offspring produced over the lifetime of an individual in the absence of density-dependent competition. The measure of recovery potential is closely related to the basic reproduction ratio in a virgin environment. For the deterministic stage model, the recovery potential, introduced in [7], is defined by using the steady-state equation as follows:
(21)gR=wARM+FA×νwJRνwJR−wJR+M+FJ.In a virgin environment, Meng et al. [7] defined the recovery potential, as a function of harvesting rates, given by the expression
(22)Recoverypotential=wARmaxM+FA×νwJRmaxνwJRmax−wJRmax+M+FJ.Meng et al. proved that a unique positive equilibrium of resource, juvenile, and adult biomass density exists when the recovery potential is larger than 1. However, extinction of the population follows when the recovery potential is smaller than 1. They derived this recovery potential, equation (22), by assuming A′=J′=0 in equations (1) and (2). However, in the stochastic approach, we will never reach an equilibrium and therefore, we will not be able to use equation (22) to find the recovery potential. Therefore, we will derive a similar expression for the recovery potential of our models in the stochastic setting. This new recovery potential expression coincides with the old definition, which we have corroborated in the Appendix; see Figures 2 and 3.Figure 2
The recovery potential is calculated by using the algorithm described in AppendixB.2. Meng et al. recovery potential (red curve) and nonequilibrium recovery potential (blue curve) w.r.t. harvesting rate for deterministic case (a, c). The mean value of the stochastic nonequilibrium recovery potential (yellow curve) and the interval of standard deviation (68%) of the stochastic nonequilibrium recovery potential (light blue curves) w.r.t. harvesting rate (b, d).
(a)(b)(c)(d)Figure 3
The recovery potential is calculated by using the algorithm described in AppendixB.2. Meng et al. recovery potential (red curve) and nonequilibrium recovery potential (blue curve) w.r.t. harvesting rate for deterministic case (a, c). The mean value of the stochastic nonequilibrium recovery potential (yellow curve) and the interval of standard deviation (68%) of the stochastic nonequilibrium recovery potential (light blue curves) w.r.t. harvesting rate (b, d).
(a)(b)(c)(d)We consider the following equations for the rates at which the biomass of juveniles and adults changes:(23)J′=wARA+wJR−vwJR−MJ−FJJ,(24)A′=vwJRJ−MA+FAA.Rearranging equation (24) yields
(25)vwJR=A′+MA+FAAJ.Note that the right-hand side does not explicitly depend onR. In what follows, our goal is to derive an expression for the recovery potential that does not include R. When this has been achieved, we can use the expectation of the recovery potential to find the probability of extinction. When wJR≠MJ+FJ, we substitute equation (25) into equation (6) yielding
(26)A′+MA+FAAJ=γwJR1−zγ.The right-hand side of this equation is an increasing continuous function with respect toR (see Appendix A.1) and hence, we can find the unique solution wJ∗R of equation (26).By substituting equation (25) and wJ∗R into equation (23), we get
(27)wAR=J′A−wJ∗RJA+A′+MA+FAAA+MJ+FJJA=J′−wJ∗R−MJ−FJJ+A′+MA+FAAA.Recall that, in this paper, we assumeF=FA=FJ. By using equations (25) and (27) and the unique solution wJ∗R, we have derived an alternative expression for the recovery potential, defined by equation (22), as
(28)RF=J′−wJ∗R−MJ−FJ+A′+MA+FAAMA+F×A′+MA+FAA′+MA+FA−wJ∗R−MJ−FJ,when wJ∗R≠MJ+F.In the casewJ∗R=MJ+F, equation (27) may be written as
(29)wAR=J′+A′+MA+FAA.By using equations (25) and (29) in equation (22), we get the recovery potential
(30)RF=J′+A′+MA+FAAMA+F×A′+MA+FAA′+MA+FA−MJ+F−MJ−FJ=J′+A′+MA+FAAMA+F.Equations (28) and (30) can be used directly in the deterministic case, but in the stochastic case, we use equations (28) and (30) to find the mean and standard deviation of the recovery potential.


### 3.5. Probability of Extinction for Stochastic Case

An important feature in population dynamics is the possibility of extinction, which has been studied in a variety of stochastic model formulations [34, 38, 39]. Population size can cause the extinction of a species through overharvesting, habitat distribution, disasters, and other factors [40–42]. One way of finding the probability of extinction is through the minimum viable population (MVP).TheMVP is an estimate of the minimum number of individuals in the population that is capable of persisting in the wild life [43, 19]. In our models, we define the probability of extinction by the probability that, during any time,
(31)A+pAJ<MVP,where pA is the survival probability for juveniles to become adults. The expression for pA is given by
(32)pA=pAF=e−M0+F/2vwRmax.See AppendixA.1 for a detailed explanation of this formula.In this paper, we present an alternative approach to investigate the probability of extinction, which is based on the results of recovery potential. The recovery potential provides an indication of the population’s probability of surviving a potential extinction caused by, e.g., environmental stochasticity or overexploitment. In other words, the biomass of an initially small population will increase in a virgin environment when the recovery potential is larger than 1 and the population will go extinct when the recovery potential is smaller than 1. We then find the probability of extinction, which is calculated by investigating the mean of recovery potential.To do this, we use initial dataR=Rmax, J=J0, and A=A0, where J0 and A0 are close to zero. We then simulate systems (1)–(3) to find the solution for two short time steps and use the central difference quotient to find an estimate for J′ and A′. The simulations are done for a range of harvesting rates F=F1,F2,⋯,FN. Since the model is stochastic, we repeat this procedure many times and use the central limit theorem to estimate the recovery potential of our models. That is, we estimate the expected recovery potential,
(33)μRF=ERF≈RF¯≕μ^RF,where RF¯ is the mean value of the simulated recovery potentials for fixed harvesting rate F. The standard deviation, σR, of the recovery potential is estimated by σ^R in a similar way.In view of the central limit theorem, the mean of the recovery potential with harvesting rateF is assumed to be a random sample from a normal distribution with cumulative distribution function, CDFx,μ,σ. The constants μ and σ represent mean and standard deviation, respectively. That is, we assume that
(34)RF¯~Nμ^RF,σ^RFand get the probability of extinction by
(35)PRF<1=CDF1,μ^RF,σ^RF.The reason for evaluating the cumulative distribution function atx=1 is that if the recovery potential is smaller than 1, then the population goes extinct; compare this with Meng et al. [7] for the deterministic setting.


## 3.1. Yield

The stage-structured models are examined by introducing a wide range of equal harvesting rates of juveniles and adults with different consequences for the yield. The yield is defined as the continuous outtake of biomass by harvesting. The juvenile and adult biomasses at equilibrium are denoted byJ∗=J∗FJ,FA and A∗=A∗FJ,FA, respectively, where FJ and FA are the proportional harvesting rates. The yield objective function is defined by
(15)Yield=FJJ∗+FAA∗.In the absence of harvesting,Ju∗=J∗0,0 and Au∗=A∗0,0 denote the juvenile and adult biomasses, respectively, at equilibrium.In the stochastic setting, the yield is calculated, using equation (15), where J∗ and A∗ are replaced by the mean biomasses of juveniles and adults.


## 3.2. Impact on Biomass and Size Structure

For the models abovementioned, we investigate the impact of harvesting on the consumer population. In addition to harvesting, the biomasses of juveniles and adults also depend on the amount of resource which in turn decreases due to the consorted foraging of all consumers on it.The impact on biomass is a measure of how much biomass of the population has decreased with respect to the harvesting rate. It is defined as(16)Impactonbiomass=1−B∗Bu∗,where B∗=J∗+A∗ and Bu∗=Ju∗+Au∗. The value of the impact on biomass is zero when no harvesting is done and increases up to one when the harvesting rate depletes the population. In some of our models, the impact on biomass temporary decreases as harvesting rate increases, which implies an increase in the population biomass; that is, these models exhibit the hydra effect. Similarly, the impact on size structure measures the relative change of biomasses of the adult population versus the juvenile population with respect to different harvesting rates. It is evaluated by the expression
(17)Impactonsizestructure=J∗J∗+A∗Ju∗Ju∗+Au∗−1−1.If the impact on size structure is positive at a certain harvesting rate, then the fraction of juveniles in the population has increased compared to the nonharvesting scenario; cf. Lundström et al. [6].Moreover, the impact on biomass of harvesting in the stochastic case is measured by the expression(18)Impactonbiomass=1−EB∗EBu∗,where EB∗=EB∗F is the steady state of the expectation of the total biomass under a harvesting rate F and EBu∗=EBu∗0.The impact on size structure for stochastic case is explored in a similar way using equation (17).


## 3.3. Resilience

Resilience is one of the important components of stability of an ecosystem. It is a measure of how fast the ecosystem recovers after population perturbation. It is strongly influenced by the types of environmental fluctuations commonly encountered by an ecosystem [37].We consider the resilience of the population, for both deterministic and stochastic cases by measuring the reciprocal of the time needed for the population to recover to the positive equilibrium given a random perturbation [6]. We denote the equilibrium of biomasses as J∗,A∗,R∗. Let κ>0; this constant scales the maximum displacement of the population from the equilibrium. A trajectory of the stage model is started from a random point uniformly distributed in the cube 0,κJ∗×0,κA∗×0,κR∗.We then find the return time as the time needed for this trajectory to be close enough to the equilibrium in the sense that(19)J∗−JtJ∗2+A∗−AtA∗2+R∗−RtR∗21/2≤εd,for some small error bound εd>0 in the deterministic case. To find the resilience, we use the definition introduced in [6]. That is, after repeating this procedure N times, the resilience is defined by
(20)Resilience=1averagevalueofthereturntimes.The higher the resilience, the smaller the risk of extinction due to random drift [6]. To estimate the resilience in the stochastic case, we find the mean value of the biomasses of the resource, the juveniles, and the adults. This is done by performing a number of simulations using the stochastic model, in which the initial values are randomly picked in the cube defined above. The value of the error bound for the stochastic case, εs, has to be larger than the one for the deterministic case according to the variance of the mean values. We then find the return time to be the time it takes for the solution to come close enough to the equilibrium in the sense of equation (19) by using the mean values and εs in place of εd. This procedure is repeated to find the average value of return times. The average return time is used to find the resilience by again using equation (20). Finally, many samples of the resilience are drawn and used to find its mean and standard deviation.


## 3.4. Recovery Potential

We consider the basic reproduction ratio which represents the average number of offspring produced over the lifetime of an individual in the absence of density-dependent competition. The measure of recovery potential is closely related to the basic reproduction ratio in a virgin environment. For the deterministic stage model, the recovery potential, introduced in [7], is defined by using the steady-state equation as follows:
(21)gR=wARM+FA×νwJRνwJR−wJR+M+FJ.In a virgin environment, Meng et al. [7] defined the recovery potential, as a function of harvesting rates, given by the expression
(22)Recoverypotential=wARmaxM+FA×νwJRmaxνwJRmax−wJRmax+M+FJ.Meng et al. proved that a unique positive equilibrium of resource, juvenile, and adult biomass density exists when the recovery potential is larger than 1. However, extinction of the population follows when the recovery potential is smaller than 1. They derived this recovery potential, equation (22), by assuming A′=J′=0 in equations (1) and (2). However, in the stochastic approach, we will never reach an equilibrium and therefore, we will not be able to use equation (22) to find the recovery potential. Therefore, we will derive a similar expression for the recovery potential of our models in the stochastic setting. This new recovery potential expression coincides with the old definition, which we have corroborated in the Appendix; see Figures 2 and 3.Figure 2
The recovery potential is calculated by using the algorithm described in AppendixB.2. Meng et al. recovery potential (red curve) and nonequilibrium recovery potential (blue curve) w.r.t. harvesting rate for deterministic case (a, c). The mean value of the stochastic nonequilibrium recovery potential (yellow curve) and the interval of standard deviation (68%) of the stochastic nonequilibrium recovery potential (light blue curves) w.r.t. harvesting rate (b, d).
(a)(b)(c)(d)Figure 3
The recovery potential is calculated by using the algorithm described in AppendixB.2. Meng et al. recovery potential (red curve) and nonequilibrium recovery potential (blue curve) w.r.t. harvesting rate for deterministic case (a, c). The mean value of the stochastic nonequilibrium recovery potential (yellow curve) and the interval of standard deviation (68%) of the stochastic nonequilibrium recovery potential (light blue curves) w.r.t. harvesting rate (b, d).
(a)(b)(c)(d)We consider the following equations for the rates at which the biomass of juveniles and adults changes:(23)J′=wARA+wJR−vwJR−MJ−FJJ,(24)A′=vwJRJ−MA+FAA.Rearranging equation (24) yields
(25)vwJR=A′+MA+FAAJ.Note that the right-hand side does not explicitly depend onR. In what follows, our goal is to derive an expression for the recovery potential that does not include R. When this has been achieved, we can use the expectation of the recovery potential to find the probability of extinction. When wJR≠MJ+FJ, we substitute equation (25) into equation (6) yielding
(26)A′+MA+FAAJ=γwJR1−zγ.The right-hand side of this equation is an increasing continuous function with respect toR (see Appendix A.1) and hence, we can find the unique solution wJ∗R of equation (26).By substituting equation (25) and wJ∗R into equation (23), we get
(27)wAR=J′A−wJ∗RJA+A′+MA+FAAA+MJ+FJJA=J′−wJ∗R−MJ−FJJ+A′+MA+FAAA.Recall that, in this paper, we assumeF=FA=FJ. By using equations (25) and (27) and the unique solution wJ∗R, we have derived an alternative expression for the recovery potential, defined by equation (22), as
(28)RF=J′−wJ∗R−MJ−FJ+A′+MA+FAAMA+F×A′+MA+FAA′+MA+FA−wJ∗R−MJ−FJ,when wJ∗R≠MJ+F.In the casewJ∗R=MJ+F, equation (27) may be written as
(29)wAR=J′+A′+MA+FAA.By using equations (25) and (29) in equation (22), we get the recovery potential
(30)RF=J′+A′+MA+FAAMA+F×A′+MA+FAA′+MA+FA−MJ+F−MJ−FJ=J′+A′+MA+FAAMA+F.Equations (28) and (30) can be used directly in the deterministic case, but in the stochastic case, we use equations (28) and (30) to find the mean and standard deviation of the recovery potential.


## 3.5. Probability of Extinction for Stochastic Case

An important feature in population dynamics is the possibility of extinction, which has been studied in a variety of stochastic model formulations [34, 38, 39]. Population size can cause the extinction of a species through overharvesting, habitat distribution, disasters, and other factors [40–42]. One way of finding the probability of extinction is through the minimum viable population (MVP).TheMVP is an estimate of the minimum number of individuals in the population that is capable of persisting in the wild life [43, 19]. In our models, we define the probability of extinction by the probability that, during any time,
(31)A+pAJ<MVP,where pA is the survival probability for juveniles to become adults. The expression for pA is given by
(32)pA=pAF=e−M0+F/2vwRmax.See AppendixA.1 for a detailed explanation of this formula.In this paper, we present an alternative approach to investigate the probability of extinction, which is based on the results of recovery potential. The recovery potential provides an indication of the population’s probability of surviving a potential extinction caused by, e.g., environmental stochasticity or overexploitment. In other words, the biomass of an initially small population will increase in a virgin environment when the recovery potential is larger than 1 and the population will go extinct when the recovery potential is smaller than 1. We then find the probability of extinction, which is calculated by investigating the mean of recovery potential.To do this, we use initial dataR=Rmax, J=J0, and A=A0, where J0 and A0 are close to zero. We then simulate systems (1)–(3) to find the solution for two short time steps and use the central difference quotient to find an estimate for J′ and A′. The simulations are done for a range of harvesting rates F=F1,F2,⋯,FN. Since the model is stochastic, we repeat this procedure many times and use the central limit theorem to estimate the recovery potential of our models. That is, we estimate the expected recovery potential,
(33)μRF=ERF≈RF¯≕μ^RF,where RF¯ is the mean value of the simulated recovery potentials for fixed harvesting rate F. The standard deviation, σR, of the recovery potential is estimated by σ^R in a similar way.In view of the central limit theorem, the mean of the recovery potential with harvesting rateF is assumed to be a random sample from a normal distribution with cumulative distribution function, CDFx,μ,σ. The constants μ and σ represent mean and standard deviation, respectively. That is, we assume that
(34)RF¯~Nμ^RF,σ^RFand get the probability of extinction by
(35)PRF<1=CDF1,μ^RF,σ^RF.The reason for evaluating the cumulative distribution function atx=1 is that if the recovery potential is smaller than 1, then the population goes extinct; compare this with Meng et al. [7] for the deterministic setting.


## 4. Simulation Results

The emergent properties are investigated in the stage-structured biomass models, systems (1)–(3), including starvation mortality. The properties we study in this section are resilience, recovery potential, and probability of extinction. These properties are compared between the different models, under a range of harvesting rates. In particular, the models are differentiated by the growth model for the resource dynamics; semichemostat growth, logistic growth, semilogistic growth, and threshold-logistic growth. The models are also differentiated by using deterministic vs. stochastic growth. In Appendix C, the steady-state biomasses, yield, impact on biomass, and impact on size structure for these growth models are investigated numerically.From here on, all parameters as well as the biomass densities for the stage-structured biomass models are considered nondimensional, using rescaling as in [6, 13]. In our simulations, we use the parameters given in Table 2.Table 2
Ecological and economic parameters. The columnUnit is the unit dimensions before rescaling into nondimensional parameters.
SymbolValueUnitInterpretationH1kg/a.u.Half-saturation constant of consumersT1Day−1Mass-specific metabolic rater1Day−1Resource turnover rateRmax2kg/a.u.Maximum resource densityϵ0.5—Efficiency of resource ingestionp0.955—Proportion of logistic vs semichemostatεd0.01—Deterministic error bound for resilienceεs0.1—Stochastic error bound for resiliencesbirth0.1cmSize at birthsmax10cmSize at maturationImax10Day−1Maximum ingestion rate per unit of biomassq0.85—Proportionality constant of ingestion ability between juveniles and adultsM00.1Day−1Natural mortality rateF—Day−1Harvesting rate of juveniles and adultsρ0.2—Volatility in resource growthR00.035—Threshold value of resource densityMVP0.009—Minimum viable populationWhen solving the stochastic models, we have simulated each model 10,000 times in order to find good estimates for the expectation and standard deviation, using the sample mean and the sample standard deviation.


### 4.1. Resilience

Resilience is increasingly used in ecology and fishery management context [6]. In simple terms, one can say that the higher the resilience value, the shorter it takes for a disturbance in the population to converge back toward its steady-state solution. We find the resilience of the population in the different resource growth models by using the mathematics in Section 3.3 together with the method and algorithms in Appendix B.1. The resource growth models investigated are semichemostat, logistic, semilogistic, and threshold-logistic resource dynamics (see Table 1), both in the deterministic and the stochastic settings.We observe that the resilience for both the deterministic models and the stochastic models (see Figures4 and 5) achieves its maximum around the same harvesting rates as the corresponding maximum values for the yield (see Figure 6). If this is not the case, one might want to aim for a slightly lower yield in order to increase the resilience; compare this with the concept of pretty good yield as defined in [6]. We notice that the resilience in the stochastic models vanishes at a lower harvesting rate than the deterministic equivalent models; this phenomenon is explained by the differences in the probability of extinctions.Figure 4
The resilience for deterministic (a, c) and stochastic (b, d) growth rates is calculated by using the algorithm described in AppendixB.1. The mean (yellow curve) and the standard deviation (68%) (light blue curves) of the resilience are shown for all models.
(a)(b)(c)(d)Figure 5
The resilience for all deterministic (a, c) and stochastic (b, d) growth rates is calculated by using the algorithm described in AppendixB.1. The mean (yellow curve) and the standard deviation (68%)(light blue curves) of the resilience are shown for all models.
(a)(b)(c)(d)Figure 6
The yield of semichemostat, logistic, semilogistic, and threshold-logistic growth rates at steady state w.r.t. harvesting rate. The grey trajectories show the yield from different simulations for the stochastic case. The yellow curve represents the yield mean value of the simulations. The light blue curves represent the interval of standard deviation (68%) of the yield. The black curve represents the yield for the deterministic case.


### 4.2. Recovery Potential

The recovery potential is the generational net biomass production (per unit of biomass) in a virgin environment. It approximates the net reproduction (expressed in biomass) of a small population introduced in an environment which is close to its maximum valueRmax; see Meng et al. [7]. We have named the recovery potential, expressed by equation (22), the Meng et al. recovery potential. As mentioned above, the Meng et al. recovery potential is evaluated in a virgin environment at equilibrium.Since our stochastic models never reach equilibrium, we reformulate the recovery potential in nonequilibrium terms through equations (28) and (30). We denote this formulation as nonequilibrium recovery potential.Figures2 and 3 show that two above formulations for the recovery potential in the deterministic case produce almost identical values. In view of this, we will use the nonequilibrium recovery potential when finding the probability of extinction (see Section 4.3).


### 4.3. The Probability of Extinction

Stability in ecological systems is important since a lack of stability promotes extinction. Under a random environment (in the deterministic setting, the probability of extinction can be evaluated using equation (31), which produces a graph with zero probability of extinction up to a certain point, and after this point, the probability of extinction is one), the lack of steady state can amplify the probability that the population goes extinct. We have two ways of finding the probability of extinction: we can either use the MVP formulation, equation (31), or we can use the method explained in Section 3.5, which we call RP formulation (here, RP stands for the recovery potential formulation used to find the probability of extinction). For MVP formulation, the number of minimum viable population of individuals is set to be 117; this is the average MVP value for fresh water fish species, given in [17].As can be seen in Figure7, the two formulations for the probability of extinction give similar results, the MVP formulation is presented in (a), (c), (e), and (g), and the RP formulation is presented in (b), (d), (f), and (h); each row corresponds to a specific growth rate model, defined in Table 1.Figure 7
Probability of extinction (PoE) with respect to the harvesting rate is calculated by using the algorithm described in AppendixB.3. The yellow curve (a, c, e, g) shows the mean of the MVP formulation of the probability of extinction, and the light blue curves represent the interval of three standard deviation (68%) of this probability of extinction. For the RP formulation of the probability of extinction (b, d, f, h), the mean is shown by the yellow curve and the three standard deviations by the light blue curves.
(a)(b)(c)(d)(e)(f)(g)(h)


## 4.1. Resilience

Resilience is increasingly used in ecology and fishery management context [6]. In simple terms, one can say that the higher the resilience value, the shorter it takes for a disturbance in the population to converge back toward its steady-state solution. We find the resilience of the population in the different resource growth models by using the mathematics in Section 3.3 together with the method and algorithms in Appendix B.1. The resource growth models investigated are semichemostat, logistic, semilogistic, and threshold-logistic resource dynamics (see Table 1), both in the deterministic and the stochastic settings.We observe that the resilience for both the deterministic models and the stochastic models (see Figures4 and 5) achieves its maximum around the same harvesting rates as the corresponding maximum values for the yield (see Figure 6). If this is not the case, one might want to aim for a slightly lower yield in order to increase the resilience; compare this with the concept of pretty good yield as defined in [6]. We notice that the resilience in the stochastic models vanishes at a lower harvesting rate than the deterministic equivalent models; this phenomenon is explained by the differences in the probability of extinctions.Figure 4
The resilience for deterministic (a, c) and stochastic (b, d) growth rates is calculated by using the algorithm described in AppendixB.1. The mean (yellow curve) and the standard deviation (68%) (light blue curves) of the resilience are shown for all models.
(a)(b)(c)(d)Figure 5
The resilience for all deterministic (a, c) and stochastic (b, d) growth rates is calculated by using the algorithm described in AppendixB.1. The mean (yellow curve) and the standard deviation (68%)(light blue curves) of the resilience are shown for all models.
(a)(b)(c)(d)Figure 6
The yield of semichemostat, logistic, semilogistic, and threshold-logistic growth rates at steady state w.r.t. harvesting rate. The grey trajectories show the yield from different simulations for the stochastic case. The yellow curve represents the yield mean value of the simulations. The light blue curves represent the interval of standard deviation (68%) of the yield. The black curve represents the yield for the deterministic case.


## 4.2. Recovery Potential

The recovery potential is the generational net biomass production (per unit of biomass) in a virgin environment. It approximates the net reproduction (expressed in biomass) of a small population introduced in an environment which is close to its maximum valueRmax; see Meng et al. [7]. We have named the recovery potential, expressed by equation (22), the Meng et al. recovery potential. As mentioned above, the Meng et al. recovery potential is evaluated in a virgin environment at equilibrium.Since our stochastic models never reach equilibrium, we reformulate the recovery potential in nonequilibrium terms through equations (28) and (30). We denote this formulation as nonequilibrium recovery potential.Figures2 and 3 show that two above formulations for the recovery potential in the deterministic case produce almost identical values. In view of this, we will use the nonequilibrium recovery potential when finding the probability of extinction (see Section 4.3).


## 4.3. The Probability of Extinction

Stability in ecological systems is important since a lack of stability promotes extinction. Under a random environment (in the deterministic setting, the probability of extinction can be evaluated using equation (31), which produces a graph with zero probability of extinction up to a certain point, and after this point, the probability of extinction is one), the lack of steady state can amplify the probability that the population goes extinct. We have two ways of finding the probability of extinction: we can either use the MVP formulation, equation (31), or we can use the method explained in Section 3.5, which we call RP formulation (here, RP stands for the recovery potential formulation used to find the probability of extinction). For MVP formulation, the number of minimum viable population of individuals is set to be 117; this is the average MVP value for fresh water fish species, given in [17].As can be seen in Figure7, the two formulations for the probability of extinction give similar results, the MVP formulation is presented in (a), (c), (e), and (g), and the RP formulation is presented in (b), (d), (f), and (h); each row corresponds to a specific growth rate model, defined in Table 1.Figure 7
Probability of extinction (PoE) with respect to the harvesting rate is calculated by using the algorithm described in AppendixB.3. The yellow curve (a, c, e, g) shows the mean of the MVP formulation of the probability of extinction, and the light blue curves represent the interval of three standard deviation (68%) of this probability of extinction. For the RP formulation of the probability of extinction (b, d, f, h), the mean is shown by the yellow curve and the three standard deviations by the light blue curves.
(a)(b)(c)(d)(e)(f)(g)(h)


## 5. Comparison Results with the Starvation Mortality Rate

In nature, animals that suffer starvation experience a higher mortality rate; see, for example, de Roos et al. [13] and Miriam and Ana [44]. Simulating our models with and without starvation, we investigate the differences in emergent properties. We include starvation in our models by changing from constant mortality rate to food-dependent mortality rate according to equations (7) and (8).In all stage-structured population resource growth rate models, the outcome of the impact on size structure is lowered when including starvation; this means that the proportion of the juvenile population increases (see Figures8 and 9). This change appears because the adult population experiences starvation at a higher resource density than the juveniles (see equations (7) and (8)), since we have set our parameter q to be less than one.Figure 8
Impact on size structure and biomass of threshold-logistic growth rate without starvation (a, c) and with starvation (b, d) at steady state w.r.t. harvesting rate. The grey trajectories show the impact on size structure and biomass from different simulations for the stochastic case. The yellow curve represents the mean value of the simulations. The light blue curves represent the interval of standard deviation (68%). The black curve represents the impact on size structure and biomass for the deterministic case. Equations (17) and (18) are used to derive the above graphs for the impact on size structure and biomass, respectively.
(a)(b)(c)(d)Figure 9
Impact on size structure and biomass of threshold-logistic growth rate without starvation (a, c) and with starvation (b, d) at steady state w.r.t. harvesting rate. The grey trajectories show the impact on size structure and biomass from different simulations for the stochastic case. The yellow curve represents the mean value of the simulations. The light blue curves represent the interval of standard deviation (68%). The black curve represents the impact on size structure and biomass for the deterministic case. Equations (17) and (18) are used to derive the above graphs for the impact on size structure and biomass, respectively.
(a)(b)(c)(d)Comparing the properties of impact on biomass and resilience shows a small variation when starvation is included or not in the models (see Figures10–12). The variations of the other properties in the simulations are minuscule between starvation and nonstarvation; i.e., steady state biomasses for resource, juvenile, and adult and yield and recovery potential as well as both formulations of the probability of extinction are not drastically changed by including starvation (see Appendix D for simulation results).Figure 10
Impact on size structure and biomass of semichemostat growth rate without starvation (a, c) and with starvation (b, d) at steady state w.r.t. harvesting rate. The grey trajectories show the impact on size structure and biomass from different simulations for the stochastic case. The yellow curve represents the mean value of the simulations. The light blue curves represent the interval of standard deviation (68%). The black curve represents the impact on size structure and biomass for the deterministic case. Equations (17) and (18) are used to derive the above graphs for the impact on size structure and biomass, respectively.
(a)(b)(c)(d)Figure 11
Impact on size structure and biomass of logistic growth rate without starvation (a, c) and with starvation (b, d) at steady state w.r.t. harvesting rate. The grey trajectories show the impact on size structure and biomass from different simulations for the stochastic case. The yellow curve represents the mean value of the simulations. The light blue curves represent the interval of standard deviation (68%). The black curve represents the impact on size structure and biomass for the deterministic case. Equations (17) and (18) are used to derive the above graphs for the impact on size structure and biomass, respectively.
(a)(b)(c)(d)Figure 12
The resilience for semichemostat, logistic, semilogistic, and threshold-logistic stochastic growth rate with nonstarvation (a, c, e, g) and starvation (b, d, f, h). The mean (yellow curve) and the standard deviation (68%) (light blue curves) of the resilience are shown for all models. The detail explanation for the simulation can be seen in Section4.1.
(a)(b)(c)(d)(e)(f)(g)(h)


## 6. Conclusions and Discussion

Models of an aquatic ecological system consisting of one fish species and its food resource were considered. Our models were two-stage-structured population models; we have investigated the impact of emergent properties of the population and resource in these models with respect to uniform harvesting rates. The models under consideration were characterized by two main features: (1) the growth rate dynamics choice for the resource and (2) whether the population experiences starvation mortality or not.The main purpose of this paper was to investigate properties of the population and the resource, both in deterministic growth rate models and the corresponding stochastic growth rate models. A major biological implication when including natural variations in the model is that the probability of extinction becomes a threat even at lower rates of harvesting, which indicates that when deterministic models are used for aquaculture, one should not harvest close to the extinction harvesting level.Many authors have investigated the behavior of stage-structured models; deterministic growth rate models have been studied by, e.g., Lundström et al. [6], Meng et al. [7], Ackleh and Jang [8], Aiello et al. [9], and de Roos et al. [13], and some stochastic growth rate models have been considered by, e.g., Rupšys et al. [45], Lv and Pitchford [46], Giet et al. [36], and Shah [34].Some of the population properties defined in the deterministic models was naturally transferable to the stochastic setting and evaluated by estimating the expectation and variance via Monte Carlo simulations. These properties included, e.g., biomasses, yield, impact on biomass, and impact on size structure. Among these properties, we found that all models under consideration shift their size structure toward juvenile individuals as harvesting is increased (see Figure13 in Appendix C). We would like to note that a shift in the size structured toward adult individuals, as harvesting was enhanced, may occur if the proportionality constant of ingestion is greater than one; see de Roos et al. [12]. This phenomenon may also occur if adult individuals cannibalize on the juveniles; see, e.g., Soudijn and de Roos [11] and Nankinga et al. [47].Figure 13
Impact on size structure of semichemostat, logistic, semilogistic, and threshold-logistic growth rate at steady state w.r.t. harvesting rate. The grey trajectories show the impact on size structure from different simulations for the stochastic case. The yellow curve represents the mean value of the simulations. The light blue curves represent the interval of standard deviation (68%). The black curve represents the impact on size structure for the deterministic case. Equation (17) is used to derive the above graphs.In contrast to the similar shifts in size structure, there is a big discrepancy between the models in that the total biomass in the logistic, semilogistic, and the threshold-logistic growth rate models is not much affected under moderate harvesting rates, whereas the total biomass in the semichemostat growth rate model declines rapidly with increasing harvesting rates (see Figure14 in Appendix C). In these three scenarios, we can in fact detect a small hydra effect; i.e., the biomass increases with increasing mortality, which is reflected by negative values of the impact on biomass at small harvesting rates in the logistic, semilogistic, and the threshold-logistic growth rate model (see Figure 14). As explained by Schröder et al. [30], the reason for this hydra effect is that moderate harvesting rates can reduce the strong competitions between individuals, leading to the higher biomass or less influence on the declination of the biomass. Since we could not detect the hydra effect in the semichemostat growth rate model, we deduced that the hydra effect depends on the rate of change for the resource growth at low levels (see Figure 1). Both the hydra effect and the discrepancy between the semichemostat and the other three models were explained by the difference of the rate of change of the resource production at low levels of the resource (see Figure 1).Figure 14
Impact on biomass of semichemostat, logistic, semilogistic, and threshold-logistic growth rate at steady state w.r.t. harvesting rate. The grey trajectories show the impact on biomass from different simulations for the stochastic case. The yellow curve represents the mean value of the simulation. The light blue curves represent the interval of standard deviation (68%). The black curve represents the impact on biomass for the deterministic case. The above graphs are calculated using equation (18).We have also compared the models with and without starvation mortality (we have chosen the same formula as [13]; see equations (7) and (8)). The main difference between including and excluding starvation mortality was manifested in the impact on biomass and impact on size structure as we discussed in Section 5. Furthermore, only minor effects were found in any of the resource growth rate models concerning the resource biomass, juvenile biomass, adult biomass, yield, resilience, recovery potential, and probability of extinction.Turning the attention to properties that cannot directly be translated from the deterministic case to the stochastic case. For example, the recovery potential defined in [7] is evaluated under the assumption that there is an equilibrium in the governing differential equations, and since this cannot be achieved with a stochastic growth rate, we have, in Section 3.4, derived an equivalent formulation in a nonequilibrium setting. Furthermore, new emergent properties become available when stochastic models are used; for example, the probability of extinction. In the literature (see, e.g., [17]), the probability of extinction can be evaluated by comparing the number of individuals of the population to the minimum viable population (MVP). We argued that the MVP is an inexact measure and hard to estimate in different environments. In Section 3.5, we proposed an alternative approach to evaluate the probability of extinction, called the RP formulation. In Section 4.3, the RP formulation was corroborated with the MVP formulation, resulting in similar shapes of the probability of extinction, but with a smoother shape in the RP formulation.The stochasticity that we introduced was focused on random growth rate for the resource. This is a natural assumption, since the resource depends on fluctuations in the environment. Another feature that cannot be naturally controlled is the harvesting. It is important to choose an optimal harvesting strategy since it affects the population survival probability, as well as the yield and other population properties. There is a variety of harvesting strategies and equipment used in the commercial fish farming industry. According to unpredicted conditions, such as weather, diseases, and climate change, the harvesting rates are stochastic in reality which is not studied in this paper. The randomness is introduced by Brownian motions; if one would like to study environmental factors including catastrophic events, other stochastic processes might be a more accurate choice. Furthermore, the simulations revealed the need of the two stages, juveniles and adults, of the fish when absolute harvesting was enforced; that is, one cannot reduce this system to a single state of biomass for the fish population, in particular, when probability of extinction or resilience is studied.The stage-structured population dynamics have been derived by bookkeeping properties for local averages of the population; in this paper, we have introduced stochasticity in the dynamics of the stage-structured population models. However, our next aim is to include the randomness in the basic assumptions of the individual state models. This will lead to a physiologically structured population model, consisting of a system of stochastic partial differential equations, from which we would like to derive the stochastic stage-structured models.



---

*Source: 3535375-2022-08-31.xml*