Abstract

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

1. Introduction

Within the framework of nonlinear discrete age-structured population models, the dynamical properties and behaviour of a great variety of species may be explored. Such species may possess an iteroparous life history which means that individuals in several age classes of the population are fertile, or they may possess a semelparous life history which is characterized by the property that only individuals of the last age class are fertile.

Regarding iteroparity, most population models (both age-structured and stage-structured) focus on cases where nonlinearities (density dependence) are built into the fecundity elements and not into the survivals. Particularly in fishery models, this has often been motivated by the assumption that most density effects are present only in the first year of life. Examples of theoretical studies which deal with nonstationary and chaotic behaviour as well as behaviour linked to concrete species may be found in [17]. In case of ergodic properties we refer to [810]; see also [11]. Another strategy is to assume constant fecundity terms and nonlinear year-to-year survival probabilities, compare [1216]. Depending on functional form of the nonlinear terms, a crude conclusion found in several of the papers referred to above is that, in case of small population densities, the models possess a stable nontrivial equilibrium where all age classes are populated. At higher population densities there may be nonstationary, periodic, and chaotic dynamics of stunning complexity.

Now, turning to the semelparous case with nonlinear fecundity element and constant survival probabilities a rather peculiar phenomenon has been detected, compare [4]. Here, in contrast to the iteroparous case, the nontrivial equilibrium tends to be unstable in large parameter regions, also in case of low population densities. Instead one finds, as time , that a cyclic state is attained where the whole population is in just one single age class at each time step. Such behaviour, which is called synchronization or SYC (single year class) dynamics, has been detected among insects; see [17]. Hoppensteadt and Keller [18] presented a model for the 17-year cicada (magicicada) which included both predation and intraspecific competition and in [19] cicada dynamics is further explored. Regarding biennials and possible SYC dynamics we refer to [20]. In [21] SYC phenomena and related MYC (multiple year class) dynamics are considered while Kon [22] discusses in a general context conditions for SYC dynamics to occur in matrix models, also compare [23]. However, if we take the opposite approach, constant fecundity term, and nonlinear year-to-year survivals probabilities, SYC dynamics appear to be a rare event; see [16]. Instead, the nontrivial equilibrium is stable whenever the population size is sufficiently small and the nonstationary dynamics has a strong resemblance of 4-cycles, either exact or approximate, the chaotic regime included.

In contrast to most of the papers quoted above, the purpose of this paper is to study the combined effect of nonlinear recruitment and nonlinear survival probabilities in semelparous population models. Several general results about possible SYC dynamics and stability properties of the nontrivial equilibrium in such problems may be obtained in [23] (two age classes) and [24] (three age classes). In case of the more general setting where an arbitrary number of age classes is considered we refer to [25, 26]. Regarding our study, we shall restrict the analysis to the case where both the fecundity and the survival elements depend on the total population and, moreover, the functional form of these elements is of Ricker type. As we prove, under such restrictions, it is possible to obtain explicit thresholds for secondary bifurcations of flip and Neimark−Sacker type and we may also prove (in some cases) that bifurcations involved are of supercritical nature. We may also investigate how the dynamics reported earlier will change as more density-dependent terms are included. For example, given SYC cycles, will the cycles persist if the strength of density dependence in the survival terms is included? If not, what kind of qualitative dynamics is it then possible to obtain? Assuming cycles of low period where all age classes are populated, does the inclusion of density dependence in the fecundity terms act in a stabilizing or destabilizing fashion? Such questions are difficult to address in the general models presented in [2426].

The paper is organized in the following way. In Section 2 we present the model, compute equilibria, and derive the th order eigenvalue equation which we need in order to perform stability and bifurcation analysis. This is followed (Sections 3 and 4) by a rigorous analysis of possible dynamic outcomes in two and three age class models, respectively. Finally, in Section 5, we unify and discuss results when the number of age classes exceeds three.

2. The Model, Fixed Points, and Stability

