Dynamics of populations with delayed density dependent birth rate regulation

Ecological Modelling 340 (2016) 64–73

Contents lists available at ScienceDirect

Ecological Modelling journal homepage: www.elsevier.com/locate/ecolmodel

Dynamics of populations with delayed density dependent birth rate regulation G.P. Neverova a,b,∗ , I.P. Yarovenko c,d , E.Ya. Frisman b a

Institute of Automation and Control Processes of the FEB RAS, 5 Radio st., 690041, Vladivostok, Russia Institute for Complex Analysis of Regional Problems of the FEB RAS, 4 Sholom-Alekhem st., 679016, Birobidzhan, Russia, Far Eastern Federal University, 8 Sukhanova st., 690950, Vladivostok, Russia, d Institute for Applied Mathematics of the FEB RAS, 7 Radio st., 690041, Vladivostok, Russia b c

a r t i c l e

i n f o

Article history: Received 19 May 2016 Received in revised form 2 September 2016 Accepted 3 September 2016 Keywords: Population dynamics Delayed density dependence Mathematical modeling Moran–Ricker model with time lag Multimodality Dynamics of insect population

a b s t r a c t This study focuses on the Moran–Ricker model with time lag for the dynamics of limited homogeneous population. The reduction of required resources is caused by their consumption by the previous generation, which corresponds to delayed density dependence. The model shows multimodality. The population can demonstrate different types of dynamics (i.e., stable, periodic, or chaotic) at the same values of demographic parameters. The Moran–Ricker model with time lag can successfully describe the population dynamics of some insect species. The point estimations of population parameters are statistically significant. The obtained point estimates that correspond to real population dynamics are located in a quasi-periodic fluctuation area and adjoin other dynamic modes. A variation of the values of demographic parameters or random fluctuations of the current population size can lead to dynamic mode change. The investigation of the possible model dynamic modes near the found point estimation of parameters is necessary to identify the dynamic modes between which transition may occur. © 2016 Elsevier B.V. All rights reserved.

1. Introduction The “golden age” of mathematical biology began in the first half of the 20th century and was marked by a burst of studies that determined the subsequent development of theoretical ecology (Lotka, 1925; Volterra, 1931; Kostitzin, 1937). The mathematical bases of these works were elegant models based on systems of differential equations. These models can describe many population phenomena observed in natural biological communities successfully, including population dynamics fluctuations, interspecies relations, competition displacement, and autoregulation processes. Owing to the works of May (1974, 1975), Shapiro (1972) and Shapiro and Luppov (1983), the mathematical population biology in the early 70 s acquired simple and impressive models based on recurrence equations. These models adequately describe the dynamics of the species with seasonal breeding. Incidentally, the recurrence equations are easily analyzed numerically using rapidly evolving computer technology. These “simple” models have a vari-

∗ Corresponding author at: Institute of Automation and Control Processes of the FEB RAS, 5 Radio st., 690041, Vladivostok, Russia. E-mail address: [email protected] (G.P. Neverova). http://dx.doi.org/10.1016/j.ecolmodel.2016.09.005 0304-3800/© 2016 Elsevier B.V. All rights reserved.

ety of dynamic modes, which evolve in a very complicated way with the change in values of model parameters. Repeated attempts to directly apply simple models based on recurrence equations to describe and forecast the dynamics of natural populations were often unproductive: model curves, while capturing the trend changes, poorly described dynamics of real populations (Dennis and Taper, 1994; Myers, 1999; Nedorezov and Sadykova, 2008; Nedorezov, 2010; Frisman et al., 2015b). The description of the population dynamics significantly improves if time lag is introduced in the equations (Turchin, 1990, 2003; Berryman and Turchin, 2001; Isaev et al., 2001; Sadykova and Nedorezov, 2013; Nedorezov and Sadykova, 2015). Delayed density-dependence can arise as a result of interspecific interaction or the adverse effect of high population size on the fecundity of the next generation (Prout and McChesney, 1985; Turchin, 1990; Williams and Liebhold, 1995). The strongest results in studies describing the dynamics of real populations have been obtained for the Moran–Ricker model and its modifications with time lag (Kendall et al., 1999; Turchin, 2003; Turchin et al., 2003; Jonzén et al., 2005; Bechtol and Kruse, 2009; Sadykova and Nedorezov, 2013; Nedorezov and Sadykova, 2015). In particular, the Moran–Ricker model with and without time lag was applied to the description and analysis of the population dynamics of red king crab off Kodiak Island (Bechtol and Kruse,

G.P. Neverova et al. / Ecological Modelling 340 (2016) 64–73

2009). Also this model with time lag has been applied to describe and analyze the dynamics of the Larch budmoth (Zeiraphera diniana Gn.) (Sadykova and Nedorezov, 2013; Nedorezov and Sadykova, 2015). As a rule, the focus of applied studies was on the estimation of qualitative and quantitative indicators that characterized the quality of correspondence between the model and the empirical datasets. Many studies (e.g. May and Oster, 1976; Skaletskaya et al., 1979) are devoted to detailed theoretical investigation of dynamical modes of the Moran-Ricker model without time lag. The central point of these studies is the conditions of stability of non-trivial equilibrium. It has been shown that stability loss occurs with increasing values of reproductive potential, and a further complication of the dynamics happens according to the Feigenbaum scenario. In particular conditions of appearance and stability of the two-year and three-year cycles for this model were obtained (Skaletskaya et al., 1979; Nedorezov, 2013). One of the most important theoretical results, namely the proof of uniqueness of the stable attractor in a Moran-Ricker model without time lag, was first obtained by M.V. Jacobson (Jacobson, 1976). The attractor type is determined by the model parameters and does not depend on the initial conditions (initial value of population size) (Skaletskaya et al., 1979). Note that intraspecific and interspecific competition is typically modelled by the Moran-Ricker model function without time lag (Maunder, 1997; Todd et al., 2004; Golinski et al., 2008; Frisman et al., 2011, 2015a). The models of local populations with a simple age structure, in which density dependent regulation is described by an exponential function, like the Moran-Ricker equation without time lag, have recently been actively investigated (Frisman et al., 2011, 2015a, 2016; Zhdanova and Frisman, 2016). These models reveal the phenomenon of multimodality, or the possibility for different stable dynamic modes (stable, periodic, chaotic) to exist under the same conditions, their transition determined by initial population size. This implies that the random variation of the current population size may lead to a change of the observed dynamic mode. These results allow to suggest that the Moran-Ricker model with a time lag should have a complex dynamic behavior. However, no theoretical study of this model with time lag has been ever published. To the best of our knowledge, no study has addressed problems of multistability and multimodality regarding this model. We set out to fill these apparent gaps. Investigation of dynamic modes of the model aims to identify and analyze the effects that may occur in natural populations due to the influence of delayed density dependent regulation. Thus, the goal of this article is detailed study of the dynamic modes of the Moran-Ricker model with time lag. We apply our results to description of insect population dynamics. The results of our analytical and numerical analysis of the model are studied in detail and compared with real population census. The general form of the Moran-Ricker model with the time lag is represented by the following equation: i=m 

