Mathematical Modeling of Population Dynamics Based on Recurrent Equations: Results and Prospects. Part I

Approaches to modeling population dynamics using discrete-time models are described in this two-part review. The development of the scientific ideas of discrete time models, from the Malthus model to modern population models that take into account many factors affecting the structure and dynamics, is presented. The most important and interesting results of recurrent equation application to biological system analysis obtained by the authors are given. In the first part of this review, the population dynamic effects that result from density-dependent regulation of population, the age and sex structures, and the influence of external factors are considered.


INTRODUCTION
The "golden age" of mathematical biology began in the first half of the 20th century, with an enchanting "splash" of works that determined for a long time the development of theoretical ecology and mathematical population genetics and provided the foundations of the synthetic theory of evolution (Haldane, 1924;Lotka, 1925;Fisher, 1930;Volterra, 1931;Wright, 1931;Kostitzin, 1937;and others). The mathematical basis of these works was the models constructed using differential equations and successfully describing many phenomena observed in biological communities (population size fluctuations, (numbers), competitive exclusion, polymorphic diversity, and many others).
In the early 1970s, mainly due to the basic studies by R. May (May, 1975) and A.P. Shapiro (Shapiro, 1972;Shapiro and Luppov, 1983), mathematical population biology was supplemented by fairly simple yet very effective models based on recurrent equations, also known as discrete-time models or maps. These models seemed quite adequate to describe the dynamics of species with a seasonal breeding period and staged development. If for mathematicians recurrent equations serve as an aid or a method to study the dynamics of complex continuous-time models, then for mathematical biologists and systemic ecologists they are an object of independent study, convenient for population dynamics simulations. The popularity of recurrent equations is determined by the fairly simple technique for constructing them (e.g., representing the life cycle in the form of a graph), available study methods (phase portrait, Lamerey's ladder, etc.) and simple interpretation of model equations and simulation results. However, it turned out that these "simple" models, even in the one-dimensional case, have a colossal variety of dynamic modes that appear in a complex manner when the parameters of the model are changed, up to chaotic dynamics. Importantly, a similar dynamic behavior in continuous-time models is noted in autonomous (with constant coefficients) systems of ordinary equations with a dimension of at least three or non-autonomous systems with a dimension of at least two.
The main results of the research of recurrent equations in the context of studying the causes and mechanisms of changes in the total number or structure of populations have been reviewed and systematized. This review does not claim to cover most of the results comprehensively. A significant part of the results presented here were obtained either by the team of the authors of this article or jointly with other researchers. Many results are known in the field of dynamical systems, biophysics, and ecology. However, this article contains a number of fairly recent results related to the development of classical models and the use of modern methods for studying dynamical systems in relation to population ecology.

DYNAMICS OF POPULATIONS WITH NONOVERLAPPING GENERATIONS
In annual plants, many insect species, as well as some fishes, amphibians, and reptiles, each individual population represents one age class, and adjacent generations of such a population do not overlap. If environmental conditions practically do not change over time, then the number of a certain generation will be determined only by the size of the previous generation. Having denoted the abundance of the nth generation as N n , we can write the following deterministic equation describing the dynamics of the size of such a population: where F is the birth rate function. The simplest form of this equation for is the Malthus model, where r is a constant. It is assumed that each individual leaves, on average, r offspring in the next generation, regardless of the size of the parental population: (1) The solution to this equation is a geometric progression with the denominator r and the initial term N 0 , which is virtually identical to the exponential growth of the population size in the absence of limiting factors.
It is well known that a long exponential increase in abundance in nature is never observed. Sooner or later, the action of limiting factors manifests itself, and r in Eq. (1) becomes a function of abundance. We set r = af(N), where f(N) is a function that describes the limitation and a is a parameter that is called the reproductive potential of the population and characterizes the rate of population growth without limitation (i.e., a is selected so that the equality f(0) = 1 holds). Now, instead of Eq. (1), we obtain (2) Such models were investigated by A.P. Shapiro (Shapiro, 1972;Shapiro and Luppov, 1983) (May, 1975). They showed that the dynamics of the population size described by Eq. (2) can be very complex if the function F(N) = aNf(N) decreases sufficiently fast at large N values (for example, faster than 1/N 2 ). The peculiarity of these models is that an increase in the reproductive potential of the population entails a loss of equilibrium stability, which is accompanied by the appearance of two-year fluctuations.
Specifying the form of the function f(N), trajectories for Eq. (2) for various initial conditions can be constructed numerically. The discrete analogue of the Verhulst model, for which f(N) = 1 -kN, and the model proposed by the ichthyologist W. Ricker (Ricker, 1954), for which f(N) = exp(-kN), are well studied. When studying Eq. (2), one usually excludes the scale parameter k and proceeds to dimensionless variables ("relative" number x = kN). In this case, the discrete analogue of the Verhulst model and the Ricker model, respectively, take the form Figure 1 shows the solutions of the Ricker model for different initial conditions: two two-year cycles with different phases of oscillations. Accordingly, such a system is bistable: according to some initial conditions, the solutions converge to one two-year cycle, and according to other initial conditions, to another ( Fig. 1). This bistability is evidenced by the results of laboratory experiments (Henson et al., 1998), showing that, at the same initial abundance and similar, though not identical conditions, two different antiphase periodic modes can be observed in populations of the flour beetle Tribolium castaneum (i.e., phase multistability is observed).
With a further increase in the reproductive potential, the solutions of the equation converge to stable four-year cycles, then to eight-year cycles, etc. The changes described in the character of population dynamics are conventionally called the first series of bifurcations, and the scenario of the formation of these cycles is called the cascade of period-doubling bifurcations, which is universal for any unimodal map (Feigenbaum universality) (Feigenbaum, 1983). It also leads to a corresponding complication of phase multistability. At even larger values of parameter a, the behavior of the population size loses any regularity and becomes chaotic. Zones of chaotic behavior are interspersed with "windows" of periodic (i.e., regular) behavior.
For some particular form of the function f(N) into Eq. (2), it is possible to construct a bifurcation diagram numerically characterizing the limiting trajectories depending on the coefficient a (Fig. 2a).   x 0 = 3 x 0 = 2 x n n A special case generalizing the Verhulst and Ricker equations is the three-parameter Hassell model (Hassell, 1975): This model makes it possible to study the population size dynamics at different intensities of ecological limitation, which is characterized by the parameter β, an increase in which leads to an increase in the rate of decline f(x) = 1/(1 + x) β at large x values. For example, at small values of β (β < 2) (i.e., for a slight limitation), nonmonotonic modes of population dynamics are not observed at any values of a, which characterizes the birth rate of the species. This is true, for example, at β = 1 (i.e., in the case when the Hassell model turns into the Beverton-Holt model). However, at large values of β (β > 3), a monotonic increase in the parameter a leads to the implementation of the first series of bifurcations and then to the appearance of a chaotic dynamics (i.e., to bifurcations analogous to those in the Ricker model).
At large fixed a values (a ≥ 30), an increase in the parameter β leads to the same changes in the population size dynamics as an increase in the parameter a at fixed β values. The difference is that, with the increase in β, bifurcations occur against the background of a decrease in the equilibrium level of abundance, whereas with an increase in a they occur against the background of its increase. Thus, the space of parameters a and β is divided into domains, each of which is characterized by its own type of dynamic behavior.
The study of the modes of dynamic behavior of Ricker, Hassell, and other models makes it possible to discover some general patterns that appear at suffi-1 (1 ) . n n n N aN kN β + = + ciently large values of the reproductive potential and the intensity of ecological limitation, when the dynamics of abundance may become chaotic or close to it (interval-periodic). If the initial population size is small, then it will increase slowly for a sufficiently long time (which may be accompanied in some generations by even small declines). This is followed by a sharp increase in abundance, accompanied in the next generation by a significant decrease in it to a value close to the initial one. However, these periodic "breaks" will not bring the population back exactly to the initial level. Therefore, despite the obvious periodic pattern of changes in abundance, complete coincidences will be found neither in abundance nor in the number of generations in the growth phase (Shlyufman et al., 2013). Such a loosely periodic behavior is characteristic of many natural populations of higher organisms, especially insects (e.g., locusts, grasshoppers, and night moths). Attempts to use one-dimensional discrete-time models to describe and predict the dynamics of certain natural populations were often unproductive: model curves, catching trends in changes in abundance, statistically poorly described the observed dynamics of real populations (Shapiro and Luppov, 1983;Dennis and Taper, 1994;Myers et al., 1999;Frisman et al., 2007Frisman et al., , 2015aNedorezov and Sadykova, 2008;Nedorezov, 2010;Tarasov et al., 2012). In some cases, the model description could be improved by introducing a delay into the equations (Nedorezov, 1986;Turchin, 1990Turchin, , 2003Berryman and Turchin, 2001;Isaev et al., 2001;Nedorezov, 2012;Sadykova and Nedorezov, 2013;Nedorezov and Sadykova, 2015). Delayed regulation occurs in natural populations as a result of interspecific interaction or in cases when a No. 1 2021 high population density adversely affects the reproduction of the next generation (Prout and McChesney, 1985;Turchin, 1990;Williams and Liebhold, 1995). The most convincing results were obtained for the Ricker (Moran-Ricker) model and its modification with delay (Kendall et al., 1999;Turchin, 2003;Turchin et al., 2003;Bechtol and Kruse, 2009;Sadykova and Nedorezov, 2013;Nedorezov and Sadykova, 2015). In particular, an attempt was made to describe and analyze the dynamics of Zeiraphera diniana Gn. on the basis of the Moran-Ricker model at different lag values (Sadykova and Nedorezov, 2013;Nedorezov and Sadykova, 2015).

DYNAMICS OF POPULATIONS WITH NONOVERLAPPING GENERATIONS AND STAGE STRUCTURE
In modeling population dynamics, in addition to the seasonality of reproduction, it is necessary to take into account its age structure, because the survival rate and the involvement of individuals of different ages in reproduction differ. Indeed, in the life cycle of any organism, several age stages can be distinguished, which are determined in certain units of time (e.g., in years) (Leslie, 1945;Caswell, 2001). In cases when the age of individuals is not known, individuals in a structured population are classified according to their stage of development rather than according to their chronological age (Lefkovitch, 1965;Caswell, 2001;Logofet and Klochkova, 2002;Logofet, 2008). In this case, the population is divided into a certain number of groups, the way of division into which is usually determined by the biological characteristics of organisms. For example, the model of dynamics of a population with nonoverlapping generations x n + 1 = ax n exp(-x n ) and staged development of individuals takes the form (3) where X n is the number of one-year-old individuals, Y n is the number of mature individuals in the nth breeding period (more precisely, by the beginning of the breeding period in the nth year), B is the reproductive potential of mature individuals, α and β are coefficients characterizing the influence of juvenile and of mature group sizes to a decrease in the birth rate, and s is the survival rate of individuals in the second year of life. It is convenient to study model (3) in the space of parameters ρ = α/(sβ) and r = sB, where ρ describes the relative contribution of the younger age group to limitation of reproduction and r characterizes the reproductive potential of individuals with allowance for the survival rate of juveniles.
It was found that the nontrivial fixed point of system (3) may lose stability according to both the period-doubling scenario and the Neimark-Sacker scenario, a discrete analog of the Poincaré-Andronov-Hopf bifurcation that leads to quasiperiodic dynamics (Kuznetsov et al., 2012). As a result, the discrete-time models, taking into account the staged structure, can describe not only long-period oscillations, which are observed in models with continuous time, but also the "noisy oscillations." Modern computational technologies make it possible to identify the areas of different dynamic modes in the model parameter space, and create the dynamic mode maps. The dynamics mode map of model (3), presented in Fig. 3 on the left-hand side shows a picture characteristic of the Neimark-Sacker bifurcation, which consists in the emergence of a region of quasiperiodic modes Q with a system of Arnold tongues immersed in it, and on the right-hand side shows a cascade of period-doubling bifurcations (Fig. 3a). Therefore, when the birth rate is limited primarily by the abundance of mature individuals, then quasiperiodic oscillations occur. When the birth rate decreases with an increase in the size of the juvenile group (i.e., at ρ > 1), then an increase in the birth rate leads to periodic fluctuations in abundance. Figure 3b shows the attraction basins of the dynamic modes of model (3) in the multistability area, where the three-year cycle is stable simultaneously with the fixed point. Under some initial conditions, the system is stabilized, whereas under others it exhibits three-year fluctuations. It should be noted that such multistability was shown earlier, in particular, for the two-dimensional Hénon map and some of its modifications (Saucedo-Solorio et al., 2002;Shrimali et al., 2008;Pisarchik and Feudel, 2014;etc.). For the Hénon map, the mechanism of the occurrence of a three-year cycle and the limiting invariant curve based on it, coexisting with a fixed point, was demonstrated (Romera et al., 2001).
Multistability within the framework of recurrent models of the dynamics of structured populations makes it possible to take a fresh look at general biological problems and notice the existing fundamental patterns that had previously remained unnoticed. In particular, multistability explains the observed differences in the dynamics of the abundance of populations of the same species dwelling under nearly identical conditions. Within the framework of a local population, the phenomenon of multistability makes it possible to explain both the appearance and disappearance of fluctuations in population size, as well as the changes in the period of the observed fluctuations. The most striking examples of changes in dynamic modes are the disappearance of cycles in populations of lemmings (Coulson and Malo, 2008;Kausrud et al., 2008;White, 2011) and some vole species (Henttonen and Wallgren, 2001;Cornulier et al., 2013). It should be understood that the dynamics model may change not only in the case of the existence of several dynamics modes at the same parameter values but also due to phase multistability, when, under the perturbation of the equation variable, the model trajectory shifts over the basins of attraction of different phases of the same dynamic mode. The longer the period of the observed oscillations, the more phases of this cycle with their basins of attraction can exist, which in a nonstationary environment will complicate the population dynamics.

DYNAMICS OF A POPULATION WITH TWO STAGES OF DEVELOPMENT
Consider a population with an age structure that by the beginning of the next breeding period can be represented by a combination of two classes: younger with abundance X n and older with abundance Y n in the nth breeding period. The breeding period ends with the appearance of newborns of the next generation. We assume that the time between two successive breeding periods is sufficient for the development of newborns to the young age. We also assume that the survival rate and birth rate of individuals of the older age class does not depend on their chronological age. This is valid for organisms with a short lifespan, which includes twothree breeding periods, as in many insects, fishes, small mammals, two-and three-year-old plants, etc.
The following system of equations that describes the dynamics of the age groups considered in adjacent generations can be written as follows: (4) where A 1 and A 2 are birth rates and B and C are the survival rates of the younger and older age classes, respectively. If A 1 , A 2 , B, and С are considered constants, then system (4) is a special case of the Lefkovitch linear matrix model (Lefkovitch, 1965;Caswell, 2001). This model is actively being developed both in theo- , n n n n n n n n n n n n n n retical and applied aspects in many studies, for example, in the studies by Logofet et al. (Logofet et al., 2006(Logofet et al., , 2017Logofet and Belova, 2007;Logofet, 2019). A matrix model aimed at describing the dynamics of the cenopopulation of Eritrichium caucasicum was proposed and calibrated (Logofet et al., 2016. However, the majority of biological species have density-dependent regulation of population growth (Lack, 1954;Gimmel'far et al., 1974;Dajo, 1975;Odum, 1975;Williamson, 1975;Svirezhev and Logofet, 1978;Boer and Reddingius, 1996;Gurney and Nisbet, 1998;Inchausti and Ginzburg, 1998;Ginzburg and Colyvan, 2004;Barraquand et al., 2017). The main manifestation of density-dependent birth rate regulation is the stress syndrome, which leads to a decrease in sexual activity and fertility of individuals, up to the resorption of part of the embryos. In particular, this is typical for species that experience strong fluctuations in abundance, for example, for lemmings and voles (Dajo, 1975;Chernyavskii and Lazutkin, 2004;Frisman et al., 2010a;Novikov et al., 2012;Krebs, 2013). In modeling the population dynamics of such species, the survival rates can be considered constant (B(X, Y) = b, C(X, Y) = c), and the birth rate, by analogy with the Ricker model, can be written as where r 1 and r 2 are the reproductive potentials of the younger and older age groups, and α and β are the coefficients characterizing the intensity of the effect of the abundances of the immature and mature age classes on the decrease in birth rate. Taking this into account and having proceeded to the dimensionless variables x = bβX and y = βY, we can transform model (4) to the form where a 1 = r 1 and a 2 = br 2 are the new designations for the coefficients that characterize the reproductive potentials and ρ = α/(bβ) describes the relative contribution of the younger age group to the birth rate limitation. Model (5) is primarily focused on describing the dynamics of populations characterized by a high maturation rate of juveniles. In such populations, mature individuals give birth to several litters per season, and the grown-up underyearlings reach maturity in the same season and begin to breed and produce offspring.
Note that a particular case of model (5) at a 1 = 0 corresponds to the situation when juveniles reach maturity and begin to participate in breeding in a year (Frisman et al., 2010c(Frisman et al., , 2011Revutskaya et al., 2016;Neverova et al., 2018).
The study of model (5) showed that a decrease in birth rate with an increase in the abundance of individuals is an effective mechanism for regulating the increase in abundance (Frisman et al., 2015bNeverova et al., 2019). However, at a high reproductive potential and high survival rates of individuals, a drop in the birth rate with an increase in abundance may lead to a loss of stability and the occurrence of fluctuations. For example, at relatively large values of a 2 or ρ (ρ ≥ 1) (i.e., in the case when the older age class primarily contributes to reproduction, and the birth rate is limited primarily by the younger age class), the loss of stability, as in the original one-dimensional Ricker model, is accompanied by the appearance of two-year fluctuations and a subsequent period-dou- x a x a y x y y x cy bling cascade. In this case, the larger the ρ value, the lower the value of the reproductive potential at which fluctuations occur.
At large a 1 or small ρ values, i.e., in the cases when underyearlings make the main contribution to reproduction, the loss of stability is accompanied by the appearance of quasiperiodic oscillations (Fig. 4a). Intensive reproduction of underyearlings and a decrease in the birth rate with an increase in the abundance of breeding individuals may lead to the emergence of complexly organized fluctuations in abundance. This mechanism apparently determines the behavior of the abundance of small mammals, such as lemmings and some voles (Chernyavskii and Lazutkin, 2004;Krebs, 2013).
In the parametric space of model (5), along with irregular dynamics, there are regions with regular dynamics-periodic oscillations (Fig. 4b). An important difference of the behavior of systems of dynamics of the populations with overlapping generations from that of the populations with nonoverlapping generations is the realization of a cascade of period-doubling bifurcations not only with an increase in the reproductive potentials of individuals but also with a change in the values of the parameter characterizing the limitation of the birth rate.
The study of age-structured models revealed their multistability. For example, a transition to a stable nontrivial fixed point or stable three-or four-year fluctuations is possible. In this case, the dynamic mode type is determined by the initial conditions, a slight variation of which may change the type of dynamics. Note that periodic changes in abundance with the predominance of three-four year cycles were

MULTISTABILITY IN POPULATION DYNAMICS MODELS WITH ALLOWANCE
FOR AGE AND SEX STRUCTURES Many authors (Caswell, and Weeks, 1986;Lindstrom and Kokko, 1998;Caswell, 2001;Shyu and Caswell, 2016;etc.) have studied discrete-time mathematical models of the dynamics of the population size with the age and sex structure. Study of the corresponding models showed (Frisman et al., 1982(Frisman et al., , 2010bCaswell and Weeks, 1986;Lindstrom and Kokko, 1998;Caswell, 2001;Revutskaya et al., 2012) that the population dynamics pattern is affected by factors such as the size of the harem, the density-dependent regulation, the competition of males for the female, and the peculiarities of the territorial distribution of females and males.
Consider a population with a seasonal breeding pattern consisting of three groups: the youngest, which includes immature individuals (P), and two older ones, represented by mature females (F) and males (M), in the nth breeding period. It is assumed that the time between two successive breeding periods is sufficient for the development of immature individuals to maturity.
It is assumed that birth rate depends on the ratio of abundances of mature females and males and can be described using the pairing function (Caswell and Weeks, 1986): where h corresponds to the average size of the harem and characterizes the type of mating relations in the population (in monogamy h = 1, in polygyny h > 1, in polyandry h < 1). In order to avoid "overestimating" the birth rate for the populations in which females give birth once per breeding period, the pairing function can be transformed to the following form (Bessa-Gomes et al., 2010): The condition for switching the pairing function corresponds to the balance of sexes in the population and has the form F = hM, where hM is the number of females that can potentially be fertilized by M males at an average harem size h.
Models with the pairing function can be used to analyze and describe the dynamics of insects, reptiles, birds, and mammals (Frisman et al., 1982;Molnar et al., 2008;Miller and Inouye, 2011;Gerber and White, 2014).
Assume that the survival rates of immature females and males by the onset of maturity do not differ and decrease linearly with increasing age class sizes, and where α and β are coefficients characterizing the intensity of the decrease in the survival rate of immature animals, which is caused by the competitive interaction between sex and age classes, respectively. We assume that, with negative values that arise at large numbers, the survival function of immature is zeroed out. From a biological point of view, this can be interpreted as the death of the offspring of the corresponding year due to the high intraspecific competition for resources.
With allowance for the assumptions described, the model can be written as a system of three recurrent equations: where a is the birth rate (the average number of offspring per pair), δ is the proportion of females among newborns, and s and v are the survival rates of mature females and males, respectively.
The stability loss of the fixed point of system (6) may proceed according to the Neimark-Saсker and period-doubling scenarios. With a further change in the model parameters, a transition to chaotic dynamics may occur. Note that the observed complex bifurcations are also characteristic of other maps, for example, for the three-dimensional Hénon map and its modifications (Gonchenko, A.S. and Gonchenko, S.V., 2016;Jiang et al., 2016;etc.). There are regions of multistability in the parametric space of model (6). The phase space of the model is split by the attraction basins of coexisting dynamic modes. In the dynamics of real populations, "jumps" over basins correspond to a change in the dynamic mode, which manifests itself either as a change in the oscillation period or as the appearance or disappearance of fluctuations. Coexisting modes arise both as a result of bifurcations and as a result of switching the pairing function.
The study of model (6) revealed complex relationships between the sex ratio and the pattern of the dynamic behavior of the population . In populations the dynamics of which can be described by model (6), transitions between different modes arise not only as a result of changes in the population parameters that determine the processes of reproduction and self-regulation but also due to a change in the principle of pair formation. The revealed multistability allows us to conclude that the variation in the population size, leading to a change in the sex ratio, complicates the population dynamics. Difficulties arise in identifying the type of dynamic mode observed, because the trajectory is formed under the influence of self-regulation processes and a change in  the pair formation type, which is caused by the difference in the abundance of sexes.

CHANGING THE DYNAMIC MODES IN POPULATIONS: INFLUENCE OF MODIFYING FACTORS
The influence of external factors on the cyclical dynamics of populations is an important field of research in mathematical ecology. In this case, not only the regular fluctuations in abundance but also clear transitions from one dynamic mode to another are analyzed. The influence of external and intrapopulation factors differs for different species of plants and animals. Moreover, it may vary significantly even for the same species depending on the natural and climatic zone in which it dwells. In this section, we consider a situation when the population dynamics is determined by both autoregulation processes and the influence of external factors.
Influence of random factors on the dynamics of a population with two age classes. Application of model (5) to the description of the population dynamics of the bank vole Myodes glareolus based on the data of long-term censuses at the Udmurt station (Zhigalskii, 2011) gave a satisfactory approximation. The model trajectory describes the trend in dynamics but does not fully fix the peaks of the bank vole population abundance (Frisman et al., 2015b). Apparently, the discrepancy between the observational and modeling data is due to the influence of external factors.
One of the main factors affecting the reproductive activity of many rodents is the reserve of food resources. The increase in the breeding activity with an increase in the abundance of food is "explosive." All other things being equal, this process can be described with an exponential function. This approach, in particular, was used successfully to describe the dependence of the increase in the birth rate of Manchurian squirrels with an increase in the yield of pine nuts (Ashichmina et al., 1985). Apparently, there are no adequate direct estimates of food abundance for the bank vole. We propose to use the Selyaninov hydrothermal coefficient (S n ) as an indicator of the abundance of food reserves, which characterizes the moisture supply of the territory. This external factor made it possible to fix the main peaks of the population size and improve the adequacy of the model. This is due to the fact that the coefficients characterizing the reproductive potentials of individuals are not constant; their values depend on the moisture supply of the territory (Fig. 5).
Further study of the model with two age classes taking into account the influence of climatic factors on reproductive processes showed that the population abundance constantly moves from one basin to another and only in certain years is into the range of parameter values at which similar dynamics are observed .
Ichthyology was one of the first areas to face the problem of sustainable use of self-reproducing resources. At the end of the 19th century, it was noted that the proportion of some fish species (sturgeon, beluga, carp, bream, roach, etc.) in catches began to decrease every year. Therefore, strategies focused on optimizing the harvesting of fish populations have been developed most (Baranov, 1918;Graham, 1935; Schaefer, 1954;Ricker, 1954;Beverton and Holt, 1957;etc.). Particular attention was given to the problems of harvest optimization for homogeneous population models (Svirezhev and Elizarov, 1972;Swan, 1975;Skaletskaya et al., 1979;Abakumov, 1993). However, a scientifically substantiated exploitation of biological resources requires a detailed study of the population structure (Svirezhev and Elizarov, 1972;Beddington and Taylor, 1973;Colley, 1979;Skaletskaya et al., 1979;Jensen, 1996;Langvatn and Loison, 1999). Firstly, this is due to the fact that population replenishment is a complex process that includes not only the birth rate but also the survival rate of immature individuals, as well as the transition of younger individuals to older age classes, etc. All these processes are affected differentially by changes in the population density and harvesting (Bergman et al., 2015). Secondly, in most cases, it is economically profitable to use selective harvesting (Colley, 1979). For example, sturgeon, whitefish, and salmon have different commercial importance depending on their age (Frisman and Last, 2005).
The best known are two harvest management strategies, which imply the density regulation of population growth. The first strategy is based on the concept of "annual harvest surplus," and the second is based on the concept of maximizing a sustained yield (Bergman et al., 2015). The study of the dynamics of a twoage exploited population showed that the maximum sustainable (optimal) yield is possible when individuals of only one age class are withdrawn (Zhdanova and Frisman, 2013;. It is well known that an optimal harvesting rate, which ensures the maximum sustained yield in a homogeneous population (i.e., without taking into account the age structure), leads to stabilization of its size (Frisman et al., 2003). At the same time, it was found that the dynamics of an exploited population with an age structure becomes more complicated if we take into account the fact that individuals of different ages affect the birth rate or survival rate with different intensity. It was shown that equilibrium harvest of a two-age population, based on the strategy of a constant optimal harvesting rate, is possible in a limited parametric space area, where the stationary solutions are stable . In the case of stability loss of the fixed point, this strategy leads to fluctuations and ceases to be both equilibrium and optimal. The system dynamics is stabilized if the harvesting strategy is based on a regular harvesting of the "surplus" over the value corresponding to the maximum reproduction of the population .
Strategies based on maximum sustained yield may lead to catastrophic consequences, up to population degeneration (Larkin, 1977;Ludwig et al., 1993;Lande et al., 1995;Finley, 2011). Irreversible changes in exploited populations can be caused not only by harvesting but also by the self-regulation processes in combination with environmental factors, leading to variations in the population growth rate and, consequently, to fluctuations (Fryxell et al., 2010) or even to a change in the dynamic mode.
Note that studies that take into account the possibility of multistability in exploited populations (Saucedo-Solorio et al., 2002;Pisarchik and Feudel, 2014;Frisman et al., 2015b) are quite rare. However, this approach makes it possible to study and take into account the possibility of changing the dynamic mode in an exploited population. In particular, it was shown that, for an exploited two-age population (both with and without sex structure), it remains possible to implement various dynamic modes at the same population parameter values depending on the initial abundance, which is characteristic of a freely developing population Revutskaya et al., 2018). Therefore, to stabilize the dynamics, the current population size should be "retained" in the basin of attraction of the mode of interest (for example, by reducing the abundance of one of the age and sex classes). However, difficulties in predicting population dynamics arise, because harvesting may shift the current abundance from one basin of attraction to another and lead to a change in the type of population dynamics.
Influence of the resource recovery rate. To study the effect of the resource recovery rate on population dynamics, it is proposed to use the Moran-Ricker model with a time lag (7) where x n is the population size at which it enters the nth breeding period and a is the reproductive potential of the population. The exponential factor characterizes the ecological limitation of the population growth, and m is the time lag value (i.e., the number of generations over which the limitations of life resources have an impact).
When m = 0, Eq. (7) is reduced to the classical Ricker model (Moran, 1950;Ricker, 1954), in which it is assumed that each generation "receives" the same amount of ecological resources required for life, regardless of the previous history of the population. This is true for those species the resource base of which can fully recover for the time between breeding periods. When m = 1, it is assumed that the resources required for the life activity of the population are able to recover fully within the time between the two breeding seasons. Equation (7)  The study of model (7) at different time lag values showed that the contribution of previous generations to the ecological limitation of reproduction processes leads to a complication of population dynamics . In general, the dependence of population dynamics on the resource recovery rate leads to the fact that, in most cases, the population exhibits quasiperiodic dynamics. Two-year fluctuations occur when the population growth is limited by densitydependent regulation in the current year.
The Moran-Ricker model at different lag values was used to describe the dynamics of insect populations . Figure 6 shows the results of describing the dynamics of populations of the pine beauty moth , the pine looper moth , and the larch bud moth (Sadykova and Nedorezov, 2013;Nedorezov and Sadykova, 2015) using model (8). Estimates of the model parameters show that, in all cases, b 1 is much greater than b 0 (ρ = b 1 /b 0 is greater than 1). Thus, the contribution of the previous generation to limiting reproduction of the population is significantly greater than the contribution of the current generation. The point estimates of the pine beauty moth and pine looper moth are located in the region of quasiperiodic oscillations, where, in the case of small fluctuations in the a and ρ values, the dynamics pattern is retained. However, the parameters may also get into the zone of one of the Arnold tongues and quasiperiodic dynamics will transit to strictly periodic ones. The point estimate of the parameters for the larch bud moth is in the periodicity window that corresponds to 18-year fluctuations. The estimate of the parameters lies near the border separating the 9-and 18-year cycles, and fluctuations in the parameter values may lead to a change in the dynamic mode .
Influence of periodic changes in environmental factors. Of special interest are periodic changes in environmental factors, which, in particular, may lead to periodic changes in the reproductive potential of populations with a seasonal reproduction pattern. This situation can be described using the Ricker equation with a periodic Malthusian parameter, which makes it possible to take into account the cyclical effects of both exogenous and endogenous factors on the population size (Shlyufman et al., 2016(Shlyufman et al., , 2017(Shlyufman et al., , 2018Frisman et al., 2019). The study showed that the Ricker equation with a Malthusian parameter varying with a period of two, along with the multistability of dynamic modes, also has phase multistability. As a result of perturbation of a variable or a phase shift in the Malthusian parameter, a phase shift in the oscillations or a change in the existing mode can be observed. The model has two different stable two-year cycles and two different series of their bifurcation according to the period-doubling scenario. Therefore, the two-year fluctuations in abundance in a periodically changing environment can be both synchronous and asynchronous to the fluctuations in the environment.

CONCLUSIONS
In modern biology, a colossal amount of time series of empirical data indicating that biological systems have a complex dynamics has been accumulated. The most striking example is the disappearance of cycles in populations of small mammals, or a change in the oscillation period. In addition, different dynamic modes can be observed in populations of the same species that dwell in the most similar conditions. All these phenomena are fully explained by the possible multistability of the dynamics of natural populations. This determines the necessity of comprehensive theoretical studies of the dynamics of nonlinear systems and dynamic chaos in ecological systems. To date, the complex irregular dynamics of physical, technical, and mathematical systems has been investigated, and features of their dynamics such as regular  and chaotic oscillations, bistability and multistability, synchronization and clustering, etc., have been discovered.
There are numerous modern studies on mathematical modeling of population dynamics. Conventionally, the dynamics of natural systems is described using continuous-time models. An important step in understanding the mechanisms that cause population fluctuations was the application of fairly simple recurrent equations to studying and modeling the population dynamics of those biological species that are characterized by a clearly expressed breeding period and staged development. The recurrence equations made it possible to explain the fluctuations in the population size under relatively constant external conditions. It was found that the periodic oscillations observed in living systems can be determined not only by external influences but also by the internal properties of the population system itself. In addition, in experimental and field studies, data are often collected regularly at certain time intervals, which also indicates the feasibility of using recurrent equations.
The first part of this review presents the methods and approaches to modeling population dynamics using discrete-time models. The evolution of scientific ideas on the way of complication of discrete-time models is described-from the Malthus model to the modern models that take into account various factors affecting the population structure and dynamics and reflect a wide range of possible dynamic modes. The dynamic effects resulting from the density-dependent regulation, the complication of age and stage structures, and the influence of external factors are considered. In the second part of this review, we will consider the dynamic effects arising in problems of simulation of the evolutionary process, as well as the influence of migration on population dynamics and the features of spatial distribution.

FUNDING
The reported study was funded by the Russian Foundation for Basic Research, project number 19-14-50326.

COMPLIANCE WITH ETHICAL STANDARDS
The authors declare that they have no conflict of interest. This article does not contain any studies involving animals or human participants performed by any of the authors.

OPEN ACCESS
This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.