# An Analysis of a Semelparous Population Model with Density-Dependent Fecundity and Density-Dependent Survival Probabilities

Authors: Arild Wikan
Journal: Journal of Applied Mathematics
 (2017)

---

## Abstract

A discrete age-structured semelparous Leslie matrix model where density dependence is included both in the fecundity and in the survival probabilities is analysed. Depending on strength of density dependence, we show in the precocious semelparous case that the nonstationary dynamics may indeed be rich, ranging from SYC (a dynamical state where the whole population is in one age class only) dynamics to cycles of low period where all age classes are populated. Quasiperiodic and chaotic dynamics have also been identified. Moreover, outside parameter regions where SYC dynamics dominates, we prove that the transfer from stability to instability goes through a supercritical Neimark−Sacker bifurcation, and it is further shown that when the population switches from possessing a precocious to a delayed semelparous life history both stability properties and the possibility of periodic dynamics become weaker.


---

## Body

## 1. Introduction

Within the framework of nonlinear discrete age-structured population models, the dynamical properties and behaviour of a great variety of species may be explored. Such species may possess an iteroparous life history which means that individuals in several age classes of the population are fertile, or they may possess a semelparous life history which is characterized by the property that only individuals of the last age class are fertile.Regarding iteroparity, most population models (both age-structured and stage-structured) focus on cases where nonlinearities (density dependence) are built into the fecundity elements and not into the survivals. Particularly in fishery models, this has often been motivated by the assumption that most density effects are present only in the first year of life. Examples of theoretical studies which deal with nonstationary and chaotic behaviour as well as behaviour linked to concrete species may be found in [1–7]. In case of ergodic properties we refer to [8–10]; see also [11]. Another strategy is to assume constant fecundity terms and nonlinear year-to-year survival probabilities, compare [12–16]. Depending on functional form of the nonlinear terms, a crude conclusion found in several of the papers referred to above is that, in case of small population densities, the models possess a stable nontrivial equilibrium where all age classes are populated. At higher population densities there may be nonstationary, periodic, and chaotic dynamics of stunning complexity.Now, turning to the semelparous case with nonlinear fecundity element and constant survival probabilities a rather peculiar phenomenon has been detected, compare [4]. Here, in contrast to the iteroparous case, the nontrivial equilibrium tends to be unstable in large parameter regions, also in case of low population densities. Instead one finds, as time t→∞, that a cyclic state is attained where the whole population is in just one single age class at each time step. Such behaviour, which is called synchronization or SYC (single year class) dynamics, has been detected among insects; see [17]. Hoppensteadt and Keller [18] presented a model for the 17-year cicada (magicicada) which included both predation and intraspecific competition and in [19] cicada dynamics is further explored. Regarding biennials and possible SYC dynamics we refer to [20]. In [21] SYC phenomena and related MYC (multiple year class) dynamics are considered while Kon [22] discusses in a general context conditions for SYC dynamics to occur in matrix models, also compare [23]. However, if we take the opposite approach, constant fecundity term, and nonlinear year-to-year survivals probabilities, SYC dynamics appear to be a rare event; see [16]. Instead, the nontrivial equilibrium is stable whenever the population size is sufficiently small and the nonstationary dynamics has a strong resemblance of 4-cycles, either exact or approximate, the chaotic regime included.In contrast to most of the papers quoted above, the purpose of this paper is to study the combined effect of nonlinear recruitment and nonlinear survival probabilities in semelparous population models. Several general results about possible SYC dynamics and stability properties of the nontrivial equilibrium in such problems may be obtained in [23] (two age classes) and [24] (three age classes). In case of the more general setting where an arbitrary number of age classes is considered we refer to [25, 26]. Regarding our study, we shall restrict the analysis to the case where both the fecundity and the survival elements depend on the total population and, moreover, the functional form of these elements is of Ricker type. As we prove, under such restrictions, it is possible to obtain explicit thresholds for secondary bifurcations of flip and Neimark−Sacker type and we may also prove (in some cases) that bifurcations involved are of supercritical nature. We may also investigate how the dynamics reported earlier will change as more density-dependent terms are included. For example, given SYC cycles, will the cycles persist if the strength of density dependence in the survival terms is included? If not, what kind of qualitative dynamics is it then possible to obtain? Assuming cycles of low period where all age classes are populated, does the inclusion of density dependence in the fecundity terms act in a stabilizing or destabilizing fashion? Such questions are difficult to address in the general models presented in [24–26].The paper is organized in the following way. In Section2 we present the model, compute equilibria, and derive the nth order eigenvalue equation which we need in order to perform stability and bifurcation analysis. This is followed (Sections 3 and 4) by a rigorous analysis of possible dynamic outcomes in two and three age class models, respectively. Finally, in Section 5, we unify and discuss results when the number of age classes exceeds three.


## 2. The Model, Fixed Points, and Stability