xn+1 = axn exp(−

bi xn−i )

(1)

i=0

where хn is the population size at time moment n, n is reproductive season number, and a is reproductive potential. The factor of i=m 

exp(−

bi xn−i )

i=0

characterizes the environmental limitation of population growth, and m is the time lag value corresponding to the number of years during which the necessary available resources constrain population development.

65

With m = 0 Eq. (1) is reduced to the basic well-studied MoranRicker model (Moran, 1950; Ricker, 1954; Skaletskaya et al., 1979): xn+1 = axn exp(−b0 xn ) · Parameter b0 is the scaling coefficient, which determines the intensity of the density-dependent limitation and indirectly characterizes the carrying capacity because if xn = 1/b0 , then the number of the next generation reaches the maximum population size of the species. Thus, each generation “receives” the same volume of necessary available resources regardless of the previous population size. This may be true for species whose carrying capacity fully recovers between breeding seasons. Otherwise, the decrease in resources owing to their consumption by previous generations should be considered. Eq. (1) considers the previous generations’ impact on population dynamics. This study examines the features of dynamic behavior of model (1) for the time lag values of 1 and 2. 2. Dynamic modes of the model with one year time lag (m = 1) The Moran–Ricker model with time lag 1 can be expressed as follows: xn+1 = axn exp(−b0 xn − b1 xn−1 ) ·

(2)

Coefficient b1 characterizes the reduction of required resources owing to their consumption by the previous generation. The decrease in required resources can be interpreted as a reduction of the reproductive potential of the current population generation. This process corresponds to delayed density dependence. Parameter b0 characterizes the intensity of the density-dependent limitation under the constraint of resources per individual. The oneyear time lag allows resources to be fully recovered during two reproductive seasons. Introduction of a new variable yn = xn−1 reduces Model (2) to a system of two recurrence equations without time lag with three parameters:



xn+1 = axn exp(−b0 xn − b1 yn ) yn+1 = xn



·

(3)

or, after substitutions b0 xn → xn and b0 yn → yn , to system xn+1 = axn exp(−xn −  · yn ), yn+1 = xn

(4)

where  = b1 /b0 · Parameter  characterizes the previous generation “participation” in the density-dependent regulation of the population birth rate. No “participation” means  = 0, and system (2) is reduced to the classical Moran–Ricker model. If  = b1 /b0 < 1 (b1 < b0 ), then the previous generation “participation” in the density-dependent regulation of the population birth rate does not exceed the “participation” of the current generation. Thus, the resources used by a population are significantly recovered between two successive periods of reproduction. The situation is inverse if = b1 /b0 > 1 (b1 > b0 ). The previous generation “participation” in the density-dependent regulation of the population birth rate is more than the “participation” of the current generation. The higher the value of parameter  ( > 1), the fewer resources are available for the next generation. System (4) has only one non-trivial invariant point x¯ = y¯ =

1 ln a, a > 0 · 1+

(5)

The solution increases monotonically with the increasing values of parameter a (for any constant value of ␳) but decreases monoton-

66

G.P. Neverova et al. / Ecological Modelling 340 (2016) 64–73

Fig. 1. Dynamic mode maps depending on the initial condition: a) initial condition х0 = y0 = 10, b) initial condition that belongs to the attraction basin of fixed point, and c) initial condition that belongs to the attraction basin of stable periodic points (3-cycle). Figures correspond to the period of observed cycles. ID stands for irregular dynamics.

Fig. 2. Dynamic mode maps with the initial condition х0 = y0 = z0 = 0.1 and different values of parameter ; Figures correspond to the period of observed cycles; ID is irregular dynamics.

G.P. Neverova et al. / Ecological Modelling 340 (2016) 64–73

67

Fig. 3. Application of the Moran–Ricker model with one year time lag to insect population dynamics a) P. flammea, b) B. piniaria, and c) Z. griseana (Sadykova and Nedorezov, 2013; Nedorezov and Sadykova, 2015).

Fig. 4. a) Dynamic modes map for the Moran–Ricker model with 1-year time lag with the point estimations for real population dynamics. Figures correspond to the period of observed cycles; ID is irregular dynamics; b) A fragment of the dynamic modes map is shown by the rectangle; c) Attraction basins of the model (4).

Fig. 5. Dynamic modes maps corresponding to the point estimations of the following species: a) B. piniaria, b) P. flammea, and c) Z. griseana: The estimations are from (Sadykova and Nedorezov, 2013; Nedorezov and Sadykova, 2015).