First we establish the model. At time we split the population into distinct nonoverlapping age classes where the total population is given by . Next, we introduce the transition matrixwhere is the average fecundity of a member of the last age class at time . denotes the survival probabilities, the (year-to-year) survival from age class to . In contrast to most papers the assumption here is that both the fecundity and the survivals are nonlinear terms. Thus, we write as where the constant and the survivals as where the constant satisfies . The parameters , may be regarded as parameters that measure the strength of density dependence. The relation between at two consecutive time steps is then expressed aswhich may also be formulated as a nonlinear map of the formBesides the trivial fixed point maps (2a) and (2b) also possess a unique nontrivial point . The latter may be expressed asThe quantities and are defined aswhere and is assumed throughout the paper in order to have a feasible equilibrium. The total equilibrium population is given asIn order to investigate stability we linearize about the fixed point. This gives birth to the th order eigenvalue equationwhere the coefficients satisfy is a locally stable hyperbolic fixed point as long as all eigenvalues of (6) are located inside the unit circle in the complex plane.

There are three ways in which may fail to be stable. It may lose its hyperbolicity when crosses the unit circle through 1 which in the general case leads to a saddle node bifurcation, alternatively through −1 which gives birth to a flip (period doubling) bifurcation, or it may fail to be hyperbolic as a pair of complex-valued eigenvalues cross the unit circle. Then a Neimark−Sacker bifurcation occurs. The Jury criteria, see the book by Murray [27], provide conditions for all eigenvalues to satisfy .

3. Two Age Classes

Let in maps (2a) and (2b). Then we have, , and . Moreover, the fixed point becomesand the eigenvalue equation may be cast in the formwhere and .

Fixed point (9) is stable whenever the Jury criteria , , and hold, that is, as long asrespectively.

There are two cases to consider: (A) the case , which means that the strength of density dependence in the fecundity is stronger than the strength of density dependence in the survival, and (B) the case .

Considering (A), it is clear from (11b) that there does not exist any stable fixed point. Moreover, since (11b) fails as an eigenvalue crosses the unit circle through −1 it is natural to search for a 2-cycle which should be stable provided is small. Evidently, such a 2-cycle must be obtained fromand here there are two possibilities:(1) which leads to the trivial 2-cycle where the unstable fixed point is the only point the cycle.(2)The points are on the form or . In this case it follows from (12) that and must satisfy the equationsand by finding from (13b) and substitute back into (13a) we arrive atGeometrically, it is now easy to see that the graph of the left hand side of (13c) and that of the right side have a unique intersection point lying in the first (positive) quadrant. In the special case we obtain and . Hence, there exists a 2-cycle on the form (SYC form)and as shown in [4] this cycle is stable in case of small ( fixed). In Figure 1(a) we show an orbit starting at which settles on the 2-cycle (14). If we continue to increase , we find that (14) goes unstable and cycles of period 2k. , are established through successive flip bifurcations. These cycles, which are all on SYC form, are stable in smaller and smaller regions as is increased. Eventually, the dynamics becomes chaotic but we emphasize that it is on SYC form also in the chaotic regime. These scenarios are demonstrated in Figures 1(b) and 1(c). Actually, we have not accounted for what happens when (or ). Here, compare (10), the eigenvalues are 1 and −1, respectively, and both the positive equilibrium and the SYC 2-cycle bifurcate forward. For proofs and details we refer to [4, 2326].

When we observe much of the same SYC dynamics as in the case. However, we may in a sense argue that an increase of acts in a stabilizing fashion. Indeed, if and , map (8) generates chaotic dynamics. On the other hand, if and ( in both cases) the outcome is a stable period 4-cycle on SYC form as shown in Figure 2.

Since the boundary of the positive cone is always invariant for semelparous Leslie matrix models of any dimension, see [26], initial conditions of the form or in the 2-dimensional case will always produce SYC dynamics. However, if and or the dynamics occurs in the vicinity (mostly as a stable 2-cycle, not on SYC form) of the unstable fixed pointwhere . Hence, does not necessarily imply SYC dynamics although it is the most likely outcome. A further discussion is postponed to Section 5.

Next, consider (B) (the case ). Our first observation is that the Jury criteria (11a) and (11b) will never be violated. Moreover, in case of sufficiently small equilibrium populations the left hand side of (11c) will be positive as well. Consequently, there exists a region in parameter space where (9) is a hyperbolic stable fixed point. However, if we increase the value of (by increasing ) such that (11c) turns into an equality, undergoes a Neimark−Sacker bifurcation, loses its hyperbolicity, and becomes unstable at the thresholdwhile the corresponding modulus 1 solutions of the eigenvalue equation (10) may be cast in the formAs is well known, bifurcations may be of both supercritical and subcritical nature. If a fixed point shall undergo a supercritical bifurcation it means that an eigenvalue (pair of eigenvalues) must cross the unit circle outwards at instability and in the Neimark−Sacker case that an attracting quasiperiodic orbit restricted to an invariant curve is created beyond instability threshold. Now, considering our bifurcation, we first express map (8), using the abbreviation , aswhere and contain second- and third-order terms of and . (Details of how (18) is derived and explicit formulas of and may be found in Appendix.)

Following Wan [28], the bifurcation will be supercritical ifis negative and that at bifurcation. The quantities in (19) are defined asDue to the complexity of and it is out of reach to compute the sign of in the most general case. However, it works in the important special case . Indeed, (18) simplifies towhere andMoreover, becomes zero so may be found aswhich is clearly negative.

Finally, at threshold , compare (11c) or (16)Hence, we have proved that the Neimark−Sacker bifurcation is of supercritical nature.

Still, assuming , let us now focus on the dynamics on the invariant curve. Whenever and is small where is the value of at threshold (for a given value of , is found from , we find a quasiperiodic orbit which fills the invariant curve with points; see Figure 3(a). In case of larger values of we observe a stable 4-cycle, comapre Figure 3(b), and as we continue to increase stable orbits of period , , are the outcomes. (Note, however, that the points in these orbits are clustered in such a way that one from an observational point of view may argue that we have almost 4-periodic dynamics in these cases as well.) The smaller is, the larger the interval becomes where the 4-periodic structure occurs. Eventually, in case of even larger fecundity values, the dynamics becomes chaotic, but even in the chaotic regime a certain kind of 4-periodicity is preserved in the sense that the chaotic attractor is divided into 4 disjoint subsets that are visited once every fourth iteration; see Figure 3(c).

The reason why we have all this 4-periodicity may be understood along the following line. Once the invariant curve is established we may regard map (8) restricted to the curve as topological equivalent to a circle map. Associated with a circle map there is a rotation number , and whenever is rational we have an orbit of finite period. If is an irrational number the orbit is quasiperiodic. Now, at bifurcation threshold the modulus 1 eigenvalues areand since is large at threshold, is located close to the imaginary axis in the left half plane. Following Guckenheimer and Holmes [29], gives asymptotic information about the rotation number; thus . This accounts for the 4-periodicity observed. Finally, note that when an orbit becomes exact periodic as a result of changing parameter , the implicit function theorem guarantees that there exists an interval around that specific parameter value where the periodicity is preserved as well.

Next, consider the case . The location of the modulus 1 eigenvalues is now given by (17) and they are not so close to the imaginary axis anymore. For “small” values of α we find the invariant curve, but the 4-periodic dynamics reported above is absent. Instead, as a result of increasing , the invariant curve becomes kinked and not topological equivalent to a circle. We may also describe the dynamics by use of the Lyapunov exponent of the orbit generated by (8) [30, 31]. In Figure 4(a) we display the values of in the range when and . Whenever , . In this interval the fixed point is stable. When the dynamics is quasiperiodic and restricted to invariant attracting curves and the corresponding Lyapunov exponent is . in the parameter window . Periodic dynamics of period 11 is the outcome here. Finally, when exceeds 32.7 and also in a tiny interval just below 30.5 we find that which means that the dynamics is chaotic. These findings are also visualized in Figure 4(b) where the dynamics is shown for selected values of . From bottom to top we recognize stable fixed points, invariant curves, kinked curves, 11-periodic dynamics, and chaos.

In case of intermediate values of , we observe a significant change of dynamics. Here, there exists a critical value , , where the third iterate of map (8) undergoes a saddle node bifurcation which results in two large amplitude 3-cycles, one stable and one unstable. Thus, in the interval we find coexistence between two stable attractors, the stable fixed point (9) and the stable 3-cycle. Consequently, the ultimate fate of an orbit depends on the initial condition but since the trapping region of the 3-cycle appears to be (much) larger, 3-periodic dynamics is the most likely outcome. Beyond there is an interval where the invariant curve established at and the stable 3-cycle coexists. At the invariant curve disappears as it is hit by the unstable 3-cycle, and if , small, the 3-cycle is the only attractor. Similar phenomena have also been found in iteroparous population models, first by Guckenheimer et al. [1], later in [4]. At even higher fecundity values the stable 3-cycle undergoes a Neimark−Sacker bifurcation which leads to three invariant curves which are visited once every third iteration. This is followed by kinked curves, periodic dynamics, and chaos through further enlargement of . Values of as well as the dynamics reported above are shown in Figures 4(c) and 4(d).

If is large but , we find that the dynamics qualitatively is quite similar to the case where is small, compare Figures 4(e) and 4(f) where and . However, there are several parameter values where the dynamics is periodic but most of these have almost no widths. The exception is the first window where the dynamics is 5-periodic. A final comment is that the larger is, , the higher is at bifurcation threshold (16), so one may argue that the strength of density dependence in the fecundity acts in a stabilizing way.

4. Three Age Classes

In this section we study the mapThe nontrivial fixed point (cf. (3)) may be expressed aswhere , , , , and .

From (6) it follows that the coefficient in the eigenvalue equationmay be written as , , and .

Equilibrium (27) is locally asymptotic stable as long as the Jury criteria , , , and hold. Hence, the stability criteria areA final observation from (28) is that when all eigenvalues , , are located on the boundary of the unit circle. Thus, if all eigenvalues move inside the unit circle when is increased, then (27) is a hyperbolic stable equilibrium in case of small. If at least one eigenvalue moves out, (27) is unstable. The Jury criteria (29a), (29b), (29c), and (29d) will help us to decide. Whenever (29c) and (29d) fail, a pair of complex valued eigenvalues will have modulus larger (or equal to) than 1 and consequently have a location outside (or on) the unit circle. If (29a) and (29b) fail, an eigenvalue crosses the unit circle through 1 or −1, respectively.

Moreover, due to the complexity of the Jury criteria we limit the discussion in the rest of this section to . Then the Jury criteria may be expressed asCriteria (30a), (30b), and (30c) obviously hold in case of small equilibrium populations . Regarding (30d) the same is true only as long as .

(A) First, consider . Then (30d) fails whenever is small but it holds for those which satisfyOn the other hand, (30b) holds only as long asTherefore, if satisfies the inequality , then is stable. Moreover, if this shall be the case, the difference must be strictly negative which is equivalent to say thatThe most transparent case to discuss is . Then , , and from (32) we find in order to ensure that . Now, if we choose a value larger than , say , we find , and finally if we obtain . Hence, there exists a parameter region where the fixed point is stable.

Turning to the dynamics, SYC dynamics is the only outcome in the region where (30d) fails. In case of small values there exists a stable 3-cycle where the points in the cycle are , , and . For higher fecundity values there are cycles of period as well as chaotic dynamics. In the intermediate region where all criteria (30a)–(30d) are valid, we find coexistence between the fixed point and chaotic SYC dynamics. Whether an orbit shall converge towards or settle on a chaotic attractor on SYC form depends on the initial condition. If we continue to increase we find that (30b) fails and as a result the fixed point will undergo a (supercritical) flip bifurcation. This gives birth to a tiny parameter region where SYC dynamics coexists with a stable 2-cycle not on SYC form. Through further enlargement of , SYC dynamics is the only outcome.

In [4] it was proved that there also exists a 3-cycle where two age classes are populated at each time. Within our framework the points in such a cycle are found to be , , and where In [21] it is proved that the cycle is unstable and further, in a more general context, it is shown in [24] that such 2-class cycles may be embedded in a heteroclinic cycle which may be attracting. Indeed, depending of the kind of competition within and between age classes, the analysis in [24] accounts for possible loops and cycles. For example, if the competition is asymmetric the bifurcating invariant loop consists of a single year class 3-cycle with a synchronous two-year class orbit that heteroclinically connect the phases of this 3-cycle. Thus, all orbits on such a heteroclinic cycle approach the single year class 3-cycle. This mechanism was also found in Bulmer’s original work on periodical cicadas. For further reading, compare [17, 24]. A final point is that when , the eigenvalue equation may be cast in the formwhere the dominant roots are . This means that independent of population size there will always be 3-cyclic dynamics in the sense that the total population at every point in the cycle equals but the structure of the points is on SYC form.

Turning to the case , it is still possible to find combinations of and such that the fixed point is stable but the parameter interval where this occurs appears to be (very) small. Indeed, assuming , the choice satisfies (32) and if it is possible to find such that as well as . However, if we change to or 0.85 it does not work. Therefore, our conclusion is that SYC dynamics dominates almost completely in the region .

(B) Next, consider . In this case all criteria (30a)–(30d) hold whenever is sufficiently small. Thus, there exists a parameter region where (27) is stable. Instability is introduced through a Neimark−Sacker bifurcation when (cf. (31a)) or equivalently whenRegarding the eigenvalue equation (28), coefficients and are always positive and at threshold (34), which is positive too. Consequently, there are no changes of signs between the coefficients in (28) which implies that there are two complex modulus 1 solutions and one real negative solution which necessarily must be and clearly . These findings allow us to express (28) asand the location of eigenvalue at threshold is given by the complex roots.

Now, scrutinizing the special case , , we find that instability threshold (34) becomesand moreover that the complex modulus 1 solutions of (35) may be expressed asFrom (36) it follows that is large at bifurcation threshold; thus the eigenvalue (37) is located in the left half plane close to the imaginary axis (in fact even closer than in the corresponding 2-age class model); therefore here too. In order to visualize the dynamics we have studied the case in somewhat more detail. In Figure 5(a), where , we show a quasiperiodic orbit which is restricted to a invariant curve. When , see Figure 5(b), an exact 4-period orbit is established through the frequency locking mechanism. In case of higher values, the fourth iterate of (24) undergoes a Neimark−Sacker bifurcation which results in four invariant curves as displayed in Figure 5(c)  . Through further enlargement of the dynamics becomes chaotic, see Figure 5(d) where .

Next, assume . Whenever is small, , we observe qualitatively the same kind of dynamics as just reported, that is, four periodicities either exact or approximate beyond instability threshold. Through further increase of the tendency towards 4-periodic dynamics becomes less pronounced and gradually disappears completely. Moreover, we experience that the fixed point may be stable for much higher fecundity values than when α is small. The nonstationary dynamics is restricted to invariant curves, not topologically equivalent to circles, compare Figure 5(e), or chaotic. The curves appear to be weakly attracting.

5. Discussion

In the previous sections we have analysed different versions of two and three-dimensional maps. We shall now consider the more general situation where there are age classes and since the parameter space is huge we will limit the discussion to the case where we have the same “strength” of survival between any two age classes. Hence, we let which means that the eigenvalue equation (6) may be cast in the formFirst, let us comment on the parameter region where . Whenever is small the fixed point is always unstable. This is proved in a general setting in [26] which also applies to the case under consideration in this paper. Next, suppose that is even and . Then left hand side of (38) may be expressed as which is clearly ≤0. On the other hand, when , then left hand side of (38) ; thus there must be a root of the equation which actually proves that the fixed point will always be unstable. When is odd the argument presented above does not hold, compare our analysis of the three age class model. However, it is possible to show that the fixed point is unstable given that is sufficiently large. To this end, assume that is odd and . Then left hand side of (38) is positive as long asand we observe that when , too, so evidently there must exist an eigenvalue where left hand side of (38) is zero. Hence, there is no stable fixed point in this case.

As the number of age classes increases there may also be cycles where more than one age class is populated at each time. Such cycles, referred to as multiple year class, MYC, cycles have already been identified in the case . In [21] where no density dependence in the survivals is assumed, MYC dynamics is further discussed and a major finding is that stability properties of such cycles on the whole are similar to the stability properties of (i.e., unstable). Consequently, assuming , we conclude that the nontrivial fixed point (3) will always be unstable except for small parameter windows when is odd. As demonstrated, such windows may be hard to find when . Regarding the dynamics, the general case is SYC dynamics. When is small, there is an attracting cycle on SYC form, as is increased we observe cycles of period , , as well as chaotic dynamics. In parameter windows where (3) is stable there is coexistence between the fixed point and chaotic dynamics on SYC form. Considering the simplest case the cycle may be written asAn increase of acts in a stabilizing fashion in the sense that the higher is, the higher value is necessary in order to generate chaotic dynamics.

From a biological point of view it is not obvious how one should interpret SYC dynamics. As far as we know, Bulmer [17] is the first who has noticed SYC dynamics in theoretical insect models. He explains its presence by saying that competition between age classes is more severe than competition within age classes. Bulmer’s argument is further strengthened by the findings in [24], see Theorem , where it is shown mathematically that strong competition within age classes in general will lead to an equilibrium where all age classes are populated (i.e., an equilibrium of form (3) in our model ) while strong competition between age classes destabilizes and promotes oscillations with missing age classes. Another argument is presented in [21]. However, it should be emphasized that the conclusion obtained there is found from a simple model.

In the remaining part of the paper we shall deal exclusively with the case . Denoting left hand side of (38) for it is immediately clear from our previous analysis that when is even the Jury criterion will always be satisfied. Hence, there is no flip bifurcation in this case. When is odd the value of at threshold is found to beand if we substitute back into (38) the constant term becomeswhere is the denumerator of (41). Moreover, since we may rule out the flip here too. Consequently, the only possible transfer from stability to instability goes through a (supercritical) Neimark−Sacker bifurcation for any .

There are two more questions which we find natural to discuss. (A) Will the 4-periodic dynamics observed when and also persist in case of larger values of and (B) what happens to the size of the stable parameter region when the number of age classes increases?

Regarding (A), assume and (the values of and with the most pronounced 4-periodic dynamics when and ). From the Jury criteria (and a lengthy calculation!) it now follows that the size of at instability threshold must satisfyThus, for given values of we may compute from (43), substitute back into (38) and solve in order to find the location of the modulus 1 eigenvalues at threshold. For example, if we find   (, ) while gives   (, ). Now, the crucial observation from calculations as above is that (as opposed to the and cases) are located in the right half of the complex plane and not particularly close to the imaginary axis. Hence, is not small, so the tendency towards 4-periodic dynamics is much weaker than in , cases, a result which has also been verified through several simulations. Therefore, our analysis clearly suggests that periodic dynamics of low period, in particular 4-periodic dynamics, is limited to species with a few number of age classes. If only the first age class is not fertile and density effects occur in the survival terms exclusively, we may find 4-periodic dynamics also for larger values of ; see [15].

Finally consider (B) and let , , and be the size of at bifurcation threshold in the and 4 cases, respectively. From (16), (31a), (31b), and (43) it now follows that , , are increasing functions of , , and moreover, that , , and . These findings are presented in Figure 6. Hence, for a given value of , which means that the value of at threshold becomes smaller as the number of age classes is increased; that is, an increase of acts as a destabilizing effect. This argument is further increased if we continue to increase . Indeed, when we have verified numerically that the fixed point is unstable in case of sufficiently small. Only quasiperiodic orbits have been identified. In [32] where a discrete stage-structured population model was analysed, the authors concluded that species who possess delayed semelparous life histories tend to be more stable than species who possess precocious semelparous life histories. Results from the analysis of age-structured models in [16] show a tendency in the opposite direction. The findings in this paper support the latter, but we feel that more work has to be done before one may conclude this issue.

Appendix

Considering map (8), , we find that the Jacobian at bifurcation threshold (11c) or (16) may be written asand as shown in the main textwhere . The corresponding eigenvector becomesFurther, define the matrix (cf. (A.3)) Then, after expanding the components of (8) up to third order, translating the bifurcation to the origin through the change of coordinates , and finally applying the transformationsmap (8) may be cast into standard form at the bifurcation as (18) where and are defined asrespectively.

Conflicts of Interest

The author declares that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The author thanks Ørjan Kristensen for computer assistance. Financial support from the publication fund of UiT The Arctic University of Norway is also gratefully acknowledged.