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

Research Article

An Analysis of Discrete Stage-Structured Prey and Prey-Predator Population Models

Arild Wikan

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 23 February 2017; Accepted 27 March 2017; Published 16 April 2017

Academic Editor: Hammad Khalil

Copyright © 2017 Arild Wikan. 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.

Discrete stage-structured prey and prey-predator models are considered. Regarding the former, we prove that the models at hand are permanent (i.e., the population will neither go extinct nor exhibit explosive oscillations) and, moreover, that the transfer from stability to nonstationary behaviour always goes through a supercritical Neimark−Sacker bifurcation. The prey model covers species that possess a wide range of different life histories. Predation pressure may both stabilize and destabilize the prey dynamics but the strength of impact is closely related to life history. Indeed, if the prey possesses a precocious semelparous life history and exhibits chaotic oscillations, it is shown that increased predation may stabilize the dynamics and also, in case of large predation pressure, transfer the population to another chaotic regime.

1. Introduction

As it is well known, discrete nonlinear age- and stage-structured population models serve as excellent tools in order to study the dynamical outcomes of various ecological populations. Such models contain species which may possess different life histories ranging from biennials to species that may live for many years. Regarding the age-structured case, examples of studies on concrete species as well as pure theoretical approaches may be obtained in [1–7]. In the stage-structured case, we refer the reader to [8–13]. A comparison of dynamic outcomes from age- and stage-structured models may be obtained in [14], and the analysis of the models that also incorporate spatial structure may be found in [15, 16].

The models referred to above may also be extended to include prey-predator interactions (see [17–22]). Some of these models are two-dimensional in the sense that neither the prey nor the predator has an internal structure, while others may, for example, be four-dimensional. Considering the latter, Wikan [21] provides an analysis of several examples where both the prey and the predator populations are divided into two age classes.

The purpose of this paper is to bring out the analysis of a case where the predator preys upon newborns only. Therefore, we have developed a three-dimensional model where we split the prey population into two separate parts by use of the same strategy as in [10], while the interaction between newborn prey and the predator is modelled in the same way as in [18]. The focus is on the stability and nonlinear behaviour. In particular, we address the question how the impact from the predator acts on prey populations that possess different life histories. The results from analysis of the prey part of the model are also provided and compared with findings obtained from other models.

The structure of the paper is as follows: In Section 2, we present and analyse the various models with respect to stability and nonstationary behaviour. Section 3 provides examples of prey models as well as prey-predator interactions while in Section 4 we unify and discuss the results.

2. Model(s), Fixed Points, and Stability

Let and be the immature and mature part of the prey population at time , respectively, and let be the predator population. and , such that and , are the fractions of the immature population and the mature (adult) population that survive from time to time . , with , is the fraction of the immature population that survives to become mature one time unit (year) later and , , is the density dependent fecundity. Depending on the species under consideration, may account for crowding effects, effects linked to shortage of food, and for some species it may also incorporate cannibalistic behaviour. Further, it is assumed that predation will take place only on the young of the part year, , of the prey population. This is accounted for by the term where measures the skill of predation. The constant , , may be interpreted as a conversion of prey into predator, or clutch parameter, the following year (see [18]). The relation between , and at two consecutive time steps may then be given as a map:

When , (1) degenerates to the prey map:

Note that (2) has a striking similarity to the general stage-structured model proposed and analysed in [10]. The difference is found in the density dependent term.