ically with the increasing values of parameter ␳ (for any constant value of a). The stability domain of the solution (5) is found based on standard method. The stability boundaries and stability domain of the solution (5) are presented in Appendix A. The following intervals can be identified for the values of parameter ␳ based on the scenario of stability loss by the non-zero solution for system (4). 1. ␳ < 1/3. The stability loss of the fixed point may happen only with the transition of one of the eigenvalues through −1 and is accompanied by a cascade of period-doubling bifurcations. Further, with increasing parameter a, dynamic mode changes according to the Feigenbaum scenario, and eventually chaos arises. This scenario is identical to the Moran-Ricker model without time lag when stability loss happens according to the Feigenbaum scenario (Skaletskaya et al., 1979). For the Moran-Ricker model period

doubling bifurcations arise when the reproductive potential value is equal to e2 . The equilibrium population size (fixed point (5)) decreases with the growth of parameter ␳, but its stability domain significantly increases. The fluctuations of population size occur at high values of reproductive potential. 2. ␳ > 1/3. Stability loss happens according to the Neimark–Sacker bifurcation scenario. The population dynamics transform into quasi-periodic fluctuations. The growth of parameter ␳ value leads to the noticeable contraction of the stability domain of population equilibrium size. The complex oscillations of population size occur at low reproductive potential. The stability domain of nontrivial solution of system (4) at ␳ ≤ 1 is not less than the stability area of the Moran–Ricker model without time lag. The transition to quasi-periodic dynamics at ␳ > 1 occurs if the value of the reproductive potential is less than e2 . The size

68

G.P. Neverova et al. / Ecological Modelling 340 (2016) 64–73

Fig. 6. The Moran–Ricker model with 2-year time lag applied to the insect population dynamics: a) E. tedella and b) Z. griseana (Schwerdtfeger, 1939, 1952; Suhovolskiy and Tarasova, 2014).

Fig. 7. Dynamic modes maps with the by point estimates. The upper row corresponds to Epinotia. tedella, and the bottom row corresponds to Z. griseana (Schwerdtfeger, 1939, 1952; Suhovolskiy and Tarasova, 2014).

fluctuations of population with delayed density dependence may be observed at low values of reproductive potential if the previous generation “participation” in the density-dependent regulation of the population birth rate begins to exceed the current generation “participation”. Codimension-two bifurcation (resonant 1:2) is found at ␳ = 1/3. Dynamic mode maps (Fig. 1) are used to analyze the possible dynamical modes of model (4). Maps are generated as follows: the limit cycle period of system (4) is calculated at each point in the plane of parameters (␳, a), and this point is painted in a specific color according to the found period. Area of irregular dynamics is denoted as ID in Fig. 1. It includes both quasi-periodic and chaotic dynamics. If stability loss occurs via the Feigenbaum scenario, then the following sequence of dynamical modes shifts 2 → 4 → 8 → 16 → 32 → .... → chaos is observed. If stability loss occurs via the Neimark-Saker scenario, then an attracting invariant curve is formed in the phase space, and a transition to quasi-periodic dynamics is observed. With increasing parameter values, stable resonant cycles emerge into an irregular dynamics area and may evolve in agreement with the NeimarkSaker scenario or with the Feigenbaum scenario. Fig. 1 shows a zone of parameters ␳ and a values with the simultaneous existence of both a stable stationary point and a stable

three-cycle (periodic point). Whether population shifts to equilibrium or to oscillations with period 3 depends on the initial condition (Fig. 1). The simultaneous existence of several stable attractors, namely, the invariant and periodic points, in the stability domain of non-trivial invariant point was originally revealed in the analysis of Henon’s model (Saucedo-Solorio et al., 2002; Pisarchik and Feudel, 2014). The cycle with period 3 emerges in the stability domain of non-trivial stationary point as a result of tangent bifurcation. The Moran-Ricker model without time lag has a cycle with period 3 in the zone of irregular dynamics. Fig. 1b and c show the green area corresponding to three-cycles that adjoins the axis of ordinates. The growth of the parameter ␳ leads to oscillations with period 3 at low values of reproductive potential. The cycles with periods 2 and 3 (2-cycle is a result of stability loss by non-zero invariant point through Feigenbaum scenario), cycle with period 3, and irregular oscillations in the instability domain (Fig. 1b and c) coexist. The particular kind of achieved dynamic mode depends on the initial condition. Therefore, the Moran–Ricker model with the time lag exhibits the phenomenon of multimodality. With the growth of parameter ␳ value, the domain of existence of a three-cycle and its bifurcations moves deeper into the stability area of fixed point. The dynamic mode changes with the further increase of parameter ␳ value into quasi-periodic dynamics. Simul-

G.P. Neverova et al. / Ecological Modelling 340 (2016) 64–73

Fig. 8. Stability domain of invariant point (5) of system (4).

taneous existence of the three limiting modes (i.e., invariant point, three-cycle, and four-cycle) is possible (Fig. 1c). Two different fourcycles emerge in the phase space of this model. One results from a period-doubling bifurcation, whereas the other is a result of tangent bifurcation. Our analyses show that quasi-periodic oscillations are observed if the present population size substantially depends on the population size of the previous year (i.e., recovery rate of required resources). The 2-cycles occur when the growth of population size is regulated by density dependence in the current year, and are possible within a narrow scope of demographic parameters. Thus, the more the reproduction of the population depends on changes in the ecological niche caused by the previous generation, the less predictable the population dynamics is. The stability domain is narrowed, and the growth of the reproductive potential leads to the emergence of quasi-periodic oscillations. However, periodic oscillations are possible if the parameter values correspond to Arnold’s tongues (i.e., resonances on the invariant curve in the area of quasi-periodic modes (Kuznetsov et al., 2012b)). 3. Dynamic modes for two years lag (m = 2) The Moran–Ricker model with time lag 2 (m = 2) is defined by the following equation: xn+1 = axn exp(−b0 xn − b1 xn−1 − b2 xn−2 ) ·

