Next Article in Journal
An Efficient GNSS Coordinate Classification Strategy with an Adaptive KNN Algorithm for Epidemic Management
Previous Article in Journal
Integrated Profitability Evaluation for a Newsboy-Type Product in Own Brand Manufacturers
Previous Article in Special Issue
Mathematical Modeling of Two Interacting Populations’ Dynamics of Onchocerciasis Disease Spread with Nonlinear Incidence Functions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Discrete-Time Model of an Exploited Population with Age and Sex Structures: Instability and the Hydra Effect

1
Institute for Complex Analysis of Regional Problems, Far East Branch, Russian Academy of Sciences, Birobidzhan 679016, Russia
2
Institute of Automation and Control Processes, Far Eastern Branch, Russian Academy of Sciences, Vladivostok 690041, Russia
*
Author to whom correspondence should be addressed.
Mathematics 2024, 12(4), 535; https://doi.org/10.3390/math12040535
Submission received: 30 December 2023 / Revised: 4 February 2024 / Accepted: 6 February 2024 / Published: 8 February 2024

Abstract

:
This study proposes a discrete-time mathematical model to investigate the impact of selective harvesting on the dynamics of a population with age and sex structures. The model assumes that the birth rate depends on the sex ratio of the population and the number of breeding pairs. The growth rate is regulated by limiting juvenile survival, where an increase in population size decreases the survival of immature individuals. We consider the following selective proportional exploitation: harvesting of juveniles and harvesting of mature males. Depending on the values of population parameters, selective harvesting can lead to the stabilization of population dynamics by dampening oscillations or the emergence and amplification of fluctuations in population size. The model reveals multistability domains in which different dynamic modes coexist, and variations in initial conditions can lead to changes in dynamic modes. Depending on the values of the population parameters, the proposed models with harvest reveal the hydra effect, indicating an increase in the equilibrium abundance of the exploited group after reproduction but before harvesting, with an increase in the harvesting rate. Selective harvesting, resulting in the hydra effect, increases the remaining population size due to reproduction and the number of harvested individuals.

1. Introduction

For numerous animal species, commercial exploitation constitutes an inherent aspect of the biological cycle intertwined with processes of reproduction and survival. Consequently, harvesting, which substantially impacts ecosystems, has emerged as a fundamental determinant of population management.
Individuals of a specific age or sex in structured populations are not uncommon subjects of commercial harvesting. From a sustainable natural resource management point of view, it is economically or ecologically beneficial to harvest only adult individuals of certain animal species, while for others, it is more suitable to collect young individuals. For example, the commercial value of sturgeon, whitefish, and salmon varies depending on their age [1]. Note that the consequences of age-specific or stage-specific harvesting have been extensively investigated [2,3,4,5,6,7,8,9,10,11,12], while the impact of harvesting individuals of different sexes on population dynamics has been less studied [13,14,15,16,17,18,19].
Selective harvesting by sex is particularly relevant for species in which females and males differ in behavior, morphology, or economic value. For instance, the extraction of adult males is quite common in commercial harvesting or trophy hunting [20,21,22,23,24,25]. The removal of mature females is often prohibited or performed simultaneously with other groups. The increased mortality of females due to harvesting may be explained, for example, by their larger size compared to that of males, especially for certain species of fish [26]. Note that in some cases, sex-specific harvesting can disrupt sex ratios, which subsequently negatively impacts population dynamics [25,27,28,29]. Many species are able to reproduce effectively with a strongly biased sex ratio. Nevertheless, even polygynous species require a minimum number of males to maintain reproduction. Surveys of saiga antelope [27] and perch [30] have revealed population size decreases likely due to the hunting of males. In the context of applied research, dynamic models of populations structured by sex and age (stage) have been successfully employed to analyze and describe the dynamics of insects, reptiles, fish, birds, and mammals [15,19,25,28,31,32,33], including exploited species [15,19,25,31]. Considering both sexes in the model has provided more realistic and biologically meaningful predictions for different species than models of the population with only stage structure [19,32,33].
Overall, hunting based on differentiation can be used to control population size and structure, which can be beneficial in the development of strategies for maintaining sustainable population growth [34,35,36,37,38,39]. Here, the study of selective harvesting in populations structured by sex and stage often focuses on developing strategies for natural resource management, maintaining the reproductive core of the population, and measures to prevent overexploitation. Additionally, two-sex models serve as valuable tools for building biological control strategies [19,40]. For example, harvesting males may lead to positive outcomes accompanied by abundance growth in populations characterized by sexually antagonistic interactions, such as the costs of excessive mating or male harassment, until female fertility becomes limited by male availability [40].
One of the expected outcomes of harvesting is a decrease in the abundance of the species being exploited. This approach is often used to manage invasive and dangerous species to reduce their population numbers to safe levels or even eliminate them [41]. However, it is important to note that the decline in the population size of valuable species due to harvesting should not lead to their extinction. In recent years, a growing body of research has identified an interesting phenomenon known as the hydra effect [41,42,43,44,45,46,47]. This phenomenon describes a situation in which population abundance grows with an increasing harvest rate and is associated with the mythical hydra that regrew new heads in place of severed heads. An increase in population size in response to higher mortality, including harvesting, is also referred to as a paradoxical increase [48] or overcompensation [41]. One of the mechanisms leading to the appearance of the hydra effect in discrete-time population models is the overcompensatory density-dependent regulation of population growth [43]. As a result, harvest decreases the current size of the exploited population, which reduces intraspecific competition and slows the further decline in abundance. This effect has been observed in homogeneous [49] and structured population models [41,47,50].
Due to the difficulty of maintaining field censuses over a long enough period to determine an accurate average density, evidence of hydra effects in populations is scarce [43]. Several theoretical and empirical studies have shown that increased mortality caused by harvesting can lead to both fluctuations and growth in population size [51,52,53]. The occurrence of the hydra effect (overcompensation) and fluctuations in population due to harvesting have been empirically confirmed in species with high reproductive potential (fecundity), such as some plants [54], invertebrates [6,51,52,55], and fish [41]. These empirical studies [6,41,51,52,54,55] have considered cases where either adults only or juveniles and adults were harvested, which resulted in overcompensation and/or instability in the population.
An example of a population of terrestrial mammals with high fecundity rates is the wild boar. Over the past decades in Europe and North America, the population size of this species has increased despite hunting. Rapid maturation and high birth rates enable these animals to increase in number in response to management by hunting and control programs [56,57]. In some cases, selective harvesting of individuals of a specific size and/or sex has been shown to lead to changes in growth rates, structure, and sex ratios [19]. With intensive exploitation, females of the wild boar population begin to reproduce at an earlier age, which increases juvenile recruitment in response to harvesting [56].
For large mammals, such as ungulates, harvesting can have direct consequences on population size, age structure, and sex ratio [34,58,59]. Indirect effects arise because changes in the size and structure of the population impact demographic processes. Phenomena such as compensation, fluctuations, and changes in the lifecycle can eventually emerge in population dynamics due to harvesting [58,60,61].
Exploitation has become an important driver of trait change in the wild, inducing selective pressures that vary in strength and direction depending on the intensity, practices, and target phenotypes of harvest [62]. Harvest-induced selection on lifecycle, morphological, and behavioral traits has been found in both fishery and hunting systems [62].
Harvesting can influence population structure and induce changes in phenotypic traits and behavior, thus indirectly impacting the population growth rate and, consequently, the dynamics and sustainable development of the population [21]. In long-lived mammals, hunting-induced selection typically affects secondary sexual traits in males, such as the antlers and horns of mule deer, elk, white-tailed deer, and bighorn sheep [63,64], as well as body mass [65].
Theoretical and empirical investigations have revealed that an increase in mortality due to harvesting can result in an earlier age at maturation, a greater probability of reproduction, and increased litter size [56,66,67]. Thus, the population response to harvest may be complicated by demographic factors, specifically variations in vital population rates for different ages, body sizes, and sexes, as well as density-dependent processes [41].
The hydra effect often occurs in discrete-time models, considering harvesting before reproduction in the population [43]. However, harvesting after reproduction, which is more common in resource management practices [68], can lead to the hydra effect in a hidden form [45]. The primary difference between models that involve harvesting before and after reproduction is the timing of the population size survey. If the population size is measured immediately after reproduction, then these two scenarios are identical [44,45].
From a biological perspective, population dynamics are the result of the combination of reproduction and harvest, and any variation in population size depends only on the timing of harvest [69,70]. The order of these processes remains unchanged in both models because harvesting follows reproduction and precedes reproduction. Consequently, the stability conditions for these two harvesting scenarios are the same, as are the population sizes after each process [70]. Depending on the time of the survey, population sizes can vary quantitatively. With measurements taken after reproduction, higher densities are obtained than after harvest [45]. Note that the stability conditions of nontrivial fixed points are Identical for discrete-time models of stage-structured population dynamics under harvesting before and after reproduction [11].
We have previously proposed and explored a discrete-time model of the dynamics of age- and sex-structured populations with [18,71] and without harvest [72,73]. This paper extends the analysis of this mathematical model by investigating the impact of selective harvesting on the dynamics of the population, depending on the sex ratio. We suggest that utilizing demographic models with sex and stage structures is a good approach to designing management strategies aimed at population control. This method allows for the study of the response of the population to exploitation by analyzing various possible scenarios of its development, thus increasing the likelihood of selecting appropriate management measures.
To take into account the contributions of both sexes to reproduction, we utilize the concept that fertility depends on the ratio of mature males to females. This can be described using mating (pair formation) functions, which have been studied by various researchers [31,74,75,76]. The specific pair formation function used can vary depending on the species being studied. Theoretical modeling has revealed that factors such as the size of the harem, population density, male competition, and territorial distribution can all significantly impact the dynamics of a two-sex population [72,73,76,77].
Thus, our study focused on the dynamic modes of population with stage and sex structure that can occur depending on the intensity and selectivity of harvesting. This problem is meaningful because birth and mortality in natural populations are hardly variable, and both are determined by species biology. As a result, population size and population structure can be managed only by harvesting. Therefore, within the framework of this research, the harvest rate can and should be considered a control parameter due to the effect of harvesting on system dynamics. This study allowed us to compare and examine the possible changes in the dynamic modes of population size depending on the intensity and selectivity of harvesting at the same values as other population parameters. In addition, such research allows us to identify and analyze the occurrence of fluctuations, the hydra effect, and extinction in populations depending on the sex ratio, density-dependent regulation of juvenile survival, and demographic parameters.

2. The Population Model with Age and Sex Structures