By adjusting the size of the parameters, it is easy to see that (2) covers species that possess different life histories. Indeed, following [10], if , the population is semelparous (i.e., reproducing only once). If , the population is iteroparous (repeated reproduction). The subclass , is often referred to as precocious semelparity which covers species with rapid development followed by only one reproduction, for example, biennials and annual plants (see [11]). Delayed semelparity occurs when and . Typical examples are periodical cicadas [11, 23], and several salmon species that live for many years before they become mature and reproduce only once. We may also divide the iteroparous case into two subclasses. The subclass , , is classified as precocious iteroparity and covers several small mammals species, among them small rodent species that start to reproduce at young age and may survive to reproduce for several years. The fourth subclass, delayed iteroparity, is characterized by , , which covers species which may live long before maturity and then survive to reproduce for many years. In this subclass, we find large mammals. Consequently, (2) may be used in order to capture the dynamics of a wide range of (prey) populations, and the role and impact of predation on newborns and young individuals are then analysed by (1).

We start by revealing some properties of (2). There are two fixed points, the trivial one and the nontrivial one:

where

and in order for (3) to be a feasible fixed point (equilibrium), we assume . Moreover, by use of stability analysis, it is straightforward to show that is stable provided . Therefore, the restriction ensures both that the origin is a repeller and that (3) is feasible.

Model (2) is said to be permanent if there exist and such that

(cf. [24]).

Theorem 1. Model (2) is permanent provided .

Proof. See Appendix A.

Considering the stability properties of the nontrivial fixed point (3), we find that the eigenvalue equation of the linearization of (2) may be expressed as

where

and according to the Jury criteria (cf. [25]), (3) is stable provided the inequalities , , and hold; that is,

Obviously, the left-hand sides, LHS, of (8a) and (8b) are positive. Therefore, there will be no transfer from stability to instability as an eigenvalue crosses the boundary of the unit circle through 1 (8a) or −1 (8b). Regarding (8c), it is clearly valid for small values of but when is increased (as a result of increasing ), LHS of (8c) eventually equals the right-hand side, RHS, and will lose its hyperbolicity through a Neimark−Sacker (Hopf) bifurcation as a pair of complex valued eigenvalues cross the unit circle. Such a bifurcation may be of supercritical or subcritical nature. In the former case, when the fixed point loses its stability, an attracting invariant curve about the unstable fixed point is established and the dynamics are restricted to that curve. In the latter case (subcritical), there is no such attracting curve. Now, considering (2), we have the following result.

Theorem 2. Consider (2) together with the fixed point given by (3). Then, for the fixed values of and , (3) will undergo a supercritical Neimark−Sacker bifurcation at the threshold

or equivalently when

Proof. See Appendix B.

Next, let us focus on the “full” prey-predator map (1). There is one obvious fixed point, namely, the trivial one . The other fixed point is in the form of

where must be found by means of numerical methods from the equation

Clearly, if , (11) implies that so one possibility is that (10) is in the form of . Depending on the value of , the same scenario persists also in the case of . Consequently, as also found in the prey-predator model analysed in [21], the interaction or skill parameter a must exceed a critical threshold in order to establish a fixed point where both species coexist. Moreover, note that implies . If , an increase of a makes larger which according to (10) leads to a reduction of and . A final observation from (10) and (11) is that a decrease of clutch parameter leads to a smaller predator equilibrium . This makes sense; the smaller the , the smaller the benefit of eating. A couple of numerical examples are presented in Tables 1 and 2. In Table 1, we have used the parameter values , , , , and which means that the prey possesses precocious semelparous life history and that the prey in the absence of the predator has a stable fixed point . In Table 2, the prey may be classified as a precocious iteroparous population. Parameter values are which means that is located at instability threshold.

Table 1

0.4
0.5
0.56
0.6
0.7
0.76
0.80
1
1.5
2
3

≈0
≈0
0.512
1.286
1.655
1.682
1.674
1.530
1.148
0.901
0.625

≈0
≈0
0.287
0.772
1.158
1.278
1.339
1.530
1.723
1.802
1.875

≈2
≈2
1.713
1.228
0.842
0.722
0.661
0.470
0.277
0.198
0.125

≈2
≈2
1.713
1.228
0.842
0.722
0.661
0.470
0.277
0.198
0.125

Table 2

1.00
1.05
1.10
1.20
1.30
1.40
1.50
1.60
2.00
3.00

≈0
1.384
1.705
1.825
1.804
1.749
1.682
1.614
1.364
0.961

≈0
1.453
1.875
2.190
2.345
2.448
2.523
2.582
2.729
2.882

3.11
1.657
1.235
0.920
0.765
0.662
0.587
0.528
0.381
0.228

1.729
0.921
0.687
0.511
0.425
0.368
0.326
0.294
0.212
0.127

Now, let us turn to stability. The eigenvalue equation of the linearization of (1) may be expressed as

where the coefficients are

As long as all eigenvalues of (12) are located within the unit circle, (10) is stable and this is ensured if the Jury criteria

hold.

Now, by use of the relation and map (1), we find that criteria (14a)–(14d) may be cast in the forms

and our first observation is that the clutch parameter drops out of the criteria. Thus, different values of correspond to different values of , and but it does not affect qualitative changes of dynamics.

Considering the criteria, assuming , we find that LHS of (15a) is an increasing function of provided

The minimal value is found to be

and the sign as well as the value of (17) depends on .

Turning to (15b), we find that LHS is a decreasing function since is decreasing and as well as is increasing. Further, observe that in case of the LHS approaches the positive expression

On the other hand, if , the LHS may be expressed as

which is negative provided is large enough.

Scrutinizing (15c), the same picture emerges. LHS is a decreasing function and when is small, LHS approaches the positive expression

When becomes large, LHS degenerates to

which is negative. Now, assuming equality in (15c), we find

and by substituting into (15b) we arrive at

which holds. Consequently, LHS of (15c) crosses the axis before LHS of (15b) so we may rule out the possibility of a period doubling bifurcation at instability threshold.

Finally, turning to (15d), when , the LHS becomes negative (notice the term in the last fraction). In the case , the LHS may be written as

This expression is positive for small values of and approaches zero when becomes large.

Thus, to summarise, if is so small that the prey in the absence of the predator possesses a stable fixed point , all LHS of (15a)–(15d) are positive. Whenever the population pressure is small, the long term dynamics is a stable fixed point . As predation pressure ay increases, the magnitude of the dominating eigenvalue(s) becomes smaller until the graph of LHS of (15a) intersects the corresponding graph of (15d) or (15c). Then, the magnitude increases and the fixed point undergoes a Neimark−Sacker bifurcation when the LHS of (15d) or (15c) becomes zero. If is so large that the prey in the absence of the predator exhibits nonstationary dynamics, (17) is negative in case of small; therefore, is unstable. However, when increases, we observe the same qualitative picture as reported above.

3. Examples

First, we concentrate on dynamics generated by the “prey map” (2) and we start with a numerical example where , , and (precocious semelparity). At threshold (9b), (or ) and the fixed point undergoes a Neimark−Sacker bifurcation. In Figure 1(a), we show the attracting invariant curve together with some initial transients in case of . When is further increased , the curve becomes kinked and not topologically equivalent to a circle anymore and in case of the curve has broken up into a diffuse cloud of points. Thus, through an enlargement of (or ), the dynamics change from stability to chaotic behaviour. The whole scenario is displayed in Figures 1(a)–1(c).

(a)

(b)

(c)

Figure 1: Dynamics generated by (2). Parameter values , , and . (a) , (b) , and (c) .

If we instead use , , and (delayed semelparity), at threshold and the corresponding fixed point is . When , we find invariant curves here too, and such a curve also persists when 0. We have not detected chaotic dynamics as in the previous case. Hence, a natural conclusion to suggest is that populations that possess delayed semelparous life histories have better stability properties than populations with precocious semelparous life histories. This is in excellent agreement with results obtained by Neubert and Caswell [10] but not with the findings in [14] where an age-structured model is analysed.

Turning to the precocious iteroparous case, exemplified by , , and , we find and at instability threshold. Beyond threshold, the behaviour is qualitatively similar to the precocious semelparous case.

Finally, considering delayed iteroparity through the example of , , and , we find from (9b) that at threshold and the corresponding fixed point becomes . As we increase beyond 45, there are invariant curves with the same kind of shapes as in the former case. The main difference between the precocious and delayed case is the size of equilibrium population at threshold. Hence, by use of the same classification as in [10], we find it natural to support the conjecture that species that possess a delayed iteroparous life history have better stability properties than populations with precocious iteroparous life histories.

Next, consider the “full” prey-predator map (1). We start our analysis by use of the same parameter values as we did when Table 1 was established (i.e., we study a case where the prey possesses a precocious semelparous life history). As already shown, implies that is stable in the absence of predation. Then, we increase (starting at , cf. Table 1) and calculate the LHS of (15a)–(15d), respectively, and the results of these calculations are presented in Figure 2. Clearly, is stable in the interval . It fails to be stable when LHS of (15d) becomes zero and a (supercritical) Neimark−Sacker bifurcation takes place. Just beyond 1.43 the nonstationary dynamics are restricted to an invariant curve, but through further enlargement of the curve becomes kinked and eventually it breaks up. This is shown in Figures 3(a)–3(c). Thus, in case of small and moderate predation stability properties are improved. On the other hand, severe predation pressure may be classified as a strong destabilizing effect.

Figure 2: The graphs of LHS of (15a), (15b), (15c), and (15d). Parameter values: , , , and . LHS of (15a) is the increasing function. LHS of (15d), (15c), and (15b) cross the -axis at 1.44, 1.61, and 1.67, respectively.

(a)

(b)

(c)

Figure 3: Dynamics generated by (1). Parameter values: , , , and . (a) , (b) , and (c) .

The results referred to above are further strengthened through our next example. In Figure 4, we present the calculations of LHS of (15a)–(15d) in case of (which in contrast to the previous example means that is unstable in case of no predation). As is shown, there is no stable fixed point in the interval (cf. (16) and (17)). Otherwise, the graphs look similar to the graphs in Figure 2. In Figure 5, we display the dynamics. In Figure 5(a), we use , and Figures 5(b) and 5(c) show the cases and .

Figure 4: LHS of (15a)–(15d) in case of . Other parameter values as in Figure 2.

(a)

(b)

(c)

Figure 5: Dynamics generated by (1). Parameter values as in Figure 3 except . (a) , (b) , and (c) .

Both examples referred to above consider species that possess a precocious semelparous life history. We have also scrutinised species with other life histories and our findings are that there is a striking resemblance between the graphs of (15a)–(15d) in these cases and Figures 2 and 4. The main difference really is that ay must be increased to a much higher value in the delayed iteroparous case in order to establish an equilibrium with coexisting species.

4. Discussion

In the previous section, we have analysed a selected number of stage-structured prey and prey-predator models. The parameter space is huge which allows us to consider prey populations that possess a wide range of different life histories. We shall now unify and discuss our findings.

First, let us comment on the prey map (2). As is shown, independent of life history, (2) possesses a stable fixed point whenever is small enough and the transfer from stability to instability occurs when the fixed point undergoes a supercritical Neimark−Sacker bifurcation at thresholds (9a) and (9b). In Figure 6, we show the value of the total equilibrium population at instability for different values of . For a given value of , the stable region is located below the corresponding curve. Hence, whenever (roughly), populations which exhibit iteroparous life histories have better stability properties than species that possess semelparous life histories. On the other hand, if , we arrive at the opposite conclusion. Moreover, since all functions are decreasing, , we may also conclude that the delayed cases appear to be more stable than the precocious ones.

Figure 6: Graphs of at bifurcation threshold in case of , 0.5, and 0.9. . When , , and 0.9 imply and .

It is tempting to compare the findings above with the results obtained from the model:

which is analysed in [10]. Note that . Thus, the difference between our model and (25) is that while we consider cases where only the mature part of population contributes to density effects, the whole population contributes in (25). This has a pronounced impact on the results. Indeed, just as (2), the fixed point of (25) is stable in case of sufficiently small and loses its hyperbolicity at

where , but in contrast to our model, becomes unstable as the dominant eigenvalue of the linearization of (25) crosses the boundary of the unit circle through −1. Consequently, when exceeds RHS of (26), a (stable) two-period orbit is established and depending on life history, further increase of may generate periodic orbits of period , , as well as chaotic dynamics. Also, note that while all stability curves in Figure 6 are decreasing (independent of ), is an increasing function whenever is large (cf. [10]). Therefore, based upon the analysis of (25), Neubert and Caswell [10] conjectured that species with precocious iteroparous life histories possess more stable dynamics than species with delayed iteroparous life histories. Due to the findings from (2), we feel that this conjecture should be modified.

Let us now focus on possible periodic dynamics. Since fails to be hyperbolic as a pair of complex valued eigenvalues cross the unit circle, we may exclude periodic dynamics of period , . Moreover, as is shown, when exceeds and is small (cf. (9b)), the dynamics is restricted to an invariant curve. On that curve, (2) is topologically equivalent to a circle map which in polar coordinates may be expressed as

(see [26]), where gives asymptotic information of the rotation number associated with the circle map. Now, at threshold (9b), the solution of (6) may be expressed as

and here we may notice that if and is small (precocious semelparity), we roughly have arg which in turn implies that the rotation number is close to . This signals approximate 4-periodic behaviour and one may even expect exact 4-period orbits through frequency locking as is further increased. Examples of such behaviour may be obtained in [5] or [27]. However, this does not occur. We have not been able to detect periodic dynamics of low period (period 3 or 4) generated by (2). Therefore, based upon our findings here, lots of simulations, and indeed also findings from other stage- and age-structured population models, it appears that if only contributes to density effects, periodic dynamics are likely to be absent.

Next, let us turn to predation. The predator population must reach a certain size in order to establish a fixed point where both species coexist. Depending on parameter values, may be stable or unstable. Moreover, starting at a low value, an increase of acts in a stabilizing fashion, but when is large, further enlargement acts in a destabilizing fashion. The “turning” point is where graphs of (15a) and (15d) intersect. This scenario takes place independent of prey life histories.

However, the “strength” of increasing ay is strongly related to life history of prey. Regarding precocious semelparous life history, we have shown that a prey population which possesses nonstationary dynamics may first be stabilized (as shown in Figure 4) and then driven to chaos through an enlargement of predation pressure. In fact (through simulations), we have also verified that a prey population which in the absence of predators exhibits chaotic oscillations may be stabilized and then brought to another chaotic regime as predation pressure is increased. Thus, it is plausible to consider predation pressure as a strong (both stabilizing and destabilizing) effect. It should also be mentioned that the results above are in excellent agreement with the findings obtained from age-structured prey-predator models with few age classes (cf. [21]). In many ways, this is expected. Indeed, if becomes small and , one may argue that stage-structured models degenerate to age-structured models. One important difference between the outcomes of discrete age- and stage-structured models in the precocious semelparous case is linked to possible periodic dynamics. In the age-structured case, it is possible to show that prey that possesses periodic dynamics of low period may force the predator to oscillate with the same kind of periodicity. Such phenomena have not been detected in our model.

The impact of predation when the prey possesses a delayed, semelparous or iteroparous, life history is somewhat different in the sense that the destabilizing effect does not seem to be very strong. It is not obvious what causes this discrepancy. In the precocious case, the value of gets smaller but not very much smaller, compared to the delayed case when a is increased. Hence, when is large, the prey is exposed to a large predator population with excellent hunting skills. This is not the situation in the delayed cases. Although an increase of leads to an increase of here too, there is a substantial reduction of predators possessing good hunting skills.

Appendix

A. Proof of Theorem 1 in the Main Text

As already shown, if , then is unstable. Moreover, the restrictions on parameters given in (2) ensure that is irreducible and that is nonnegative for all . Consequently, (2) is forward invariant. We must show that there exists a compact set such that for all there exists a satisfying for all . To this end, let be a constant satisfying . Then, from (2) (using difference equation notation),

and by induction

Then, there exists such that for

Further, in case if , we also have

and once again (induction) we find that for

Finally, take and ; then, for , , , and we are done. This proof is based upon Kon et al. (2004).

B. Proof of Theorem 2 in the Main Text

The proof consists of two parts. First we show that the eigenvalues really leave the unit circle at threshold. To this end, note that (6) implies that

hence evaluated at threshold (9b)

Thus, the eigenvalues leave the unit circle.

Next, note that the Jacobian of (2) (evaluated at threshold) may be written as

with associated complex valued modulus 1 eigenvalues

where .

Define the matrix

whose columns are the real and imaginary parts, respectively, of the eigenvectors associated with the eigenvalues.

Then, after expanding the first component of (2) up to the third order, applying the change of coordinates (in order to transform the bifurcation to the origin) together with the transformations

map (4) may be cast into “standard form” as

where

and .

Now, by use of Theorem  3.5.2 in Guckenheimer and Holmes (1990), the bifurcation will be of supercritical nature if the quantity a defined through

is negative. For the problem at hand,

This yields

which is negative and we conclude that the bifurcation is supercritical.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

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

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.

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.

J. M. Cushing, B. Dennis, R. A. Desharnais, and R. F. Costantino, “An interdisciplinary approach to understanding nonlinear ecological dynamics,” Ecological Modelling, vol. 92, no. 2-3, pp. 111–119, 1996.

R. F. Costantino, R. A. Desharnais, J. M. Cushing, and B. Dennis, “Chaotic dynamics in an insect population,” Science, vol. 275, no. 5298, pp. 389–391, 1997.

M. G. Neubert and H. Caswell, “Density-dependent vital rates and their population dynamic consequences,” Journal of Mathematical Biology, vol. 41, no. 2, pp. 103–121, 2000.

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.

A. Wikan and A. Eide, “An analysis of a nonlinear stage-structured cannibalism model with application to the northeast arctic cod stock,” Bulletin of Mathematical Biology, vol. 66, no. 6, pp. 1685–1704, 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.

A. Wikan, “Age or stage structure? A comparison of dynamic outcomes from discrete age- and stage-structured population models,” Bulletin of Mathematical Biology, vol. 74, no. 6, pp. 1354–1378, 2012.

M. Gyllenberg, G. Söderbacka, 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.

A. Wikan, “On reserves, stability and the maximum sustainable yield problem,” Journal of Mathematics and Statistics, vol. 9, no. 4, pp. 325–333, 2013.

J. R. Beddington, C. A. Free, and J. H. Lawton, “Dynamic complexity in predator-prey models framed in difference equations,” Nature, vol. 255, no. 5503, pp. 58–60, 1975.

A. Hastings, “Age-dependent predation is not a simple process. II. Wolves, ungulates, and a discrete time model for predation on juveniles with a stabilizing tail,” Theoretical Population Biology, vol. 26, no. 2, pp. 271–282, 1984.

H. A. Lauwerier and J. A. Metz, “Hopf bifurcation in host-parasitoid models,” IMA Journal of Mathematics Applied in Medicine and Biology, vol. 3, no. 3, pp. 191–209, 1986.

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.

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

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.

J. D. Murray, Mathematical Biology: Spatial Models and Biomedical Applications, Biomathematics, Springer, Berlin, Germany, 3rd edition, 1993.

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

A. Wikan, “Dynamic consequences of reproductive delay in Leslie matrix models with nonlinear survival probabilities,” Mathematical Biosciences, vol. 146, no. 1, pp. 37–62, 1997.