(6)

By analogy with the previous case, the change of variables yn = xn−1 , zn = xn−2 , b0 xn → xn , b0 yn → yn , and b0 zn → zn transforms the Model (6) into a system of three equations without any time lag:

⎧ x = axn exp(−xn − yn − zn ) ⎪ ⎨ n+1 ⎪ ⎩

yn+1 = xn

,

(7)

zn+1 = yn

where  = b1 /b0 and  = b2 /b0 · Model (7) has only one non-trivial fixed point x¯ = y¯ = z¯ =

1 ln a, a > 0 · 1++

Analysis of the stability of this solution, as well as a graph of its stability domain, are presented in Appendix B. The following value ranges of parameter  can be defined according to the type of stability domain deformation possible bifurcations.

69

1. 0 <  ≤ 0.25. There are two possible scenarios of stability loss when  = 0: the emergence of two-cycle when  is small, and the emergence and destruction of the invariant curve when  is large. Stability domain significantly expands due to the increase in parameter  value. Convergence of the population size to a stable fixed point is possible for a high value of the reproductive potential. 2. 0.25 <  ≤ 0.5 . Growth of parameter  leads to the contraction of the stability domain. The stability area for  = 0.5 is similar to one for  = 0. The types of bifurcations remain. 3. 0.5 <  < 1 . The stability domain continues to contract with increasing values of the parameter. The range of the parameter , where loss of stability through the cascade of period-doubling bifurcations is possible, narrows. There exists such value of the parameter  that the stability area is no less than the stability area for the Moran-Ricker model without time lag. 4.  ≥ 1 . The increase of both parameters  and  leads to the narrowing of the stability area, which is less than the size of the stability domain of the basic Moran-Ricker model. Moreover, the stability loss of fixed point cannot happen by Feigenbaum scenario. The stability loss and emergence of quasi-periodic oscillations can be observed at low values of reproductive potential.

Fig. 2 illustrates that the increase in the coefficient  values leads to a non-monotonic change in the area of stability domain of the non-trivial fixed point. The expansion is followed by a quick contraction. The three-cycle (i.e., the result of a tangent bifurcation) occurs at high values of reproductive potential. A three-year cycle (i.e., periodic points) that coexists with the non-zero fixed point is replaced by cycles with periods 4 and 5. With the increasing value of parameter , four and five-cycles occur with higher values of reproductive potential. Fig. 2 shows that model (7) has two different cycles with period 4 as well as the model with m = 1. Both four-cycles can coexist at the same values of demographic parameters (i.e.,  = 0.1 and  = 0.2, Fig. 2). Special attention should be given to the map of dynamic modes located at the center of Fig. 2 ( = 0.5). Here, we observe effect described in (Kuznetsov et al., 2012a): a small part of the boundary, which corresponds to stability loss in accordance with the Feigenbaum scenario, leads to the hard transition through the multiplier −1, rather than to period doubling bifurcation. Given the further increase of the parameter  ( ≥ 1), the stability domain of the non-trivial invariant point decreases, and fluctuations are possible only in the area of irregular dynamics. The area of stability region is substantially less than for the classical Moran-Ricker model.

4. Application of the Moran–Ricker model with time lag to insect population dynamics As mentioned before, the Moran–Ricker model with the time lag has been successfully used to describe the population dynamics of several insect species (e.g. Sadykova and Nedorezov, 2013; Nedorezov and Sadykova, 2015). In these studies researchers have focused on the quality of correspondence between model trajectory and empirical dynamics, while omitting detailed analysis of possible behavior and stability. We next apply the results of our analysis to models of real insect populations. We used the time series on insect dynamics that can be freely downloaded online from the Global Population Dynamics Database (GPDD). We applied the Moran–Ricker model with different time lags to describe the dynamics of the Larch budmoth (Z. Griseana) (Schwerdtfeger, 1939; Schwerdtfeger 1952; Suhovolskiy and Tarasova, 2014), Spruce needleminer (Epinotia tedella) (Data set

70

G.P. Neverova et al. / Ecological Modelling 340 (2016) 64–73

Fig. 9. Stability domain of system (7) fixed point in parameter space (a, ) with different values of .

20631), Pine beauty (Panolis flammea) (Data set 10065), and Pine looper (Bupalus piniaria) (Data set 2728). We also employed some results that apply the Moran–Ricker model to the description of Larch budmoth (Zeiraphera griseana) population dynamics (GPDD: Data set 1525). The results are presented in the papers of Sadykova and Nedorezov (2013) and Nedorezov and Sadykova (2015). Model (1) was used to describe the time series of population size for the species listed. Traditional methods (Bard, 1974; Kremer, 2004; Draper and Smith, 2014) were adopted to analyze the correspondence of empirical time series and model trajectories. For every considered modification of the Moran–Ricker model (m = 0, 2), we estimate values of model parameters and initial condition (x0 m − i , i = 0,¯m) minimizing sums of squared deviations between time series and model trajectories. The estimated model parameters following the hypotheses were checked with regard to the equivalence of the average of residuals to zero (Student t-test), the normality of deviations (Кolmogorov test), and the absence of serial correlation in sequences of residuals (Durbin–Watson test). Moreover, the average standard error, coefficient of determination, and adjusted coefficient of determination were calculated by considering the length of time series and the number of estimated model parameters. The Akaike information criterion was calculated for each version of the model (Akaike, 1974). Among all models with different time lags, preference was given to the model with the lowest Akaike criterion value. Table 1 presents the obtained estimates of the model parameters and initial population size for the model with the lowest Akaike criterion value. The estimates of statistical indicators and criteria show that the Moran–Ricker model with time lag satisfactorily describes the insects’ population dynamics. The behavior of real and model trajectories for different values of time lag is presented in Fig. 3a and b. The trajectory of Moran-Ricker model without time lag describes transition to a stable state that does not correspond to the real dynamics. However, this same model with time lag 1 describes the trend population dynamics and catches the key outbreaks. The descriptions of the actual dynamics by the models with time lags 1 and 2 are close. Due to the value of the Akaike information criterion of model with time lag 1 being less than for model with time lag 2, preference was given to the model with time lag 1. It should be