First we establish the model. At timet we split the population xt into n distinct nonoverlapping age classes x-t=(x1,…,xn)T where the total population is given by x=x1+⋯+xn. Next, we introduce the transition matrix(1)A=00⋯0fp1000p20⋱00pn-10,where f is the average fecundity of a member of the last age class at time t. pi denotes the survival probabilities, the (year-to-year) survival from age class i-1 to i. In contrast to most papers the assumption here is that both the fecundity and the survivals are nonlinear terms. Thus, we write f as f=Fexp⁡(-αx) where the constant F>1 and the survivals pi as pi=Pexp⁡(-βix) where the constant P satisfies 0<P≤1. The parameters α≥0, βi≥0 may be regarded as parameters that measure the strength of density dependence. The relation between x- at two consecutive time steps is then expressed as(2a)x-t+1=Ax-twhich may also be formulated as a nonlinear map of the form(2b)f:Rn⟶Rnx-⟶Ax-.Besides the trivial fixed point x~=(0,…,0)T maps (2a) and (2b) also possess a unique nontrivial point x-∗=(x1∗,…,xn∗)T. The latter may be expressed as(3)x1∗,…,xi∗,…,xn∗=1Dx∗,…,Pi-1aβ1+⋯+βi-1Dx∗,…,Pn-1aβ1+⋯+βn-1Dx∗.The quantities a and D are defined as(4)a=R0-1/α+β1+⋯+βn-1,D=1+Paβ1+P2aβ1+β2+⋯+Pn-1aβ1+⋯+βn-1,where R0=FPn-1 and R0>1 is assumed throughout the paper in order to have a feasible equilibrium. The total equilibrium population x∗ is given as(5)x∗=1α+∑i=1n-1βiln⁡R0.In order to investigate stability we linearize about the fixed point. This gives birth to the nth order eigenvalue equation(6)λn+a1λn-1+a2λn-2+⋯+an-1λ+an=0,where the coefficients a1,…,an satisfy(7)a1a2a3⋮an-1an=αβ1β2βn-2βn-1βn-1αβ1βn-3βn-2βn-2βn-1αβn-4βn-3⋮⋱⋮β2β3β4αβ1β1β2β3βn-1αx1∗x2∗x3∗⋮xn-1∗xn∗+000⋮0-1.x-∗ is a locally stable hyperbolic fixed point as long as all eigenvalues λ of (6) are located inside the unit circle in the complex plane.There are three ways in whichx-∗ may fail to be stable. It may lose its hyperbolicity when λ crosses the unit circle through 1 which in the general case leads to a saddle node bifurcation, alternatively through −1 which gives birth to a flip (period doubling) bifurcation, or it may fail to be hyperbolic as a pair of complex-valued eigenvalues cross the unit circle. Then a Neimark−Sacker bifurcation occurs. The Jury criteria, see the book by Murray [27], provide conditions for all eigenvalues to satisfy |λ|<1.


## 3. Two Age Classes

