Discrete Dynamics in Nature and Society
Volume 2017 (2017), Article ID 1964286, 11 pages
https://doi.org/10.1155/2017/1964286

Research Article

Nonstationary and Chaotic Dynamics in Age-Structured Population Models

Arild Wikan and Ørjan Kristensen

School of Business and Economics, The Arctic University of Norway, Campus Harstad, Harstad, Norway

Correspondence should be addressed to Arild Wikan; arild.wikan@uit.no

Received 31 March 2017; Revised 10 May 2017; Accepted 18 May 2017; Published 6 August 2017

Academic Editor: Abdelalim Elsadany

Copyright © 2017 Arild Wikan and Ørjan Kristensen. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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 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.

2. The Model, Fixed Points, and Stability

First, we establish the model. At time , we split the population into distinct nonoverlapping age classes , 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

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 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 and , 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 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 ( “very roughly”) and . 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 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 ; see [13, 36] and the overcompensatory Ricker relation, , 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 ( gives the Beverton and Holt relation while results in the Ricker case). Also, observe that is written as a product of a constant term and a density dependent part, for example, . A final comment is that . 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) possesses a unique nontrivial fixed point of hyperbolic type. 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 . 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 is small, the Jury criteria (see the books [38, 39]) give conditions for to be stable. In cases where is large, the criteria become complicated and not applicable unless in very special situations. It may be a severe problem to decide if all satisfy when is large. If we increase a parameter ( or ) such that an eigenvalue leaves the unit circle, the fixed point becomes unstable. The threshold where loses its hyperbolicity (i.e., where ) 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 , , or = ±, where . If , 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]. results in a flip or periodic doubling bifurcation. Hence, assuming that the bifurcation is of supercritical nature, an attracting period orbit is established when the fixed point goes unstable. When a pair of complex modulus 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 that , and . Then (2) is a linear map of the form

where is nonnegative and primitive. On the other hand, if , and , map (2) becomes

where 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 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 which is a simple root of the characteristic equation . Moreover, the eigenvalue is strictly greater than the magnitude of any other eigenvalue, for . The eigenvector corresponding to is real and strictly positive. may not be the only positive eigenvalue, but if there are others, they do not have nonnegative eigenvectors.(2)If is irreducible but imprimitive (cyclic) with index of imprimitivity , there exists a real eigenvalue which is a simple root of with associated eigenvector . The eigenvalues satisfy for , but there are complex eigenvalues equal in magnitude to whose values are λ ₁.

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, is real and positive, , 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 and the stable population structure . Thus implies an exponential increasing population and an exponential decreasing population, where we in all cases have the stable age distribution . Turning to map (4), we find in a similar way by use of part 2 of the Perron-Frobenius theorem that

Thus, opposed to the findings from (6), it now follows from (7) that is not stable in the sense that an initial population not proportional to 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.

(a)

(b)

Figure 1: (a) Dynamics generated by map (3). (b) Dynamics generated by map (4).

The dominant eigenvalue (the growth rate) of is , while the stable age distribution is given by the associated eigenvector . The eigenvalues of are , , and . Consequently, 3-cycle behaviour and no stable age distribution is the outcome.

3.2. Example 2 (Density Dependent Fecundity)

Assume , where is a constant larger than unity and . Then model (2) possesses a unique nontrivial fixed point

where , and . denotes the inverse of the function and . 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 (the Ricker case). Then and (10) may be written as

Now, rewrite (12) as , where . Clearly the equation has roots located inside the unit circle . On the boundary, we have

and whenever . Thus according to Rouche’s theorem, the equations and have the same number of zeros located inside . Hence guarantees a stable fixed point.

Next, assume (the Beverton and Holt case). Then and (10) takes the form

Now, by use of exactly the same kind of arguments as in the Ricker case, we find that stability is achieved if 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 . We shall now explore this further.

To this end, consider a two-dimensional version of (1) and (2), where 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 , the fixed point fails to be hyperbolic at the threshold or equivalently when and undergoes a supercritical flip bifurcation as proved in [16]. Hence, just beyond threshold an attracting period orbit is established. Through further increase of (which is done by increasing and fixing ) we observe stable orbits of period Eventually the dynamics becomes chaotic as displayed in Figure 2(a). Note that the attractor is divided into disjoint subsets (branches) which are visited once every fourth iteration. In case of higher values, these branches merge together.