noted that, if a model with time lag m satisfactorily describes real data, then quality descriptions of real data by the model can either improve or leave it unchanged with further increase of time lag. Fig. 3 is supplemented by the results of Larch budmoth (Zeiraphera griseana) population dynamics (GPDD: Data set 1525), which were reported by Sadykova and Nedorezov (2013) and Nedorezov and Sadykova (2015) (Fig. 3c). The following estimates for the Moran-Ricker model with 1 year time lag were previously obtained for the population of Larch budmoth (Sadykova and Nedorezov, 2013; Nedorezov and Sadykova, 2015): a = 8.2996, b0 = 2.205 · 10−3 , b1 = 0.02214, and  = b1 /b0 ≈ 10.05. The coefficients’ values correspond to an asymptotic mode close to a cycle with period 9. This result shows the good fit of the time series (Sadykova and Nedorezov, 2013; Nedorezov and Sadykova, 2015). The model parameter estimations show that b1 value is greater than b0 (parameter  = b1 /b0 is greater than 1). The previous generation “participation” in the density-dependent regulation of the population birth rate is significantly greater than the “participation” of the current generation. Fig. 4 presents a dynamic modes map for the model (4) with m = 1 in the parameter space (a, ). The point estimates of parameters for the insect species (including the point estimate from (Sadykova and Nedorezov, 2013; Nedorezov and Sadykova, 2015)) are also presented. Fig. 4a shows that the point estimations for P. flammea and B. piniaria are located in the area of quasi-periodic dynamics. The type of dynamic mode does not change with small fluctuations in the values of parameters a and ␳. The transition from quasi-periodic to periodic dynamics is possible because the parameter values may fall into one of Arnold’s tongues. However, this scenario is unlikely because of the narrowness of the tongues. The point estimation of the parameters for Z. griseana (GPDD: Data set 1525) is in Arnold’s tongues of quasi-periodic dynamics area and corresponds to oscillations with period 18. The 18-cycle is the result of the period doubling bifurcation of the nine-cycle. The doubled elements of the cycle are close to one another. The point estimation is placed near the border that separates the cycles with periods 9 and 18. A small fluctuation in the values of the parameters may lead to a change in the dynamic mode.

G.P. Neverova et al. / Ecological Modelling 340 (2016) 64–73

71

Table 1 Values of statistical indicators and criteria characterizing the correspondence of empirical time series and model trajectories. Species

Pine beauty (Panolis flammea)

Pine looper or Bordered white (Buupalus piniaria)

Larch budmoth, (Zeiraphera griseana)

Spruce needleminer, (Epinotia tedella)

Data source

GPDD: Data set 10065:

GPDD: Data set 2728

GPDD: Data set 20631

Time lag that corresponds to the lowest Akaike criterion value Estimations of model parameters

m=1

m=1

(Schwerdtfeger, 1939, 1952; Suhovolskiy and Tarasova, 2014) m=2

m=2

a = 5.1, b0 = 2.35, b1 = 7.52, x0 = 1.2 · 10−3 , y0 = 0.06

a = 3.809, b0 = 0.015, b1 = 0.147, x0 = 20.27, y0 = 9.835

 = 3.2

 = 9.8

0.852 0.79

0.829 0.7

a = 6.28, b0 = 7.14 · 10−3 , b1 = 9.2 · 10−4 , b2 = 0.022, x0 = 311, y0 = 85, z0 = 27  = 0.129  = 3.082 0.948 0.927

a = 6.86, b0 = 2.42 · 10−4 , b1 = 9.3 · 10−6 , b2 = 0.0017, x0 = 928, y0 = 652, z0 = 718  = 0.038  = 6.964 0.98 0.967

0.159 < 2.86 t < t␣ = 0.01 (accepted) 0.8 < 1.63  <  ␣=0.01 (accepted) d* < d < d* 1.27 < 2.1 < 2.73 (accepted)

0.134 < 3.012 t < t␣ = 0.01 (accepted) 0.54 < 1.63  <  ␣ = 0.01 (accepted) d* < d < d* 1.25 < 1.875 < 2.75 (accepted)

0.8 < 2.77 t < t␣ = 0.01 (accepted) 1.59 < 1.63  <  ␣ = 0.01 (accepted) d* < d < d* 1.41 < 2.2 < 2.59 (accepted)

0.218 < 2.9 t < t␣= 0.01 (accepted) 0.472 < 1.63 ␭ <  ␣ = 0.01 (accepted) d* < d < d* 1.44 < 1.85 < 2.56 (accepted)

Values of parameters  and  Determination coefficient Adjusted determination coefficient Student t-test with 1% significance level

Кolmogorov test with 1% significance level Durbin–Watson test with 1% significance level

