Nonstationary and Chaotic Dynamics in Age-Structured Population Models

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-dimensionalmodels. Therefore amore 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. Thepaper is structured as follows: in Section 2, 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.


Introduction
As it is well known (cf.[1][2][3][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][7][8][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][14][15][16][17][18][19][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][24][25], prey-predator systems [26][27][28][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][34][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 Section 2, 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.

The Model, Fixed Points, and Stability
First, we establish the model.At time , we split the population  into  distinct nonoverlapping age classes   = ( 1, , . . .,  , )  , where  = ∑   is the total population.Next, we introduce the Leslie matrix where   is the average fecundity of a member of the th age class at time , while   denotes the survival probability of age class .Then the relation between  at two consecutive time steps is given by the map  : R  → R   → . ( 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   > 0 for all  ≤ , we say that the species has an iteroparous life history and matrix  is classified as a nonnegative and primitive matrix.On the other hand, if   = 0,  <  and   > 0, the species possess a semelparous life history and in this case we say that  is an irreducible but imprimitive (cyclic) matrix.It is also customary to divide both the iteroparous and semelparous cases into two subclasses.If  is small and   > 0 for all , 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   = 0,  <  ( ≈ /2 "very roughly") and   > 0,  ≤  ≤ .Humans and large mammals fall into this category.Regarding semelparity, whenever the number of age classes  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 constant   =   ,   =   , 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,   =   , while the fecundities are density dependent.In the striped bass fishery model [12], it was assumed that   =   (), where  = ∑     is a weighted sum of age classes.Assuming two age classes only, cannibalism may be accounted for by the assumption  1 =  1 ( 2 ) which suggests that individuals from the second age class prey upon members from the first age class.A third possibility is   =   () or   =   () which means that the elements in (1) depend on the total population .Two of the most frequently used fecundity functions are the compensatory Beverton and Holt relation   () =   (1 + ) −1 ; see [13,36] and the overcompensatory Ricker relation,   () =   exp(−), 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   () =   (1 − ) 1/ ( = −1 gives the Beverton and Holt relation while  → 0 results in the Ricker case).Also, observe that   () is written as a product of a constant term   and a density dependent part, for example, exp(−).A final comment is that    () ≤ 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.Assuming   =   () and (or)   =   (), model (2) 2), the eigenvalues of the linearization may leave the unit circle through  = 1,  = −1, or  = ± exp(), where  = 2/,  = 0, 1, . . .,  − 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].
In order to reveal the dynamics generated by ( 3) and (4), we assume in both cases solutions of the form   =    and apply the Perron-Frobenius theorem which may be formulated as follows.Theorem 1 (Perron-Frobenius).
(1) If  is a positive or nonnegative and primitive, then there exists a real eigenvalue  1 > 0 which is a simple root of the characteristic equation | − | = 0.Moreover, the eigenvalue is strictly greater than the magnitude of any other eigenvalue, The eigenvector  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.
Proofs may be obtained in [39,41,42].Now, considering map (3), it follows directly from the assumption   =    and part 1 of the Perron-Frobenius theorem that   may be expressed as Here,   are numbers,  1 is real and positive,  1 > |  |,  ̸ = 1, and the eigenvalues in (5) are numbered in order of decreasing magnitude.Moreover, Consequently, the long-term dynamics of the population is described by the growth rate  1 and the stable population structure  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  1 .
Turning to map (4), we find in a similar way by use of part 2 of the Perron-Frobenius theorem that lim Thus, opposed to the findings from (6), it now follows from (7) that  1 is not stable in the sense that an initial population not proportional to  1 will not converge to it.Instead the limit in ( 7) is periodic with period .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 respectively.

Example 2 (Density Dependent Fecundity).
Assume   =   f(),  = 1, . . ., , where   is a constant larger than unity and f () ≤ 0,   =   , 0 <   < 1,  = 1, . . .,  − 1.Then model (2) possesses a unique nontrivial fixed point where  1 = 1,   =  1  2 ⋅ . . .⋅  −1 ,  = ∑  =1   , and  * = f−1 (1/).f−1 denotes the inverse of the function f and  = ∑  =1     .Following [43], the eigenvalue equation may be cast in the form where the negative parameter  is expressed as Note that ( 10) is valid for a wide range of fecundity functions including all members of the Deriso Schnute family.First, assume   =  exp (−) (the Ricker case).Then  * = ∑   and (10) may be written as Now, rewrite (12) as ()+ℎ() = 0, where () =   .Clearly the equation () = 0 has  roots located inside the unit circle .On the boundary, we have Now, by use of exactly the same kind of arguments as in the Ricker case, we find that stability is achieved if  * > −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  * exceeds 2. We shall now explore this further.
To this end, consider a two-dimensional version of (1) and (2), where  1 =  2 =  exp (−) and   = , that is, the map The fixed point is easily found to be and by use of the Jury criteria (cf.[38,39]), it is straightforward to show that ( 16) is stable provided Assuming 0 <  < 1/2, the fixed point fails to be hyperbolic at the threshold  * = 2(1 + ) −1 or equivalently when  = (1 + ) −1 exp(2/(1 + )) 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  * (which is done by increasing  and fixing ) we observe stable orbits of period 2  ,  = 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  values, these branches merge together.Consider the latter case, 1/2 <  < 1, ( * 1 ,  * 2 ) is the only stable attractor as long as  * is small, but when  * exceeds a certain threshold   which is lower than bifurcation threshold   given by (17b), the third iterate  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  (which has 7 fixed points) may be expressed as −( − + 1 ) ( −  +  1 )} .

