# Nonstationary and Chaotic Dynamics in Age-Structured Population Models

Authors: Arild Wikan; Ørjan Kristensen
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

The dynamics from nonlinear discrete age-structured population models is under consideration. Focus is on bifurcations, as well as nonstationary and chaotic dynamics. We also explore different mechanisms which may lead to periodic phenomena. Some new results are also presented, in particular from models where both fecundity and survival terms contain nonlinear elements.


---

## Body

## 1. Introduction

As it is well known (cf. [1–4]), simple one-dimensional maps of biological relevance can exhibit an extraordinary rich dynamical behaviour ranging from stable fixed points to chaotic oscillations. The review paper [5] provides an excellent summary of results from the papers quoted above.However, there are several important biological factors which are not possible to include in one-dimensional models. Therefore a more realistic approach is to apply age-structured models. In the discrete case, such matrix models (often referred to as Leslie matrix models) were independently developed in the 1940s (see [6–9]), but perhaps somewhat strange not adapted by ecologists until the 1970s (cf. [10] for a counterexample). The papers referred to above mainly consider linear models. In the 1970s and later, it has been more and more common to consider nonlinear (density dependent) terms too, as accounted for in the classical papers [11, 12], as well as in [13–20]. Some of these papers scrutinize the dynamics of concrete species, and others deal with theoretical aspects exclusively.Additionally it should be emphasized that matrix models may serve as basic tools for related problems like migration [21, 22], harvest [23–25], prey-predator systems [26–29], ergodicity [30, 31], and permanence [32]. Finally, we will also want to stress that several papers have been published recently where continuous systems of fractional order together with their discretizations have been analyzed and compared from different perspectives. Excellent examples may be obtained in [33–35]. The purpose of this paper is, by way of examples, to give a thorough description of different nonlinear phenomena which may occur in models quoted above. Focus is on bifurcations and nonstationary behaviour and parts of the content have a certain review component.The paper is structured as follows: in Section2, we present the model and analyze equilibria and stability. In Section 3, we provide several examples of nonstationary and chaotic dynamics while we in the last section state some concluding remarks.


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

First, we establish the model. At timet, we split the population x into n distinct nonoverlapping age classes x¯t=(x1,t,…,xn,t)T, where x=∑xi is the total population. Next, we introduce the Leslie matrix(1)A=f1f2⋯fnp10⋯00⋱⋱000pn-10,where fi is the average fecundity of a member of the ith age class at time t, while pi denotes the survival probability of age class i. Then the relation between x¯ at two consecutive time steps is given by the map(2)f:Rn⟶Rnx¯⟶Ax¯.The rationale behind model (2) is that sexual maturity is linked to age or that other properties than age are irrelevant; alternatively, if such relevant properties exist, they must be highly correlated with age.Map (2) covers species who possess a wide range of different life histories. Indeed, if fi>0 for all i≤n, we say that the species has an iteroparous life history and matrix A is classified as a nonnegative and primitive matrix. On the other hand, if fi=0,i<n and fn>0, the species possess a semelparous life history and in this case we say that A is an irreducible but imprimitive (cyclic) matrix. It is also customary to divide both the iteroparous and semelparous cases into two subclasses. If n is small and fi>0 for all i, one says that a species possesses a precocious iteroparous life history. Typical examples of species in this subclass are small mammals and small rodents. The delayed iteroparous subclass is characterized by fi=0,i<k (k≈n/2 “very roughly”) and fi>0,k≤i≤n. Humans and large mammals fall into this category. Regarding semelparity, whenever the number of age classes n is small, we classify the species as having a precocious semelparous life history. Examples may be found among annual plants, biennials, and insects. In the last subclass, delayed semelparity, we find species who may live for many years and then reproduce only once. Excellent examples are cicadas and several salmon species.If all matrix elements are constantfi=Fi, pi=Pi, model (2) is linear (constant terms are indicated by use of capital letters). In fishery models, it is often assumed that the survival probabilities are constants, pi=Pi, while the fecundities are density dependent. In the striped bass fishery model [12], it was assumed that fi=fi(y), where y=∑αixi is a weighted sum of age classes. Assuming two age classes only, cannibalism may be accounted for by the assumption p1=p1(x2) which suggests that individuals from the second age class prey upon members from the first age class. A third possibility is fi=fi(x) or pi=pi(x) which means that the elements in (1) depend on the total population x. Two of the most frequently used fecundity functions are the compensatory Beverton and Holt relation fi(x)=Fi(1+βx)-1; see [13, 36] and the overcompensatory Ricker relation, fi(x)=Fiexp⁡(-βx), first proposed by [37] and later applied in [11, 14, 16]. Note that both fecundity functions referred to above are members of the Deriso Schnute family fi(x)=Fi(1-γβx)1/γ (γ=-1 gives the Beverton and Holt relation while γ→0 results in the Ricker case). Also, observe that fi(x) is written as a product of a constant term Fi and a density dependent part, for example, exp⁡(-βx). A final comment is that fi′(x)≤0. Hence, except for Allee effects, a wide range of different biological scenarios may be investigated through model (2).Next, let us turn to the analysis. Assumingfi=fi(x) and (or) pi=pi(x), model (2) possesses a unique nontrivial fixed point (x1∗,…,xn∗) of hyperbolic type. x∗=x1∗+⋯+xn∗ is the total equilibrium population. Stability analysis is performed by linearizing about the fixed point which in turn results in an eigenvalue equation of degree n. The fixed point is locally asymptotically stable as long as all eigenvalues λ are located inside the unit circle in the complex plane. Provided the number of age classes n is small, the Jury criteria (see the books [38, 39]) give conditions for (x1∗,…,xn∗) to be stable. In cases where n is large, the criteria become complicated and not applicable unless in very special situations. It may be a severe problem to decide if all λi satisfy λi<1 when n is large. If we increase a parameter (F or P) such that an eigenvalue λ leaves the unit circle, the fixed point becomes unstable. The threshold where (x1∗,…,xn∗) loses its hyperbolicity (i.e., where λ=1) is termed instability or bifurcation threshold. Moreover, the location where an eigenvalue leaves the unit circle has crucial impact on the dynamics just beyond instability threshold. Considering the multidimensional model (2), the eigenvalues of the linearization may leave the unit circle through λ=1, λ=-1, or λ = ±exp⁡(iθ), where θ=2πk/n,k=0,1,…,n-1. If λ=1, the general case is that the fixed point will undergo a saddle node bifurcation at threshold. The other possibilities are the pitchfork and the transcritical bifurcations; for details, see [40]. λ=-1 results in a flip or periodic doubling bifurcation. Hence, assuming that the bifurcation is of supercritical nature, an attracting period 2 orbit is established when the fixed point goes unstable. When a pair of complex modulus 1 eigenvalues leave the unit circle, the fixed point will go through a Neimark Sacker (Hopf) bifurcation. Then, in the supercritical case, assuming we are outside the strongly resonant cases where λ is third or fourth root of unity, an invariant attracting closed curve is established around the unstable fixed point. Bifurcations of subcritical nature appear to be rare events in population models like (2), but an excellent example may be obtained in the prey-predator model discussed in [26].