The dynamic modes map shows that Arnold’s tongues are superimposed on one another with large values of reproductive potential and low values of the coefficient ␳. The enlarged fragment (Fig. 4b) demonstrates stripes of nine-cycle on top of ten-cycle and its subsequent bifurcation (i.e., cycles with periods 20 and 40, etc., and irregular dynamics). Nine-cycle coexists with the above dynamic modes. Therefore, the realization of cycle with period 9 or companion mode (i.e., 10-cycle or results of its bifurcations) depends on the initial conditions. Fig. 6c shows the attraction basins that demonstrate the coexistence of cycles with periods 9 and 10. The type of dynamic mode determines both population parameters’ values and initial population size if the point estimation that corresponds to the real dynamics is located in the multistability area. A small fluctuation of the current population size may be considered as a modification of initial conditions caused by random influence of external factors. Therefore, if the point estimation that corresponds to the real dynamics is located in the multistability area, then the type of dynamic mode determines both population parameter values and current population size. Additional maps of dynamic modes in the parameter space (Fig. 5) were constructed for the calculated point estimation. Fig. 5a and b show that the point estimations of P. flammea and B. piniaria are located in the area of quasi-periodic dynamics. This dynamic mode is not disturbed when slight fluctuations are observed in the parameter values. The map of dynamic modes of Z. griseana population model contains many different coterminous dynamic modes. Thus, dynamic modes are sensitive to the variation of values of coefficients characterizing the effect of density-dependent factors on the reproduction process. The behavior of the real and model trajectories that correspond to the dynamics of E. tedella and Z. griseana (Schwerdtfeger, 1939, 1952; Suhovolskiy and Tarasova, 2014) are shown in Fig. 6. The trajectory of the Moran-Ricker model without time lag for E. tedella describes the transition to a stable equilibrium that does not correspond to real dynamics. The trajectory of the Moran-Ricker model without time lag for Z. griseana corresponds to irregular fluctuations weakly associated with the dynamics of real population. The

best results were obtained for the model with 2-year time lag. The trajectories corresponding to this model describe the main features of the population dynamics for these species. The values of statistical indicators and criteria that correspond to obtained parameter estimations (Table 1) show the Moran–Ricker model with 2-year time lag satisfactorily describes the dynamics of E. tedella and Z. griseana. The model parameter estimations show that b2 is significantly greater than b0 and b1 . The “participation” of the generation born two years ago in the density-dependent regulation of the population birth rate is greater than the “participations” of the generations of current and the previous years. The results of point estimations are supplemented by the dynamic mode maps in parametric spaces (a, ), (b1 , b0 ), (b2 , b0 ) and (b2 , b1 ) (Fig. 7). The point estimation of the demographic parameters of the Epinotia tedella population is located in quasi-periodic dynamics area. Moreover, almost the entire region is permeated by a large number of areas of cycles with different periods. The variation of the parameters b2 and b0 can lead to a dynamic mode shift between quasi-periodic dynamics and a whole range of “Arnold’s tongues.” Consequently, a real dynamics trajectory can be a set of transients from one mode to another. The point estimate of population parameters of Zeiraphera griseana that inhabit the Alps (Schwerdtfeger, 1939, 1952; Suhovolskiy and Tarasova, 2014) is located in the quasi-periodic dynamics area near “Arnold’s tongues,” which correspond to a 19year cycle. The elements of the 19-cycle are located on the invariant curve that corresponds to the quasi-periodic oscillations. This estimate is most sensitive to the value variation of parameter b2 , which characterizes the intensity of density dependent regulation with a two-year delay. Thus, the analysis of model dynamic modes at point estimation facilitates understanding of possible scenarios of a population development with the demographic parameter variation of a population.

72

G.P. Neverova et al. / Ecological Modelling 340 (2016) 64–73

5. Conclusions

 = e±i : a = exp

Accounting for delay factors due to previous generation contribution in environmental limit in birth rate leads to a significantly higher complexity of behavior of mathematical models for local population dynamics. The model shows the phenomenon of multimodality, which represents the possibility for different dynamic modes (i.e., equilibrium, periodic oscillations, or chaotic fluctuations) to exist under the same parameters’ values. The transition of these modes is determined by initial population size. Our analysis suggests that the dependence of the population dynamics on the recovery rate of available resources leads to quasi-periodic dynamics in most cases. Two-year fluctuations are possible in a narrow range of demographic parameter values. These fluctuations occur with limited population growth by density-dependent regulation in the current year. The Moran–Ricker model with time lag was applied to describe the dynamics of the Z. griseana, E. tedella, P. flammea, and B. piniaria. The obtained estimates of population parameters are statistically significant. The obtained point estimates corresponding to real population dynamics are located in a quasi-periodic fluctuations area that adjoins other dynamic modes. A variation of value demographic parameters (e.g., as a result of evolution or global climate change) can lead to a dynamic mode change. Moreover, the dynamic mode shift in the case of several dynamic modes existence may be due to current population size changes caused by the influence of modifying factors. Thus, to analyze and describe the dynamics of real populations, it is necessary to check calculated point estimate models for adequacy and compliance to the real object. Further investigation of possible model dynamic modes near this point makes it possible to identify dynamic modes between which transition may occur and to describe their possible transformations. Acknowledgement



The boundaries of stability area of the fixed point of system (7) are found based on the characteristic equation of the system (7) 3 − S2 + H − J = 0 (S,H and J are three invariants of the Jacobian matrix) and correspond to specific bifurcation transitions (Kuznetsov, Sedova, 2012): 1. Transcritical bifurcation,  = 1 (the dominant root is a real and is equal to 1): H = S + J − 1; 2. Period-doubling bifurcation,  = −1 (the dominant root is a real and equal to −1): H = −S − J − 1; 3. Neimark-Sacker bifurcation,  = e±i (the real root modulus is less than 1, and the absolute values of the complex conjugate roots are equal to 1): H = SJ − J 2 + 1. Characteristic polynomial coefficients of system (7) can be easily expressed through the model (7) parameters: ln a − 1 −  −   ln a  ln a ,H= J=− · 1++ 1++ 1++

Therefore, the boundaries of the stability domain of the stationary point of model (7) are defined by the following conditions.  = 1 : a = 1,

The standard method of finding the stability domain is based on the following theorem. The solutions of the characteristic equation 2 − S + J = 0 belong to the circle || < 1 if and only if

 = −1 : a = exp −2

(8)

(11)