(a)

(b)

(c)

(d)

(e)

Figure 2: (a) Dynamics generated by map (15). Parameter values . (b) Coexisting attractors. . (c) Convergence towards the invariant curve. and . (d) Convergence towards the -cycle. and . (e) Chaotic dynamics generated by map (15). .

Consider the latter case, 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 -cycles, one stable and one unstable, are created. This may be justified along the following line: First, note that map (which has fixed points) may be expressed as

Then we compute the linearization , where 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 into and compute the eigenvalues and . The results of such calculations are and from which we conclude that a saddle node bifurcation takes place. For example, if , we find which implies that and and consequently that threshold may be expressed as . Thus there exists an interval , where there are two stable attractors, the stable -cycle and the stable fixed point (16). At ( in case of ), 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 -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 -cycle. The scenario referred to above persists in an interval . At ( when ), the invariant curve is hit by the branches of the unstable -cycle and disappears. Therefore, whenever and is small, the stable -cycle is the only attractor. If we continue to increase (by increasing ), successive flip bifurcations are the outcomes, creating orbits of period . 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

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 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 , 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 , so . Consequently, just beyond threshold (22), the dynamics is an almost -periodic orbit restricted to the curve. As we continue to increase we find through frequency locking an exact -periodic orbit. This is shown in Figure 3(a) where . As we continue to increase , periodic orbits of period are created. In Figure 3(b), we display an orbit of period . Eventually the dynamics becomes chaotic. In Figure 3(c), we show a chaotic attractor but clearly there is a tendency towards -periodic dynamics in the chaotic regime as well.

(a)

(b)

(c)

Figure 3: (a) A -periodic orbit generated by map (19). , . (b) A -periodic orbit generated by map (19). Parameter values as in (a), except . (c) Chaotic dynamics generated by map (19). The attractor is divided into subsets which are visited once every fourth iteration. Parameter values as in (a), except .

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 , is found to be fourth root of unity at bifurcations threshold 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 -periodic orbit with infinitesimal small amplitude. From map (26), one finds that when the fixed point becomes unstable, a large amplitude -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 , , and . 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 , , 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 . 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 of (30) which actually proves that fixed point (29) will always be unstable.

Still assuming , let us scrutinize the case 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 , it is natural to search for a -cycle which should be stable provided is small. Evidently, such a -cycle must be obtained from

and here there are two possibilities:

(1) which leads to the trivial -cycle where the unstable fixed point is the only point in the cycle.(2)The points are of the form or . 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 , we obtain , . Hence there exists a -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 which converges to the -cycle (36). If we increase we find that (36) goes unstable and cycles of period , 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).

(a)

(b)

(c)

Figure 4: (a) Convergence to a -cycle of SYC form. . (b) A -cycle of SYC form. . (c) Chaotic SYC dynamics. .

We close this example by making a few comments of the remaining case . Still assuming , if , then we are essentially back in Example 3. Hence, the dynamics has a strong resemblance of -cycles, either exact or approximate. Whenever , 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 -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 . In order to scrutinize the dynamics in somewhat more detail we have also computed the Lyapunov exponent for -values between and (cf. Figure 5(b)). As long as , which means that the fixed point is stable. in the interval , which corresponds to quasiperiodic orbits restricted to invariant curves. in the parameter window . Here we find periodic dynamics of period . Finally when exceeds and in a tiny interval just below we find that , 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 ().

(a)

(b)

(c)

Figure 5: (a) The invariant curve has become kinked. , . (b) Values of Lyapunov exponents for -values in the interval (). (c) Dynamics for selected values of between and . Stable fixed points, invariant curves, kinked curves, periodic dynamics, and chaos, respectively.

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 , , 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 -cycle. There is also another interval where we observe coexistence between an invariant curve and the -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. When is odd, the transfer from stability to instability always goes through a (supercritical) flip bifurcation (), but when is even a Neimark Sacker bifurcation occurs when is large (when the Neimark Sacker interval is ). 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 -periodic dynamics (exact or approximated) found when 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 . 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 . Moreover, the periodic dynamics found when persists also when (cf. the simpler model presented in [39]), but the route to chaos is different from the case. When , only quasiperiodic dynamics seems to occur beyond instability threshold.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The publication charges for this article have been covered by a grant from the publication fund of UiT The Arctic University of Norway.