## 3. Examples

### 3.1. Example 1 (The Linear Case)

Suppose thatfi=Fi,i=1,…,n, and pi=Pi,i=1,…,n-1. Then (2) is a linear map of the form(3)x¯⟶AIx¯,where AI is nonnegative and primitive. On the other hand, if fi=0,i<n,fn=Fn, and pi=Pi,i=1,…,n-1, map (2) becomes(4)x¯⟶ASx¯,where AS is an irreducible and imprimitive (cyclic) matrix.In order to reveal the dynamics generated by (3) and (4), we assume in both cases solutions of the form x¯t=λtw¯ and apply the Perron-Frobenius theorem which may be formulated as follows.Theorem 1 (Perron-Frobenius).
(1)
IfA is a positive or nonnegative and primitive, then there exists a real eigenvalue λ1>0 which is a simple root of the characteristic equation A-λI=0. Moreover, the eigenvalue is strictly greater than the magnitude of any other eigenvalue, λ1>λi for i≠1. The eigenvector w¯1 corresponding to λ1 is real and strictly positive. λ1 may not be the only positive eigenvalue, but if there are others, they do not have nonnegative eigenvectors.
(2)
IfA is irreducible but imprimitive (cyclic) with index of imprimitivity d+1, there exists a real eigenvalue λ1>0 which is a simple root of A-λI=0 with associated eigenvector w¯1>0¯. The eigenvalues λi satisfy λ1≥λi for i≠1, but there are d complex eigenvalues equal in magnitude to λ1 whose values are λ1exp⁡(2kπi/(d+1)),k=1,…,d.Proofs may be obtained in [39, 41, 42].Now, considering map (3), it follows directly from the assumption x¯t=λtw¯ and part 1 of the Perron-Frobenius theorem that x¯t may be expressed as(5)x¯t=c1λ1tw¯1+c2λ2tw¯2+⋯+cnλntw¯n.Here, ci are numbers, λ1 is real and positive, λ1>λi,i≠1, and the eigenvalues in (5) are numbered in order of decreasing magnitude. Moreover,(6)limt→∞x¯tλ1t=c1w¯1.Consequently, the long-term dynamics of the population is described by the growth rate λ1 and the stable population structure w¯1. Thus λ1>1 implies an exponential increasing population and 0≤λ1<1 an exponential decreasing population, where we in all cases have the stable age distribution w¯1. Turning to map (4), we find in a similar way by use of part 2 of the Perron-Frobenius theorem that(7)limt→∞x¯tλ1t=c1w¯1+∑k=2ncke2kπ/nitw¯k.Thus, opposed to the findings from (6), it now follows from (7) that w¯1 is not stable in the sense that an initial population not proportional to w¯1 will not converge to it. Instead the limit in (7) is periodic with period n. These are the most important cases when model (2) is linear. In Figure 1, we show the dynamics generated by maps (3) and (4) in the cases (8)AI=1120.30000.20,AS=0061/20001/30,respectively.Figure 1
(a) Dynamics generated by map (3). (b) Dynamics generated by map (4).
(a)
(b)The dominant eigenvalue (the growth rate) ofAI is λ=1.3, while the stable age distribution is given by the associated eigenvector w¯1=(0.97,0.22,0.03)T. The eigenvalues of AS are λ1=1, λ2,3=exp(±2πi/3), and AS3=I. Consequently, 3-cycle behaviour and no stable age distribution is the outcome.


### 3.2. Example 2 (Density Dependent Fecundity)