Appendix B. Stability domain of model (7)

Appendix A. Stability domain of model (4)

|S| − 1 < J < 1,

·

The stability domain of the invariant point (5) in parameter space (, a) is depicted in Fig. 8. Boundary (9) separates the stability domains of non-trivial and trivial invariant points. If a < 1, then system (4) has a unique zero equilibrium point, which is globally stable (i.e., population degenerates). The trivial (zero) solution loses its stability, and a new stable non-zero (nontrivial) solution (5) appears if parameters change and the line a = 1 is crossed over when a < 1. Thus, the system (4) shows the population extinction as well as the classical model of Moran-Ricker for a <1 and any value of ␳. Boundaries (10) and (11) have an intersection point in the plane of parameters (, a). The point coordinates are (1/3, e4 ). The abscissa =1/3 separates two scenarios of stability loss.

S=−

This work was partially supported by the Russian Foundation for Basic Research (project no. 15-29-02658 ofi m).

 + 1



 = e±i : a = exp

1++ −−1

(12) ;

(13)



2 + 2 + 5 2 − 4) 1 (1 +  + )(− −  ± 2 2  −

(14)

where S and J are trace and determinant of the Jacobian matrix (Shapiro and Luppov, 1983; Kuznetsov and Sedova, 2012). Inequalities (8) are defined in the plane (S, J) as a “triangle of stability.” The boundaries of these inequalities are given by the following lines.

Stability area is formed by surfaces (12)–(14). Stability domain of system (11) fixed point in parameter space (a, ) with different values of  is depicted in Fig. 9.

1) 1 − S + J = 0, where one of the eigenvalues ∗ is equal to 1; 2) 1 + S + J = 0, where one of the eigenvalues ∗ is equal to −1; 3) J = 1, where the eigenvalues are complex numbers 1 2 = 1. The segment (−2 < J < 2) limits the “stability triangle.” The eigenvalues also conjugate the following: 1 = exp(i) and 2 = exp(−i).

References

Thus, S = 1 − In a/(1 + ), J =  · lna/1 +  therefore, the boundaries of the stability domain of the stationary point (5) for the model (4) are defined by the following conditions.  = 1 : a = 1,  = −1 : a = exp

 2( + 1) 1−

(9) ,

(10)

Akaike, H., 1974. A new look at the statistical model identification Automatic Control. IEEE Trans. Autom. Control 19 (6), 716–723. Bard, Y., 1974. Nonlinear Parameter Estimation. Academic Press Inc., New York. Bechtol, W.R., Kruse, G.H., 2009. Analysis of a stock-Recruit relationship for red king crab off kodiak island, Alaska. Mar. Coast. Fisher.: Dyn. Manag. Ecosyst. Sci. 1, 29––44. Berryman, A.A., Turchin, P., 2001. Identifying the density-dependent structure underlying ecological time series. Oikos 92, 265–270. Dennis, B., Taper, M.L., 1994. Density dependence in time series observations of natural populations: estimation and testing. Ecol. Monogr. 64 (2), 205–224. Draper, N.R., Smith, H., 2014. Applied Regression Analysis, 3rd ed. Wiley and Sons Inc., New York. Frisman, E.Ya., Neverova, G.P., Revutskaya, O.L., 2011. Complex dynamics of the population with a simple age structure. Ecol. Modell. 222, 1943–1950. Frisman, E.Ya., Neverova, G.P., Kulakov, M.P., Zhigalskii, O.A., 2015a. Multimode phenomenon in the population dynamics of animals with short live cycles. Dokl. Biol. Sci. 460, 42–47.

G.P. Neverova et al. / Ecological Modelling 340 (2016) 64–73 Frisman, E.Ya., Revutskaya, O.L., Neverova, G.P., 2015b. Basic trends of game mammal population dynamics in the Russian Middle Amur river area: the observation and simulation results. Sib. J. For. Sci. 3, 103–114 (in Russian). Frisman, E.Y., Neverova, G.P., Kulako, M.P., 2016. Change of dynamic regimes in the population of species with short life cycles: results of an analytical and numerical study. Ecol. Complex., http://dx.doi.org/10.1016/j.ecocom.2016.02. 001. Golinski, M., Bauch, C., Anand, M., 2008. The effects of endogenous ecological memory on population stability and resilience in a variable environment. Ecol. Modell. 212 (3), 334–341. Isaev, A.S., Hlebopros, R.G., Nedorezov, L.V., Kiselev, V.V., Kondakov, Yu.P., Suhovolskiy, V.G., 2001. Populyatsionnaya Dinamika Lesnyih Nasekomyih. Nauka, Moskva (in Russian). Jacobson, M.V., 1976. Properties of parameter family of dynamical systems x→A·x·exp(-x). Successes Math. Sci. 31 (2), 239–240 (in Russian). Jonzén, N., Pople, A.R., Grigg, G.C., Possingham, H.P., 2005. Of sheep and rain: large-scale population dynamics of the red kangaroo. J. Anim. Ecol. 74 (1), 22–30. Kendall, B.E., Briggs, C.J., Murdoch, W.W., Turchin, P., Ellner, S.P., McCauley, E., Nisbet, R.M., Wood, S.N., 1999. Why do population cycle? A synthesis of statistical and mechanistic modeling approaches. Ecology 80, 1789–1805. Kostitzin, V.A., 1937. La Biologie Mathematique. A. Colin, Paris. Kremer, N.S., 2004. Probability Theory and Mathematical Statistics. Publishing house Unity- Dana, Moscow (in Russian). Kuznetsov, A.P., Sedova, J.V., 2012. Bifurcations of three- and four-dimensional maps: universal properties. Rus. J. Nonlin. Dyn. 20 (5), 461–471 (in Russian). Kuznetsov, A.P., Kuznetsov, S.P., Pozdnyakov, M.V., Sedova, J.V., 2012a. Universal two-dimensional map and its radiophysical realization. Rus. J. Nonlin. Dyn. 8 (3), 461–471 (in Russian). Kuznetsov, A.P., Savin, A.V., Sedova, Yu.V., Tyuryukina, L.V., 2012b. Bifurkatsii otobrazheniy. Izdatelskiy tsentr (Nauka), Saratov (in Russian). Lotka, A.J., 1925. Elements of Physical Biology. Williams and Wilkins Company, Baltimore, Maryland. Maunder, M.N., 1997. Investigation of density dependence in salmon spawner–egg relationships using queueing theory. Ecol. Modell. 104 (2), 189–197. May, R.M., Oster, G.F., 1976. Bifurcations and dynamic complexity in simple ecological models. Am. Nat. 110 (974), 573–599. May, R.M., 1974. Biological populations with non-overlapping generations: stable points, stable cycles, and chaos. Science 186, 645–647. May, R.M., 1975. Biological population obeying difference equations: stable points, stable cycles, and chaos. J. Theor. Biol. 51 (2), 511–524. Moran, P.A.P., 1950. Some remarks on animal population dynamics. Biometrics 6 (3), 250–258. Myers, R.A., Bowen, K.G., Barrowman, N.J., 1999. Maximum reproductive rate of fish at low population sizes. Can. J. Fish. Aquat.Sci. 56 (12), 2404––2419. Nedorezov, L.V., Sadykova, D.L., 2008. Green oak leaf roller moth dynamics: an application of discrete time mathematical models. Ecol. Modell. 212 (1), 162–170. Nedorezov, L.V., Sadykova, D.L., 2015. Dynamics of larch bud moth populations: application of Moran-Ricker models with time lag. Ecol. Modell. 297, 26–32. Nedorezov, L.V., 2010. Analysis of pine looper population dynamics using discrete time mathematical models. Math. Biol. Bioinf. 5 (2), 114–123 (in Russian).