References

A. N. Sarkovskii, “Coexistence of cycles of a continuous map of a line into itself,” Ukrainian Mathematical Journal, vol. 16, pp. 61–71, 1964.

T. Y. Li and J. A. Yorke, “Period three implies chaos,” The American Mathematical Monthly, vol. 82, no. 10, pp. 985–992, 1975.

R. M. May, “Simple mathematical models with very complicated dynamics,” Nature, vol. 261, no. 5560, pp. 459–467, 1976.

R. M. May and G. F. Oster, “Bifurcations and Dynamic Complexity in Simple Ecological Models,” The American Naturalist, vol. 110, no. 974, pp. 573–599, 1976.

H. Thunberg, “Periodicity versus chaos in one-dimensional dynamics,” SIAM Review, vol. 43, no. 1, pp. 3–30, 2001.

H. Bernardelli, “Population waves,” Journal of Burma Research Society, vol. 31, pp. 1–18, 1941.

E. G. Lewis, “On the generation and growth of a population,” Sankhya: The Indian Journal of Statistics, vol. 6, pp. 93–96, 1942.

P. H. Leslie, “On the use of matrices in certain population mathematics,” Biometrika, vol. 33, pp. 183–212, 1945.

P. H. Leslie, “Some further notes on the use of matrices in population mathematics,” Biometrika, vol. 35, pp. 213–245, 1948.

L. Pennycuick, “A computer model of the oxford great tit population,” Journal of Theoretical Biology, vol. 22, no. 3, pp. 381–400, 1969.

J. Guckenheimer, G. Oster, and A. Ipaktchi, “The dynamics of density dependent population models,” Journal of Mathematical Biology, vol. 4, no. 2, pp. 101–147, 1977.

S. A. Levin and C. P. Goodyear, “Analysis of an age-structured fishery model,” Journal of Mathematical Biology, vol. 9, no. 3, pp. 245–274, 1980.

J. A. Silva and T. G. Hallam, “Compensation and stability in nonlinear matrix models,” Mathematical Biosciences, vol. 110, no. 1, pp. 67–101, 1992.

J. A. Silva and T. G. Hallam, “Effects of delay, truncations and density dependence in reproduction schedules on stability of nonlinear Leslie matrix models,” Journal of Mathematical Biology, vol. 31, no. 4, pp. 367–395, 1993.

A. Wikan and E. Mjølhus, “Periodicity of 4 in age-structured population models with density dependence,” Journal of Theoretical Biology, vol. 173, no. 2, pp. 109–119, 1995.

A. Wikan and E. Mjølhus, “Overcompensatory recruitment and generation delay in discrete age-structured population models,” Journal of Mathematical Biology, vol. 35, no. 2, pp. 195–239, 1996.

J. M. Cushing, “Nonlinear semelparous Leslie models,” Mathematical Biosciences and Engineering. MBE, vol. 3, no. 1, pp. 17–36, 2006.

E. Mjølhus, A. Wikan, and T. Solberg, “On synchronization in semelparous populations,” Journal of Mathematical Biology, vol. 50, no. 1, pp. 1–21, 2005.

H. Behncke, “Periodical cicadas,” Journal of Mathematical Biology, vol. 40, no. 5, pp. 413–431, 2000.

I. Ugarcovici and H. Weiss, “Chaotic dynamics of a nonlinear density dependent population model,” Nonlinearity, vol. 17, no. 5, pp. 1689–1711, 2004.

M. Gyllenberg, G. Soderbacka, and S. Ericsson, “Does migration stabilize local population dynamics? Analysis of a discrete metapopulation model,” Mathematical Biosciences, vol. 118, no. 1, pp. 25–49, 1993.

F. Lutscher and M. A. Lewis, “Spatially-explicit matrix models,” Journal of Mathematical Biology, vol. 48, no. 3, pp. 293–324, 2004.

K. L. Cooke, R. Elderkin, and M. Witten, “Harvesting procedures with management policy in iterative density-dependent population models,” Natural Resource Modeling, vol. 2, no. 3, pp. 383–420, 1988.