Letn=2 in maps (2a) and (2b). Then we have(8)x1,x2⟶Fe-αxx2,Pe-β1xx1.R0=FP, R0>1, and x∗=(α+β1)-1ln⁡R0. Moreover, the fixed point becomes(9)x1∗,x2∗=11+Paβ1x∗,Paβ11+Paβ1x∗and the eigenvalue equation may be cast in the form(10)λ2+a1λ+a2=0,where a1=αx1∗+β1x2∗ and a2=β1x1∗+αx2∗-1.Fixed point (9) is stable whenever the Jury criteria 1+a1+a2>0, 1-a1+a2>0, and 1-|a2|>0 hold, that is, as long as(11a)α+β1x∗>0,(11b)β1-α1-Paβ1x1∗>0,(11c)2-β1x1∗+αx2∗>0,respectively.There are two cases to consider: (A) the caseα>β1, which means that the strength of density dependence in the fecundity is stronger than the strength of density dependence in the survival, and (B) the case β1>α.Considering (A), it is clear from (11b) that there does not exist any stable fixed point. Moreover, since (11b) fails as an eigenvalue crosses the unit circle through −1 it is natural to search for a 2-cycle which should be stable provided x∗ is small. Evidently, such a 2-cycle must be obtained from(12)x1,t+2=R0e-β1x1,t+x2,te-αx1,t+1+x2.t+1x1,tx2,t+2=R0e-αx1,t+x2,te-β1x1,t+1+x2,t+1x2,tand here there are two possibilities:(1)
xt=xt+1 which leads to the trivial 2-cycle where the unstable fixed point x1∗,x2∗ is the only point the cycle.(2)
The points are on the form(A,0) or (0,B). In this case it follows from (12) that A and B must satisfy the equations(13a)A-1Pe-β1AB=0(13b)β1A-αB+ln⁡R0=0and by finding B from (13b) and substitute back into (13a) we arrive at(13c)αPeβ1A=β1+ln⁡R0A.Geometrically, it is now easy to see that the graph of the left hand side of (13c) and that of the right side have a unique intersection point lying in the first (positive) quadrant. In the special case β1=0 we obtain A=(αP)-1ln⁡R0 and B=α-1ln⁡R0. Hence, there exists a 2-cycle on the form (SYC form)(14)1αPln⁡R0,0,0,1αln⁡R0and as shown in [4] this cycle is stable in case of F small (P fixed). In Figure 1(a) we show an orbit starting at (x10,x20)≠(0,0) which settles on the 2-cycle (14). If we continue to increase F, we find that (14) goes unstable and cycles of period 2k. k=2,3,…, are established through successive flip bifurcations. These cycles, which are all on SYC form, are stable in smaller and smaller regions as F is increased. Eventually, the dynamics becomes chaotic but we emphasize that it is on SYC form also in the chaotic regime. These scenarios are demonstrated in Figures 1(b) and 1(c). Actually, we have not accounted for what happens when R0=1 (or x∗=0). Here, compare (10), the eigenvalues are 1 and −1, respectively, and both the positive equilibrium and the SYC 2-cycle bifurcate forward. For proofs and details we refer to [4, 23–26].Figure 1
(a) An orbit approaching the SYC 2-cycle(P-1ln⁡R0.0), (0,ln⁡R0). Parameter values (α,β)=(1,0), (P,F)=(0.5,10). (b) A 4-cycle on SYC form. (α,β)=(1,0), (P,F)=(0.5,15). (c) Chaotic SYC dynamics. (α,β)=(1,0), (P,F)=(0.5,25).
(a)
(b)
(c)When0<β1<α we observe much of the same SYC dynamics as in the β1=0 case. However, we may in a sense argue that an increase of β1 acts in a stabilizing fashion. Indeed, if β1=0 and F=25, map (8) generates chaotic dynamics. On the other hand, if β1=0.8 and F=25 (P=0.5 in both cases) the outcome is a stable period 4-cycle on SYC form as shown in Figure 2.Figure 2
4-periodic SYC dynamics.(α,β)=(1,0.8) and (P,F)=(0.5,25).Since the boundary of the positive cone is always invariant for semelparous Leslie matrix models of any dimension, see [26], initial conditions of the form (v,0) or (0,w) in the 2-dimensional case will always produce SYC dynamics. However, if β1→α and (x10,x2,0)≠(v,0) or (0,w) the dynamics occurs in the vicinity (mostly as a stable 2-cycle, not on SYC form) of the unstable fixed point(15)x1∗,x2∗=FF+Px∗,PF+Px∗,where x∗=(2α)-1ln⁡R0. Hence, 0<β1<α does not necessarily imply SYC dynamics although it is the most likely outcome. A further discussion is postponed to Section 5.Next, consider (B) (the case0<α<β1). Our first observation is that the Jury criteria (11a) and (11b) will never be violated. Moreover, in case of sufficiently small equilibrium populations x1∗,x2∗ the left hand side of (11c) will be positive as well. Consequently, there exists a region in parameter space where (9) is a hyperbolic stable fixed point. However, if we increase the value of x∗ (by increasing F) such that (11c) turns into an equality, (x1∗,x2∗) undergoes a Neimark−Sacker bifurcation, loses its hyperbolicity, and becomes unstable at the threshold(16)ln⁡R0=α+β121+Paβ1β1+αPaβ1while the corresponding modulus 1 solutions of the eigenvalue equation (10) may be cast in the form(17)λ=-α+β1Paβ1β1+αPaβ1±β12-α21-P2a2β1β1+αPaβ1i.As is well known, bifurcations may be of both supercritical and subcritical nature. If a fixed point shall undergo a supercritical bifurcation it means that an eigenvalue (pair of eigenvalues) λ must cross the unit circle outwards at instability and in the Neimark−Sacker case that an attracting quasiperiodic orbit restricted to an invariant curve is created beyond instability threshold. Now, considering our bifurcation, we first express map (8), using the abbreviation b=β12-α21-P2a2β1, as(18)uv⟶-α+β1Paβ1β1+αPaβ1-bβ1+αPaβ1bβ1+αPaβ1-α+β1Paβ1β1+αPaβ1uv+fu,vgu,v,where f(u,v) and g(u,v) contain second- and third-order terms of u and v. (Details of how (18) is derived and explicit formulas of f(u,v) and g(u,v) may be found in Appendix.)Following Wan [28], the bifurcation will be supercritical if(19)γ=-Re⁡1-2λλ¯21-λξ11ξ20-12ξ112-ξ022+Re⁡λ¯ξ21is negative and that dλ/dF>0 at bifurcation. The quantities in (19) are defined as(20)ξ11=14fuu+fvv+iguu+gvvξ20=18fuu-fvv+2guv+iguu-gvv-2fuvξ02=18fuu-fvv-2guv+iguu-gvv+2fuvξ21=116fuuu+fuvv+guuv+gvvv+iguuu+guvv-fuuv-fvvv.Due to the complexity of f(u,v) and g(u,v) it is out of reach to compute the sign of γ in the most general case. However, it works in the important special case α=0. Indeed, (18) simplifies to(21)uv⟶-1F-bFbF-1Fuv+fu,vgu,v,where b=β1F2-1 and(22)fu,v=β1buv-12Fb2uv2+1bF2b2β1v3,gu,v=β1bFfu,v.Moreover, ξ11 becomes zero so γ may be found as(23)γ=-ξ022+Re⁡λ-ξ21=-β1216β12+F2-1which is clearly negative.Finally, at thresholdx2∗=2/F, compare (11c) or (16)(24)dλdF=F+22FF+1>0.Hence, we have proved that the Neimark−Sacker bifurcation is of supercritical nature.Still, assumingα=0, let us now focus on the dynamics on the invariant curve. Whenever F>FT and |F-FT| is small where FT is the value of F at threshold (for a given value of P, FT is found from P=(FT)-1exp⁡2(1+FT)/FT), we find a quasiperiodic orbit which fills the invariant curve with points; see Figure 3(a). In case of larger values of F we observe a stable 4-cycle, comapre Figure 3(b), and as we continue to increase F stable orbits of period 4·2k, k=1,2,…, are the outcomes. (Note, however, that the points in these orbits are clustered in such a way that one from an observational point of view may argue that we have almost 4-periodic dynamics in these cases as well.) The smaller P is, the larger the interval (Fmin,Fmax) becomes where the 4-periodic structure occurs. Eventually, in case of even larger fecundity values, the dynamics becomes chaotic, but even in the chaotic regime a certain kind of 4-periodicity is preserved in the sense that the chaotic attractor is divided into 4 disjoint subsets that are visited once every fourth iteration; see Figure 3(c).Figure 3
(a) An invariant curve.(α,β)=(0,1) and (P,F)=(0.5,20). (b) 4-periodic dynamics generated by map (8). (α,β)=(0,1) and (P,F)=(0.5,25). (c) Chaotic dynamics generated by (8). (α,β)=(0,1) and (P,F)=(0.5,32).
(a)
(b)
(c)The reason why we have all this 4-periodicity may be understood along the following line. Once the invariant curve is established we may regard map (8) (α=0) restricted to the curve as topological equivalent to a circle map. Associated with a circle map there is a rotation number σ, and whenever σ is rational we have an orbit of finite period. If σ is an irrational number the orbit is quasiperiodic. Now, at bifurcation threshold the modulus 1 eigenvalues are(25)λ=-1F±1FF2-1iand since F is large at threshold, λ is located close to the imaginary axis in the left half plane. Following Guckenheimer and Holmes [29], arg⁡λ gives asymptotic information about the rotation number; thus σ≈1/4. This accounts for the 4-periodicity observed. Finally, note that when an orbit becomes exact periodic as a result of changing parameter F, the implicit function theorem guarantees that there exists an interval around that specific parameter value where the periodicity is preserved as well.Next, consider the case0<α<β. The location of the modulus 1 eigenvalues is now given by (17) and they are not so close to the imaginary axis anymore. For “small” values of α we find the invariant curve, but the 4-periodic dynamics reported above (α=0) is absent. Instead, as a result of increasing F, the invariant curve becomes kinked and not topological equivalent to a circle. We may also describe the dynamics by use of the Lyapunov exponent L of the orbit generated by (8) [30, 31]. In Figure 4(a) we display the values of L in the F range 10<F<50 when (α,β)=(0.1,1.0) and P=0.5. Whenever 10<F<20.341, L<0. In this interval the fixed point is stable. When 20.341≤F<29.3 the dynamics is quasiperiodic and restricted to invariant attracting curves and the corresponding Lyapunov exponent is L=0. L<0 in the parameter window 30.5<F<32.7. Periodic dynamics of period 11 is the outcome here. Finally, when F exceeds 32.7 and also in a tiny interval just below 30.5 we find that L>0 which means that the dynamics is chaotic. These findings are also visualized in Figure 4(b) where the dynamics is shown for selected values of F. From bottom to top we recognize stable fixed points, invariant curves, kinked curves, 11-periodic dynamics, and chaos.Figure 4
Values of Lyapunov exponents and the dynamics generated by (8) when P=0.5. Parameter values: (a) and (b) (α,β)=(0.1,1.0), (c) and (d) (α,β)=(0.5,1.0), and (e) and (f) (α,β)=(0.8,1.0). (L<0 corresponds to stable periodic orbits, L=0 corresponds to quasiperiodic orbits restricted to attracting invariant curves, and L>0 corresponds to chaotic dynamics.)
(a)
(b)
(c)
(d)
(e)
(f)In case of intermediate values ofα, (α,β)=(0.5,10) we observe a significant change of dynamics. Here, there exists a critical F value FC, FC<FT, where the third iterate of map (8) undergoes a saddle node bifurcation which results in two large amplitude 3-cycles, one stable and one unstable. Thus, in the interval FC<F<FT we find coexistence between two stable attractors, the stable fixed point (9) and the stable 3-cycle. Consequently, the ultimate fate of an orbit depends on the initial condition but since the trapping region of the 3-cycle appears to be (much) larger, 3-periodic dynamics is the most likely outcome. Beyond FT there is an interval FT<F<FK where the invariant curve established at FT and the stable 3-cycle coexists. At F=FK the invariant curve disappears as it is hit by the unstable 3-cycle, and if F>FK, |F-FK| small, the 3-cycle is the only attractor. Similar phenomena have also been found in iteroparous population models, first by Guckenheimer et al. [1], later in [4]. At even higher fecundity values the stable 3-cycle undergoes a Neimark−Sacker bifurcation which leads to three invariant curves which are visited once every third iteration. This is followed by kinked curves, periodic dynamics, and chaos through further enlargement of F. Values of L as well as the dynamics reported above are shown in Figures 4(c) and 4(d).Ifα is large but α<β, we find that the dynamics qualitatively is quite similar to the case where α is small, compare Figures 4(e) and 4(f) where (α,β)=(0.8,1.0) and P=0,5. However, there are several parameter values where the dynamics is periodic but most of these have almost no widths. The exception is the first window where the dynamics is 5-periodic. A final comment is that the larger α is, (α<β), the higher F is at bifurcation threshold (16), so one may argue that the strength of density dependence in the fecundity (α<β) acts in a stabilizing way.


