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:
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
[
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):
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 . Therefore, the coefficients rF and rM are and , 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.,
. 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]:
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,
, then the number of pairs formed corresponding to the possible abundance of fertilized females will equal the number of mature females
. Consequently, the fertility values of female
rF and male
rM are as follows:
When
, 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.,
. In this case, the fertility values of female
rF and male
rM are as follows:
Thus, the fertility functions of females and males can be rewritten as follows:
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: , . Here, αi (αi ≥ 0), βi (βi ≥ 0), and γi (γi ≥ 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: . 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
. 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:
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:
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:
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
f ≤
hm, then
and the coordinates of the non-zero fixed point are determined by the following formulas:
where
,
. Solution (7) has biological meaning when the number of age classes is positive, that is,
.
If
f >
hm, then
and the coordinates of the non-zero fixed point are as follows:
and are positive with
, where
, and
.
The boundaries of the stability area of fixed points in the system (6) are found based on the characteristic polynomial.
where
S,
H, and
D are three invariants of the Jacobian matrix of the system (6).
Let us rewrite Equation (9) in the form
, where
J is the Jacobian matrix of system (6) evaluated at the fixed point
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:
where
,
, and
.
The coefficients of the polynomial
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:
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 at ;
- (2)
The period-doubling bifurcation line (PD) is at ;
- (3)
The Neimark–Sacker bifurcation line (NS) is at and .
The coefficients of the polynomial (9) for system (6) at
f ≤
hm take the following values:
The coefficients of the polynomial (9) for system (6) at
f >
hm are as follows:
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 and corresponds to Equality 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:
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 , , and , 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 , 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., , and the number of pairs formed corresponds to the number of mature females, which is , 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:
The coordinates of the non-zero fixed point of the model (10) after reproduction but before harvesting are as follows:
In the case of a shortage of males in the breeding population , 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, , and the coordinates of the non-zero fixed point of the model (10) are found as follows:
The solution after harvesting is as follows:
The solution after reproduction but before harvesting is as follows:
Non-zero positive values of the equilibrium numbers of the system (10) are possible only if , where is the critical value and equals for solutions (11)–(12) or for solutions (13)–(14). If the harvest rate up is higher than , then the population will go extinct.
The equality
corresponds to the balance of sexes in a stable population. Substituting the expressions for the equilibrium numbers of the model (10) into
, we can obtain the threshold value of the parameter
, providing the balance of sexes:
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
, then (i) the number of males for breeding is sufficient or exceeds the required number because
; (ii) the offspring abundance is determined only by the number of mature females, and (iii) the equilibrium numbers are determined by the expressions
. However, with
, the population has a shortage of males because
, and some females cannot participate in reproduction; thus, the equilibrium numbers of population groups are
.
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:
The coefficients for solutions (13) and (14) are as follows:
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
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,
.
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
, 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
, where
,
. After detecting periodicity, the corresponding pixel on the diagram is colored according to the period obtained. If
, then we determine such dynamics to be irregular. Additionally, we check the condition
at each iteration step, as the survival function of juveniles
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 . 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
and is found as follows:
where
.
The maximums of the equilibrium numbers of mature females and males can be found in the formulas and and are both achieved when the harvest rate is the same as the maximum pre-exploitation abundance of juveniles . 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
.
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
(
is the solution of equation
), 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 <
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
, the juvenile abundance after reproduction can reach
, while the post-harvest abundance of
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
, then the hydra effect occurs. Consequently, the hydra effect is hidden if the abundance of
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
, 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
(
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
, 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
, 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:
where the harvest rate
um corresponds to the fraction of captured mature males.
4.1. Fixed Points of Model (16) and Their Stability
If , then and the stationary population size is determined by the following formulas:
where
.
If
, then
and the solution after harvesting is as follows:
The solution after reproduction but before harvesting is as follows:
Fixed points (17)–(18) are positive and exist with
. Solutions (19)–(20) exist at
, where
is the critical value. If the harvest rate
um is higher than
, 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:
When , the sex ratio considering harem size h is . 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 , then the number of mature females is greater than the maximum number of pairs: . 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., , 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 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 , 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:
The coefficients for solutions (19) and (20) are
Let us consider the change in the stability domain of the system (16) depending on the value of the harvest rate
um at
, which corresponds to an unexploited population with a sufficient number of males (
Figure 8a,c). With
um >
, 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,
, 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:
(
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
, 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
(
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
, an increase in
um up to
narrows the stability area concerning values of the parameter of intraspecific competition (
) (
Figure 8a). If
um >
, then the sex ratio in the population is
, and further growth of
um expands the stability domain concerning
(
Figure 8b). If initially, with
, 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 .
With values of the system (16) parameters corresponding to a non-exploited population with a sufficient number of males, i.e.,
, the hydra effect area does not depend on the harvest rate
um and is located between the lines:
and
, and
(
Figure 9a–d). In
Figure 9a–d, the graphs
and
correspond to the lines
TC(1) and
H(1), respectively. Note that the curves
and
intersect at
. At the same time, the maximum equilibrium numbers of juveniles and mature females, which are calculated using the formulas
and
, are achieved at values greater than
and coincide with the lower boundary of the hydra’s area. Note that the boundary
is a condition for the existence of a nontrivial solutions (17) and (18) of the system (16). Therefore, if
, then the equilibrium numbers of juveniles
and mature females
always increase with the growth of harvest rate of males
um.
When
, 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
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
, 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
(
Figure 9). At the same time, the maximum equilibrium population of mature females is found from
, which corresponds to the curve
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
, is achieved at biologically meaningless values of the system (16) parameters. Thus, with
, 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:
, 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
, 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
increases with a higher harvest rate of males (
um) if the population has a sufficient number of males, i.e.,
(
Figure 11b–d). Conversely, with
, 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
always decreases. However, the equilibrium number of males after reproduction but before harvesting
can reveal the hydra effect.
With the values of the model parameters from hydra’s domain, if the sex ratio is skewed toward males
, then an increase in the harvest rate of males (
um) up to the threshold value
leads to a growth of male numbers after reproduction
(
Figure 11a).
In the case of a non-high birth rate (
a), the pre-exploitation number of males
rises to the level of
, after which it begins to decline at
um >
and
. Such behavior is depicted in
Figure 11b for
a = 2.5. In this case, the threshold value
coincides with the harvest rate
, and the maximum equilibrium pre-exploitation number is achieved at
um =
. 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
.
With a higher birth rate
a, the maximum equilibrium number of females is achieved at
um >
, corresponding to
, as shown in
Figure 11c for
a = 3.5. As before, the hydra effect is observed at the harvest rate value of
um =
.
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
, with a higher harvest rate of males, the pre-harvesting equilibrium number
increases slowly. After surpassing the threshold value
, the sex ratio is skewed toward females:
. Moreover, the parameter values continue to fall within the hydra domain, and pre-harvest abundance
continues to increase with growth in
um. When
um surpasses
, the equilibrium male number after reproduction begins to decrease. Therefore, the pre-harvest number of males increases to
at
with
and then reaches
at
with
.
4.3. Instability and Multistability
A further growth of the reproductive potential (
a) destabilizes the dynamics of the system (16), when
, 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
, that is, the number of formed pairs in the population corresponds to the number of females, and irregular dynamics emerge with
. 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
, 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
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
and then at
(
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.,
, 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
, 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].