Assumefi=Fif^(x),i=1,…,n, where Fi is a constant larger than unity and f^′(x)≤0,pi=Pi,0<Pi<1,i=1,…,n-1. Then model (2) possesses a unique nontrivial fixed point(9)x1∗,…,xn∗=L1Kx∗,…,LnKx∗,where L1=1,Li=P1P2·…·Pi-1,K=∑i=1nLi, and x∗=f^-1(1/R). f^-1 denotes the inverse of the function f^ and R=∑i=1nLiFi. Following [43], the eigenvalue equation may be cast in the form(10)λn-∑i=1nLifi∗+θλn-i=0,where the negative parameter θ is expressed as(11)θ=f^′x∗x1∗∑i=1nLiFi.Note that (10) is valid for a wide range of fecundity functions including all members of the Deriso Schnute family.First, assumefi=Fexp(-x) (the Ricker case). Then x∗=∑FLi and (10) may be written as(12)λn-1x∑i=1nFLi1-x∗λn-i=0.Now, rewrite (12) as g(λ)+h(λ)=0, where g(λ)=λn. Clearly the equation g(λ)=0 has n roots located inside the unit circle C. On the boundary, we have (13)hλ=-1K1-x∗∑i=1nLiλn-i≤L1K1-x∗+⋯+LnK1-x∗≤1-x∗and whenever x∗<2,1-x∗<g(λ)=1. Thus according to Rouche’s theorem, the equations g(λ)+h(λ)=0 and g(λ)=0 have the same number of zeros located inside C. Hence x∗=ln⁡(FLi)<2 guarantees a stable fixed point.Next, assumefi=F(1+x)-1 (the Beverton and Holt case). Then x∗=∑FLi-1 and (10) takes the form(14)λn-1K11+x∗∑i=1nLiλn-1=0.Now, by use of exactly the same kind of arguments as in the Ricker case, we find that stability is achieved if x∗>-2 which of course is always satisfied. Consequently, by use of the Beverton and Holt recruitment function, the only possible dynamics is a stable fixed point. Regarding the Ricker case, the dynamics may be nonstationary and possibly chaotic if x∗ exceeds 2. We shall now explore this further.To this end, consider a two-dimensional version of (1) and (2), where f1=f2=Fexp(-x) and pi=P, that is, the map(15)x1,x2⟶Fe-xx1+Fe-xx2,Px1.The fixed point is easily found to be(16)x1∗,x2∗=11+Px∗,P1+Px∗and by use of the Jury criteria (cf. [38, 39]), it is straightforward to show that (16) is stable provided (17a)x∗<21+P0<P<12,(17b)x∗<1+2PP12<P<1.Assuming 0<P<1/2, the fixed point fails to be hyperbolic at the threshold x∗=2(1+P)-1 or equivalently when F=(1+P)-1exp⁡(2/(1+P)) and undergoes a supercritical flip bifurcation as proved in [16]. Hence, just beyond threshold an attracting period 2 orbit is established. Through further increase of x∗ (which is done by increasing F and fixing P) we observe stable orbits of period 2k,k=1,2,…. Eventually the dynamics becomes chaotic as displayed in Figure 2(a). Note that the attractor is divided into 4 disjoint subsets (branches) which are visited once every fourth iteration. In case of higher F values, these branches merge together.Figure 2
(a) Dynamics generated by map (15). Parameter values (P,F)=(0.2,58.5). (b) Coexisting attractors. (P,F)=(0.9,11.83). (c) Convergence towards the invariant curve. (P,F)=(0.9,11.83) and (x1,0,x2,0)=(2.0,1.78825). (d) Convergence towards the 3-cycle. (P,F)=(0.9,11.83) and (x1,0,x2,0)=(2.0,1.78830). (e) Chaotic dynamics generated by map (15). (P,F)=(0.9,18.5).
(a)
(b)
(c)
(d)
(e)Consider the latter case,1/2<P<1,(x1∗,x2∗) is the only stable attractor as long as x∗ is small, but when x∗ exceeds a certain threshold xA which is lower than bifurcation threshold xB given by (17b), the third iterate g of map (15) undergoes a saddle node bifurcation. Hence, two 3-cycles, one stable and one unstable, are created. This may be justified along the following line: First, note that map g (which has 7 fixed points) may be expressed as(18)x1,x2⟶Fe-FeFe-xx+Px1Fe-xx+Px1+PFe-xx·Fe-Fe-xx+Px1Fe-xx+Px1+PFe-xx,PFe-Fe-xx+Px1Fe-xx+Px1.Then we compute the linearization Dg(x¯1,x¯2), where (x¯1,x¯2) is a fixed point of (18) and find for fixed values of P the value of F where the 3-cycle attractor disappears. Next, we use this F and substitute a corresponding fixed point (x¯1,x¯2) into Dg(x¯1,x¯2) and compute the eigenvalues μ1 and μ2. The results of such calculations are μ1≈1 and 0<μ2<1 from which we conclude that a saddle node bifurcation takes place. For example, if P=0.9, we find F=10.036 which implies that μ1=1.0005 and μ2=0.48 and consequently that threshold xA may be expressed as xA=ln⁡(10.036·1.9)≈2.94. Thus there exists an interval xA<x∗<xB, where there are two stable attractors, the stable 3-cycle and the stable fixed point (16). At xB=P-1(1+2P) (xB=3.11 in case of P=0.9), the fixed point experiences a supercritical Neimark Sacker bifurcation. Consequently, just above threshold, we find coexistence between an attracting invariant curve where the dynamics is quasiperiodic and a stable large amplitude 3-cycle. This is displayed in Figure 2(b). In Figure 2(c), we see the behaviour of points which belong to the trapping region of the invariant curve, and in Figure 2(d) the initial point belongs to the trapping region of the 3-cycle. The scenario referred to above persists in an interval xB<x∗<xC. At xC (xC=3.143 when P=0.9), the invariant curve is hit by the branches of the unstable 3-cycle and disappears. Therefore, whenever x∗>xC and x∗-xC is small, the stable 3-cycle is the only attractor. If we continue to increase x∗ (by increasing F), successive flip bifurcations are the outcomes, creating orbits of period 3·2k. Eventually the dynamics becomes chaotic here too. This is shown in Figure 2(e). Thus to summarize, depending on parameter values, map (15) may generate dynamics of stunning complexity, all generic bifurcations (saddle node, flip, and Neimark Sacker) occur, and there are coexisting attractors as well as chaotic dynamics.


### 3.3. Example 3 (Density Dependent Survival)