## 4. Three Age Classes

In this section we study the map(26)f:R3⟶R3x1,x2,x3⟶Fe-αxx3,Pe-β1xx1,Pe-β2xx2.The nontrivial fixed point (cf. (3)) may be expressed as(27)x1∗,x2∗,x3∗=1Dx∗,Paβ1Dx∗,P2aβ1+β2Dx∗,where a=R0-(1/(α+β1+β2)), D=1+Paβ1+P2aβ1+β2, R0=FP2, R0>1, and x∗=(α+β1+β2)-1ln⁡R0.From (6) it follows that the coefficient in the eigenvalue equation(28)λ3+a1λ2+a2λ+a3=0may be written as a1=αx1∗+β1x2∗+β2x3∗, a2=αx2∗+β2x1∗+β1x3∗, and a3=αx3∗+β1x1∗+β2x2∗-1.Equilibrium (27) is locally asymptotic stable as long as the Jury criteria 1+a1+a2+a3>0, 1-a1+a2-a3>0, 1-a3>0, and 1-a32-a2-a2a1>0 hold. Hence, the stability criteria are(29a)α+β1+β2x∗>0(29b)2-α+β1-β2x1∗-β1+β2-αx2∗-α-β1+β2x3∗>0(29c)2-β1x1∗+β2x2∗+αx3∗>0(29d)2β1-α-β2x1∗+2β2-α-β1x2∗+2α-β1-β2x3∗+αβ1-β12x1∗2+β1β2-β22x2∗2+αβ2-α2x3∗2+β12+αβ2-2β1β2x1∗x2∗+α2+β1β2-2αβ1x1∗x3∗+αβ1+β22-2αβ2x2∗x3∗>0.A final observation from (28) is that when x∗=0 all eigenvalues λ=exp⁡(2kπ/3)i, k=0,1,2, are located on the boundary of the unit circle. Thus, if all eigenvalues move inside the unit circle when x∗ is increased, then (27) is a hyperbolic stable equilibrium in case of x∗ small. If at least one eigenvalue moves out, (27) is unstable. The Jury criteria (29a), (29b), (29c), and (29d) will help us to decide. Whenever (29c) and (29d) fail, a pair of complex valued eigenvalues will have modulus larger (or equal to) than 1 and consequently have a location outside (or on) the unit circle. If (29a) and (29b) fail, an eigenvalue crosses the unit circle through 1 or −1, respectively.Moreover, due to the complexity of the Jury criteria we limit the discussion in the rest of this section toβ1=β2=β. Then the Jury criteria may be expressed as(30a)α+2βx∗>0(30b)2-αx1∗-x2∗+x3∗-2βx2∗>0(30c)2-αx3∗-βx1∗+x2∗>0(30d)β-αx1∗-x3∗1-αx3∗+βx1∗+x2∗+x2∗-x3∗>0.Criteria (30a), (30b), and (30c) obviously hold in case of small equilibrium populations x∗. Regarding (30d) the same is true only as long as β>α.(A) First, considerβ<α. Then (30d) fails whenever x∗ is small but it holds for those x∗ which satisfy(31a)x∗>1+Paβ-2P2a2βD1-P2a2βαP2a2β+β1+Paβ=u.On the other hand, (30b) holds only as long as(31b)x∗<2Dα1-Paβ+P2a2β+2βPaβ=v.Therefore, if x∗ satisfies the inequality u<x∗<v, then (x1∗,x2∗,x3∗) is stable. Moreover, if this shall be the case, the difference u-v must be strictly negative which is equivalent to say that(32)α-2β+2β-α10β-13α23α-2β<Paβ<1.The most transparent case to discuss is β=0. Then u=(1+P-2P2)(1+P+P2)(P2-P4)-1, v=2(1+P+P2)(1-P+P2)-1, and from (32) we find t1=(1+13)/6<P<1 in order to ensure that u-v<0. Now, if we choose a P value larger than t1, say P=0.9, we find u=4.93, v=5.96 and finally if 170.92<F<476.64 we obtain u<x∗<v. Hence, there exists a parameter region where the fixed point is stable.Turning to the dynamics, SYC dynamics is the only outcome in the region where (30d) fails. In case of small F values there exists a stable 3-cycle where the points in the cycle are (α-1P-2x∗,0,0), (0,α-1P-1x∗,0), and (0,0,α-1x∗). For higher fecundity values there are cycles of period 3·2k as well as chaotic dynamics. In the intermediate region where all criteria (30a)–(30d) are valid, we find coexistence between the fixed point and chaotic SYC dynamics. Whether an orbit shall converge towards x1∗,x2∗,x3∗ or settle on a chaotic attractor on SYC form depends on the initial condition. If we continue to increase F we find that (30b) fails and as a result the fixed point will undergo a (supercritical) flip bifurcation. This gives birth to a tiny parameter region where SYC dynamics coexists with a stable 2-cycle not on SYC form. Through further enlargement of F, SYC dynamics is the only outcome.In [4] it was proved that there also exists a 3-cycle where two age classes are populated at each time. Within our framework the points in such a cycle are found to be (α-1P-1(1-P2)z∗,α-1P(1-P)z∗,0), (0,α-1(1-P2)z∗,α-1P2(1-P)z∗), and (α-1(1-P)z∗,0,α-1P(1-P2)z∗) where z∗=(1-P2)-1x∗. In [21] it is proved that the cycle is unstable and further, in a more general context, it is shown in [24] that such 2-class cycles may be embedded in a heteroclinic cycle which may be attracting. Indeed, depending of the kind of competition within and between age classes, the analysis in [24] accounts for possible loops and cycles. For example, if the competition is asymmetric the bifurcating invariant loop consists of a single year class 3-cycle with a synchronous two-year class orbit that heteroclinically connect the phases of this 3-cycle. Thus, all orbits on such a heteroclinic cycle approach the single year class 3-cycle. This mechanism was also found in Bulmer’s original work on periodical cicadas. For further reading, compare [17, 24]. A final point is that when P=1, the eigenvalue equation may be cast in the form(33)λ3+13ln⁡R0λ2+13ln⁡R0λ+13ln⁡R0-1=0,where the dominant roots are exp⁡±(2π/3)i. This means that independent of population size there will always be 3-cyclic dynamics in the sense that the total population at every point in the cycle equals x∗ but the structure of the points is on SYC form.Turning to the case0<β<α, it is still possible to find combinations of P and F such that the fixed point is stable but the parameter interval where this occurs appears to be (very) small. Indeed, assuming β=α/100, the choice Paβ=0.8 satisfies (32) and if P=0.845 it is possible to find F such that Paβ=0.8 as well as u<x∗<v. However, if we change to P=0.84 or 0.85 it does not work. Therefore, our conclusion is that SYC dynamics dominates almost completely in the region 0<β<α.(B) Next, considerβ>α. In this case all criteria (30a)–(30d) hold whenever x∗ is sufficiently small. Thus, there exists a parameter region where (27) is stable. Instability is introduced through a Neimark−Sacker bifurcation when x∗=u (cf. (31a)) or equivalently when(34)P=1Fexp⁡α+2β1+Paβ-2P2a2β1+Paβ+P2a2β1-P2a2βαP2a2β+β1+Paβ.Regarding the eigenvalue equation (28), coefficients a1 and a2 are always positive and at threshold (34), a3=Paβ(1+Paβ)-1 which is positive too. Consequently, there are no changes of signs between the coefficients in (28) which implies that there are two complex modulus 1 solutions and one real negative solution which necessarily must be λ=-Paβ(1-Paβ)-1 and clearly -1<λ<0. These findings allow us to express (28) as(35)λ2+α1+Paβ-P2a2β+2βP2a2βαP2a2β+β1+Paβλ+1λ+Paβ1+Paβ=0and the location of eigenvalue at threshold is given by the complex roots.Now, scrutinizing the special caseα=0, β>0, we find that instability threshold (34) becomes(36)P=1Fexp⁡2F+F+1F+F-2F+FF-1and moreover that the complex modulus 1 solutions of (35) may be expressed as(37)λ=-F-1FF-1±1-F-12FF-12i.From (36) it follows that F is large at bifurcation threshold; thus the eigenvalue (37) is located in the left half plane close to the imaginary axis (in fact even closer than in the corresponding 2-age class model); therefore arg⁡λ≈π/2 here too. In order to visualize the dynamics we have studied the case P=0.7 in somewhat more detail. In Figure 5(a), where F=24, we show a quasiperiodic orbit which is restricted to a invariant curve. When F=30, see Figure 5(b), an exact 4-period orbit is established through the frequency locking mechanism. In case of higher F values, the fourth iterate of (24) undergoes a Neimark−Sacker bifurcation which results in four invariant curves as displayed in Figure 5(c)  (F=39). Through further enlargement of F the dynamics becomes chaotic, see Figure 5(d) where F=43.Figure 5
(a) A quasiperiodic orbit which settles on an invariant curve.(α,β)=(0,1) and (P,F)=(0.7,24). (b) Convergence towards a 4-cycle. (α,β)=(0,1) and (P,F)=(0.7,30). (c) 4 invariant curves which are visited once every fourth iteration. (α,β)=(0,1) and (P,F)=(0.7,39). (d) Chaotic dynamics. (α,β)=(0,1) and (P,F)=(0.7,43). (e) Dynamics in the case (α,β)=(0.9,1) and (P,F)=(0.7,95).
(a)
(b)
(c)
(d)
(e)Next, assumeα>0. Whenever α is small, α≪β, we observe qualitatively the same kind of dynamics as just reported, that is, four periodicities either exact or approximate beyond instability threshold. Through further increase of α the tendency towards 4-periodic dynamics becomes less pronounced and gradually disappears completely. Moreover, we experience that the fixed point may be stable for much higher fecundity values than when α is small. The nonstationary dynamics is restricted to invariant curves, not topologically equivalent to circles, compare Figure 5(e), or chaotic. The curves appear to be weakly attracting.