73

Nedorezov, L.V., 2013. About an approach to population periodic dynamics analysis (on an example of larch bud moth fluctuations). Popul. Dyn.: Anal. Modell. Forecast. 2 (1), 23–37. Pisarchik, A.N., Feudel, U., 2014. Control of multistability. Phys. Rep. 540, 167––218. Prout, T., McChesney, F., 1985. Competition among immatures affects their adult fertility: population dynamics. Am. Nat. 126 (4), 521–558. Ricker, W.E., 1954. Stock and recruitment. J. Fish. Res board Can. 11 (5), 559–623. Sadykova, D.L., Nedorezov, L.V., 2013. Larch bud moth dynamics: can we explain periodicity of population fluctuations by the time lag dependence in birth rate? Popul. Dyn.: Anal. Modell. Forecast. 2 (4), 154––181. Saucedo-Solorio, J.M., Pisarchik, A.N., Aboites, V., 2002. Shift of critical points in the parametrically modulated Hénon map with coexisting attractors. Phys. Lett. A304, 21–29. Schwerdtfeger, F., 1939. Uber Die Kritische Eizahl Und Eiparasitierung Beim Kiefernspanner, Bupalus Piniarius L. MFF, Hannover. Schwerdtfeger, F., 1952. Untersuchunger uber der eisen bestand fon kiefernspannern (Bupalus piniarius L.), forleule (Panolis flammea schiff.) und kiefernswarmer (Hyloicus pinastri L.). zeitschrift furangew. Entomology 34 (2), 216–283. Shapiro, A.P., Luppov, S.P., 1983. Rekurrentnye Uravneniya V Teorii Populyatsionnoy Biologii. Nauka, Moskva (in Russian). Shapiro, A.P., 1972. K voprosu o tsiklah v vozvratnyih posledovatelnostyah. Upravlenie i informatsiya, Vip. 3. DVNTs AN SSSR, Vladivostok (in Russian). Skaletskaya, E.I., Frisman, E.Ya., Shapiro, A.P., 1979. Diskretnyie modeli dinamiki chislennosti populyatsiy i optimizatsiya promyisla. Nauka, Moskva (in Russian). Suhovolskiy, V.G., Tarasova, O.V., 2014. Prognozyi chislennosti nasekomyih-vrediteley, riski povrezhdeniya derevev i prinyatie resheniy v zadachah lesozaschityi//Mezhdunarodnyiy seminar «Matematicheskie modeli v teoreticheskoy ekologii i zemledelii» (Poluektovskie chteniya), posvyaschennyiy pamyati professora Ratmira Aleksandrovicha Poluektova. Tezisyi dokladov, 99–102 (in Russian). Todd, C.R., Nicol, S.J., Koehn, J.D., 2004. Density-dependence uncertainty in population models for the conservation management of trout cod, Maccullochella macquariensis. Ecol. Modell. 171 (4), 359–380. Turchin, P., Wood, S.N., Ellner, S.P., Kendall, B.E., Murdoch, W.W., Fischlin, A., Casas, J., McCauley, E., Briggs, C.J., 2003. Dynamical effects of plant quality and parasitism on population cycles of larch budmoth. Ecology 84 (5), 1207––1214. Turchin, P., 1990. Rarity of density dependence or population regulation with lags? Nature 344, 660–663. Turchin, P., 2003. Complex Population Dynamics: A Theoretical/Empirical Synthesis. Princeton University Press, Princeton. Volterra, V., Lecons sur la theorie mathematique de la lutte pour la vie, 1931, Paris. Williams, D.W., Liebhold, A., 1995. Detection of delayed density dependence: effects of autocorrelation in an exogenous factor. Ecology 76 (3), 1005–1008. Zhdanova, O.L., Frisman, E.Ya., 2016. Manifestation of multimodality in a simple ecological-Genetic model of population evolution. Rus. J. Genet. 52 (8), 870–878.