Now suppose density dependent survival and consider the map(19)x1,x2⟶Fx2,P1-γβx1/γx1.The fixed point may be expressed as(20)x1∗,x2∗=F1+F1β1-FP-γ,11+F1β1-FP-γand under the assumption(21)γ>γC=-F21+Fwe prove in [44] that the fixed point (20) goes through a supercritical Neimark Sacker bifurcation at the threshold(22)P=1F1+γ21+FF1/γ.Hence, beyond threshold an attracting invariant curve is established. We shall now describe the dynamics on such a curve. First, observe that, on the curve, map (19) is topological equivalent to a circle map which actually means that the map does nothing but move or rotate points around the curve. Following [40], associated with a circle map, there is also a rotation number σ which in the context under consideration here may be expressed as(23)σ=c+bdaFH-F,where c=arg⁡λ and FH is the F-value at bifurcation threshold (22) for a fixed value of P. σ may be irrational or rational. In the former case, it is customary to call the orbit of a point quasiperiodic, and by performing a sufficiently large number of iterations the orbit will cover the entire curve as already exemplified in Figure 2(b). If σ=1/n, the dynamic outcome is a periodic orbit of period n. Moreover, the implicit function theorem guarantees that if σ is rational for a value F, there is also an open interval around the parameter where the periodicity is maintained. Actually, on several occasions we have experienced that such intervals may be large. Now, considering (19) we find at threshold (22) that the eigenvalues become(24)λ=-1F±1-1F2iand from (22) it follows that F is “large.” Therefore, λ is located close to the imaginary axis in the left half plane which again implies that c=arg⁡λ≈π/2, so σ≈1/4. Consequently, just beyond threshold (22), the dynamics is an almost 4-periodic orbit restricted to the curve. As we continue to increase F we find through frequency locking an exact 4-periodic orbit. This is shown in Figure 3(a) where γ=-1/10. As we continue to increase F, periodic orbits of period 4·2k are created. In Figure 3(b), we display an orbit of period 16. Eventually the dynamics becomes chaotic. In Figure 3(c), we show a chaotic attractor but clearly there is a tendency towards 4-periodic dynamics in the chaotic regime as well.Figure 3
(a) A4-periodic orbit generated by map (19). γ=-0.1, (P,F)=(0.6,30). (b) A 16-periodic orbit generated by map (19). Parameter values as in (a), except F=40. (c) Chaotic dynamics generated by map (19). The attractor is divided into 4 subsets which are visited once every fourth iteration. Parameter values as in (a), except F=43.
(a)
(b)
(c)Through Examples 2 and 3, we have discussed several mechanisms which may lead to periodic dynamics, both of even and odd periods. However, there are still more possibilities. Indeed, in the map(25)x1,x2⟶Fx2,Pe-x1x1,where 0<P≤1,PF>1, λ is found to be fourth root of unity at bifurcations threshold F=P-1exp⁡(2) and in the map(26)x1,x2⟶Fx1+x2e-x1+x2,x1λ is the third root of unity at threshold. These cases are referred to as the strong resonant cases. They are more difficult to analyze since the location of eigenvalues at threshold violates the assumptions that are necessary in order to prove that an invariant curve beyond bifurcation threshold is established in the Neimark Sacker case.In short, regarding (25), there is no invariant curve beyond threshold. The fixed point bifurcates directly to a 4-periodic orbit with infinitesimal small amplitude. From map (26), one finds that when the fixed point becomes unstable, a large amplitude 3-cycle is the only attractor beyond threshold.


### 3.4. Example 4 (Both Density Dependent Fecundity and Survival)

In our last example, we consider a semelparous population model where we assume density dependence both in the fecundity and in the survival terms. Hence, in the Leslie matrix (1), we assume fi=0,i<n, fn=Fexp(-αx), and pi=Pexp(-βx),i=1,…,n-1. α and β may be regarded as parameters which measure the “strength” of density dependence. If α>β the strength of density dependence is larger in the fecundity than in the survivals.By definingR0, a, and D as (27)R0=FPn-1,FPn-1>1,a=R0-1/α+n-1β,D=∑i=0n-1Paβi,we may express the total equilibrium population as(28)x∗=1α+∑i=1n-1βiln⁡R0while the fixed point may be written as(29)x1∗,…,xi∗,…,xn∗=1Dx∗,…,Paβi-1Dx∗,…,Paβn-1Dx∗.The eigenvalue equation may be cast in the form(30)λn+a1λn-1+a2λn-2+⋯+an=0,where the coefficients are (31)a1=αx1∗+β∑i=2nxi∗,a2=αx2∗+β∑i=1i≠2nxi∗,an=αxn∗+β∑i=1n-1xi∗-1.First let us comment on the parameter region where α≥β. Suppose that n is even and λ=-1. Then LHS (left hand side) of (30) may be expressed as (32)x1∗-x2∗β-α+x3∗-x4∗β-α+⋯+xn-1∗-xn∗β-αwhich clearly is ≤0. On the other hand, when λ→-∞, LHS of (30) →+∞; thus there must be a root λ~<-1 of (30) which actually proves that fixed point (29) will always be unstable.Still assumingα≥β, let us scrutinize the case n=2 in somewhat more detail. In order for the fixed point to be stable, it follows from the Jury criteria that the inequalities(33a)α+βx∗>0,(33b)β-αx1∗-x2∗>0,(33c)2-βx1∗+αx2∗>0must hold, but as already accounted for, (33b) fails when α≥β. Moreover, since (33b) is violated 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(34)x1,t+2=R0e-βx1,t+x2,te-αx1,t+1+x2,t+1x1,t,x2,t+2=R0e-αx1,t+x2,te-βx1,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 in the cycle.(2)
The points are of the form(A,0) or (0,B). In this case, it follows from (34) that A and B must satisfy the equations(35)A-1Pe-βAB=0,βA-αB+ln⁡R0=0.In general, system (35) must be solved by means of numerical methods, but in the special case β=0, we obtain A=(αP)-1ln⁡R0, B=α-1ln⁡R0. Hence there exists a 2-cycle where only one age class is populated at each time, namely, the cycle(36)1αPln⁡R0,0,0,1αln⁡R0.Dynamics where only one age class is populated at each time as above is referred to as SYC (single year class) dynamics or synchronization (cf. [18, 45]). A cycle like (36) is said to be of SYC form. In Figure 4(a), we show an orbit starting at (x1,0,x2,0)≠(0,0) which converges to the 2-cycle (36). If we increase F we find that (36) goes unstable and cycles of period 2k,k=2,3,…, are established. These cycles, which are all of SYC form, are stable in smaller and smaller regions as F is increased. Eventually, the dynamics becomes chaotic but we emphasize that it is still of SYC form. These scenarios are demonstrated in Figures 4(b) and 4(c).Figure 4
(a) Convergence to a2-cycle of SYC form. (P,F)=(0.5,10). (b) A 4-cycle of SYC form. (P,F)=(0.5,15). (c) Chaotic SYC dynamics. (P,F)=(0.5,40).
(a)
(b)
(c)We close this example by making a few comments of the remaining caseβ>α. Still assuming n=2, if α=0, then we are essentially back in Example 3. Hence, the dynamics has a strong resemblance of 4-cycles, either exact or approximate. Whenever α>0, but β>α, the eigenvalues at threshold (which are found by use of (30) and (33c)) are not so close to the imaginary axis as the eigenvalues given by (24). This affects the dynamics beyond threshold. For “small” values of α, we find the invariant curve but there is no tendency towards 4-periodic dynamics. Instead we observe that when F is increased, the invariant curve becomes kinked and not topological equivalent to a circle; see Figure 5(a) where (α,β)=(0.1,1.0). In order to scrutinize the dynamics in somewhat more detail we have also computed the Lyapunov exponent L for F-values between 10 and 50 (cf. Figure 5(b)). As long as 10<F<20.341,L<0, which means that the fixed point is stable. L=0 in the interval 20.341≤F<29.3, which corresponds to quasiperiodic orbits restricted to invariant curves. L<0 in the parameter window 30.5<F<32.7. Here we find periodic dynamics of period 11. Finally when F exceeds 32.7 and in a tiny interval just below F=30.5 we find that L>0, which means that the dynamics is chaotic. These findings are visualized in Figure 5(c) too, where the dynamics is shown in case of specific values of F (10≤F≤35).Figure 5
(a) The invariant curve has become kinked.(α,β)=(0.1,1.0), (P,F)=(0.5,30). (b) Values of Lyapunov exponents for F-values in the interval [10,50] (P=0.5). (c) Dynamics for selected values of F between 10 and 35. Stable fixed points, invariant curves, kinked curves, periodic dynamics, and chaos, respectively.
(a)
(b)
(c)Whenα is close to β, we detect qualitatively much of the same behaviour but since the value of F at threshold in this case is larger than in the former case, one may argue that an increase of strength in the fecundity acts in a stabilizing fashion. Finally, in case of intermediate values of α, α≈β/2, there is a significant change of dynamics. Through the same mechanism as we described in Example 2, there exists an F-interval where there is coexistence between the fixed point and a large amplitude 3-cycle. There is also another interval where we observe coexistence between an invariant curve and the 3-cycle as exemplified in Figure 2(b).