## 5. Discussion

In the previous sections we have analysed different versions of two and three-dimensional maps. We shall now consider the more general situation where there aren age classes and since the parameter space is huge we will limit the discussion to the case where we have the same “strength” of survival between any two age classes. Hence, we let β=β1=⋯=βn-1 which means that the eigenvalue equation (6) may be cast in the form(38)λn+αx1∗+β∑i=2nxi∗λn-1+αx2∗+β∑i=1i≠2nxi∗λn-2+αx3∗+β∑i=1i≠3nxi∗λn-3+⋯+αxn-1∗+β∑i=1i≠n-1nxi∗λ+αxn∗+β∑i=1n-1xi∗-1=0.First, let us comment on the parameter region where α≥β. Whenever x∗ is small the fixed point is always unstable. This is proved in a general setting in [26] which also applies to the case under consideration in this paper. Next, suppose that n is even and λ=-1. Then left hand side of (38) may be expressed as x1∗-x2∗β-α+x3∗-x4∗β-α+⋯+xn-1∗-xn∗β-α which is clearly ≤0. On the other hand, when λ→-∞, then left hand side of (38) →+∞; thus there must be a root λ~<-1 of the equation which actually proves that the fixed point (x1∗,…,xn∗) will always be unstable. When n is odd the argument presented above does not hold, compare our analysis of the three age class model. However, it is possible to show that the fixed point is unstable given that x∗ is sufficiently large. To this end, assume that n is odd and λ=-1. Then left hand side of (38) is positive as long as(39)x∗>2∑j=0n-1Paβjα∑j=0n-1-1jPaβjand we observe that when λ→-∞, λn→-∞ too, so evidently there must exist an eigenvalue λ~<-1 where left hand side of (38) is zero. Hence, there is no stable fixed point in this case.As the number of age classes increases there may also be cycles where more than one age class is populated at each time. Such cycles, referred to as multiple year class, MYC, cycles have already been identified in then=3 case (β=0). In [21] where no density dependence in the survivals is assumed, MYC dynamics is further discussed and a major finding is that stability properties of such cycles on the whole are similar to the stability properties of (x1∗,…,xn∗) (i.e., unstable). Consequently, assuming α>β, we conclude that the nontrivial fixed point (3) will always be unstable except for small parameter windows when n is odd. As demonstrated, such windows may be hard to find when 0<β<α. Regarding the dynamics, the general case is SYC dynamics. When x∗ is small, there is an attracting n cycle on SYC form, as x∗ is increased we observe cycles of period 2kn, k=1,2,…, as well as chaotic dynamics. In parameter windows where (3) is stable there is coexistence between the fixed point and chaotic dynamics on SYC form. Considering the simplest case (α,β)=(1,0) the n cycle may be written as(40)P-n-1x∗,0,…,0,…,0,…,0,P-i-1x∗,0,…,0,…,0,…,x∗.An increase of β(β<α) acts in a stabilizing fashion in the sense that the higher β is, the higher F value is necessary in order to generate chaotic dynamics.From a biological point of view it is not obvious how one should interpret SYC dynamics. As far as we know, Bulmer [17] is the first who has noticed SYC dynamics in theoretical insect models. He explains its presence by saying that competition between age classes is more severe than competition within age classes. Bulmer’s argument is further strengthened by the findings in [24], see Theorem 9, where it is shown mathematically that strong competition within age classes in general will lead to an equilibrium where all age classes are populated (i.e., an equilibrium of form (3) in our model ) while strong competition between age classes destabilizes and promotes oscillations with missing age classes. Another argument is presented in [21]. However, it should be emphasized that the conclusion obtained there is found from a simple model.In the remaining part of the paper we shall deal exclusively with the caseβ>α. Denoting left hand side of (38) for P(λ) it is immediately clear from our previous analysis that when n is even the Jury criterion -1nP-1>0 will always be satisfied. Hence, there is no flip bifurcation in this case. When n is odd the value of x∗ at threshold -1nP-1=0 is found to be(41)x∗=2D2βPaβ∑j=0n-3/2Paβ2j+α∑j=0n-1-1jPaβjand if we substitute back into (38) the constant term becomes(42)an=2β∑j=0n-3/2Paβ2j-α∑j=0b-1-1jPaβjK,where K is the denumerator of (41). Moreover, since an>1 we may rule out the flip here too. Consequently, the only possible transfer from stability to instability goes through a (supercritical) Neimark−Sacker bifurcation for any n≥2.There are two more questions which we find natural to discuss. (A) Will the 4-periodic dynamics observed whenn=2 and n=3 also persist in case of larger values of n and (B) what happens to the size of the stable parameter region when the number of age classes increases?Regarding (A), assumen=4 and (α,β)=(0,1) (the values of α and β with the most pronounced 4-periodic dynamics when n=2 and n=3). From the Jury criteria (and a lengthy calculation!) it now follows that the size of x∗ at instability threshold must satisfy(43)D-P3D-P3P+D-1x∗2+D1+P2+4P4-2D1+2Px∗+4PD2=0.Thus, for given values of P we may compute x∗ from (43), substitute back into (38) and solve in order to find the location of the modulus 1 eigenvalues λ at threshold. For example, if P=0.6 we find λ1,2=0.111928±0.993716i  (λ3=0.203731, λ4=0.915659) while P=0.2 gives λ1,2=0.110914±0.99383i  (λ3=0.605291, λ4=0.916205). Now, the crucial observation from calculations as above is that λ1,2 (as opposed to the n=2 and n=3 cases) are located in the right half of the complex plane and not particularly close to the imaginary axis. Hence, arg⁡λ1,2-π/2 is not small, so the tendency towards 4-periodic dynamics is much weaker than in n=2, n=3 cases, a result which has also been verified through several simulations. Therefore, our analysis clearly suggests that periodic dynamics of low period, in particular 4-periodic dynamics, is limited to species with a few number of age classes. If only the first age class is not fertile and density effects occur in the survival terms exclusively, we may find 4-periodic dynamics also for larger values of n; see [15].Finally consider (B) and letx(2)∗(P), x(3)∗(P), and x(4)∗(P) be the size of x∗ at bifurcation threshold in the n=2,3 and 4 cases, respectively. From (16), (31a), (31b), and (43) it now follows that x(i)∗(P), i=2,3,4, are increasing functions of P, 0<P<1, and moreover, that 2<x(2)∗(P)<4, 1<x(3)∗(P)<9/4, and 0<x(4)∗(P)<4/3. These findings are presented in Figure 6. Hence, for a given value of P, x(4)∗(P)<x(3)∗(P)<x(2)∗(P) which means that the value of x∗ at threshold becomes smaller as the number of age classes is increased; that is, an increase of n acts as a destabilizing effect. This argument is further increased if we continue to increase n. Indeed, when n=5 we have verified numerically that the fixed point is unstable in case of P sufficiently small. Only quasiperiodic orbits have been identified. In [32] where a discrete stage-structured population model was analysed, the authors concluded that species who possess delayed semelparous life histories tend to be more stable than species who possess precocious semelparous life histories. Results from the analysis of age-structured models in [16] show a tendency in the opposite direction. The findings in this paper support the latter, but we feel that more work has to be done before one may conclude this issue.Figure 6
The curvesx(i)∗ at threshold. From bottom to top: x(4)∗, x(3)∗, and x(2)∗.



---

*Source: 8934295-2017-12-17.xml*