(18)
Then we compute the linearization ( 1 ,  2 ), where ( 1 ,  2 ) is a fixed point of ( 18) and find for fixed values of  the value of  where the 3-cycle attractor disappears.Next, we use this  and substitute a corresponding fixed point ( 1 ,  2 ) into ( 1 ,  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  = 0.9, we find  = 10.036 which implies that  1 = 1.0005 and  2 = 0.48 and consequently that threshold   may be expressed as   = ln (10.036 ⋅ 1.9) ≈ 2.94.Thus there exists an interval   <  * <   , where there are two stable attractors, the stable 3-cycle and the stable fixed point (16).At   =  −1 (1 + 2) (  = 3.11 in case of  = 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   <  * <   .At   (  = 3.143 when  = 0.9), the invariant curve is hit by the branches of the unstable 3-cycle and disappears.Therefore, whenever  * >   and | * −  | is small, the stable 3-cycle is the only attractor.If we continue to increase  * (by increasing ), successive flip bifurcations are the outcomes, creating orbits of period 3 ⋅ 2  .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.

Example 3 (Density Dependent Survival
).Now suppose density dependent survival and consider the map The fixed point may be expressed as and under the assumption we prove in [44] that the fixed point ( 20) goes through a supercritical Neimark Sacker bifurcation at the threshold 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 where  = arg  and   is the -value at bifurcation threshold (22) for a fixed value of . 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/, the dynamic outcome is a periodic orbit of period .Moreover, the implicit function theorem guarantees that if  is rational for a value , 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 and from (22) it follows that  is "large."Therefore,  is located close to the imaginary axis in the left half plane which again implies that  = 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  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 , periodic orbits of period 4 ⋅ 2  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.
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 where 0 <  ≤ 1,  > 1,  is found to be fourth root of unity at bifurcations threshold  =  −1 exp(2) and in the map 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 4periodic 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.

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   = 0,  < ,   =  exp (−), and   = exp (−),  = 1, . . .,  − 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 defining  0 , , and  as we may express the total equilibrium population as while the fixed point may be written as The eigenvalue equation may be cast in the form where the coefficients are First let us comment on the parameter region where  ≥ .Suppose that  is even and  = −1.Then LHS (left hand side) of ( 30) may be expressed as 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  = 2 in somewhat more detail.In order for the fixed point to be stable, it follows from the Jury criteria that the inequalities must 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  * is small.Evidently, such a 2-cycle must be obtained from and here there are two possibilities: (1)   =  +1 which leads to the trivial 2-cycle where the unstable fixed point ( * 1 ,  * 2 ) is the only point in the cycle.
(2) The points are of the form (, 0) or (0, ).In this case, it follows from (34) that  and  must satisfy the equations In general, system (35) must be solved by means of numerical methods, but in the special case  = 0, we obtain  = () −1 ln  0 ,  =  −1 ln  0 .Hence there exists a 2-cycle where only one age class is populated at each time, namely, the cycle 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 ( 1,0 ,  2,0 ) ̸ = (0, 0) which converges to the 2-cycle (36).If we increase  we find that (36) goes unstable and cycles of period 2  ,  = 2, 3, . .., are established.These cycles, which are all of SYC form, are stable in smaller and smaller regions as  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).
We close this example by making a few comments of the remaining case  > .Still assuming  = 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  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  for -values between 10 and 50 (cf.Figure 5(b)).As long as 10 <  < 20.341,  < 0, which means that the fixed point is stable. = 0 in the interval 20.341 ≤  < 29.3, which corresponds to quasiperiodic orbits restricted to invariant curves. < 0 in the parameter window 30.5 <  < 32.7.Here we find periodic dynamics of period 11.Finally when  exceeds 32.7 and in a tiny interval just below  = 30.5 we find that  > 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  (10 ≤  ≤ 35).
When  is close to , we detect qualitatively much of the same behaviour but since the value of  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 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).

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.When  is odd, the transfer from stability to instability always goes through a (supercritical) flip bifurcation (0 <  < 1), but when  is even a Neimark Sacker bifurcation occurs when  is large (when  = 4 the Neimark Sacker interval is 0.62 <  < 1).Moreover, when  is large,  * () is an increasing function which clearly suggests that an increase of  acts stabilizing to the dynamics.For further reading, see [12,16].Turning to Example 3, the 4-periodic dynamics (exact or approximated) found when  = 2 takes over when  becomes large.Formally, this may be proved through an asymptotic argument where one shows that the dominant eigenvalues of the linearization of the -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  is odd, it is proved in a forthcoming paper that SYC dynamics is the only possibility in cases of small and large equilibrium populations  * .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  acts in a destabilizing fashion.Actually, if  is sufficiently small, it is not possible to achieve a stable fixed point at all if  = 5.Moreover, the periodic dynamics found when  = 2 persists also when  = 3 (cf.the simpler model presented in [39]), but the route to chaos is different from the  = 2 case.When  ≥ 4, only quasiperiodic dynamics seems to occur beyond instability threshold.