## 3.1. Example 1 (The Linear Case)

Suppose thatfi=Fi,i=1,…,n, and pi=Pi,i=1,…,n-1. Then (2) is a linear map of the form(3)x¯⟶AIx¯,where AI is nonnegative and primitive. On the other hand, if fi=0,i<n,fn=Fn, and pi=Pi,i=1,…,n-1, map (2) becomes(4)x¯⟶ASx¯,where AS is an irreducible and imprimitive (cyclic) matrix.In order to reveal the dynamics generated by (3) and (4), we assume in both cases solutions of the form x¯t=λtw¯ and apply the Perron-Frobenius theorem which may be formulated as follows.Theorem 1 (Perron-Frobenius).
(1)
IfA is a positive or nonnegative and primitive, then there exists a real eigenvalue λ1>0 which is a simple root of the characteristic equation A-λI=0. Moreover, the eigenvalue is strictly greater than the magnitude of any other eigenvalue, λ1>λi for i≠1. The eigenvector w¯1 corresponding to λ1 is real and strictly positive. λ1 may not be the only positive eigenvalue, but if there are others, they do not have nonnegative eigenvectors.
(2)
IfA is irreducible but imprimitive (cyclic) with index of imprimitivity d+1, there exists a real eigenvalue λ1>0 which is a simple root of A-λI=0 with associated eigenvector w¯1>0¯. The eigenvalues λi satisfy λ1≥λi for i≠1, but there are d complex eigenvalues equal in magnitude to λ1 whose values are λ1exp⁡(2kπi/(d+1)),k=1,…,d.Proofs may be obtained in [39, 41, 42].Now, considering map (3), it follows directly from the assumption x¯t=λtw¯ and part 1 of the Perron-Frobenius theorem that x¯t may be expressed as(5)x¯t=c1λ1tw¯1+c2λ2tw¯2+⋯+cnλntw¯n.Here, ci are numbers, λ1 is real and positive, λ1>λi,i≠1, and the eigenvalues in (5) are numbered in order of decreasing magnitude. Moreover,(6)limt→∞x¯tλ1t=c1w¯1.Consequently, the long-term dynamics of the population is described by the growth rate λ1 and the stable population structure w¯1. Thus λ1>1 implies an exponential increasing population and 0≤λ1<1 an exponential decreasing population, where we in all cases have the stable age distribution w¯1. Turning to map (4), we find in a similar way by use of part 2 of the Perron-Frobenius theorem that(7)limt→∞x¯tλ1t=c1w¯1+∑k=2ncke2kπ/nitw¯k.Thus, opposed to the findings from (6), it now follows from (7) that w¯1 is not stable in the sense that an initial population not proportional to w¯1 will not converge to it. Instead the limit in (7) is periodic with period n. These are the most important cases when model (2) is linear. In Figure 1, we show the dynamics generated by maps (3) and (4) in the cases (8)AI=1120.30000.20,AS=0061/20001/30,respectively.Figure 1
(a) Dynamics generated by map (3). (b) Dynamics generated by map (4).
(a)
(b)The dominant eigenvalue (the growth rate) ofAI is λ=1.3, while the stable age distribution is given by the associated eigenvector w¯1=(0.97,0.22,0.03)T. The eigenvalues of AS are λ1=1, λ2,3=exp(±2πi/3), and AS3=I. Consequently, 3-cycle behaviour and no stable age distribution is the outcome.


## 3.2. Example 2 (Density Dependent Fecundity)

