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.
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 probability . 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 Section 2, 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 at birth, and the maximum size of individuals is denoted by . Both stages forage for the shared resource, , and the metabolic requirement, , 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 and , 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, , and adults, , which forages on a shared resource . 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 (), adult biomass (), and food resource () which are
Here, and are the net biomass production per unit of body mass of juveniles and adults, respectively. Furthermore, and 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 assume uniform harvesting rates, i.e., equal proportional harvesting rates 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 . The maturation rate, , is the resource-dependent rate at which juveniles mature and become adults. The net biomass production rates for juveniles and adults are given by where where is the efficiency coefficient for the assimilated ingested resource. Here, is the half-saturation constant of consumers and is the maximum ingestion rate for juveniles. The mass-specific metabolic rate is denoted by . The ratio between the adult and juvenile feeding rates is denoted by . The factor 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 given in [12]. We set to reduce the number of parameters. In the stage-structured models, the juvenile maturation rate is given by where . In equations (1)–(6), the constant is measured in biomass per unit of volume, while the constants , , and as well as the mortality rates and are expressed per unit of time. All parameters as well as the biomass densities , , and 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]) where 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 . 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 either semichemostat 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.
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 a threshold-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 Figure 1, 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 where is the resource turnover rate and 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:
Note that in the threshold-logistic growth model, the function means that the resource available for consumers is the resource density above the threshold . When using the semilogistic growth rate model, we use the parameter to denote the proportion of the semichemostat growth and for the proportion of logistic growth.
2.2.2. Stochastic Setting
In this section, the resource dynamics described in Section 2.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 where 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
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 Table 1 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].
Semilogistic: the resource dynamics in this model is a linear combination of the stochastic semichemostat and the logistic growth rates. Given a constant , we define the semilogistic dynamics as where and are independent Brownian motions.
When consumers are introduced to these models, the rate of change of the available resource biomasses, , 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 by and , respectively, where and are the proportional harvesting rates. The yield objective function is defined by
In the absence of harvesting, and denote the juvenile and adult biomasses, respectively, at equilibrium.
In the stochastic setting, the yield is calculated, using equation (15), where and 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 where and . 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
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 where is the steady state of the expectation of the total biomass under a harvesting rate and .
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 . Let ; 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 .
We then find the return time as the time needed for this trajectory to be close enough to the equilibrium in the sense that for some small error bound in the deterministic case. To find the resilience, we use the definition introduced in [6]. That is, after repeating this procedure times, the resilience is defined by
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, , 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 in place of . 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:
In a virgin environment, Meng et al. [7] defined the recovery potential, as a function of harvesting rates, given by the expression
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 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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
We consider the following equations for the rates at which the biomass of juveniles and adults changes:
Rearranging equation (24) yields
Note that the right-hand side does not explicitly depend on . In what follows, our goal is to derive an expression for the recovery potential that does not include . When this has been achieved, we can use the expectation of the recovery potential to find the probability of extinction. When , we substitute equation (25) into equation (6) yielding
The right-hand side of this equation is an increasing continuous function with respect to (see Appendix A.1) and hence, we can find the unique solution of equation (26).
By substituting equation (25) and into equation (23), we get
Recall that, in this paper, we assume . By using equations (25) and (27) and the unique solution , we have derived an alternative expression for the recovery potential, defined by equation (22), as when .
In the case , equation (27) may be written as
By using equations (25) and (29) in equation (22), we get the recovery potential
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 ().
The 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, where is the survival probability for juveniles to become adults. The expression for is given by
See Appendix A.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 data , , and , where and 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 and . The simulations are done for a range of harvesting rates . 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, where is the mean value of the simulated recovery potentials for fixed harvesting rate . The standard deviation, , of the recovery potential is estimated by in a similar way.
In view of the central limit theorem, the mean of the recovery potential with harvesting rate is assumed to be a random sample from a normal distribution with cumulative distribution function, . The constants and represent mean and standard deviation, respectively. That is, we assume that and get the probability of extinction by
The reason for evaluating the cumulative distribution function at 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.
When 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 Figures 4 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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
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 value ; 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.
Figures 2 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 Figure 7, 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.
(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 Figures 8 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 to be less than one.
(a)
(b)
(c)
(d)
(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 Figures 10–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).
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(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 Figure 13 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].
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 Figure 14 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).
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 (). 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.
Appendix
A. Mathematical Derivations
A.1. Uniqueness Proof of the Solution
In Section 3.4 we stated that the right-hand side of equation (6) is an increasing continuous function with respect to . Since is a strictly increasing function when it is positive, it has an inverse function. Hence, for simplicity, we use the substitution and . Equation (6) then becomes
Our goal is to show that this function is injective. Since this function is continuous for all (this function is actually not defined at but the limit exists and equals zero) and differentiable for all such that , we study its derivatives. We find the derivative of to check if is an increasing function.
We see that the first factor in the above differential equation is positive; thus, it is enough to show that the second factor is positive; we rewrite
Then, is zero when (remember that is a constant); we are left to show that is positive when .
Setting gives
We must show that when and . The derivative of is
Observe that, since , we get . We now denote
It is now sufficient to show that when and when . This follows from the fact that the derivative of is
The calculations above in fact show that the function is bijective, and therefore, the solution is unique, which completes the proof.
A.2. Deriving
In this section, we derive an estimate for the probability of juveniles becoming adults, . The current amount of juveniles which do not reproduce obeys the equation for which we get the solution as
Let be the average time to become adult (assuming an equal spread of ages in the juvenile population), i.e., , where is the juvenile maturation rate. The total death rate for the juvenile population is .
The expected number of adults from current population is proportional to the biomass
When this population is close to extinction, i.e., close to zero, the resource is close to its maximum and since the population is small, we have zero starvation mortality; hence, we assume and . Thus, we get from equation (A.9)
We will be using this in the calculations for the MVP formulation of the probability of extinction, because we are only interested in knowing when the expected population goes extinct, and in this circumstance, the population will be close to zero.
B. Algorithms
Algorithms of our simulations are presented in this section, both in deterministic and stochastic settings.
B.1. Algorithm of Resilience
For the resilience, the initial values in each simulation are chosen at random for the resource, the juvenile, and the adult biomass in the cube. We then run the simulations, and for each trajectory, we determine the return time by using inequality (19). The resilience is measured by taking the reciprocal of the average value of the return times over a large number of simulations by the use of equation (20).
B.1.1. Algorithm of Resilience in Deterministic Setting
(i)The reciprocal of the mean value of return time is taken over the number of simulations with each random initial values for the deterministic case(ii)We find the return times when the solution trajectory is close enough to the equilibrium in the sense of Equation (19) by the mean value of 50 such simulations
B.1.2. Algorithm of Resilience in Stochastic Setting
(i)The reciprocal of the mean value of return times is taken over the number of simulations. The mean value of 20 simulations is used in equation (19) to find the value of return times by using the same random initial values(ii)We then find the average value of the return times after repeating the above procedure 20 times by using equation (20) when the solution trajectory is sufficiently close to the equilibrium in the sense of equation (19)(iii)Finally, we repeat the above procedure 15 times to find the mean and standard deviation of the resilience
B.2. Algorithm of Recovery Potential
For the recovery potential, we use the solutions of the juvenile, adult, and resource through the number of simulations. For the deterministic Meng et al. recovery potential, equation (22) is used. For the nonequilibrium recovery potential, equations (28) and (30) are utilized.
B.2.1. Algorithm of Nonequilibrium Recovery Potential in Deterministic Setting
(i)We use initial data , , and , where and are close to zero(ii)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 and (iii)Finally, equations (28) and (30) are used to find the nonequilibrium recovery potential
B.2.2. Algorithm of Nonequilibrium Recovery Potential in Stochastic Setting
(i)We use initial data , , and , where and are close to zero(ii)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 and (iii)We repeat this procedure 10,000 times to find the nonequilibrium recovery potential for each simulation by using equations (28) and (30)(iv)Finally, the mean and standard deviation of the nonequilibrium recovery potential are calculated
B.3. Algorithm of Extinction Probability
B.3.1. Algorithm of Extinction Probability in Stochastic Setting
The algorithm of extinction probability and the numerical values in the calculations are presented as follows: (i)Each model in systems (1)–(3) is run for 10,000 simulations(ii)For the MVP formulation of the probability of extinction, we have used , which we deduced from our simulations and the value of , which is based on the assumptions that (1) the simulated lake contains an expected value of 10,000 adult individuals at steady state without harvesting and (2) the steady state of the expected density (weight/volume) in our simulations(iii)For the RP formulation of the probability of extinction, we find the mean value of the nonequilibrium recovery potential for 20 simulations and call this stochastic variable . We then repeat this procedure 500 times and follow the explanations given in Section 3.5 to find the probability of extinction
C. Further Simulation Results
C.1. Stage-Structured Biomass Model with Resource Dynamics
For completeness, in this appendix, we study the steady-state biomasses, yield, impact on biomass, and size structure with starvation mortality for all growth rates.
C.2. Stage-Structured Biomass Dynamics and Yield
Our consumer-resource models are based on the models derived by de Roos et al. [13], which are reliable approximations of a fully size-structured population models. They showed that stage-structured population models formulated in this way incorporate key individual life-history processes. In this section, we investigate the model with starvation mortality for four types of resource dynamics which depend on the different harvesting rates for the stochastic case and include the deterministic model as a base case.
The stochastic models are the same as the deterministic models, with the exception that a random perturbation is added to the resource growth rate (see Section 2.2.2). We investigate the mean value and standard deviation of all our findings when the solution has reached a steady state according to equation (19).
We see that the resource density is increased with the harvesting rate in Figure 15 as the population becomes overexploited. In addition, when the biomasses of the juveniles (Figure 16) and the adults (Figure 17) are decreased, the population also goes extinct (Figure 7). We see that the resource reaches after all individuals die out, due to the harvesting strategies.
In Figure 16, for all growth models, the juvenile biomass first increases with harvesting rate as they can synthesize more proteins at metabolic costs close to the theoretical minimum compared to the adults, due to the difference in ingestion rates, which compensates the reduction of reproduction of offspring. At higher harvesting rates, the increment in food cannot compensate the loss of reproduction of offspring any longer and the juvenile biomass will therefore start to decrease. When the harvesting rate becomes too high, both the juvenile and the adult populations go extinct (see Figures 16 and 17).
In Figure 17, we see that the adult biomass decreases, for all models, as the harvesting rate increases.
Figure 6 represents how the yield changes with respect to harvesting rates. The yield for semichemostat first increases faster than the yield for logistic, semilogistic, and threshold-logistic as harvesting rate increases. The yield decreases as the population approaches the MVP; the yield for semichemostat, logistic, semilogistic, and threshold-logistic growth approaches zero as the harvesting rate becomes too large. The logistic, semilogistic, and threshold-logistic models reach the maximum yield closer to the point of extinction, than in the semichemostat model. In our models, the harvesting rate is deterministic, but in real life, this is not realistic and might also cause extinction (see Conclusions and Discussion for a deeper discussion).
C.3. Impact on Biomass and Size Structure
The population size structure is completely determined by the distribution of biomass between juveniles and adults, as compared to the distribution of the same stages at steady state without any harvesting. We investigate changes in population biomass and size structured in response to harvesting of juveniles and adults at equal rates. In the deterministic model by Meng et al. [7], harvesting juveniles and adults equally always leads to an increase in the percentage of juvenile biomass in the population, but in other deterministic stage-structured models, this percentage might decrease, e.g., in cannibalistic models. We follow these ideas and hence study the consequences of harvesting through the impact measures which influence biomass and size structure for deterministic and stochastic cases.
The mean value of the impact on size structure in stochastic simulations is always lower than the corresponding deterministic simulation (see Figure 13). We find that the mean value of the logistic model is relatively close to the deterministic solution, whereas in the other three growth rate models, the mean value of the stochastic case is significantly more than one standard deviation below the corresponding deterministic growth rate. The reason why the stochastic growth rate models on average has a lower impact on size structure compared to the deterministic models can be understood by a careful examination of how the net biomass production fluctuates under stochastic growth rates. The net biomass production rates, defined by equations (4) and (5), are strictly positive at steady state in the deterministic models; otherwise, we would not get any flux of biomass between the two stages; and therefore, the population would die out due to background mortality, hence not a steady state. In the stochastic setting, on the one hand, starvation will occur more often when the harvesting rate is zero, compared to higher harvesting rates. This is due to the fact that the population is larger, forcing the food demand to be higher, which causes the mean value to be higher than the corresponding deterministic value. On the other hand, when the harvesting rate is larger, thus, the population demands less of the food resource, and consequently, the stochasticity will not cause the values of the net biomass production to become negative and therefore not cause starvation. But when the resource is above the mean value, the net biomass production for the juveniles will increase more than the corresponding adult term (because we have chosen the proportionality constant of ingestion ability ), giving a higher rate of maturation of juveniles to adults, than the rate of reproduction. This difference shifts the mean value of the fraction , in equation (17) toward a lower value.
The impact on size structure in all our models is positive; the main reason for this is the choice of . A more careful exploration of the shape of impact on size structure can be found in de Roos and Persson [28].
Figure 14 shows that the impact on biomass, in the semichemostat case, increases almost linearly from zero up to one, whereas in the logistic, semilogistic, and threshold-logistic case, the impact on biomass resembles an exponential growth curve. The relative flatness in the logistic, semilogistic, and the threshold-logistic cases of the impact on biomass implies that the population size is hardly affected for low harvesting rates; this is further discussed in Conclusions and Discussion. Furthermore, when the harvesting rate is absent, the impact on biomass and size structure will not be affected. In Figures 13 and 14, this condition is not included when we simulate the impact on size structure and biomass.
It is also worth noting that, in contrast to the intuitive feeling that a population must decrease when harvesting is introduced, in the stochastic logistic, semilogistic, and threshold-logistic cases, the mean value of the impact on biomass may be negative. This implies that the total population biomass can increase under small harvesting rates in comparison with the steady-state solution without any harvesting. This phenomenon has also been observed by Schröder et al. [30] but does not occur in our deterministic model.
C.4. Resilience and Recovery Potential
The emergent properties in the stage-structured biomass models which are resilience and recovery potential for the semilogistic and threshold-logistic growth rates are presented in this section under the deterministic and stochastic settings.
In Figures 2 and 3, we estimate the recovery potential in stochastic settings, using the nonequilibrium recovery potential. The Meng et al. recovery potential in the deterministic equivalence is plotted in the same figure to show that the stochastic recovery potentials do not diverge much from the deterministic recovery potentials.
D. The Impact of Starvation in the Models
The simulations of the emergent properties for the stage-structured population models are investigated with and without starvation. The changes in the general outcome for the resilience and impact on size structure and biomass are presented when the model is simulated to compare the results with and without starvation.
Data Availability
We do not have any underlying data supporting the results of our study.
Disclosure
This manuscript is an extended version of the conference paper [48].
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by International Science Programme (ISP) in collaboration with South-East Asia Mathematical Network (SEAMaN). The authors thank Masood Aryapoor for his valuable feedback to finalize this paper. The authors are grateful to Lennart Persson for his insightful comments about growth rate models in our manuscript.