S. Tang and L. Chen, “The effect of seasonal harvesting on stage-structured population models,” Journal of Mathematical Biology, vol. 48, no. 4, pp. 357–374, 2004.

A. Wikan, “On the interplay between cannibalism and harvest in stage-structured population models,” Journal of Marine Biology, vol. 2015, Article ID 580520, 8 pages, 2015.

M. G. Neubert and M. Kot, “The subcritical collapse of predator populations in discrete-time predator-prey models,” Mathematical Biosciences, vol. 110, no. 1, pp. 45–66, 1992.

A. Wikan, “From chaos to chaos. An analysis of a discrete age-structured prey-predator model,” Journal of Mathematical Biology, vol. 43, no. 6, pp. 471–500, 2001.

S. A. Gourley and Y. Kuang, “A stage structured predator-prey model and its dependence on maturation delay and death rate,” Journal of Mathematical Biology, vol. 49, no. 2, pp. 188–200, 2004.

A. Wikan, “An analysis of discrete stage-structured prey and prey-predator population models,” Discrete Dynamics in Nature and Society. An International Multidisciplinary Research and Review Journal, Article ID 9475854, 11 pages, 2017.

J. M. Cushing, “A strong ergodic theorem for some nonlinear matrix models for the dynamics of structured populations,” Natural Resource Modeling, vol. 3, no. 3, pp. 331–357, 1989.

K. M. Crowe, “A nonlinear ergodic theorem for discrete systems,” Journal of Mathematical Biology, vol. 32, no. 3, pp. 179–191, 1994.

R. Kon, Y. Saito, and Y. Takeuchi, “Permanence of single-species stage-structured models,” Journal of Mathematical Biology, vol. 48, no. 5, pp. 515–528, 2004.

A. A. Elsadany and A. E. Matouk, “Dynamical behaviors of fractional-order Lotka-Volterra predator-prey model and its discretization,” Journal of Applied Mathematics and Computing, vol. 49, no. 1-2, pp. 269–283, 2015.

A. E. Matouk, A. A. Elsadany, E. Ahmed, and H. N. Agiza, “Dynamical behavior of fractional-order Hastings-Powell food chain model and its discretization,” Communications in Nonlinear Science and Numerical Simulation, vol. 27, pp. 153–167, 2015.

A. E. Matouk and A. A. Elsadany, “Dynamical analysis, stabilization and discretization of a chaotic fractional-order GLV model,” Nonlinear Dynamics. An International Journal of Nonlinear Dynamics and Chaos in Engineering Systems, vol. 85, no. 3, pp. 1597–1612, 2016.

D. L. Deangelis, L. J. Svoboda, S. W. Christensen, and D. S. Vaughan, “Stability and return times of Leslie matrices with density-dependent survival: applications to fish populations,” Ecological Modelling, vol. 8, no. C, pp. 149–163, 1980.

E. W. Ricker, “Stock and recruitment,” Journal of the Fisheries Research Board of Canada, vol. 11, no. 5, pp. 559–623, 1954.

J. D. Murray, Mathematical Biology: Spatial Models and Biomedical Applications, Springer, New York, NY, USA, 3rd edition, 1993.

A. Wikan, “Discrete dynamical systems,” Bookboon.com (Ventus publishing Aps), 2013. ISBN 978-87-403-0327-8.

J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcation of Vector Fields, Springer, New York, NY, USA, 1983.

R. A. Horn and C. R. Johnson, Matrix Analysis, Cambridge University Press, Cambridge, UK, 1985.

J. C. Frauenthal, “Analysis of age-structure models,” in Mathematical Ecology of Biomathematics, T. G. Hallam and S. A. Levin, Eds., vol. 17 of Biomathematics, pp. 117–147, Springer, Berlin, Germany, 1986.

A. Wikan, “Dynamical consequences of harvest in discrete age-structured population models,” Journal of Mathematical Biology, vol. 49, no. 1, pp. 35–55, 2004.

A. Wikan, “Four-periodicity in Leslie matrix models with density dependent survival probabilities,” Theoretical Population Biology, vol. 53, no. 2, pp. 85–97, 1998.

N. V. Davydova, O. Diekmann, and S. A. van Gils, “Year class coexistence or competitive exclusion for strict biennials?” Journal of Mathematical Biology, vol. 46, no. 2, pp. 95–131, 2003.