Assumefi=Fif^(x),i=1,…,n, where Fi is a constant larger than unity and f^′(x)≤0,pi=Pi,0<Pi<1,i=1,…,n-1. Then model (2) possesses a unique nontrivial fixed point(9)x1∗,…,xn∗=L1Kx∗,…,LnKx∗,where L1=1,Li=P1P2·…·Pi-1,K=∑i=1nLi, and x∗=f^-1(1/R). f^-1 denotes the inverse of the function f^ and R=∑i=1nLiFi. Following [43], the eigenvalue equation may be cast in the form(10)λn-∑i=1nLifi∗+θλn-i=0,where the negative parameter θ is expressed as(11)θ=f^′x∗x1∗∑i=1nLiFi.Note that (10) is valid for a wide range of fecundity functions including all members of the Deriso Schnute family.First, assumefi=Fexp(-x) (the Ricker case). Then x∗=∑FLi and (10) may be written as(12)λn-1x∑i=1nFLi1-x∗λn-i=0.Now, rewrite (12) as g(λ)+h(λ)=0, where g(λ)=λn. Clearly the equation g(λ)=0 has n roots located inside the unit circle C. On the boundary, we have (13)hλ=-1K1-x∗∑i=1nLiλn-i≤L1K1-x∗+⋯+LnK1-x∗≤1-x∗and whenever x∗<2,1-x∗<g(λ)=1. Thus according to Rouche’s theorem, the equations g(λ)+h(λ)=0 and g(λ)=0 have the same number of zeros located inside C. Hence x∗=ln⁡(FLi)<2 guarantees a stable fixed point.Next, assumefi=F(1+x)-1 (the Beverton and Holt case). Then x∗=∑FLi-1 and (10) takes the form(14)λn-1K11+x∗∑i=1nLiλn-1=0.Now, by use of exactly the same kind of arguments as in the Ricker case, we find that stability is achieved if x∗>-2 which of course is always satisfied. Consequently, by use of the Beverton and Holt recruitment function, the only possible dynamics is a stable fixed point. Regarding the Ricker case, the dynamics may be nonstationary and possibly chaotic if x∗ exceeds 2. We shall now explore this further.To this end, consider a two-dimensional version of (1) and (2), where f1=f2=Fexp(-x) and pi=P, that is, the map(15)x1,x2⟶Fe-xx1+Fe-xx2,Px1.The fixed point is easily found to be(16)x1∗,x2∗=11+Px∗,P1+Px∗and by use of the Jury criteria (cf. [38, 39]), it is straightforward to show that (16) is stable provided (17a)x∗<21+P0<P<12,(17b)x∗<1+2PP12<P<1.Assuming 0<P<1/2, the fixed point fails to be hyperbolic at the threshold x∗=2(1+P)-1 or equivalently when F=(1+P)-1exp⁡(2/(1+P)) and undergoes a supercritical flip bifurcation as proved in [16]. Hence, just beyond threshold an attracting period 2 orbit is established. Through further increase of x∗ (which is done by increasing F and fixing P) we observe stable orbits of period 2k,k=1,2,…. Eventually the dynamics becomes chaotic as displayed in Figure 2(a). Note that the attractor is divided into 4 disjoint subsets (branches) which are visited once every fourth iteration. In case of higher F values, these branches merge together.Figure 2
(a) Dynamics generated by map (15). Parameter values (P,F)=(0.2,58.5). (b) Coexisting attractors. (P,F)=(0.9,11.83). (c) Convergence towards the invariant curve. (P,F)=(0.9,11.83) and (x1,0,x2,0)=(2.0,1.78825). (d) Convergence towards the 3-cycle. (P,F)=(0.9,11.83) and (x1,0,x2,0)=(2.0,1.78830). (e) Chaotic dynamics generated by map (15). (P,F)=(0.9,18.5).
(a)
(b)
(c)
(d)
(e)Consider the latter case,1/2<P<1,(x1∗,x2∗) is the only stable attractor as long as x∗ is small, but when x∗ exceeds a certain threshold xA which is lower than bifurcation threshold xB given by (17b), the third iterate g of map (15) undergoes a saddle node bifurcation. Hence, two 3-cycles, one stable and one unstable, are created. This may be justified along the following line: First, note that map g (which has 7 fixed points) may be expressed as(18)x1,x2⟶Fe-FeFe-xx+Px1Fe-xx+Px1+PFe-xx·Fe-Fe-xx+Px1Fe-xx+Px1+PFe-xx,PFe-Fe-xx+Px1Fe-xx+Px1.Then we compute the linearization Dg(x¯1,x¯2), where (x¯1,x¯2) is a fixed point of (18) and find for fixed values of P the value of F where the 3-cycle attractor disappears. Next, we use this F and substitute a corresponding fixed point (x¯1,x¯2) into Dg(x¯1,x¯2) and compute the eigenvalues μ1 and μ2. The results of such calculations are μ1≈1 and 0<μ2<1 from which we conclude that a saddle node bifurcation takes place. For example, if P=0.9, we find F=10.036 which implies that μ1=1.0005 and μ2=0.48 and consequently that threshold xA may be expressed as xA=ln⁡(10.036·1.9)≈2.94. Thus there exists an interval xA<x∗<xB, where there are two stable attractors, the stable 3-cycle and the stable fixed point (16). At xB=P-1(1+2P) (xB=3.11 in case of P=0.9), the fixed point experiences a supercritical Neimark Sacker bifurcation. Consequently, just above threshold, we find coexistence between an attracting invariant curve where the dynamics is quasiperiodic and a stable large amplitude 3-cycle. This is displayed in Figure 2(b). In Figure 2(c), we see the behaviour of points which belong to the trapping region of the invariant curve, and in Figure 2(d) the initial point belongs to the trapping region of the 3-cycle. The scenario referred to above persists in an interval xB<x∗<xC. At xC (xC=3.143 when P=0.9), the invariant curve is hit by the branches of the unstable 3-cycle and disappears. Therefore, whenever x∗>xC and x∗-xC is small, the stable 3-cycle is the only attractor. If we continue to increase x∗ (by increasing F), successive flip bifurcations are the outcomes, creating orbits of period 3·2k. Eventually the dynamics becomes chaotic here too. This is shown in Figure 2(e). Thus to summarize, depending on parameter values, map (15) may generate dynamics of stunning complexity, all generic bifurcations (saddle node, flip, and Neimark Sacker) occur, and there are coexisting attractors as well as chaotic dynamics.


## 3.3. Example 3 (Density Dependent Survival)

Now suppose density dependent survival and consider the map(19)x1,x2⟶Fx2,P1-γβx1/γx1.The fixed point may be expressed as(20)x1∗,x2∗=F1+F1β1-FP-γ,11+F1β1-FP-γand under the assumption(21)γ>γC=-F21+Fwe prove in [44] that the fixed point (20) goes through a supercritical Neimark Sacker bifurcation at the threshold(22)P=1F1+γ21+FF1/γ.Hence, beyond threshold an attracting invariant curve is established. We shall now describe the dynamics on such a curve. First, observe that, on the curve, map (19) is topological equivalent to a circle map which actually means that the map does nothing but move or rotate points around the curve. Following [40], associated with a circle map, there is also a rotation number σ which in the context under consideration here may be expressed as(23)σ=c+bdaFH-F,where c=arg⁡λ and FH is the F-value at bifurcation threshold (22) for a fixed value of P. σ may be irrational or rational. In the former case, it is customary to call the orbit of a point quasiperiodic, and by performing a sufficiently large number of iterations the orbit will cover the entire curve as already exemplified in Figure 2(b). If σ=1/n, the dynamic outcome is a periodic orbit of period n. Moreover, the implicit function theorem guarantees that if σ is rational for a value F, there is also an open interval around the parameter where the periodicity is maintained. Actually, on several occasions we have experienced that such intervals may be large. Now, considering (19) we find at threshold (22) that the eigenvalues become(24)λ=-1F±1-1F2iand from (22) it follows that F is “large.” Therefore, λ is located close to the imaginary axis in the left half plane which again implies that c=arg⁡λ≈π/2, so σ≈1/4. Consequently, just beyond threshold (22), the dynamics is an almost 4-periodic orbit restricted to the curve. As we continue to increase F we find through frequency locking an exact 4-periodic orbit. This is shown in Figure 3(a) where γ=-1/10. As we continue to increase F, periodic orbits of period 4·2k are created. In Figure 3(b), we display an orbit of period 16. Eventually the dynamics becomes chaotic. In Figure 3(c), we show a chaotic attractor but clearly there is a tendency towards 4-periodic dynamics in the chaotic regime as well.Figure 3
(a) A4-periodic orbit generated by map (19). γ=-0.1, (P,F)=(0.6,30). (b) A 16-periodic orbit generated by map (19). Parameter values as in (a), except F=40. (c) Chaotic dynamics generated by map (19). The attractor is divided into 4 subsets which are visited once every fourth iteration. Parameter values as in (a), except F=43.
(a)
(b)
(c)Through Examples 2 and 3, we have discussed several mechanisms which may lead to periodic dynamics, both of even and odd periods. However, there are still more possibilities. Indeed, in the map(25)x1,x2⟶Fx2,Pe-x1x1,where 0<P≤1,PF>1, λ is found to be fourth root of unity at bifurcations threshold F=P-1exp⁡(2) and in the map(26)x1,x2⟶Fx1+x2e-x1+x2,x1λ is the third root of unity at threshold. These cases are referred to as the strong resonant cases. They are more difficult to analyze since the location of eigenvalues at threshold violates the assumptions that are necessary in order to prove that an invariant curve beyond bifurcation threshold is established in the Neimark Sacker case.In short, regarding (25), there is no invariant curve beyond threshold. The fixed point bifurcates directly to a 4-periodic orbit with infinitesimal small amplitude. From map (26), one finds that when the fixed point becomes unstable, a large amplitude 3-cycle is the only attractor beyond threshold.


## 3.4. Example 4 (Both Density Dependent Fecundity and Survival)

In our last example, we consider a semelparous population model where we assume density dependence both in the fecundity and in the survival terms. Hence, in the Leslie matrix (1), we assume fi=0,i<n, fn=Fexp(-αx), and pi=Pexp(-βx),i=1,…,n-1. α and β may be regarded as parameters which measure the “strength” of density dependence. If α>β the strength of density dependence is larger in the fecundity than in the survivals.By definingR0, a, and D as (27)R0=FPn-1,FPn-1>1,a=R0-1/α+n-1β,D=∑i=0n-1Paβi,we may express the total equilibrium population as(28)x∗=1α+∑i=1n-1βiln⁡R0while the fixed point may be written as(29)x1∗,…,xi∗,…,xn∗=1Dx∗,…,Paβi-1Dx∗,…,Paβn-1Dx∗.The eigenvalue equation may be cast in the form(30)λn+a1λn-1+a2λn-2+⋯+an=0,where the coefficients are (31)a1=αx1∗+β∑i=2nxi∗,a2=αx2∗+β∑i=1i≠2nxi∗,an=αxn∗+β∑i=1n-1xi∗-1.First let us comment on the parameter region where α≥β. Suppose that n is even and λ=-1. Then LHS (left hand side) of (30) may be expressed as (32)x1∗-x2∗β-α+x3∗-x4∗β-α+⋯+xn-1∗-xn∗β-αwhich clearly is ≤0. On the other hand, when λ→-∞, LHS of (30) →+∞; thus there must be a root λ~<-1 of (30) which actually proves that fixed point (29) will always be unstable.Still assumingα≥β, let us scrutinize the case n=2 in somewhat more detail. In order for the fixed point to be stable, it follows from the Jury criteria that the inequalities(33a)α+βx∗>0,(33b)β-αx1∗-x2∗>0,(33c)2-βx1∗+αx2∗>0must hold, but as already accounted for, (33b) fails when α≥β. Moreover, since (33b) is violated 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(34)x1,t+2=R0e-βx1,t+x2,te-αx1,t+1+x2,t+1x1,t,x2,t+2=R0e-αx1,t+x2,te-βx1,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 in the cycle.(2)
The points are of the form(A,0) or (0,B). In this case, it follows from (34) that A and B must satisfy the equations(35)A-1Pe-βAB=0,βA-αB+ln⁡R0=0.In general, system (35) must be solved by means of numerical methods, but in the special case β=0, we obtain A=(αP)-1ln⁡R0, B=α-1ln⁡R0. Hence there exists a 2-cycle where only one age class is populated at each time, namely, the cycle(36)1αPln⁡R0,0,0,1αln⁡R0.Dynamics where only one age class is populated at each time as above is referred to as SYC (single year class) dynamics or synchronization (cf. [18, 45]). A cycle like (36) is said to be of SYC form. In Figure 4(a), we show an orbit starting at (x1,0,x2,0)≠(0,0) which converges to the 2-cycle (36). If we increase F we find that (36) goes unstable and cycles of period 2k,k=2,3,…, are established. These cycles, which are all of SYC form, are stable in smaller and smaller regions as F is increased. Eventually, the dynamics becomes chaotic but we emphasize that it is still of SYC form. These scenarios are demonstrated in Figures 4(b) and 4(c).Figure 4
(a) Convergence to a2-cycle of SYC form. (P,F)=(0.5,10). (b) A 4-cycle of SYC form. (P,F)=(0.5,15). (c) Chaotic SYC dynamics. (P,F)=(0.5,40).
(a)
(b)
(c)We close this example by making a few comments of the remaining caseβ>α. Still assuming n=2, if α=0, then we are essentially back in Example 3. Hence, the dynamics has a strong resemblance of 4-cycles, either exact or approximate. Whenever α>0, but β>α, the eigenvalues at threshold (which are found by use of (30) and (33c)) are not so close to the imaginary axis as the eigenvalues given by (24). This affects the dynamics beyond threshold. For “small” values of α, we find the invariant curve but there is no tendency towards 4-periodic dynamics. Instead we observe that when F is increased, the invariant curve becomes kinked and not topological equivalent to a circle; see Figure 5(a) where (α,β)=(0.1,1.0). In order to scrutinize the dynamics in somewhat more detail we have also computed the Lyapunov exponent L for F-values between 10 and 50 (cf. Figure 5(b)). As long as 10<F<20.341,L<0, which means that the fixed point is stable. L=0 in the interval 20.341≤F<29.3, which corresponds to quasiperiodic orbits restricted to invariant curves. L<0 in the parameter window 30.5<F<32.7. Here we find periodic dynamics of period 11. Finally when F exceeds 32.7 and in a tiny interval just below F=30.5 we find that L>0, which means that the dynamics is chaotic. These findings are visualized in Figure 5(c) too, where the dynamics is shown in case of specific values of F (10≤F≤35).Figure 5
(a) The invariant curve has become kinked.(α,β)=(0.1,1.0), (P,F)=(0.5,30). (b) Values of Lyapunov exponents for F-values in the interval [10,50] (P=0.5). (c) Dynamics for selected values of F between 10 and 35. Stable fixed points, invariant curves, kinked curves, periodic dynamics, and chaos, respectively.
(a)
(b)
(c)Whenα is close to β, we detect qualitatively much of the same behaviour but since the value of F at threshold in this case is larger than in the former case, one may argue that an increase of strength in the fecundity acts in a stabilizing fashion. Finally, in case of intermediate values of α, α≈β/2, there is a significant change of dynamics. Through the same mechanism as we described in Example 2, there exists an F-interval where there is coexistence between the fixed point and a large amplitude 3-cycle. There is also another interval where we observe coexistence between an invariant curve and the 3-cycle as exemplified in Figure 2(b).


## 4. Summary

As it is clear from the examples, nonlinear age-structured population models are excellent tools in order to study the dynamical outcomes of biological populations. Depending on location of density dependent elements and parameter values, the dynamics may vary from stable fixed point to complicated chaotic behaviour. We close the paper by adding a few more results. In Example 2 (the Ricker case), we have to distinguish between even and odd age classes. Whenn is odd, the transfer from stability to instability always goes through a (supercritical) flip bifurcation (0<P<1), but when n is even a Neimark Sacker bifurcation occurs when P is large (when n=4 the Neimark Sacker interval is 0.62<P<1). Moreover, when n is large, x∗(P) is an increasing function which clearly suggests that an increase of n acts stabilizing to the dynamics. For further reading, see [12, 16]. Turning to Example 3, the 4-periodic dynamics (exact or approximated) found when n=2 takes over when n becomes large. Formally, this may be proved through an asymptotic argument where one shows that the dominant eigenvalues of the linearization of the n-dimensional version of map (19) cross the unit circle very close to the imaginary axis at the threshold. For details, see [44]. Regarding Example 4, we want to stress the following: assuming α>β and n is odd, it is proved in a forthcoming paper that SYC dynamics is the only possibility in cases of small and large equilibrium populations x∗. However in case of intermediate values there may exist a parameter region where the fixed point is stable. Depending on initial conditions, an orbit may converge to the fixed point or settle on a chaotic attractor of SYC form. This is demonstrated in [18] when β=0. Therefore, by combining the results above with the findings presented in Example 4, we find it natural to conclude that as long as the strength of the density dependence in the fecundity is larger than in the survivals, SYC dynamics is indeed the most likely dynamic outcome. Finally, when β>α, we find in contrast to the iteroparous model discussed in Example 2 that an increase of the number of age classes n acts in a destabilizing fashion. Actually, if P is sufficiently small, it is not possible to achieve a stable fixed point at all if n=5. Moreover, the periodic dynamics found when n=2 persists also when n=3 (cf. the simpler model presented in [39]), but the route to chaos is different from the n=2 case. When n≥4, only quasiperiodic dynamics seems to occur beyond instability threshold.



---

*Source: 1964286-2017-08-06.xml*