We assume that a population with a seasonal breeding cycle can be represented by three groups: immature individuals (juveniles), mature females, and mature males, denoted as P (P ≥ 0), F (F ≥ 0), and M (M ≥ 0), respectively. The time step of the model is the interval between successive breeding seasons. We assume that the number of newborn individuals depends on the abundance of mature females and mature males in the population at the previous time step. Note that we segregate young individuals by sex when they have matured and replenish the adult part of the population.
A two-sex model in which fecundity depends on male and female abundances [74,76,77] can be written as follows:
( P n + 1 F n + 1 M n + 1 ) = ( 0 r F r M δ w 1 s 0 ( 1 δ ) w 2 0 v ) ( P n F n M n )
where n is the number of breeding seasons, rF (rF > 0) and rM (rM > 0) represent the fecundity values of females and males, respectively; δ (0 ≤ δ ≤ 1) is the proportion of newborn females; w1 (0 ≤ w1 < 1) and w2 (0 ≤ w2 < 1) are the survival rates of immature individuals; and s (0 ≤ s < 1) and v (0 ≤ v < 1) are the survival rates of mature males and females, respectively.
In a two-sex model, the fecundity values for male rM and female rF are derived from a total birth function B = B ( F n , M n ) [74,76]. The total number of offspring P in the year (n + 1) is the sum of individuals of both sexes, taking into account the fertility coefficients of the females (rF) and the males (rM):
P n + 1 = B ( F n , M n ) = r F F n + r M M n .
Let us assume that the primary sex ratio in zygotes after fertilization is 1:1, meaning that each zygote has one parent of each sex, such that r F F n = r M M n . Therefore, the coefficients rF and rM are r F = B ( F n , M n ) / ( 2 F n ) and r M = B ( F n , M n ) / ( 2 M n ) , respectively, where the factor of 1/2 prevents double counting of offspring from both males and females.
On the other hand, the number of newborns Pn is determined by the number of mating pairs C participating in reproduction and by the average number of offspring (a) per pair, i.e., P n + 1 = B ( F n , M n ) = a C ( F n , M n ) . We assume that the number of formed pairs C depends on the ratio of the numbers of females and males in the population, and can be described using the modified harmonic mean mating function [78]:
C ( F n ,   M n ) = min ( F n ,   2 F n M n / ( F n / h + M n ) ) ,
which determines the number of fertilized females, allows us to avoid overestimation of the birth rate for populations whose females produce offspring once per breeding season. The parameter h corresponds to the average size of the harem and characterizes the following mating relationships in the population: monogamy with h = 1, polygyny with h > 1, and polyandry at 0 < h < 1.
The condition for switching this function corresponds to the balance of sexes in the population and is expressed as Fn = hMn, where hMn is the number of females that can be fertilized by males whose number is M and the average harem size is h.
In other words, if the abundance of males is sufficient, that is, M n F n / h , then the number of pairs formed corresponding to the possible abundance of fertilized females will equal the number of mature females min ( F n ,   2 F n M n / ( F n / h + M n ) ) = F n . Consequently, the fertility values of female rF and male rM are as follows:
r F = B ( F n , M n ) 2 F n = a F n 2 F n = a 2 ,   r M = B ( F n , M n ) 2 M n = a F n 2 M n .
When M n < F n / h , which indicates a shortage of males in the breeding population, the number of pairs formed corresponding to the number of fertilized females is calculated as the harmonic mean of the number of males and females, taking into account harem size h, i.e., min ( F n ,   2 F n M n / ( F n / h + M n ) ) = 2 F n M n / ( F n / h + M n ) . In this case, the fertility values of female rF and male rM are as follows:
r F = B ( F n , M n ) 2 F n = a M n F n / h + M n ,   r M = B ( F n , M n ) 2 M n = a F n F n / h + M n .
Thus, the fertility functions of females and males can be rewritten as follows:
r F = a min ( 1 2 , M n F n / h + M n ) ,   r M = a min ( F n 2 M n , F n F n / h + M n ) .
Note that functional dependencies (3) were previously used in a two-sex matrix model to describe the fertility functions of the population of the California sea lion, Zalophus californianus [32].
We assume that density-dependent factors influence population dynamics and that the survival rates of immature females (w1) and males (w2) are the most sensitive parameters to population size. To describe the density-dependent regulation of juvenile survival, we use a discrete analog of the Verhulst equation, which takes into account self-limiting processes and competitive interactions among age classes: w 1 = 1 α 1 P β 1 F γ 1 M , w 2 = 1 α 2 P β 2 F γ 2 M . Here, αii ≥ 0), βii ≥ 0), and γii ≥ 0) are coefficients that characterize the intensities of declining juvenile survival due to competition between juveniles, adult females, and adult males, respectively. In this study, we assume that the impacts of adult females and males are the same, γi = βi and that the survival of juvenile females and males does not differ, that is, w1 = w2 = w; as a result, α1 = α2 = α and β1 = β2 = β. Therefore, survival rates linearly depend on the progeny and abundances of mature males and females: w 1 = w 2 = 1 α P β ( F + M ) . The coefficients α (α ≥ 0) and β (β ≥ 0) characterize the intensity of decreasing juvenile survival due to ecological limitations caused by competition for resources.
The linear dependence of young survival on the number of age-sex groups is convenient for analytical investigation, but it can lead to the loss of biological meaning in the model at large numbers, as in such cases, juvenile survival can be negative. Therefore, to maintain the meaning of the model, it is necessary to satisfy the condition 0 w 1 = w 2 = 1 α P β ( F + M ) < 1 . Our previous studies [72,73] discussed conditions for maintaining the biological meaning of model (6) in detail, depending on the parameters and initial conditions.
The matrix form of model (1) with the fecundity functions of the female rF and the male rM (3) and the survival of the immature females w1 and males w2 is as follows:
( P n + 1 F n + 1 M n + 1 ) = ( 0 a min ( 1 2 , M n F n / h + M n ) a min ( F n 2 M n , F n F n / h + M n ) δ ( 1 α P n β ( F n + M n ) ) s 0 ( 1 δ ) ( 1 α P n β ( F n + M n ) ) 0 v ) ( P n F n M n )
Thus, the number of individuals in each sex and age group of a population without harvesting can be described by a system of three recurrent equations proposed in [51] and has the following form in scalar form:
{ P n + 1 = a min ( F n ,   2 F n M n / ( F n / h + M n ) )   F n + 1 = δ ( 1 α P n β ( F n + M n ) ) P n + s F n M n + 1 = ( 1 δ ) ( 1 α P n β ( F n + M n ) ) P n + v M n ,
where n is the number of breeding seasons.
The following substitutions of p = αP, f = αF, and m = αM transform model (5) to a simpler form:
{ p n + 1 = a min ( f n ,   2 f n m n / ( f n / h + m n ) )   f n + 1 = δ ( 1 p n ρ ( f n + m n ) ) p n + s f n m n + 1 = ( 1 δ ) ( 1 p n ρ ( f n + m n ) ) p n + v m n .
Here, p, f, and m describe the relative abundances or densities of the corresponding age and sex groups, respectively; parameter ρ = β/α characterizes the relative contribution of mature individuals to the limitation of juvenile survival.
Stationary numbers corresponding to the coordinates of the fixed point of system (4) or (5) depend on the sex ratio. If fhm, then min ( f ,   2 f m / ( f / h + m ) ) = f and the coordinates of the non-zero fixed point are determined by the following formulas:
p ¯ ( 1 ) = ( 1 v ) ( s + a δ 1 ) ρ ( ( s v ) δ + 1 s ) + a δ ( 1 v ) ,   f ¯ ( 1 ) = 1 a p ¯ ( 1 ) ,   m ¯ ( 1 ) = ( 1 δ ) ( 1 s ) a δ ( 1 v ) p ¯ ( 1 ) ,
where a ,     δ 0 , v 1 . Solution (7) has biological meaning when the number of age classes is positive, that is, a > ( 1 s ) / δ .
If f > hm, then min ( f ,   2 f m / ( f / h + m ) ) = 2 f m / ( f / h + m ) and the coordinates of the non-zero fixed point are as follows:
p ¯ ( 2 ) = ( 1 s ) ( 1 v ) ( ( 1 δ ) h ( s + 2 a δ 1 ) δ ( 1 v ) ) ( ( s v ) δ + 1 s ) ( ( 1 s ) ( 1 δ ) h + ( 1 v ) δ ) ρ + 2 a h δ ( 1 δ ) ( 1 s ) ( 1 v ) , f ¯ ( 2 ) = ( 1 s ) ( 1 δ ) h + ( 1 v ) δ ( 1 s ) 2 a h ( 1 δ ) p ¯ ( 2 ) ,   m ¯ ( 2 ) = ( 1 s ) ( 1 δ ) h + ( 1 v ) δ ( 1 v ) 2 a h δ p ¯ ( 2 )
and are positive with a > ( δ ( 1 v ) + h ( 1 δ ) ( 1 s ) ) / ( 2 h δ ( 1 δ ) ) , where a ,   h ,   δ 0 , and s ,     v ,     δ 1 .
The boundaries of the stability area of fixed points in the system (6) are found based on the characteristic polynomial.
λ 3 + S λ 2 + H λ + D = 0 ,
where S, H, and D are three invariants of the Jacobian matrix of the system (6).
Let us rewrite Equation (9) in the form det ( J λ I ) = 0 , where J is the Jacobian matrix of system (6) evaluated at the fixed point ( p ¯ ( i ) , f ¯ ( i ) , m ¯ ( i ) ) at i = 1, 2, λ are the eigenvalues of matrix J, and I is the identity matrix. The matrix J consists of elements, which are the partial derivatives of the right-hand side of the system (6), and is as follows:
J = [ g p g f g m q p q f q m b p b f b m ] ( p ¯ ( i ) , f ¯ ( i ) , m ¯ ( i ) ) ,
where g = a min ( f ,   2 f m / ( f / h + m ) )   , q = δ ( 1 p ρ ( f + m ) ) p n + s f , and b   = ( 1 δ ) ( 1 p ρ ( f + m ) ) p + v m .
The coefficients of the polynomial Δ ( λ ) = | λ I J | = λ 3 + S λ 2 + H λ + D are expressed in terms of the invariants of the Jacobian matrix J and have the following meanings: S is the trace of the Jacobian matrix; H is the sum of the principal second-order minors of the Jacobian matrix J; and D is the determinant of the Jacobian matrix J. Therefore, S, H, and D are determined using the following formulas:
S = t r a c e ( J ) = g p q f b m , H = | g p g f q p q f | + | g p g m b p b m | + | q f q m b f b m | , D = det ( J ) .
Based on these coefficients, we can find stability area boundaries, each of which is a hypersurface that corresponds to specific bifurcations of codimension 1. The conditions for different bifurcation lines have the following form [79]:
(1)
The transcritical bifurcation line (TC) is H = D 1 S at λ = 1 ;
(2)
The period-doubling bifurcation line (PD) is H = D 1 + S at λ = 1 ;
(3)
The Neimark–Sacker bifurcation line (NS) is H = S D D 2 + 1 at λ = e ± i ϕ and | λ | = 1 .
The coefficients of the polynomial (9) for system (6) at fhm take the following values:
S ( 1 ) = ρ p ¯ ( 1 ) ( s + v ) , H ( 1 ) = δ a ( 2 p ¯ ( 1 ) + ρ f ¯ ( 1 ) + ρ m ¯ ( 1 ) 1 ) ρ p ¯ ( 1 ) ( ( 1 δ ) s + δ v ) + s v , D ( 1 ) = δ a v ( 2 p ¯ ( 1 ) + ρ f ¯ ( 1 ) + ρ m ¯ ( 1 ) 1 ) .
The coefficients of the polynomial (9) for system (6) at f > hm are as follows:
S ( 2 ) = ρ p ¯ ( 2 ) ( s + v ) , H ( 2 ) = 2 a h ( ( 1 δ ) ( f ¯ ( 2 ) ) 2 + δ h ( m ¯ ( 2 ) ) 2 ) ( 2 p ¯ ( 2 ) + ρ f ¯ ( 2 ) + ρ m ¯ ( 2 ) 1 ) ρ δ v p ¯ ( 2 ) ( f ¯ ( 2 ) + h m ¯ ( 2 ) ) 2 ( f ¯ + h m ¯ ) 2 + s ( v ρ ( 1 δ ) p ¯ ( 2 ) ) , D ( 2 ) = 2 a h ( ( 1 δ ) s ( f ¯ ( 2 ) ) 2 + δ h v ( m ¯ ( 2 ) ) 2 ) ( 2 p ¯ ( 2 ) + ρ f ¯ ( 2 ) + ρ m ¯ ( 2 ) 1 ) ( f ¯ ( 2 ) + h m ¯ ( 2 ) ) 2 .
The right-hand side of the system (6) is nonsmooth, which is the main feature of the model (6) that does not affect the stability analysis of its fixed points. The parameter space of model (6) is devided by hsw where h s w = δ ( 1 v ) / ( ( 1 s ) ( 1 δ ) ) and corresponds to Equality f = h m taking into account equilibrium numbers (7) and (8). Crossing value of hsw leads to a switching of function (2) corresponding to the number of pairs formed. As a result, the stability domain of system (6) in the ( ρ , a) plane is determined by bifurcation lines that correspond to either fixed point (7) or (8).
A study on the stability of system (6) was presented in [73], which showed that fixed points of system (6) lose stability via both the Neimark–Sacker scenario, leading to the emergence of quasiperiodics, and the period-doubling bifurcation, giving rise to regular oscillations in population size.

3. Selective Harvest of Immature Individuals

Let us consider the scenario of an undifferentiated harvest of juveniles in a population. Such selective harvesting is rational due to the economic value of the harvested species and/or the preservation of population sustainability. Examples of juvenile-specific exploitation due to its economic value include Greenland seals (Phoca groenlandica) [80], northern fur seals (Callorhinus ursinus L.) [25], and sables (Martes zibellina) [81]. Moreover, the harvesting of young wild boars [82], deer [83], and red deer (Cervus elaphus) [38] has been shown to maintain population sustainability. In general, the differentiated harvest of juveniles leads to the preservation of the reproductive core of the population, which consists of mature females and large breeding males [83].
We assume that a constant share of the juveniles is harvested after the breeding season. The number of captured individuals is proportional to the immature group size. Generally, model (6) with juvenile harvest can be expressed as follows:
{ p n + 1 = a min ( f n ,   2 f n m n / ( f n / h + m n ) )   ( 1 u p ) f n + 1 = δ ( 1 p n ρ ( f n + m n ) ) p n + s f n m n + 1 = ( 1 δ ) ( 1 p n ρ ( f n + m n ) ) p n + v m n .

3.1. Fixed Points of Model (10)

We distinguished the sizes of the population groups immediately after reproduction but before harvesting, and denoted them as p ˜ , f ˜ , and m ˜ , and the number of individuals remaining after harvesting, which are p, f, and m, respectively. The sizes of the juvenile group before and after harvesting are related by the equation p = p ˜ ( 1 u p ) , where the harvest rate up is a fraction of the captured juveniles. The sizes of the unexploited sex groups in the population before and after harvesting remain unchanged.
The stationary population sizes corresponding to the coordinates of the fixed point of the system (10) depend on the sex ratio. When the number of males in the stable population is sufficient, i.e., m ¯ f ¯ / h , and the number of pairs formed corresponds to the number of mature females, which is min ( f ,   2 f m / ( f / h + m ) ) = f , the stationary sizes of the stage and sex groups are determined by the following formulas:
The coordinates of the non-zero fixed point of the model (10) after harvesting are as follows:
p ¯ ( 1 ) = ( 1 v ) ( s + a δ ( 1 u p ) 1 ) ρ ( ( s v ) δ + 1 s ) + a δ ( 1 v ) ( 1 u p ) ,   f ¯ ( 1 ) = 1 a ( 1 u p ) p ¯ ( 1 ) ,   m ¯ ( 1 ) = ( 1 δ ) ( 1 s ) a δ ( 1 v ) ( 1 u p ) p ¯ ( 1 ) ;
The coordinates of the non-zero fixed point of the model (10) after reproduction but before harvesting are as follows:
p ˜ ¯ ( 1 ) = p ¯ ( 1 ) / ( 1 u p ) ,   f ˜ ¯ ( 1 ) = f ¯ ( 1 ) ,   m ˜ ¯ ( 1 ) = m ¯ ( 1 ) .
In the case of a shortage of males in the breeding population m ¯ < f ¯ / h , the number of pairs formed is calculated as the harmonic mean of the number of males and females, taking into account the harem size h, that is, min ( f ,   2 f m / ( f / h + m ) ) = 2 f m / ( f / h + m ) , and the coordinates of the non-zero fixed point of the model (10) are found as follows:
The solution after harvesting is as follows:
p ¯ ( 2 ) = ( 1 s ) ( 1 v ) ( ( 1 δ ) h ( s + 2 a δ ( 1 u p ) 1 ) δ ( 1 v ) ) ( ( s v ) δ + 1 s ) ( ( 1 s ) ( 1 δ ) h + ( 1 v ) δ ) ρ + 2 a h δ ( 1 δ ) ( 1 s ) ( 1 v ) ( 1 u p ) , f ¯ ( 2 ) = ( 1 s ) ( 1 δ ) h + ( 1 v ) δ ( 1 s ) 2 a h ( 1 δ ) ( 1 u p ) p ¯ ( 2 ) ,   m ¯ ( 2 ) = ( 1 s ) ( 1 δ ) h + ( 1 v ) δ ( 1 v ) 2 a h δ ( 1 u p ) p ¯ ( 2 ) ;
The solution after reproduction but before harvesting is as follows:
p ˜ ¯ ( 2 ) = p ¯ ( 2 ) / ( 1 u p ) , f ˜ ¯ ( 2 ) = f ¯ ( 2 ) ,   m ˜ ¯ ( 2 ) = m ¯ ( 2 ) .
Non-zero positive values of the equilibrium numbers of the system (10) are possible only if 0 u p < u p c r , where u p c r is the critical value and equals u p c r = 1 ( 1 s ) / a δ for solutions (11)–(12) or u p c r = 1 ( 1 s ) / 2 a δ ( 1 v ) / 2 a h ( 1 δ ) for solutions (13)–(14). If the harvest rate up is higher than u p c r , then the population will go extinct.
The equality f ¯ = h m ¯ corresponds to the balance of sexes in a stable population. Substituting the expressions for the equilibrium numbers of the model (10) into f ¯ = h m ¯ , we can obtain the threshold value of the parameter h s w , providing the balance of sexes:
h s w = δ ( 1 v ) / ( ( 1 s ) ( 1 δ ) ) ,
which coincides with the condition to switch the pair-formation function in the model (6) without harvest. Consequently, if the model parameters satisfy the inequality h h s w , then (i) the number of males for breeding is sufficient or exceeds the required number because m ¯ f ¯ / h ; (ii) the offspring abundance is determined only by the number of mature females, and (iii) the equilibrium numbers are determined by the expressions ( p ¯ ( 1 ) , f ¯ ( 1 ) , m ¯ ( 1 ) ) . However, with h < h s w , the population has a shortage of males because m ¯ < f ¯ / h , and some females cannot participate in reproduction; thus, the equilibrium numbers of population groups are ( p ¯ ( 2 ) , f ¯ ( 2 ) , m ¯ ( 2 ) ) .

3.2. Stability of Fixed Points of Model (10)

To find the boundaries of the stability area for fixed points of the system (10), we use the coefficients of the characteristic polynomial based on the methodology presented in Section 2. The coefficients of the characteristic polynomials (9) for solutions (11) and (12) coincide and are as follows:
D ( 1 ) = a δ v ( 1 u 1 ) ( 1 ρ ( f ¯ ( 1 ) + m ¯ ( 1 ) ) + 2 p ¯ ( 1 ) , H ( 1 ) = a δ ( 1 u 1 ) ( ρ ( f ¯ ( 1 ) + m ¯ ( 1 ) ) + 2 p ¯ ( 1 ) 1 ) + ρ p ¯ ( 1 ) ( δ s δ v s ) + s v ) , S ( 1 ) = p ¯ ( 1 ) ρ s v .
The coefficients for solutions (13) and (14) are as follows:
D ( 2 ) = 2 a h ( 1 u 1 ) ( 1 ρ ( f ¯ ( 2 ) + m ¯ ( 2 ) ) 2 p ¯ ( 2 ) ) ( δ h v ( m ¯ ( 2 ) ) 2 + s ( f ¯ ( 2 ) ) 2 ( 1 δ ) ) / ( f ¯ ( 2 ) + h m ¯ ( 2 ) ) 2 , H ( 2 ) = ( 2 a h ρ ( 1 δ ) ( 1 u 1 ) ( f ¯ ( 2 ) ) 3 + ( f ¯ ( 2 ) ) 2 ( 2 a h ( 1 δ ) ( 1 u 1 ) ( 2 p ¯ ( 2 ) + ρ m ¯ ( 2 ) 1 ) ρ p ¯ ( 2 ) ( s ( 1 δ ) + + δ v ) + s v ) 2 h f ¯ ( 2 ) m ¯ ( 2 ) ( ρ p ¯ ( 2 ) ( s ( 1 δ ) + δ v ) s v ) + 2 a δ h 2 ρ ( 1 u 1 ) ( m ¯ ( 2 ) ) 3 + + h 2 ( m ¯ ( 2 ) ) 2 ( 2 a δ ( 1 u 1 ) ( 2 p ¯ ( 2 ) + ρ f ¯ ( 2 ) 1 ) ρ p ¯ ( 2 ) ( s ( 1 δ ) + δ v ) + s v ) ) / ( f ¯ ( 2 ) + h m ¯ ( 2 ) ) 2 , S ( 2 ) = p ¯ ( 2 ) ρ s v .
In the parameter space of model (10), passing through hyperplane (15) leads to the switching of the pair-formation function. As a result, the bifurcation lines for solutions (11) or (13) bound the stability area of the system (10) (Figure 1a). Crossing the NS boundary, with an increasing birth rate, results in the quasiperiodic dynamics of the population and its stage groups. The transition through the PD line gives rise to stable oscillations due to the cascade of period-doubling bifurcations.
Analysis of the locations of the stability areas of solutions (11) and (13) for different sex ratios shows that at small ρ values corresponding to the low contribution of mature individuals in competition, the largest stability range with the numerical superiority of males m ¯ f ¯ / h is observed (Figure 1a). With increasing ρ , the stability area of the fixed point decreases for a population with a numerical superiority of males and expands for those with a shortage of males; that is, m ¯ < f ¯ / h .
As shown in Figure 1b,c, an increase in up with any sex ratio shifts the bifurcation lines NS and PD toward higher values of parameter a, which corresponds to an expansion of the stability area of the non-trivial fixed point of the system (10). This means that population dynamics stabilize with an increase in the harvest rate of juveniles.
Simultaneously, the non-trivial fixed point becomes unstable in accordance with the Neimark–Sacker scenario at low ρ values when crossing the NS line and in accordance with the period-doubling scenario at higher ρ values when crossing the PD line (Figure 1 and Figure 2).
An increase in the birth rate (a) leads to a narrowing of the stability region for the non-trivial fixed point at any sex ratio and causes its displacement toward higher harvest rates of juveniles up (Figure 2a,b).
An increase in the survival coefficient of females (s) or males (v) expands the stability area and shifts it toward higher values of intraspecific competition ( ρ ) (Figure 2c,d).
A higher harvest rate can result in a change in the dynamic modes of the population. Accordingly, if the estimated population parameter values are within the area of quasiperiodic or regular oscillations near the stability domain boundary, then the introduction of harvesting or an increase in intensity can dampen oscillations and stabilize population dynamics. Additionally, note that the bifurcation line PD has a convex form for certain parameter values (Figure 1 and Figure 2). Hence, an increase in the harvest rate of juveniles can lead to a shift in the observed dynamic mode, namely, periodic oscillations. Moreover, the increase in harvesting pressure can excite oscillations in a stable population.

3.3. Multistability of Dynamics

We investigate the dynamic modes of model (10) by means of dynamic mode maps. To generate these maps, we use a scanning method. The size of each image is 500 × 500 pixels. For example, in Figure 3a, for each pair of parameter values ρ and a at the nodes of a uniform grid that covers the area { ( ρ , a )   |   0 ρ 20 ,     0 a 10 } , we computed 10,000 iterations of the mapping. Subsequently, we analyzed the nature of the established dynamics after the transient process. We use the results of the last k = 500 steps to determine a cycle period T, which satisfies the inequality sup | x t x t + T | < 10 6 , where t = 1 , 2 , , k , T = 1 , 2 , , k / 2 . After detecting periodicity, the corresponding pixel on the diagram is colored according to the period obtained. If T > k / 2 , then we determine such dynamics to be irregular. Additionally, we check the condition 1 p n ρ ( f n + m n ) > 0 at each iteration step, as the survival function of juveniles w 1 = w 2 = 1 p n ρ ( f n + m n ) can become negative at large population sizes. If this condition is violated at a given step, we reset the survival function to zero at that step.
From a biological point of view, zero survival can be interpreted as complete offspring mortality for a given year due to high intraspecific competition for life resources among age classes caused by overpopulation. For instance, density-dependent regulation of juvenile survival, leading to extinction due to overcrowding, has been observed in the population of the flour beetle (Tribolium confusum), whose dynamics were studied in laboratory experiments [84]. The higher the initial density of the eggs, the fewer individuals survive to the adult stage. Furthermore, a situation may arise where the population becomes extinct because developing beetles consume all available food before any of them reach maturity [84].
Figure 3a shows the stability loss of a non-trivial fixed point either through a period-doubling bifurcation or via the Neimark–Sacker scenario. Moreover, there are wide areas of irregular dynamics because juvenile survival is a piecewise-linear function during numerical experiments, and newborn group size is determined by a minimum function with switching between different pair-formation functions depending on the sex ratio. For example, Figure 3a shows such a region of irregular dynamics within the 2-cycle area. Note that at the boundaries of this region, there are “islands” with a cascade of period-doubling bifurcations.
A numerical investigation of the model (10) also reveals multistability areas, where different dynamic modes coexist with the same values of model parameters. The map of dynamic modes in Figure 3 shows the multistability area, where a stable state coexists with a 3-cycle. This follows from the fact that the period-3 Arnold tongue overlies the stability domain of the non-trivial fixed point (Figure 3). Note that the 3-cycle arises due to a tangential bifurcation and coexists with a stable fixed point and dynamic modes emerging due to the stability loss of this fixed point.
The dynamic mode that will attract depends on the initial values of the age-sex group sizes of the population. Consequently, variations in initial conditions, for example, due to the influence of environmental factors or a changing harvest rate, can lead to either one dynamic mode or another, which will result in altering the nature of population dynamics.
The initial condition space can be divided by the attraction basins of several dynamic modes due to multistability. As shown in Figure 4a, different stable cycles can coexist. For example, a stable fixed point and a 3-periodic fixed point can exist together (Figure 4a), and cycles with 2 and 3 periods divide the initial condition space (Figure 4b). Generally, 3-cycle and dynamic modes, which arise from the loss of stability of the non-trivial fixed point (Figure 4a) or 2-cycle (Figure 4b), are attractive simultaneously.
As seen in Figure 4, a higher harvest rate up narrows the 3-cycle attraction basins, which is accompanied by the capture of the entire space of initial conditions by the stable fixed point. In addition, as the dynamic modes change, quasi-periodic oscillations transit to equilibrium (Figure 4a) or a 2-cycle (Figure 4b).
The cascade of period-doubling bifurcations, which is typical of Verhulst’s law, is not observed in Figure 4b; the 2-cycle loses stability via the Neimark–Sacker scenario. The bifurcation diagrams in Figure 5a,c allow us to study in more detail the transitions between dynamic modes by increasing the bifurcation parameter. As shown in Figure 5a, under certain initial conditions, a decrease in the harvest rate up leads to stability loss of the 2-cycle, resulting in quasiperiodic dynamics of the system (10). In the system (10) phase space, closed invariant curves emerge around each cycle element (Figure 5b). Under other initial conditions, a 3-cycle arises and bifurcates via a period-doubling scenario as the values of up decrease. Figure 5c demonstrates that the 3-cycle transits to a 6-cycle, which is then broken down with the emergence of six invariant curves.
Variations in population size and/or population parameters, caused by environmental and human-induced factors, can lead to a transition between basins of attraction. This transition in the dynamics of natural populations with age and sex structure is accompanied by a change in the period of oscillation or the emergence or disappearance of fluctuations. Additionally, a variation in the harvest rate can also cause a shift in the dynamic mode, such as swinging or dampening oscillations.

3.4. The Hydra Effect

The hydra effect refers to an increase in the equilibrium abundance of an exploited group of individuals in a population with a higher harvesting rate, i.e., mortality due to exploitation. Since this study considers harvesting after reproduction and the equilibrium population size is the number of individuals after harvesting, the hydra effect is concealed. To determine this effect, we analyze changes in population size postharvest and immediately after reproduction but before harvesting.
For system (10), the existence domain of the hydra effect is bounded by the surface H(k) of the maximum values of the equilibrium abundances of the exploited group of individuals in the population after reproduction. H(k) can be determined from the equation ( p ˜ ¯ ( k ) ) / u p = 0 . Consequently, the area where the hydra effect occurs is located within the stability domain of the equilibrium solutions of the model (10).
The value of the harvest rate at which the pre-harvest equilibrium abundance of juveniles reaches its maximum is determined from the equation ( p ˜ ¯ ( k ) ) / u p = 0 and is found as follows:
u H ( 1 ) = ( 1 v ) ( s + a δ 1 ) L a ( 1 v ) δ   at   m ¯ f ¯ / h ; u H ( 2 ) = ( 1 s ) ( 1 v ) ( ( 1 δ ) ( s + 2 a δ 1 ) h δ ( 1 v ) ) ( ( 1 δ ) ( 1 s ) h + δ ( 1 v ) ) L 2 a ( 1 v ) ( 1 s ) ( 1 δ ) h δ   at   m ¯ < f ¯ / h ,
where L = ( 1 s ) ( 1 v ) ( ρ ( s v ) δ + ( 1 s ) ( ρ v + 1 ) ) .
The maximums of the equilibrium numbers of mature females and males can be found in the formulas ( f ¯ ( k ) ) / u p = 0 and ( m ¯ ( k ) ) / u p = 0 and are both achieved when the harvest rate u H ( k ) is the same as the maximum pre-exploitation abundance of juveniles p ˜ ¯ H ( k ) . Consequently, an increase in the numbers of females and males results in a hydra effect on the pre-harvest abundance of juveniles.
An example of the existence area of the hydra effect is depicted in Figure 6a, which shows a decrease in the hydra domain with an increasing harvest rate of immature individuals in the parameter plane ( u p , ρ ) . Figure 5b,c illustrates the change in the sizes of the population groups with increasing up at different values of the parameter of intraspecific competition ρ .
When ρ > 4 a h δ ( 1 s ) ( 1 v ) ( 1 δ ) ( ( 1 δ ) ( s + a δ + 1 ) h + δ ( 1 v ) ) ( 1 δ s + s + δ v ) ( ( 1 δ ) ( 1 s a δ ) h + δ ( 1 v ) ) 2 = ρ ( 2 ) * ( ρ ( 2 ) * is the solution of equation u H ( 2 ) = 0 ), as seen in Figure 6a, where ρ = 10, an increase in the harvest rate of juveniles up results in a decrease in the numbers of all age and sex groups in the population (Figure 6b).
If the estimated population parameters are within the hydra effect area, such as ρ = 1.2 < ρ ( 2 ) * in Figure 6a, then the introduction of harvesting or a higher harvest intensity will increase the abundance of all population groups (Figure 6c). Figure 6c shows that with an increase in the harvest rate up to a level of u H ( 2 ) , the juvenile abundance after reproduction can reach p ˜ ¯ H ( k ) , while the post-harvest abundance of p ¯ decreases for any value of up. Therefore, if we take into account the offspring number when considering the equilibrium abundance of the exploited class after reproduction p ˜ ¯ , then the hydra effect occurs. Consequently, the hydra effect is hidden if the abundance of p ¯ is counted after harvesting.
The hydra effect can be observed in the domain of unstable solutions of the system (10). By calculating the average population size under exploitation, we can see this effect in the instability domain of the nontrivial solution (Figure 6d). As the harvesting rate up increases, the average group size of immature individuals after reproduction also increases. At u p = u * , the inverse Neimark–Sacker bifurcation occurs in system (10), and further growth of up stabilizes the population dynamics. As a result, the graphs of the equilibrium post-reproduction number of juveniles, the asymptotic abundance of juveniles, and the average value coincide and show that the hydra effect is observed up to the value u H ( 2 ) (Figure 6d).
Figure 7 illustrates the changes in the area of occurrence of the hydra effect at different sex ratios with increasing harvest rates of juveniles up (Figure 7a,b) and birth rate a (Figure 7c,d). As shown, in the parameter plane ( ρ , a ) , an increase in up at different sex ratios displaces the hydra effect area toward higher values of birth rate (Figure 7a,b). Similarly, in the parameter plane ( u p , ρ ) , a higher birth rate (a) results in a shift of the hydra domain toward higher values of up (Figure 7c,d).

4. Selective Harvest of Mature Males

The hunting of game species is usually non-selective when individuals of different sexes and ages are harvested. However, during a hunting season, only adult males can be allowed to yield, for example, the following species in Russia: noble deer, spotted deer, roe deer, European moose, Siberian moose, elk, and hind [24]. As a rule, hunting large mammals involves harvesting mature individuals; for instance, in the Selous Game Reserve (Tanzania, Africa), tourists are allowed to shoot a large number of game species [85]. Recently, trophy hunting has become increasingly popular, where animals with impressive antlers, tusks, or skulls are hunted to obtain derivatives. As a result, non-selective harvesting transforms into selective harvesting [24]. Preference is given to mature individuals with large body sizes and remarkable features of such “trophy” species as the Saiga antelope (Saiga tatarica), elephants, moose, and other animal species [15,65].
In the case of selective proportional harvesting of mature males after the breeding season, model (6) takes the following form:
{ p n + 1 = a min ( f n ,   2 f n m n / ( f n / h + m n ) )   f n + 1 = δ ( 1 p n ρ ( f n + m n ) ) p n + s f n m n + 1 = ( ( 1 δ ) ( 1 p n ρ ( f n + m n ) ) p n + v m n ) ( 1 u m ) ,
where the harvest rate um corresponds to the fraction of captured mature males.

4.1. Fixed Points of Model (16) and Their Stability

If m ¯ f ¯ / h , then min ( f ,   2 f m / ( f / h + m ) ) = f and the stationary population size is determined by the following formulas:
  • The coordinates of the non-zero fixed point of the model (16) after harvesting are as follows:
p ¯ ( 1 ) = ( 1 v ( 1 u m ) ) ( s + a δ 1 ) ( ( ( 1 s + v ) u m + s v ) δ + ( 1 s ) ( 1 u m ) ) ρ + a δ ( 1 v ( 1 u m ) ) , f ¯ ( 1 ) = p ¯ ( 1 ) / a ,   m ¯ ( 2 ) = ( 1 s ) ( 1 δ ) h + ( 1 v ) δ ( 1 v ) 2 a h δ ( 1 u p ) p ¯ ( 2 ) ;
  • The coordinates of non-zero fixed point after reproduction but before harvesting are as follows:
p ˜ ¯ ( 1 ) = p ¯ ( 1 ) , f ˜ ¯ ( 1 ) = f ¯ ( 1 ) ,   m ˜ ¯ ( 1 ) = m ¯ ( 1 ) / ( 1 u m ) ,
where a > ( 1 s ) / δ .
If m ¯ < f ¯ / h , then min ( f ,   2 f m / ( f / h + m ) ) = 2 f m / ( f / h + m ) and the solution after harvesting is as follows:
p ¯ ( 2 ) = ( 1 s ) ( 1 v ( 1 u m ) ) ( ( 1 δ ) ( 1 u m ) h ( s + 2 a δ 1 ) δ ( 1 v ( 1 u m ) ) ) h ( 1 δ ) ( 1 s ) ( 1 u m ) ( G ρ + 2 a δ ( 1 v ( 1 u m ) ) ) + ρ δ G ( 1 v ( 1 u m ) ) , f ¯ ( 2 ) = ( 1 s ) ( 1 δ ) ( 1 u m ) h + ( 1 v ( 1 u m ) ) δ 2 a h ( 1 s ) ( 1 δ ) ( 1 u m ) p ¯ ( 2 ) , m ¯ ( 2 ) = ( 1 s ) ( 1 δ ) ( 1 u m ) h + ( 1 v ( 1 u m ) ) δ 2 a h δ ( 1 v ( 1 u m ) ) p ¯ ( 2 ) ;
The solution after reproduction but before harvesting is as follows:
p ˜ ¯ ( 2 ) = p ¯ ( 2 ) , f ˜ ¯ ( 2 ) = f ¯ ( 2 ) ,   m ˜ ¯ ( 2 ) = m ¯ ( 2 ) / ( 1 u m ) .
Fixed points (17)–(18) are positive and exist with a > ( 1 s ) / δ . Solutions (19)–(20) exist at 0 u m < u m c r , where u m c r = 1 δ / ( δ v h ( δ 1 ) ( s + 2 a δ 1 ) is the critical value. If the harvest rate um is higher than u m c r , then the population will go extinct. The condition for switching the pair-formation function for model (16), determining the threshold value of the parameter u, is as follows:
u m s w = 1 δ / ( ( 1 δ ) ( 1 s ) h + δ v ) .
When u m < u m s w , the sex ratio considering harem size h is m ¯ f ¯ / h . This means that the abundance of offspring is determined by the number of mature females. In this case, Equations (17)–(18) can be used to find the equilibrium population size. If u m u m s w , then the number of mature females is greater than the maximum number of pairs: f ¯ > h m ¯ . In this case, Equations (19)–(20) can be used to determine the equilibrium population size. A high harvest of mature males can cause a shortage of males, meaning that some mature females will not be able to produce offspring.
Note that if an abundance of males is plentiful in an unexploited population, i.e., m ¯ f ¯ / h , and the sex ratio is skewed toward males, then the introduction of exploitation of males with values of harvest rate um higher than the threshold level u m s w will cause a shift in the sex ratio. If a population without harvesting has a sex ratio that is skewed toward mature females due to a shortage of males m ¯ < f ¯ / h , then an increase in the harvest rate of males um will not change the sex ratio.
The boundaries of the stability area of fixed points in the system (16) are found based on the coefficients of the characteristic polynomial (9). The coefficients for solutions (17) and (18) have the following form:
D ( 1 ) = a δ v ( 1 u 3 ) ( 1 ρ ( f ¯ ( 1 ) + m ¯ ( 1 ) ) 2 p ¯ ( 1 ) , H ( 1 ) = a δ ( ρ ( f ¯ ( 1 ) + m ¯ ( 1 ) ) 1 ) + p ¯ ( 1 ) ( 2 a δ ρ ( 1 u 3 ) ( s δ ( v + s ) ) + s v ( 1 u 3 ) , S ( 1 ) = ρ p ¯ ( 1 ) ( 1 u 3 ( 1 δ ) ) v ( 1 u 3 ) s .
The coefficients for solutions (19) and (20) are
D ( 2 ) = 2 a h ( 1 u 3 ) ( 1 ρ ( f ¯ ( 2 ) + m ¯ ( 2 ) ) 2 p ¯ ( 2 ) ) ( δ h v ( m ¯ ( 2 ) ) 2 + s ( f ¯ ( 2 ) ) 2 ( 1 δ ) ) / ( f ¯ ( 2 ) + h m ¯ ( 2 ) ) 2 , H ( 2 ) = ( 2 a h ρ ( 1 δ ) ( 1 u 3 ) ( f ¯ ( 2 ) ) 3 + ( 1 u 3 ) ( f ¯ ( 2 ) ) 2 ( 2 a h ( 1 δ ) ( 2 p ¯ ( 2 ) + ρ m ¯ ( 2 ) 1 ) ρ p ¯ ( 2 ) ( s ( 1 δ ) + δ v ) + s v ) + 2 h f ¯ ( 2 ) m ¯ ( 2 ) ( 1 u 3 ) ( s v ρ p ¯ ( 2 ) ( s ( 1 δ ) + δ v ) ) + 2 a δ h 2 ρ ( m ¯ ( 2 ) ) 3 + + h 2 ( m ¯ ( 2 ) ) 2 ( ( 1 u 3 ) ( ρ p ¯ ( 2 ) ( s ( 1 δ ) + δ v ) + s v ) + 2 a δ ( 2 p ¯ ( 2 ) + ρ f ¯ ( 2 ) 1 ) ) ) / ( f ¯ ( 2 ) + h m ¯ ( 2 ) ) 2 , S ( 2 ) = ρ p ¯ ( 2 ) ( 1 u 3 ( 1 δ ) ) v ( 1 u 3 ) s .
Let us consider the change in the stability domain of the system (16) depending on the value of the harvest rate um at m ¯ f ¯ / h , which corresponds to an unexploited population with a sufficient number of males (Figure 8a,c). With um > u m s w , an increase in the value of the coefficient um leads to a change in the sex ratio and, accordingly, the stability area of the model (16) (Figure 8b,d).
If the offspring number is only determined by the number of females, m ¯ f ¯ / h , then the introduction of male harvesting and an increase in its intensity contract the stability area in order to the birth rate (a); the Neimark–Sacker bifurcation lines (NS) shift toward smaller values of a (Figure 8a). Consequently, when an abundance of males is plentiful, a higher harvest rate can result in the emergence and amplification of quasiperiodic oscillations. With an increase in the value of the parameter um, crossing the threshold value defined by relation (21) means that the number of mature females in the population exceeds the number of forming pairs: f ¯ > h m ¯ (Figure 8b). Simultaneously, the stability region of the nontrivial fixed point of the system (16) expands for the same values of the birth rate a. Therefore, an increase in the intensity of male harvesting um in a population with a sex ratio skewed toward females can dampen quasiperiodic oscillations.
A variation in the parameter of intraspecific competition, ρ , can cause changes in the stability domain for any sex ratio in the population that are not necessarily consistent. When m ¯ f ¯ / h , there are population parameter values at which an increase in the harvest rate um can move the period-doubling bifurcation line PD to either lower values (Figure 8a) or higher values of ρ (Figure 8c). The same is true for the case of f ¯ > h m ¯ (Figure 8b,d). Therefore, an increase in male harvest intensity for any sex ratio can result in the emergence, amplification, or dampening of periodic oscillations, depending on the values of system (16) parameters.
Note that, with m ¯ f ¯ / h , an increase in um up to u m s w narrows the stability area concerning values of the parameter of intraspecific competition ( ρ ) (Figure 8a). If um > u m s w , then the sex ratio in the population is f ¯ > h m ¯ , and further growth of um expands the stability domain concerning ρ (Figure 8b). If initially, with m ¯ f ¯ / h , a higher um expands the stability area (Figure 8c), then crossing the threshold value and subsequent growth in um lead to the stability area construction concerning values of ρ (Figure 8d).

4.2. The Hydra Effect

The functional dependence determining the boundary of the hydra effect area for system (16) H(k) is found from the equation ( m ˜ ¯ ( k ) ) / u m = 0 .
With values of the system (16) parameters corresponding to a non-exploited population with a sufficient number of males, i.e., m ¯ f ¯ / h , the hydra effect area does not depend on the harvest rate um and is located between the lines: a 1 H ( 1 ) = ( 1 s ) / δ and a 2 H ( 1 ) = ρ ( ( s v 1 ) δ + 1 s ) / v δ , and a 2 H ( 1 ) > a 1 H ( 1 ) (Figure 9a–d). In Figure 9a–d, the graphs a = a 1 H ( 1 ) and a = a 2 H ( 1 ) correspond to the lines TC(1) and H(1), respectively. Note that the curves a 2 H ( 1 ) and a 1 H ( 1 ) intersect at ρ * = v ( 1 s ) / ( 1 s δ ( 1 s + v ) ) . At the same time, the maximum equilibrium numbers of juveniles and mature females, which are calculated using the formulas ( f ¯ ( 1 ) ) / u m = 0 and ( p ¯ ( 1 ) ) / u m = 0 , are achieved at values greater than a = ( 1 s ) / δ and coincide with the lower boundary of the hydra’s area. Note that the boundary a 1 H ( 1 ) = ( 1 s ) / δ is a condition for the existence of a nontrivial solutions (17) and (18) of the system (16). Therefore, if m ¯ f ¯ / h , then the equilibrium numbers of juveniles p ¯ ( 1 ) and mature females f ¯ ( 1 ) always increase with the growth of harvest rate of males um.
When m ¯ f ¯ / h , a decrease in the survival rate of mature males v (Figure 9b) and the proportion of newborn females δ (Figure 9d) leads to the expansion of the hydra effect area and a shift of this area toward lower values of the intraspecific competition parameter ρ . Conversely, an increase in the survival rate of mature females s (Figure 8c) contracts the hydra domain.
A further increase in the harvest rate um beyond the threshold level u m s w leads to a change in the sex ratio (Figure 9). Thus, for higher values of um, there exists a hydra effect area for a fixed point (20). If the population has a shortage of males m ¯ < f ¯ / h , then the area where the hydra effect occurs depends on the value of um. The boundary of the hydra effect area is defined by the surface H(2), which is determined by the equation ( m ˜ ¯ ( 2 ) ) / u m = 0 (Figure 9). At the same time, the maximum equilibrium population of mature females is found from ( f ¯ ( 2 ) ) / u m = 0 , which corresponds to the curve M f 2 in Figure 9. Note that the highest values of the equilibrium number of females occur at lower values of reproductive potential (a) compared to the maximum number of males. The local maximum of the equilibrium juvenile abundance, determined by ( p ¯ ( 2 ) ) / u m = 0 , is achieved at biologically meaningless values of the system (16) parameters. Thus, with m ¯ < f ¯ / h , the stationary juvenile number always decreases as the harvest rate of mature males um increases.
If the sex ratio is skewed toward mature females in the unexploited population: m ¯ < f ¯ / h , then an increase in the harvest rate of mature males um does not alter the sex ratio. Figure 10 shows the changes in the hydra effect area due to variations in the population parameters and an increase in the harvest rate um.
With m ¯ < f ¯ / h , an increase in the survival coefficients of females (s) and males (v), the proportion of newborn females (δ), and the male harvest rate (um) leads to a shift in the hydra area toward higher birth rates (a) and intraspecific competition ( ρ ) (Figure 10).
Additionally, Figure 11 shows the changes in the hydra effect area and the equilibrium numbers of age-sex groups in the population with increasing um for different sex ratios. The equilibrium number of immature individuals p ¯ increases with a higher harvest rate of males (um) if the population has a sufficient number of males, i.e., m ¯ f ¯ / h (Figure 11b–d). Conversely, with m ¯ < f ¯ / h , the number of juveniles decreases. A growth in the number of mature females can be observed with increasing um for any sex ratio in the population. The equilibrium post-harvest number of males   m ¯ always decreases. However, the equilibrium number of males after reproduction but before harvesting m ˜ ¯ can reveal the hydra effect.
With the values of the model parameters from hydra’s domain, if the sex ratio is skewed toward males m ¯ f ¯ / h , then an increase in the harvest rate of males (um) up to the threshold value u m s w leads to a growth of male numbers after reproduction m ˜ (Figure 11a).
In the case of a non-high birth rate (a), the pre-exploitation number of males m ˜ ¯ rises to the level of m ˜ ¯ H ( 1 ) , after which it begins to decline at um > u m s w and m ¯ < f ¯ / h . Such behavior is depicted in Figure 11b for a = 2.5. In this case, the threshold value u m s w coincides with the harvest rate u H ( 1 ) , and the maximum equilibrium pre-exploitation number is achieved at um = u m s w . Simultaneously, with an increase in the pre-exploitation male number, the numbers of juveniles and females also grow with increasing values of um up to u m s w .
With a higher birth rate a, the maximum equilibrium number of females is achieved at um > u m s w , corresponding to m ¯ < f ¯ / h , as shown in Figure 11c for a = 3.5. As before, the hydra effect is observed at the harvest rate value of um = u m s w .
A further increase in parameter a leads to the emergence of the hydra effect on the number of males for any sex ratio in the population; for example, a = 6.5 in Figure 11d. As can be seen, at m ¯ f ¯ / h , with a higher harvest rate of males, the pre-harvesting equilibrium number m ˜ increases slowly. After surpassing the threshold value u m s w , the sex ratio is skewed toward females: m ¯ < f ¯ / h . Moreover, the parameter values continue to fall within the hydra domain, and pre-harvest abundance m ˜ continues to increase with growth in um. When um surpasses u m = u H ( 2 ) , the equilibrium male number after reproduction begins to decrease. Therefore, the pre-harvest number of males increases to m ˜ ¯ H ( 1 ) at u m s w = u H ( 1 ) with m ¯ f ¯ / h and then reaches m ˜ ¯ H ( 2 ) at u m = u H ( 2 ) with m ¯ < f ¯ / h .

4.3. Instability and Multistability

A further growth of the reproductive potential (a) destabilizes the dynamics of the system (16), when m ¯ < f ¯ / h , for example, the case of a = 7.6 in Figure 11a. Here, the hydra effect occurs within a limited area with stable fixed points of the system (16). Figure 12a,b depict changes in equilibrium and average numbers of males after reproduction with a bifurcation diagram showing the asymptotic dynamic mode of male abundance for specific values of the parameters. As shown, in the instability area, the average number of males can both increase and decrease with the growing harvest rate um (Figure 12a,b). Consequently, in the domain of unstable solutions of the system (16), the hydra effect has a non-monotonic nature. In general, a growth of the average number of males with an increase in the harvest rate um allows us to conclude the existence of the hydra effect in the instability area of the system (16).
Note that different dynamic modes are observed under different initial conditions (Figure 12a,b), which indicates multistability in the system (16). Thus, one set of initial values gives equilibrium dynamics of population size (Figure 12a), while the other one results in the occurrence of a 3-cycle (Figure 12b). Figure 12c shows the basin of attraction of these dynamic modes, where points B1 and B2 are the initial conditions for bifurcation diagrams depicted in Figure 12a,b.
Note that in the model under consideration multistability arises due to the tangent bifurcation caused by the nonlinearity of the system (16) (Figure 12) and due to switching in the pair formation function. Note that switching the pair formation function significantly complicates the observed dynamics of the system (16). For example, in the switching region, the coexistence of three dynamic modes is possible: a fixed point, a 3-cycle, and irregular dynamics (Figure 12d). These fixed and periodic points occur when m ¯ f ¯ / h , that is, the number of formed pairs in the population corresponds to the number of females, and irregular dynamics emerge with m ¯ < f ¯ / h . Consequently, if the model parameter values are within the switching condition (21) area, then the existing dynamic modes overlap, and the observed dynamics depend on the current abundance of the population.
Figure 12e–g shows basins of attraction of coexisting dynamic modes at m ¯ < f ¯ / h , where irregular dynamics are possible. Note that with a decrease in the harvest rate um, the following series of bifurcations of fixed point with 3-cycle are observed. The fixed point (1-cycle) bifurcates via the Neimark–Sacker scenario, which leads to the emergence of irregular dynamics. The cycle with period 3 loses stability through the cascade of period-doubling bifurcations; thus, the 3-cycle is replaced by 6-cycle, 12-cycle, 26-cycle, etc. As a result, irregular dynamics (ID) can coexist with 26-cycle (Figure 12e) or 12-cycle (Figure 12f), also 1- and 3-cycles can simultaneously attract (Figure 12g).
The dynamic modes maps in Figure 13a,b illustrate the change in the dynamic modes in the parameter plane (um, a) due to the variation in the initial conditions, which correspond to points B1 and B2 in Figure 12c. Figure 13 shows that there exist multistability areas; for example, resonant cycles with different periods overlap, and 1- and 3-cycles coexist. In these domains, the initial conditions determine which dynamic mode will be attractive.
A further increase in the birth rate a destabilizes the dynamics of the system (16) for any sex ratio in the population, for example, a = 8.1 in Figure 11a and Figure 13c. Similar to the previous case, the hydra effect occurs within a limited region where the fixed points of the system (16) are stable. Figure 13c illustrates the dynamics of the male number after reproduction with an increasing harvest rate um. Initially, the pre-harvest number of males m ˜ is stable and grows slowly. Then, the non-trivial fixed point of the system (16) loses stability via the Neimark–Sacker bifurcation (Figure 11a), and irregular oscillations in population size occur, initially at m ¯ f ¯ / h and then at m ¯ < f ¯ / h (Figure 13a). Figure 13c also demonstrates the possibility of multistability within this range. A further increase in um stabilizes the dynamics of the system (16).

5. Discussion

The research conducted in this paper demonstrates complex relationships between the sex ratio, harvesting, and the nature of population dynamics.
The stability loss of fixed points of systems (10) and (16) is shown to occur via the Neimark–Sacker bifurcation, leading to the emergence of quasiperiodic dynamics and the period-doubling bifurcation, resulting in oscillations with finite periods. In most cases, an increase in the harvest rate stabilizes dynamics by dampening oscillations. However, shifts in dynamics due to changes in sex ratios and non-monotonic behavior of the stability domain boundaries of systems (10) and (16) with variations in the values of population parameters can lead to unexpected scenarios of population development.
Specifically, increasing the intensity of harvesting immature individuals up can both dampen and induce oscillations depending on the reproductive potential a (Figure 1 and Figure 2). In the case corresponding to a non-exploited population with a sufficient number of males, i.e., m ¯ f ¯ / h , the introduction of exploitation of males with harvest rate um can lead to the emergence and amplification of quasiperiodic oscillations (Figure 8a,c). Conversely, in populations with a shortage of males m ¯ < f ¯ / h , such selective harvesting dampens quasiperiodic oscillations (Figure 8b,d). However, it is important to note that depending on the value parameters of the system (16) for the population with any sex ratio, an increase in the harvest rate of males can result in the emergence and amplification or suppression of periodic oscillations (Figure 8). The possibility of oscillation emergence or stabilization with increasing harvesting intensity has been demonstrated previously in models of the dynamics of ecologically limited populations structured by stage or age, e.g. [11,50].
Moreover, models (10) and (16) of exploited population dynamics reveal the possibility of achieving various dynamic modes under different initial population sizes due to multistability, as in model (6) without harvesting [73]. Variations in initial conditions may lead to the realization of either one dynamic mode or another (Figure 3, Figure 4 and Figure 5, Figure 12 and Figure 13). Multistability arises due to the tangent bifurcation caused by the system’s nonlinearity and switching in the pair formation function (Figure 12 and Figure 13). Therefore, slight variations in the current population size, leading to changes in the sex ratio, can complicate population dynamics due to multistability. Simultaneously, challenges arise in identifying the type of observed dynamic mode, as self-regulation processes and the current sex ratio in the population impact the behavior of the trajectory.
The harvest of juveniles is shown to not significantly change the sex ratio of mature individuals in the population (Figure 6). For an unexploited population at any sex ratio with the offspring abundance determined by the number of females, the introduced harvest of juveniles does not alter the sex ratio but can decrease or increase the equilibrium numbers of females and males (Figure 6b and c, respectively). Evidently, the number of immature individuals decreases with increasing harvest rate. Moreover, under specific values of parameters, with an increase in harvest rate, the number of juveniles after harvest can increase by the next breeding season before harvest. Note that the growth of the juvenile number after reproduction and the numbers of mature females and males occur simultaneously; as a result, their highest numbers correspond to the same harvest rate. Consequently, the selective harvest of immature individuals can lead to the emergence of the hydra effect on the equilibrium number of juveniles in the population after reproduction but before harvesting (Figure 6c).
It is generally accepted that the occurrence of the hydra effect requires harvesting to precede reproduction [43,48]. This paper, as the study [45], demonstrates that the hydra effect can occur during harvesting after reproduction. When considering population size immediately after harvesting, the hydra effect is hidden.
At low reproductive potential, the hydra effect arises at low values of the intrapopulation competition coefficient and the harvest rate of juveniles (Figure 7). With an increase in population reproductive potential, the hydra effect area within the stability domain of equilibrium solutions expands with respect to intraspecific competition parameters and the harvest rate of immature individuals. A higher birth rate increases the intensity of density-dependent regulation and leads to the emergence of oscillations. Consequently, in this case, the hydra area for stable equilibrium numbers narrows and shifts toward higher harvest rates (Figure 7).
In contrast to the exploitation of juveniles, the harvest of mature males alters the sex ratio among breeding individuals. In the case of a non-exploited population with a sufficient number of mature males, which corresponds to a common situation in nature, the introduction of male harvesting and an increase in its intensity initially led to a balanced sex ratio, followed by a shortage of males in the population (Figure 8). However, if there is an unexploited population with a sex ratio skewed toward females, then harvesting males will lead to an even greater deficit of males. As a result, the abundance of all age-sex groups in the population can significantly decrease due to excessive harvesting. There are known cases of dramatic declines in reproduction and population collapse in various species due to the intensive hunting of males [27,30].
Similarly to the juvenile harvest, the exploitation of males demonstrates the hydra effect in the exploited group numbers after reproduction but before exploitation (Figure 9, Figure 10 and Figure 11). The equilibrium pre-harvest number of individuals in the exploited population group increases with a higher harvest rate due to overcompensatory density-dependent regulation of juvenile survival. The individuals removed can be considered to be in excess abundance, the harvest of which decreases competitive pressure and leads to growth in population size. This effect is observed in the population with any sex ratio. In the case of mature male exploitation, the hydra effect occurs at relatively high birth rates, even if the survival rates of females and males are relatively low. Moreover, with a higher harvesting rate, the population size increases with a higher birth rate. Conversely, in a population with a sufficient number of males, an increase in the harvesting mortality of mature males can result in the growth of the pre-harvesting number at relatively low birth rates (Figure 9).
Moreover, the numerical investigation of the model proposed in this study revealed the non-monotonic nature of the hydra effect in the instability domain (Figure 12a,b and Figure 13c). Note that the increase in the equilibrium number of a stable population due to harvesting is driven by overcompensatory density regulation [43,48]. However, in cases of unstable dynamics, the hydra effect can occur due to changes in the amplitude and nature of population fluctuations [43]. A higher harvesting mortality often dampens oscillations, which leads to a higher average population size. However, the analysis carried out in this study shows that an increase in the harvesting rate does not always suppress oscillations, so it is necessary to consider the type and amplitude of fluctuations.
The existence of the hydra effect and its influence on the older age class of the population have previously been investigated for overcompensatory models with proportional exploitation [41,50] and threshold exploitation [47]. Such studies have been developed since adult individual harvest is the most promising from a fishing point of view [47]. Note that the existence of excess abundance in fish populations and its harvesting have been actively discussed [43].
The emergence of the hydra effect during differentiated harvesting can be used for the sustainable management of renewable resources [45]. Specifically, the hydra effect leads to an increase in the numbers of both the exploited group of the population accounting for offspring after reproduction and harvested individuals. However, under specific parameter values, the value of the harvest rate, corresponding to the maximum equilibrium number of the exploited group, is near a critical harvest rate, leading to population extinction. Consequently, the organization of harvesting requires a balanced approach, especially when population parameters are unknown.
The removal of individuals can be used to combat biological invasions or pests. In this case, management measures can have the opposite effect when the pest number increases due to the hydra effect [41,45,55]. For example, in a coastal marine ecosystem, an eradication program resulted in stage-specific overcompensation and a 30-fold, single-year increase in the population of the invasive European green crab (Carcinus maenas) [55]. One hypothesis for this overpopulation is that the removal of adults reduces the cannibalism of young individuals and leads to an overcompensatory replenishment in the population of this species. Cannibalism of young individuals by adult individuals is crucial for density-dependent regulation of population replenishment [55].

6. Conclusions

Thus, depending on the values of population parameters, the selective harvest of individuals from a specific age-sex group (juveniles or mature males) can result in either the damping of fluctuations and stabilization of population dynamics or the emergence of oscillations in abundance. The proposed models reveal multistability domains in which different dynamic modes coexist. In these domains, the initial conditions determine which dynamic mode will be achieved. Notably, multistability arises due to bifurcations caused by the nonlinearity of the systems and due to switching in the pair formation function. Therefore, slight variations in the current population size, leading to changes in the sex ratio, complicate population behavior and can result in a shift in the observed dynamic mode. Depending on the values of the population parameters, the harvesting scenarios considered in this paper reveal the hydra effects, which correspond to an increase in the equilibrium abundance of the harvested age-sex group after reproduction with an increase in the harvest rate. Selective harvesting, causing the hydra effect, leads to an increase in the abundance of the remaining population after reproduction but before harvest and to a growth in the number of harvested individuals.

Author Contributions

Conceptualization, E.F. and O.R.; investigation, O.R.; writing—original draft preparation, O.R. and G.N.; writing—review and editing, G.N. and E.F. All authors have read and agreed to the published version of the manuscript.

Funding

This work was carried out within the framework of the state targets of the Institute for Complex Analysis of Regional Problems of the Far Eastern Branch of the Russian Academy of Sciences.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Frisman, E.Y.; Last, E.V. Nonlinear effects on population dynamics related to age structure and fishery impact. Biol. Bull. 2005, 32, 425–437. [Google Scholar] [CrossRef]
  2. Skaletskaya, E.I.; Frisman, E.Y.; Shapiro, A.P. Discrete Models of Population Dynamics and Harvest Optimization; Nauka: Moscow, Russia, 1979. (In Russian) [Google Scholar]
  3. Abakumov, A.I. Management and Optimization in Models of Harvested Populations; Dal’nauka: Vladivostok, Russia, 1993. (In Russian) [Google Scholar]
  4. Tyutyunov, Y.; Senina, I.; Jost, C.; Arditi, R. Risk assessment of the harvested pike-perch population of the Azov sea. Ecol. Model. 2002, 149, 297–311. [Google Scholar] [CrossRef]
  5. Freckleton, R.P.; Silva Matos, D.M.; Bovi, M.L.A.; Watkinson, A.R. Predicting the impacts of harvesting using structured population models: The importance of density-dependence and timing of harvest for a tropical palm tree. J. Appl. Ecol. 2003, 40, 846–858. [Google Scholar] [CrossRef]
  6. Cameron, T.C.; Benton, T.G. Stage-structured harvesting and its effects: An empirical investigation using soil mites. J. Anim. Ecol. 2004, 73, 996–1006. [Google Scholar] [CrossRef]
  7. Hunter, C.M.; Caswell, H. Selective harvest of sooty shearwater chicks: Effects on population dynamics and sustainability. J. Anim. Ecol. 2005, 74, 589–600. [Google Scholar] [CrossRef]
  8. Wikström, A.; Ripa, J.; Jonzén, N. The role of harvesting in age-structured populations: Disentangling dynamic and age truncation effects. Theor. Popul. Biol. 2012, 82, 348–354. [Google Scholar] [CrossRef] [PubMed]
  9. Wallace, K.; Leslie, A.; Coulson, T. Re-evaluating the effect of harvesting regimes on Nile crocodiles using an integral projection model. J. Anim. Ecol. 2013, 82, 155–165. [Google Scholar] [CrossRef] [PubMed]
  10. Abakumov, A.I.; Izrailsky, Y.G. The harvesting effect on a fish population. Math. Biol. Bioinform. 2016, 11, 191–204. [Google Scholar] [CrossRef]
  11. Neverova, G.P.; Abakumov, A.I.; Yarovenko, I.P.; Frisman, E.Y. Mode change in the dynamics of exploited limited population with age structure. Nonlinear Dyn. 2018, 94, 827–844. [Google Scholar] [CrossRef]
  12. Tyutyunov, Y.V.; Titova, L.I.; Senina, I.N.; Dashkevich, L.V. Quasi-extinction risk assessment practices for harvested fish species based on a long-term forecast modelling of population dynamics. Proc. South. Res. Cent. Russ. Acad. Sci. 2020, 8, 181–198. (In Russian) [Google Scholar] [CrossRef]
  13. Tyutyunov, Y.; Arditi, R.; Büttiker, B.; Dombrovsky, Y.; Staub, E. Modelling fluctuations and optimal harvesting in perch populations. Ecol. Model. 1993, 69, 19–42. [Google Scholar] [CrossRef]
  14. Lindström, J. Harvesting and sex differences in demography. Wildl. Biol. 1998, 4, 213–221. [Google Scholar] [CrossRef]
  15. Snyder, K.T.; Freidenfelds, N.A.; Miller, T.E. Consequences of sex-selective harvesting and harvest refuges in experimental metapopulations. Oikos 2014, 123, 309–314. [Google Scholar] [CrossRef]
  16. Tahvonen, O.; Kumpula, J.; Pekkarinen, A.-J. Optimal harvesting of an age-structured, two-sex herbivore–plant system. Ecol. Model. 2014, 272, 348–361. [Google Scholar] [CrossRef]
  17. Shyu, E.; Caswell, H. Mating, births, and transitions: A flexible two-sex matrix model for evolutionary demography. Popul. Ecol. 2018, 60, 21–36. [Google Scholar] [CrossRef] [PubMed]
  18. Revutskaya, O.L.; Neverova, G.P.; Frisman, E.Y. Influence of Harvest on the Dynamics of Populations with Age and Sex Structures. Math. Biol. Bioinform. 2018, 13, 270–289. [Google Scholar] [CrossRef]
  19. Stubberud, M.W.; Vindenes, Y.; Vøllestad, L.A.; Winfield, I.J.; Stenseth, N.C.; Langangen, Ø. Effects of size- and sex-selective harvesting: An integral projection model approach. Ecol. Evol. 2019, 9, 12556–12570. [Google Scholar] [CrossRef]
  20. Ginsberg, J.R.; Milner-Gulland, E.J. Sex-biased harvesting and population dynamics in ungulates: Implications for conservation and sustainable use. Conserv. Biol. 1994, 8, 157–166. [Google Scholar] [CrossRef]
  21. Whitman, K.; Starfield, A.M.; Quadling, H.S.; Packer, C. Sustainable trophy hunting of African lions. Nature 2004, 428, 175–178. [Google Scholar] [CrossRef] [PubMed]
  22. Nilsen, E.B.; Solberg, E.J. Patterns of hunting mortality in Norwegian moose (Alces alces) populations. Eur. J. Wildl. Res. 2006, 52, 153–163. [Google Scholar] [CrossRef]
  23. Milner, J.M.; Nilsen, E.B.; Andreassen, H.P. Demographic side effects of selective hunting in ungulates and carnivores. Conserv. Biol. 2007, 21, 36–47. [Google Scholar] [CrossRef]
  24. Danilkin, A.A. Biological Basis of Hunting Trophy Business; Association of Scientific Publications KMK: Moscow, Russia, 2010; 150p. (In Russian) [Google Scholar]
  25. Zhdanova, O.; Kuzin, A.; Frisman, E. The Harvest Effect on Dynamics of Northern Fur Seal Population: Mathematical Modeling and Data Analysis Results. Mathematics 2022, 10, 3067. [Google Scholar] [CrossRef]
  26. Gerasimov, Y.V. (Ed.) Fishes of the Rybinsk Reservoir: Population Dynamics and Ecology; I.D. Papanin Institute for Biology of Inland Waters RAS, Filigran: Yaroslavl, Russia, 2015. (In Russian) [Google Scholar]
  27. Milner-Gulland, E.J.; Bukreeva, O.M.; Coulson, T.; Lushchekina, A.A.; Kholodova, M.V.; Bekenov, A.B.; Grachev, I.A. Reproductive collapse in saiga antelope harems. Nature 2003, 422, 135. [Google Scholar] [CrossRef] [PubMed]
  28. Miller, T.E.X.; Inouye, B.D. Confronting two-sex demographic models with data. Ecology 2011, 92, 2141–2151. [Google Scholar] [CrossRef]
  29. Miller, T.E.X.; Inouye, B.D. Sex and stochasticity affect range expansion of experimental invasions. Ecol. Lett. 2013, 16, 354–361. [Google Scholar] [CrossRef]
  30. Langangen, Ø.; Edeline, E.; Ohlberger, J.; Winfield, I.J.; Fletcher, J.M.; James, J.B.; Stenseth, N.C.; Vøllestad, L.A. Six decades of pike and perch population dynamics in Windermere. Fish. Res. 2011, 109, 131–139. [Google Scholar] [CrossRef]
  31. Frisman, E.Y.; Skaletskaya, E.I.; Kuzyn, A.E. A mathematical model of the population dynamics of a local northern fur-seal herd. Ecol. Model. 1982, 16, 151–172. [Google Scholar] [CrossRef]
  32. Gerber, L.R.; White, E.R. Two-sex matrix models in assessing population viability: When do male dynamics matter? J. Appl. Ecol. 2014, 51, 270–278. [Google Scholar] [CrossRef]
  33. Eberhart-Phillips, L.J.; Küpper, C.; Miller, T.E.X.; Cruz-López, M.; Maher, K.H.; dos Remedios, N.; Stoffel, M.A.; Hoffman, J.I.; Kruge, O.; Székely, T. Sex-specific early survival drives adult sex ratio bias in snowy plovers and impacts mating system and population growth. Proc. Natl. Acad. Sci. USA 2017, 56, 201620043. [Google Scholar] [CrossRef]
  34. Getz, W.M.; Haight, R.G. Population Harvesting: Demographic Models of Fish, Forest, and Animal Resources. In Monographs in Population Biology; Princeton University Press: Princeton, NJ, USA, 1989; Volume 27. [Google Scholar]
  35. Clark, C.W. Mathematical Bio-Economics: The Optimal Management of Renewable Resources, 2nd ed.; Wiley: New York, NY, USA, 1990. [Google Scholar]
  36. Lande, R.; Engen, S.; Saether, B.E. Optimal harvesting, economic discounting and extinction risk in fluctuating populations. Nature 1994, 372, 88. [Google Scholar] [CrossRef]
  37. Fryxell, J.M.; Packer, C.; McCann, K.; Solberg, E.J.; Sæther, B.E. Resource management cycles and the sustainability of harvested wildlife populations. Science 2010, 328, 903–906. [Google Scholar] [CrossRef]
  38. Milner, J.M.; Bonenfant, C.; Mysterud, A. Hunting Bambi—Evaluating the basis for selective harvesting of juveniles. Eur. J. Wildl. Res. 2011, 57, 565–574. [Google Scholar] [CrossRef]
  39. Abakumov, A.; Izrailsky, Y. Optimal Harvest Problem for Fish Population—Structural Stabilization. Mathematics 2022, 10, 986. [Google Scholar] [CrossRef]
  40. Rankin, D.J.; Kokko, H. Do males matter? The role of males in population dynamics. Oikos 2007, 116, 335–348. [Google Scholar] [CrossRef]
  41. Zipkin, E.F.; Kraft, C.E.; Cooch, E.G.; Sullivan, P.J. When can efforts to control nuisance and invasive species backfire? Ecol. Appl. 2009, 19, 1585–1595. [Google Scholar] [CrossRef] [PubMed]
  42. Abrams, P.A.; Matsuda, H. The effect of adaptive change in prey on the dynamics of an exploited predator population. Can. J. Fish. Aquat. Sci. 2005, 62, 758–766. [Google Scholar] [CrossRef]
  43. Abrams, P.A. When does greater mortality increase population size? The long history and diverse mechanisms underlying the hydra effect. Ecol. Lett. 2009, 12, 462–474. [Google Scholar] [CrossRef] [PubMed]
  44. Liz, E.; Ruiz-Herrera, A. The hydra effect, bubbles, and chaos in a simple discrete population model with constant effort harvesting. J. Math. Biol. 2012, 65, 997–1016. [Google Scholar] [CrossRef] [PubMed]
  45. Hilker, F.M.; Liz, E. Harvesting, census timing and “hidden” hydra effects. Ecol. Complex. 2013, 14, 95–107. [Google Scholar] [CrossRef]
  46. Schröder, A.; van Leeuwen, A.; Cameron, T.C. When less is more: Positive population level effects of mortality. Trends Ecol. Evol. 2014, 29, 614–624. [Google Scholar] [CrossRef]
  47. Liz, E.; Sovrano, E. Stability, bifurcations and hydra effects in a stage-structured population model with threshold harvesting. Commun. Nonlinear Sci. Numer. Simul. 2022, 109, 106280. [Google Scholar] [CrossRef]
  48. Seno, H. A paradox in discrete single species population dynamics with harvesting/thinning. Math. Biosci. 2008, 214, 63–69. [Google Scholar] [CrossRef] [PubMed]
  49. Ricker, W. Stock and recruitment. J. Fish. Board Can. 1954, 11, 559–623. [Google Scholar] [CrossRef]
  50. Liz, E.; Pilarczyk, P. Global dynamics in a stage-structured discrete-time population model with harvesting. J. Theor. Biol. 2012, 297, 148–165. [Google Scholar] [CrossRef]
  51. Costantino, R.F.; Cushing, J.M.; Dennis, B.; Desharnais, R.A. Experimentally-induced transitions in the dynamic behavior of insect populations. Nature 1995, 375, 227–230. [Google Scholar] [CrossRef]
  52. Costantino, R.F.; Desharnais, R.A.; Cushing, J.M.; Dennis, B. Chaotic dynamics in an insect population. Science 1997, 275, 389–391. [Google Scholar] [CrossRef]
  53. Cushing, J.M.; Costantino, R.F.; Dennis, B.; Desharnais, R.A.; Henson, S.M. Nonlinear population dynamics: Models, experiments and data. J. Theor. Biol. 1998, 194, jt980736. [Google Scholar] [CrossRef]
  54. Buckley, Y.M.; Hinz, H.L.; Matthies, D.; Rees, M. Interactions between density-dependent processes population dynamics and control of an invasive plant species, Tripleurospermum perforatum (scentless chamomile). Ecol. Lett. 2001, 4, 551–558. [Google Scholar] [CrossRef]
  55. Grosholz, E.; Ashton, G.; Bradley, M.; Brown, C.; Ceballos-Osuna, L.; Chang, A.; de Rivera, C.; Gonzalez, J.; Heineke, M.; Marraffini, M.; et al. Stage-specific overcompensation, the hydra effect, and the failure to eradicate an invasive predator. Proc. Natl. Acad. Sci. USA 2021, 118, e2003955118. [Google Scholar] [CrossRef] [PubMed]
  56. Servanty, S.; Gaillard, J.-M.; Ronchi, F.; Focardi, S.; Baubet, E.; Gimenez, O. Influence of harvesting pressure on demographic tactics: Implications for wildlife management. J. Appl. Ecol. 2011, 48, 835–843. [Google Scholar] [CrossRef]
  57. Treichler, J.W.; VerCauteren, K.C.; Taylor, C.R.; Beasley, J.C. Changes in wild pig (Sus scrofa) relative abundance, crop damage, and environmental impacts in response to control efforts. Pest Manag. Sci. 2023, 79, 4765–4773. [Google Scholar] [CrossRef] [PubMed]
  58. Solberg, E.J.; Saether, B.-E.; Strand, O.; Loison, A. Dynamics of a harvested moose population in a variable environment. J. Anim. Ecol. 1999, 68, 186–204. [Google Scholar] [CrossRef]
  59. Strand, O.; Nilsen, E.B.; Solberg, E.J.; Linnell, J.C.D. Can management regulate the population size of wild reindeer (Rangifer tarandus) through harvest? Can. J. Zool. 2012, 90, 163–171. [Google Scholar] [CrossRef]
  60. Skogland, T. Density-dependence in a fluctuating wild reindeer herd; maternal vs. offspring effects. Oecologia 1989, 84, 442–450. [Google Scholar] [CrossRef]
  61. Langvatn, R.; Loison, A. Consequences of harvesting on age structure, sex ratio and population dynamics of red deer Cervus elaphus in central Norway. Wildl. Biol. 1999, 5, 213–223. [Google Scholar] [CrossRef]
  62. de Walle, J.V.; Pelletier, F.; Zedrosser, A.; Swenson, J.E.; Jenouvrier, S.; Bischof, R. The interplay between hunting rate, hunting selectivity, and reproductive strategies shapes population dynamics of a large carnivore. Evol. Appl. 2021, 14, 2414–2432. [Google Scholar] [CrossRef]
  63. Strickland, B.K.; Demarais, S.; Castle, L.E.; Lipe, J.W.; Lunceford, W.H.; Jacobson, H.A.; Frels, D.; Miller, K.V. Effects of selective-harvest strategies on white-tailed deer antler size. Wildl. Soc. Bull. 2001, 29, 509–520. [Google Scholar]
  64. Pigeon, G.; Festa-Bianchet, M.; Coltman, D.W.; Pelletier, F. Intense selective hunting leads to artificial evolution in horn size. Evol. Appl. 2016, 9, 521–530. [Google Scholar] [CrossRef]
  65. Tenhumberg, B.; Tyre, A.J.; Pople, A.R.; Possingham, H.P. Do harvest refuges buffer kangaroos against evolutionary responses to selective harvesting? Ecology 2004, 85, 2003–2017. [Google Scholar] [CrossRef]
  66. Olsen, E.M.; Heino, M.; Lilly, G.R.; Morgan, M.J.J.; Brattey, J.; Ernande, B.; Dieckmann, U. Maturation trends indicative of rapid evolution preceded the collapse of northern cod. Nature 2004, 428, 932–935. [Google Scholar] [CrossRef]
  67. Proaktor, G.; Coulson, T.; Milner-Gulland, E.J. Evolutionary responses to harvesting in ungulates. J. Anim. Ecol. 2007, 76, 669–678. [Google Scholar] [CrossRef]
  68. Kokko, H. Optimal and suboptimal use of compensatory responses to harvesting: Timing of hunting as an example. Wildl. Biol. 2001, 7, 141–150. [Google Scholar] [CrossRef]
  69. Cid, B.; Hilker, F.M.; Liz, E. Harvest timing and its population dynamic consequences in a discrete single-species model. Math. Biosci. 2014, 248, 78–87. [Google Scholar] [CrossRef]
  70. Åström, M.; Lundberg, P.; Lundberg, S. Population dynamics with sequential density-dependencies. Oikos 1996, 75, 174–181. [Google Scholar] [CrossRef]
  71. Revutskaya, O.L.; Frisman, E.Y. Harvesting impact on population dynamics with age and sex structure: Optimal harvesting and the hydra effect. Comput. Res. Model. 2022, 14, 1107–1130. [Google Scholar] [CrossRef]
  72. Revutskaya, O.; Neverova, G.; Frisman, E. Complex Dynamic Modes in a Two-Sex Age-Structured Population Model. In Models of the Ecological Hierarchy: From Molecules to the Ecosphere; Jordán, F., Jørgensen, S.E., Eds.; Elsevier: Amsterdam, The Netherlands, 2012; pp. 149–162. [Google Scholar]
  73. Revutskaya, O.L.; Kulakov, M.P.; Neverova, G.P.; Frisman, E.Y. Changing of the dynamics modes in populations with age and sex structure. Dokl. Biol. Sci. 2017, 477, 239–243. [Google Scholar] [CrossRef]
  74. Caswell, H.; Weeks, D.E. Two-sex models: Chaos, extinction, and other dynamic consequences of sex. Am. Nat. 1986, 128, 707–735. [Google Scholar] [CrossRef]
  75. Hadeler, K.P.; Waldstätter, R.; Wörz-Busekros, A. Models for Pair Formation in Bisexual Populations. J. Math. Biol. 1988, 26, 635–649. [Google Scholar] [CrossRef] [PubMed]
  76. Caswell, H. Matrix Population Models: Construction, Analysis, and Interpretation, 2nd ed.; Sinauer Associates Inc.: Sunderland, UK, 2001; 722p. [Google Scholar]
  77. Lindström, J.; Kokko, H. Sexual reproduction and population dynamics: The role of polygyny and demographic sex differences. Proc. R. Soc. Lond. B 1998, 265, 483–488. [Google Scholar] [CrossRef] [PubMed]
  78. Bessa-Gomes, C.; Legendre, S.; Clobert, J. Discrete two-sex models of population dynamics: On modeling the mating function. Acta Oecologica 2010, 36, 439–445. [Google Scholar] [CrossRef]
  79. Kuznetsov, A.P.; Sedova, J.V. Bifurcations of three- and four-dimensional maps: Universal properties. Izv. VUZ. Appl. Nonlinear Dyn. 2012, 20, 26–43. (In Russian) [Google Scholar] [CrossRef]
  80. Conrad, J.; Bjørndal, T. A Bioeconomic Model of the Harp Seal in the Northwest Atlantic. Land Econ. 1991, 67, 158–171. [Google Scholar] [CrossRef]
  81. Monakhov, V.G. Hunting Selectivity and Its Influence on the Structure of Sable Populations in the Cis-Ural Region. Russ. J. Ecol. 2018, 49, 434–441. [Google Scholar] [CrossRef]
  82. Asinovsky, A.I. Principles of wild boar shooting as a means of preserving the species and evaluating the trophy. Bull. Assoc. Rosokhotrybolovsoyuz 2011, 2, 80–100. (In Russian) [Google Scholar]
  83. Danilkin, A.A. Deer (Cervidae). In Mammals of Russia and Adjacent Regions; GEOS: Moscow, Russia, 1999; 552p. (In Russian) [Google Scholar]
  84. Bigon, M.; Zarper, D.; Townsend, K. Ecology. Individuals, Populations, and Communities; Mir: Moscow, Russia, 1989; Volume 1, 667p. (In Russian) [Google Scholar]
  85. Caro, T.M.; Young, C.R.; Cauldwell, A.E.; Brown, D.D.E. Animal breeding systems and big game hunting: Models and application. Biol. Conserv. 2009, 142, 909–929. [Google Scholar] [CrossRef]
Figure 1. The stability areas of the system (10) on the plane ( ρ , a) with various values of other parameters. (a) The dotted lines show the boundaries of the stability domain of solution (11); the solid lines correspond to those for solution (13). Changes in the stability regions of solutions (13) (b) and (11) (c) at different values of harvest rate up.
Figure 1. The stability areas of the system (10) on the plane ( ρ , a) with various values of other parameters. (a) The dotted lines show the boundaries of the stability domain of solution (11); the solid lines correspond to those for solution (13). Changes in the stability regions of solutions (13) (b) and (11) (c) at different values of harvest rate up.
Mathematics 12 00535 g001
Figure 2. Stability regions of solutions (13) (a,c) and (11) (b,d) of the system (10) on the plane (up, ρ ) with s = 0.8, v = 0.5, δ = 0.5 (a,b), a = 5, h = 2, δ = 0.5 (c,d). TC, NS, and PD are the lines of transcritical, Neimark–Sacker, and period-doubling bifurcations, respectively.
Figure 2. Stability regions of solutions (13) (a,c) and (11) (b,d) of the system (10) on the plane (up, ρ ) with s = 0.8, v = 0.5, δ = 0.5 (a,b), a = 5, h = 2, δ = 0.5 (c,d). TC, NS, and PD are the lines of transcritical, Neimark–Sacker, and period-doubling bifurcations, respectively.
Mathematics 12 00535 g002
Figure 3. Dynamic mode maps of the model (10). The figures correspond to the period of observed cycles. ID stands for irregular dynamics. l0 corresponds to population extinction. The white lines TC, NS, and PD are the transcritical, Neimark–Sacker, and period-doubling bifurcations, respectively.
Figure 3. Dynamic mode maps of the model (10). The figures correspond to the period of observed cycles. ID stands for irregular dynamics. l0 corresponds to population extinction. The white lines TC, NS, and PD are the transcritical, Neimark–Sacker, and period-doubling bifurcations, respectively.
Mathematics 12 00535 g003
Figure 4. Basins of attraction for the coexisting dynamic modes of the model (10) at s = 0.8, v = 0.5, δ = 0.5, and h = 2 with variation in the harvest rate up value. The figures correspond to the periods of the observed cycles. Q stands for quasi-periodic dynamics.
Figure 4. Basins of attraction for the coexisting dynamic modes of the model (10) at s = 0.8, v = 0.5, δ = 0.5, and h = 2 with variation in the harvest rate up value. The figures correspond to the periods of the observed cycles. Q stands for quasi-periodic dynamics.
Mathematics 12 00535 g004
Figure 5. Bifurcation diagrams of the dynamic variable p with respect to the parameter up at various values of the initial conditions (a,c). (b) Phase portrait of system (10) with s = 0.8, v = 0.5, δ = 0.5, h = 2, a = 6.75, and ρ = 10.
Figure 5. Bifurcation diagrams of the dynamic variable p with respect to the parameter up at various values of the initial conditions (a,c). (b) Phase portrait of system (10) with s = 0.8, v = 0.5, δ = 0.5, h = 2, a = 6.75, and ρ = 10.
Mathematics 12 00535 g005
Figure 6. (a) In the plane (up, ρ ), the stability area of the fixed point of system (10) (shading) and the existence domain of the hydra effect (hatching). TC, NS, and PD are the lines of transcritical, Neimark–Sacker, and period-doubling bifurcations, respectively. H is the line of maximum values of the pre-harvest abundance of the population. (b,c) Dependences of the equilibrium numbers of the system (10) on the harvest rate of juveniles up. The curves p ˜ ¯ , p ¯ ,   f ¯   , and   m ¯ correspond to the numbers of immature individuals before and after harvesting, mature females, and mature males, respectively. (d) Bifurcation diagrams of the dynamic variable p ˜ n corresponding to the number of juveniles after reproduction (blue dots). The green line is the equilibrium number of juveniles after reproduction p ˜ ¯ . The red line corresponds to the average number of juveniles. u p c r is the critical value, if the harvest rate up is higher than u p c r , then the population will go extinct. u * is the bifurcation value at which the inverse Neimark–Sacker bifurcation occurs.
Figure 6. (a) In the plane (up, ρ ), the stability area of the fixed point of system (10) (shading) and the existence domain of the hydra effect (hatching). TC, NS, and PD are the lines of transcritical, Neimark–Sacker, and period-doubling bifurcations, respectively. H is the line of maximum values of the pre-harvest abundance of the population. (b,c) Dependences of the equilibrium numbers of the system (10) on the harvest rate of juveniles up. The curves p ˜ ¯ , p ¯ ,   f ¯   , and   m ¯ correspond to the numbers of immature individuals before and after harvesting, mature females, and mature males, respectively. (d) Bifurcation diagrams of the dynamic variable p ˜ n corresponding to the number of juveniles after reproduction (blue dots). The green line is the equilibrium number of juveniles after reproduction p ˜ ¯ . The red line corresponds to the average number of juveniles. u p c r is the critical value, if the harvest rate up is higher than u p c r , then the population will go extinct. u * is the bifurcation value at which the inverse Neimark–Sacker bifurcation occurs.
Mathematics 12 00535 g006
Figure 7. The existence domains of the hydra effect are highlighted by shading at up = 0.5 (a,b) and a = 4 (c,d) and hatching at up = 0.1 (a,b) and a = 10 (c,d) with s = 0.8, h = 2, and δ = 0.5. TC, NS, and PD are the lines of transcritical, Neimark–Sacker, and period-doubling bifurcations, respectively. H is the line of the maximum pre-harvest abundance of the population.
Figure 7. The existence domains of the hydra effect are highlighted by shading at up = 0.5 (a,b) and a = 4 (c,d) and hatching at up = 0.1 (a,b) and a = 10 (c,d) with s = 0.8, h = 2, and δ = 0.5. TC, NS, and PD are the lines of transcritical, Neimark–Sacker, and period-doubling bifurcations, respectively. H is the line of the maximum pre-harvest abundance of the population.
Mathematics 12 00535 g007
Figure 8. Stability regions of fixed points (17) (a,c) and (19) (b,d) in the ( ρ , a) plane at h = 2 and δ = 0.5 with variation in the harvest rate of mature males um. TC, NS, and PD are the lines of transcritical, Neimark–Sacker, and period-doubling bifurcations, respectively.
Figure 8. Stability regions of fixed points (17) (a,c) and (19) (b,d) in the ( ρ , a) plane at h = 2 and δ = 0.5 with variation in the harvest rate of mature males um. TC, NS, and PD are the lines of transcritical, Neimark–Sacker, and period-doubling bifurcations, respectively.
Mathematics 12 00535 g008
Figure 9. The domains of the hydra effect are highlighted by shading at um = 0.1 for the fixed point (18) and hatching at um = 0.6 for the fixed point (20) in the parameter plane ( ρ , a). TC, NS, and PD are the lines of transcritical, Neimark–Sacker, and period-doubling bifurcations, respectively. H is the line of the maximum values of pre-harvest abundance of mature males. Mf is the line of the maximum equilibrium values of the abundance of mature females.
Figure 9. The domains of the hydra effect are highlighted by shading at um = 0.1 for the fixed point (18) and hatching at um = 0.6 for the fixed point (20) in the parameter plane ( ρ , a). TC, NS, and PD are the lines of transcritical, Neimark–Sacker, and period-doubling bifurcations, respectively. H is the line of the maximum values of pre-harvest abundance of mature males. Mf is the line of the maximum equilibrium values of the abundance of mature females.
Mathematics 12 00535 g009
Figure 10. Changes in the domains of the hydra effect for fixed point (20) with variations in the values of population parameters and an increase in harvest rate um. The domains are highlighted in shading at um = 0.1 (a,b) and um = 0.7 (c,d) and hatching at um = 0.6 (a,b) and um = 0.8 (c,d). TC, NS, and PD are the lines of transcritical, Neimark–Sacker, and period-doubling bifurcations, respectively. H is the line of the maximum values of pre-harvest abundance of mature males.
Figure 10. Changes in the domains of the hydra effect for fixed point (20) with variations in the values of population parameters and an increase in harvest rate um. The domains are highlighted in shading at um = 0.1 (a,b) and um = 0.7 (c,d) and hatching at um = 0.6 (a,b) and um = 0.8 (c,d). TC, NS, and PD are the lines of transcritical, Neimark–Sacker, and period-doubling bifurcations, respectively. H is the line of the maximum values of pre-harvest abundance of mature males.
Mathematics 12 00535 g010
Figure 11. (a) The stability domains of the fixed point of the system (16) and the hydra effect areas are highlighted by shading and hatching, respectively, in the parameter plane (um, a). TC and NS are the lines of transcritical and Neimark–Sacker bifurcations, respectively. H is the line of the maximum values of pre-harvest abundance of mature males. Mf is the line of the maximum values of the abundance of mature females. (bd) Dependences of the equilibrium numbers of the system (16) on the harvest rate of mature males (um) with variations in parameter a.  p ¯ ,   f ¯   , m ˜ ¯ , and   m ¯ are the stationary numbers for immature individuals, mature females, and mature males before and after harvesting, respectively. The values of u m s w separate regions Ω I ( m ¯ f ¯ / h ) and Ω II ( m ¯ < f ¯ / h ). uc is the critical value; if the harvest rate um is higher than u m c r , then the population will go extinct.
Figure 11. (a) The stability domains of the fixed point of the system (16) and the hydra effect areas are highlighted by shading and hatching, respectively, in the parameter plane (um, a). TC and NS are the lines of transcritical and Neimark–Sacker bifurcations, respectively. H is the line of the maximum values of pre-harvest abundance of mature males. Mf is the line of the maximum values of the abundance of mature females. (bd) Dependences of the equilibrium numbers of the system (16) on the harvest rate of mature males (um) with variations in parameter a.  p ¯ ,   f ¯   , m ˜ ¯ , and   m ¯ are the stationary numbers for immature individuals, mature females, and mature males before and after harvesting, respectively. The values of u m s w separate regions Ω I ( m ¯ f ¯ / h ) and Ω II ( m ¯ < f ¯ / h ). uc is the critical value; if the harvest rate um is higher than u m c r , then the population will go extinct.
Mathematics 12 00535 g011
Figure 12. (a,b) Bifurcation diagrams of the dynamic variable m ˜ n corresponding to the number of mature males after reproduction (blue dots) with respect to the parameter um under different initial conditions with m0 = 0.025. The green line is the equilibrium number of males after reproduction m ˜ ¯ . The red line corresponds to the average number of males. The parameter values are s = 0.4, v = 0.2, δ = 0.5, h = 2, a = 7.6, and ρ = 4.8. The values of u m s w separate regions Ω I ( m ¯ f ¯ / h ) and Ω II ( m ¯ < f ¯ / h ). u m c r is the critical value, if the harvest rate um is higher than u m c r , then the population will go extinct. (cg) The basins of attraction for the coexisting dynamic modes of the model (16) at m0 = 0.025. The figures correspond to the periods of observed cycles. ID stands for irregular dynamics.
Figure 12. (a,b) Bifurcation diagrams of the dynamic variable m ˜ n corresponding to the number of mature males after reproduction (blue dots) with respect to the parameter um under different initial conditions with m0 = 0.025. The green line is the equilibrium number of males after reproduction m ˜ ¯ . The red line corresponds to the average number of males. The parameter values are s = 0.4, v = 0.2, δ = 0.5, h = 2, a = 7.6, and ρ = 4.8. The values of u m s w separate regions Ω I ( m ¯ f ¯ / h ) and Ω II ( m ¯ < f ¯ / h ). u m c r is the critical value, if the harvest rate um is higher than u m c r , then the population will go extinct. (cg) The basins of attraction for the coexisting dynamic modes of the model (16) at m0 = 0.025. The figures correspond to the periods of observed cycles. ID stands for irregular dynamics.
Mathematics 12 00535 g012
Figure 13. (a,b) Dynamic mode maps of model (16) for s = 0.4, v = 0.2, δ = 0.5, h = 2, ρ = 4.8, and m0 = 0.025. The figures correspond to the period of observed cycles. 10 corresponds to the extinction of the population. ID stands for irregular dynamics. The white lines TC, NS, and PD are the transcritical, Neimark–Sacker, and period-doubling bifurcations, respectively. (c) Bifurcation diagrams of the dynamic variable m ˜ n corresponding to the number of mature males after reproduction (blue dots) with respect to the parameter um at m0 = 0.025. The green line is the equilibrium number of males after reproduction m ˜ ¯ . The red line corresponds to the average number of males. The values of u m s w separate regions Ω I ( m ¯ f ¯ / h ) and Ω II ( m ¯ < f ¯ / h ). u m c r is the critical value, if the harvest rate um is higher than u m c r , then the population will go extinct.
Figure 13. (a,b) Dynamic mode maps of model (16) for s = 0.4, v = 0.2, δ = 0.5, h = 2, ρ = 4.8, and m0 = 0.025. The figures correspond to the period of observed cycles. 10 corresponds to the extinction of the population. ID stands for irregular dynamics. The white lines TC, NS, and PD are the transcritical, Neimark–Sacker, and period-doubling bifurcations, respectively. (c) Bifurcation diagrams of the dynamic variable m ˜ n corresponding to the number of mature males after reproduction (blue dots) with respect to the parameter um at m0 = 0.025. The green line is the equilibrium number of males after reproduction m ˜ ¯ . The red line corresponds to the average number of males. The values of u m s w separate regions Ω I ( m ¯ f ¯ / h ) and Ω II ( m ¯ < f ¯ / h ). u m c r is the critical value, if the harvest rate um is higher than u m c r , then the population will go extinct.
Mathematics 12 00535 g013
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Revutskaya, O.; Neverova, G.; Frisman, E. Discrete-Time Model of an Exploited Population with Age and Sex Structures: Instability and the Hydra Effect. Mathematics 2024, 12, 535. https://doi.org/10.3390/math12040535

AMA Style

Revutskaya O, Neverova G, Frisman E. Discrete-Time Model of an Exploited Population with Age and Sex Structures: Instability and the Hydra Effect. Mathematics. 2024; 12(4):535. https://doi.org/10.3390/math12040535

Chicago/Turabian Style

Revutskaya, Oksana, Galina Neverova, and Efim Frisman. 2024. "Discrete-Time Model of an Exploited Population with Age and Sex Structures: Instability and the Hydra Effect" Mathematics 12, no. 4: 535. https://doi.org/10.3390/math12040535

APA Style

Revutskaya, O., Neverova, G., & Frisman, E. (2024). Discrete-Time Model of an Exploited Population with Age and Sex Structures: Instability and the Hydra Effect. Mathematics, 12(4), 535. https://doi.org/10.3390/math12040535

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop