The winner takes it all: how semelparous insects can become periodical

The aim of this short note is to give a simple explanation for the remarkable periodicity of Magicicada species, which appear as adults only every 13 or 17 years, depending on the region. We show that a combination of two types of density dependence may drive, for large classes of initial conditions, all but 1 year class to extinction. Competition for food leads to negative density dependence in the form of a uniform (i.e., affecting all age classes in the same way) reduction of the survival probability. Satiation of predators leads to positive density dependence within the reproducing age class. The analysis focuses on the full life cycle map derived by iteration of a semelparous Leslie matrix.


Periodical insects: a conundrum
According to Bulmer (1977), "An insect population is said to be periodical if the life cycle has a fixed length of k years (k > 1) and if the adults do not appear every year but only every kth year." He then provides several examples, one by quoting Lloyd and Dybas (1966) Periodical cicadas, Magicicada spp., have the longest life cycles known for insects. In any one population, all but a tiny fraction are the same age. The nymphs suck juices from the roots of deciduous forest trees in eastern United States. Mature nymphs finally emerge from the ground, become adults, mate, lay their eggs, and die within the same few weeks of every seventeenth (or, in the South, every thirteenth) year. Not one species does this, but three. There are three separate and distinct species that occur together over most of the range of periodical cicadas and, wherever the species coexist, they are invariably synchronized. In different regions, different "broods" of periodical cicadas may be out of synchrony by several years, but the species (in a given region) never are. Finally, the same three species-the same as nearly as anyone can tell by looking at them or listening to the songs of their males-exist both as 17 and as 13-year periodical cicadas! ; Note that the three species are now considered to be seven species, each with their own period of either 13 or 17 years (Marshall 2008)] In the wake of Bulmer's pioneering paper a rich literature on discrete time models for semelparous organisms developed, see (Behncke 2000;Webb 2001;Diekmann et al. 2005;Cushing 2009;Cushing and Henson 2012;Kon 2012Kon , 2017Cushing 2015;Wikan 2017) and the references given there. The population splits into year classes according to the year of birth (equivalently: the year of reproduction) counted modulo k. As a year class is reproductively isolated from the other year classes, it forms a population by itself. Accordingly, competition may lead to exclusion.
Despite this rich literature, the spatio-temporal pattern displayed by periodical cicada population dynamics is, we think, still an enigma. The interest of one of the authors in the subject was rekindled while refereeing an early version of (Machta et al. 2019). The model introduced and studied in that paper is characterized by -uniform (i.e., age-class unspecific) negative density dependence during development; -positive density dependence for reproduction.
The first of these captures competition for food, the second satiation of predators to eat adults in a given year (Williams et al. 1993).
The paper (Machta et al. 2019) considers the limit k → ∞ and employs a continuous time description of the competition for resources (whence the word "hybrid" in the title). The main result establishes the instability of any steady state with more than 1 year class present.
Here we focus on the "full life cycle" map, i.e., the k-times iterated Leslie matrix, cf. (Davydova et al. 2003. The full life cycle map is represented by a diagonal matrix. If competition is uniform, in the sense that in any year any of the survival probabilities is the product of a fixed age-class specific factor and a year-specific factor that is the same for all age-classes, the diagonal elements all have the same overall survival factor. In other words, negative density dependence results over a full life cycle for all year classes in exactly the same multiplication factor. Each diagonal element has an additional reproduction factor. Positive density dependence causes this factor to be highest for the year class that had, in the year that it reproduced, the highest density. In this way, a numerical advantage is reinforced and asymptotically there remains, "generically", only 1 year class: all the others go extinct. See Fig. 1 for an example. Note that this holds independently of the resulting single year class (SYC) dynamics in the sense that the winners may exhibit ultimately steady state, periodic or even chaotic dynamics. Also note that it depends on the initial condition which one of the k year classes will win the competition. Domains of attraction (or, dominance) are separated by sets of special initial conditions for which several year classes coexist forever (e.g., steady coexistence states and their stable manifolds). In principle these separatrices may have an intricate structure [e.g., see (Hofbauer et al. 2004)].
Our main result provides an explicit condition that guarantees that the year class having highest age at the census time will win. Please note that there are k year classes and that each is as 'strong' as any other one. The intrinsic symmetry is broken when we try to predict which year class is going to dominate on the basis of the abundances at a particular time, simply since the various year classes are at that time in a different phase of their life cycle. The main result yields, by iteration, an implicit partial characterization of the domain in the state space such that the year class that has age j, with 1 ≤ j < k, will win. The implicit character makes these less informative, yet we should keep in mind that the condition described in our main result is supplemented by k − 1 implicit conditions for victory of the other year classes. This raises the fundamental question: is the union of the domains of domination/attraction of the k year classes everywhere dense? In Fig. 2 we provide numerical evidence that the answer is 'yes' for at least some parameter values. In the "Appendix" we show that there are parameter values for which the answer is 'no'.
The upshot is that the combination of uniform negative density dependence and concentrated (in one point of the life cycle) positive density dependence, as assumed in (Machta et al. 2019), leads, as a rule, to exclusion of all but at most 1 year class. This does not prove, of course, that this is the mechanism underlying the observed phenomenon of SYC dynamics in Magicicada. But in the spirit of Occam's Razor, this probably currently offers the simplest, and hence most plausible, potential explanation.
Anyhow, as far as we are aware, we provide below the first analytical demonstration of 1 year class driving all other year classes to extinction without severe restrictions on the initial conditions and without any condition on the resulting SYC dynamics, except for boundedness.

Model formulation
Let k be an integer bigger than one. Let N(t) be the k-vector such that N j (t) is the sub-population density of individuals that at the census moment in year t have age j. In view of linear algebra we take j = 1, 2, . . . , k, even though biologically j = 0, 1, . . . , k − 1 would perhaps make more sense if, as we assume, the census moment is in the early autumn, so after reproduction took place. We assume that (1) where is a "semelparous" Leslie matrix (here the adjective "semelparous" indicates that in the first row only the last element is non-zero, reflecting that individuals reproduce at age k and then die). The dependence of the matrix elements h i on N captures density dependence. We assume that where Here the σ i are survival probabilities under "standard" conditions and the (uniform, i.e., i-independent) factor π(N) reflects the reduction of the survival probability as a result of crowding and competition for food. The last age group too needs to survive the winter before they can emerge as adults, so the population density of emerging adults equals σ k π(N)N k . The function β is a composite model ingredient. It describes the number of offspring at the census moment per emerging adult. The monotonicity of β reflects that the per adult capita risk of falling victim to predation decreases with adult population density. This is the effect of predators becoming satiated (or even over-satiated) when prey density is very high, and that the time interval between emergence and death of the adults is too short for a numerical increase in the number of predators (mainly birds).

The main results
So far we did not specify any property of π that justifies the description "negative" density dependence. But now we require that population densities remain bounded and implicitly this is an assumption concerning the function π .
Hypothesis 1 (dissipativity (Hale 1988)) There exists R > 0 such that is forward invariant and attracts all orbits.
By "attracts all orbits" we mean that for every . We adopt Hypothesis 1 for the rest of the paper. In Sect. 4 below we provide easily verifiable conditions on π and β that guarantee that Hypothesis 1 is satisfied.
At the end of Sect. 1 we already observed that (due to symmetry, see Davydova et al. 2005) for a bit more detail) any year class is a candidate for becoming the sole inhabiter of the world. It turns out to be relatively easy to describe a set of initial conditions such that the year class that has age k at the initial time will outcompete all the other year classes.
The underlying idea is the following. During a full life cycle, a year class has highest density at the first census after birth and lowest density at the last census before reproduction, simply since inbetween the density is reduced by mortality. Likewise, for a steady state with all year classes present, abundance decreases with age. (Possibly the same holds for any solution of period k with several year classes present, but we do not know.) The age classes at a particular census, on the other hand, reflect the relative abundances of the year classes, albeit in a manner that does not allow general straightforward meaningful comparison. However, if the highest age class has highest abundance, this clearly shows that the corresponding year class is dominant and on the way to becoming the winner. As shown in the following theorem, this test can be refined by correcting for the predictable mortality (as captured by the σ 's) before reaching the highest age.

Theorem 1 The set
is forward invariant under the full life cycle map, and if N(0) ∈ Ω k and for given is a strictly decreasing function of m. This function has limit zero for m → ∞ in case lim sup m→∞ N k (mk) > 0.
. Left: full dynamics of all cohorts; middle: the dynamics of the full life cycle map, for all cohorts, shown from the year at which the orbit enters Ω 3 ; right: dynamics generated by the full life cycle map, with Ω 3 in black, with red points of the orbit outside Ω 3 and blue those inside. We observe convergence to a single year class fixed point inside Ω 3 that lies on the N 3 -axis. Note that the choice σ 1 = σ 2 = σ 3 = 1 is no restriction, as we can scale the variables N k with these parameters, as indicated in Sect. 4. Do note, however, that stronger mortality increases the size of Ω k (color figure online) The set Ω 3 is depicted in Fig. 1. Note that all year classes go extinct if N k (mk) → 0 for m → ∞. Hypothesis 1 rules out the possibility that both N i (mk) and N k (mk) grow beyond any bound for m → ∞ while their ratio goes to zero. So we conclude that for initial conditions belonging to Ω k only the k-th year class can possibly persist.
Proof In order to avoid overburdening the reader with notational detail, we focus on the proof of the case k = 3. The proof for general k does not require new arguments.
From (1) it follows that The product of the three Leslie matrices is a diagonal matrix. We claim that the three elements of the diagonal of this matrix have a factor in common and are in fact of the form β(·)θ 1 with the argument of β being θ 1 N 1 (0) for the first element, θ 2 N 2 (0) with θ 2 := σ 2 σ 3 π(N(0))π(N(1)) for the second element, and θ 3 N 3 (0) with θ 3 := σ 3 π(N(0)) for the third element.
To verify the claim, we elaborate (10) from right to left, using π 1 to denote π(N(0)), π 2 to denote π(N(1)) and π 3 to denote π(N(2)). All factors in the products in the vectors below follow the order of events over time.
Noting that θ 1 = σ 1 σ 2 σ 3 π 1 π 2 π 3 , inspection reveals that the claim is correct. It follows that and we conclude that, since β is strictly increasing, Hence, when N 1 (0) > 0, Essentially the same argumentation yields that, when N(0) ∈ Ω 3 and N 2 (0) > 0, then It follows that Ω 3 is forward invariant under the full life cycle map and that for i = 1, 2 the sequences are decreasing in m, strictly when positive. It remains to be shown that either their limit is zero or all year classes go extinct.
Theorem 1 does neither restrict nor predict the resulting SYC dynamics. But does it allow for a large class of initial conditions? The answer, of course, heavily depends on what one means by "large". Certainly the set is open. Moreover, if a given initial condition does not belong to Ω k , one can apply (1) once and check whether N(1) belongs to Ω k (if it does, the year class that had label k − 1 in year 0 is the winner). By applying (1) repeatedly one can thus check for any of the year classes whether or not Theorem 1 guarantees that they will win the competition. The example of a steady state with all year classes present clearly shows that this procedure may fail to point out a winner, for the simple reason that there might not be a winner. The next question then is: how exceptional is it that several year classes persist?
It is tempting to conjecture that a fixed point of the full life cycle map, with more than 1 year class having positive density, is necessarily unstable, as was found to be the case for the model considered by Machta et al. (2019). However, in the "Appendix" we show that such a fixed point can in fact be stable if the negative density dependence incorporated in π is stronger than the positive density dependence incorporated in β, in the sense that an increase of the oldest age group N k can lead to a decrease of π(N)N k and thus to a decrease of the youngest age group in the next year. So without further restrictions on the class of models, it is not guaranteed that generically there will be a single winner. The restriction that for all non-negative N 1 , . . . , N k−1 the map N k → π(N)N k is increasing seems both meaningful and promising. See the "Appendix" for an example.
As illustrated in Fig. 2, numerical tests indicate that the occurrence of a single winner is in fact rather common. In the same figure we also evaluate the performance of the following Fig. 2 The long-term dynamics as determined numerically, illustrating that for almost all initial conditions, a single year class emerges as the "winner". Top row: the winning cohort number, as a function of initial values (N 1 (0), N 2 (0), N 3 (0)), shown in slices in the N 1 -direction (left) and in the N 3 -direction (right). Bottom row: the cohort number that is predicted to win using the Prediction Tool in the text. This tool is good, but not perfect, in predicting which cohort will eventually rule the world. Model ingredients are β(x) = max{0, 10(1 − 0.1/x)}, and π(N) = e −0.7(N 1 +N 2 +N 3 ) , and σ 1 = σ 2 = σ 3 = 1 Prediction Tool For an arbitrary initial condition N(0), find the index j that maximizes the diagonal elements in the full life cycle map, i.e., The prediction is that the year class having age j in year 0 will drive all other year classes to extinction. It appears that the tool works well, but is not infallible. We now present a theorem in which we incorporate more restrictions, but also describe more precisely the ultimate dynamics. The full life cycle map acts on R k . Any subspace characterised by k − 1 components being zero is forward invariant. The restriction to such a forward invariant subspace corresponds to a map from R to R that we shall call a SYC full life cycle map. The k maps are different but equivalent, see (Davydova 2004, Chapter V;Davydova et al. 2003, Section 7). To stay in line with Theorem 1 we work primarily with the map corresponding to N 1 = N 2 = · · · = N k−1 = 0.

Theorem 2 If a single year class fixed point or periodic orbit of the SYC full life cycle map is linearly stable with respect to that map, so with respect to within year class perturbations, it is automatically also linearly stable with respect to perturbations that involve small introductions of other year classes.
Proof We again give a proof for the case k = 3. The proof for general k uses no new arguments.
(Of course, all the π i are functions of N 1 ,N 2 , and N 3 .) We first consider a SYC fixed point and choose the phase such that only the third component of the fixed point vector is positive. Then necessarily M 3,3 = θ 1 β(θ 3 N 3 ) = 1 when evaluated in the fixed point. N 3 is itself a fixed point of the SYC map and is, by assumption, stable as such.
The linearisation of the full life cycle map in the fixed point is represented by a matrix as well, of course, L, say. Since the first two components of the fixed point are zero, L 1,2 = L 1,3 = L 2,3 = L 2,1 = 0, making L a lower-diagonal matrix, with eigenvalues equal to the diagonal elements. Direct computation shows that L must be of the form The diagonal element in the third and last row is complicated, but is equal to the linearisation of the SYC full life cycle map in the fixed point. Hence, since we have assumed that fixed point to be linearly stable, |L 3,3 | < 1. So all eigenvalues of L are less than one in absolute value and the fixed point of the three dimensional full life cycle map is linearly stable.
For an orbit of period p, we consider the p times iterated full life cycle map. The structure is exactly the same as before, the only difference is that diagonal elements have p factors θ 1 β. The Jacobi matrix is again lower-diagonal, and of the form By the same arguments as before, θ 1 ( p − 1) · · · θ 1 (0)β p (0) < 1, and |K 3,3 | < 1 by the assumption that the p-periodic orbit is linearly stable. So exactly as before we reach the conclusion that all eigenvalues are less than one in absolute value.

Model ingredients
To obtain more insight which choices of model ingredients π and β ensure persistence of populations and induce dissipative dynamics as defined in Hypothesis 1, we scale the variables in the following way. Let and define functions β new and π new according to In terms of these new variables, and dropping the superscripts for convenience, the map is seen to be simplified to Let us now choose the following form for π : So π(N) is determined by first computing a weighted population size E and next applying a scalar map p that assigns to E a component of the survival probability, i.e., a number between zero and one. Assume c 3 > 0, and define Let us now consider the SYC dynamics. To facilitate the description of that case, we introduce Dropping again the "new" superscripts, the SYC dynamics are given by To ensure boundedness of the population consisting of just 1 year class, the graph of the SYC full life cycle map, should be below the 45 • line for large values of N 3 . Our model ingredients β and p must be chosen accordingly. Figure 3 gives an impression of the form of the graph of the SYC full life cycle map, for different choices of the model ingredients.
The derivative at the trivial fixed point N 3 = 0 corresponds to multiplication by If this quantity exceeds 1, the population cannot go extinct, whereas there is an Allee effect when this quantity is less than one and yet part of the graph lies above the 45 • line (see the green and orange curves in Fig. 3 for examples). Some possible choices for p and β include This choice for β is strictly increasing only where it takes strictly positive values and equals zero in a neighbourhood of x = 0, so there is certainly an Allee effect. It corresponds to the total population size of adults being reduced at a constant rate during a fixed time window. Alternatively, we can solve, with P, the predator density, as a parameter, for a fixed period of time. This leads to an implicitly defined function β. For bifurcation diagrams of Single Year Class Maps, see Chapter V of (Davydova 2004). Concerning dissipativity (Hypothesis 1), let for N ∈ R 3 Then (15) implies that If for some choice of R > 0 and > 0 the inequality holds for all |N| ≥ R, then the set is forward invariant and attracts all orbits. Indeed, any orbit starting outside the set {N ∈ R 3 + : |N| ≤ R} has to enter this set and once inside this set we can use (18) and the fact that π(N) ≤ 1 to conclude that, if β(∞) > 1, it may again leave the ball with radius R but not the ball with radius β(∞)R.

Discussion
The discrete time population dynamics of semelparous species with one reproduction opportunity per year, is described by a special kind of Leslie matrices, viz. those for which in the first (= reproduction) row only the last element is non-zero. This reflects that an individual born in a certain year reproduces, if at all, exactly k years later, where k is the length of the life cycle. So the population splits into k reproductively isolated year classes. Mathematically this shows up in the fact that the k times iterated matrix is diagonal, so fully reducible.
Although they are reproductively isolated, the year classes do interact by competition for resources. When resource availability is constant in time, each year class does as well or bad as any other year class. But when resource availability fluctuates in time, the phase of the life cycle relative to resource peaks and troughs can be decisive for success or failure. Thus one year class can drive other year classes to extinction by inducing, for instance, periodic environmental conditions. The work of Davydova and co-workers (Davydova 2004;Davydova et al. 2003Davydova et al. , 2005Diekmann et al. 2005) focused on this phenomenon.
The present paper is inspired by Machta et al. (2019) and analyses a model such that, over a full life cycle, competition for food is neutral with respect to the year class distinction. But the model includes a second form of density dependence: it assumes that, due to predator satiation effects, per capita reproduction success of adults is larger when the cohort of adults is larger. So individuals belonging to a large cohort get much offspring while the negative impact of cohort size on survival probabilities is affecting individuals of all other year classes equally strongly, when measured over a full life cycle. Thus the 'strongest' year class becomes even stronger and eventually drives all other year classes to extinction and SYC (single year class) dynamics (Mjølhus et al. 2005;Davydova et al. 2003) gets established. The idea that a combination of predator satiation and intraspecific competition among larvae gives rise to single year class dynamics goes at least back to Hoppensteadt and Keller (1976) and Bulmer (1977), and has been demonstrated in many of the models in the literature by way of numerical experiments. The assumption of uniform competition introduced in Machta et al. (2019) allows, as we have shown above, to actually prove that SYC dynamics results for large classes of initial conditions. In our opinion, the precise nature of the relationship between mechanisms and phenomena is better revealed by a proof than by simulations.
To explain, in the context of the model, that several species coexist in synchrony would probably require the assumption that both the functions π and the functions β for the various species are proportional. So this is still rather puzzling.
We now briefly consider some of the ecological evidence supporting the two chief modelling ingredients.
The long larval stages associated with periodical insects likely result from development constrained by certain abiotic factors such as low temperature, poor food availability, and large adult body size (Danks 1992;Hellövaara et al. 1994). Magicicada larvae feed underground on xylem found in plant roots. Xylem is a transport liquid in plants and is poor in nutrients. The competition for this shared food resource thus likely affects all cohorts feeding on them. Cicada nymphs have been shown to be uniformly spatially distributed, suggesting that cohorts do indeed compete for xylem (Williams and Simon 1995).
Predator satiation is well-documented for periodical cicadas (Williams and Simon 1995, and many references therein). The first cicadas to emerge still face a high predation pressure, with as much as 40% eaten within the first days. As numbers increase dramatically in the days following the start of the event, with up to 3,5 million adult cicadas per hectare, per capita predation pressure drops practically to zero (Williams et al. 1993). The predators, mainly birds such as cuckoos, woodpeckers, jays and other larger birds, do not increase much in number during the year of the outbreak, but several show increased population sizes in the 1-3 years to follow (Koenig and Liebhold 2005). It has also been shown that for several of these species, predator population sizes are in fact lower right before an outbreak event, thus decreasing predation pressure and allowing the insects to benefit even more from their numerical dominance (Koenig and Liebhold 2013).
In this paper we have focussed purely on the problem of elimination of year classes and how a single cohort arises from a starting situation with multiple cohorts. In our modelling framework we have not allowed for so-called stragglers, insects that have a longer developmental time and thus emerge in the year after an outbreak. This has been investigated recently using predominantly numerical simulations by Blackwood et al. (2018).
The enigma of periodical insect species is not confined to Magicicada, but is found also in several other insect orders, notably among Xestia moths, and in some wellknown beetle species such as Melonontha cockchafer beetles (Hellövaara et al. 1994). Life spans may vary between species, and even between populations of the same species. Magicicada are exceptional, however, in their life spans, which are either 13 or 17 years, and are among the longest of all insect life spans.
There are several intriguing open theoretical questions regarding periodical insects. The adult insects, having developed under ground over a period of 13 or 17 years, all emerge from their burrows within an extremely short time span, usually between 7 and 10 days (Williams and Simon 1995, and references therein). Just how they synchronise so precisely remains unclear (Williams and Simon 1995), although temperature cues have been suggested.
In some outbreak years, as much as half of all individuals fail to develop into adults, and eclose the next year (White and Lloyd 1979). This has been associated with extremely poor nutrition conditions. Differences in the developmental rate of individual larvae are apparently common, with late-developing nymphs 'catching up' before emerging as adults with the rest. To incorporate such phenomena would require developing a physiologically structured population model, with developmental rate directly influenced by food availability and competition (de Roos and Persson 2013).
classes. In order that simple computations suffice to reach this conclusion, we choose k = 2, σ 1 = σ 2 = 1, So we focus our attention on If (N 1 ,N 2 ) is a steady state, we should havē The trivial steady state (0, 0) is locally stable but, due to positive density dependence incorporated in β, this does not exclude that nontrivial steady states exist. The equation has for two positive solutions, one withN 2 < 1 3 and one withN 2 > 1 3 . Each of these yields, when combined with (21), a nontrivial steady state. The corresponding linearized system is given by x 2 (t + 1) = e −N 2 x 1 (t) −N 2 x 2 (t).
So if we think in terms of bifurcations that happen when the parameter β 0 is increased, the scenario is as follows: a) at β 0 = 3e a saddle-node bifurcation of nontrivual steady states happens, but both steady states are unstable b) the larger, with respect to both components, of the two steady states undergoes a period-doubling bifurcation at β 0 = e 3 ; for e 3 < β 0 < 2 3 e 9/2 this steady state is stable. c) at β 0 = 2 3 e 9/2 this steady state loses stability in a Neimark-Sacker bifurcation. For completeness, let us have a brief look at SYC dynamics. If N 1 (t) = 0 we obtain with stable trivial steady state and two nontrivial steady states for β 0 > 2e arising by saddle-node bifurcation at β 0 = 2e with valueN 2 = 1 2 . The linearized recursion is so we see that the larger of the two is stable for 2e < β 0 < 2 3 e 3 , losing its stability by period doubling at the upper boundary of this window. It appears that the bifurcation diagrams of the "all year class" steady state and of the "single year class" steady state are in no way related to each other.
It is easy to repeat the bifurcation analysis with π in (19) replaced by π(N) = 1 1 + N 2 .
It turns out that in this case the "two year class" steady state is unstable for all parameter values.
In Fig. 4 we visualize the outcome of numerical experiments with k = 3 and model ingredients such that 2 year classes can coexist in a stable fixed point of the full life cycle map. By symmetry there are three such fixed points. It appears that the domains of attraction of these three fixed points are rather small. As a rule, there is an ultimate winner. Fig. 4 The long-term dynamics of the full life cycle map, with model ingredients set to the k = 3 analogue of (19). The figures show the winning cohort (1, 2 or 3), or whether a mixed steady state was reached, as a function of initial values (N 1 (0), N 2 (0), N 3 (0)). They are shown in slices in the N 1 -direction (left) and in the N 3 -direction (right). Mixed steady states appear on the boundaries of the basins of attraction for each year class, and involve the year classes on either side of these boundaries. Note that SYC dynamics still predominates