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

• Odo Diekmann
• Robert Planqué
Open Access
Article

## Abstract

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.

## Keywords

Age-structured population models Nonlinear dynamics Nonlinear Leslie matrix models Single year class dynamics Semelparous insects

37N25 92B25

## 1 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! [Pp. 133–134; 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 2012, 2017; Cushing 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\rightarrow \infty$$ 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, 2005). 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 \le 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.

## 2 Model formulation

Let k be an integer bigger than one. Let $$\varvec{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, \ldots ,k$$, even though biologically $$j= 0,1,\ldots ,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
\begin{aligned} \varvec{N}(t+1) = L(\varvec{h}(\varvec{N}(t)) \varvec{N}(t). \end{aligned}
(1)
where
\begin{aligned} L(\varvec{h}) = \begin{pmatrix} 0 &{} 0 &{} \cdots &{} 0 &{} h_k\\ h_1 &{} 0 &{} &{} 0 &{}0\\ 0 &{} h_2 &{} &{} 0 &{}0\\ \vdots &{} \vdots &{} \ddots &{} \vdots &{} \vdots \\ 0 &{} 0 &{} \cdots &{} h_{k-1} &{} 0 \end{pmatrix} \end{aligned}
(2)
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 $$\varvec{N}$$ captures density dependence. We assume that
\begin{aligned} h_i(\varvec{N})&= \sigma _i \pi (\varvec{N}) \quad \text { for } i=1,\ldots ,k-1, \end{aligned}
(3)
\begin{aligned} h_k(\varvec{N})&= \sigma _k \pi (\varvec{N})\beta (\sigma _k \pi (\varvec{N}) N_k). \end{aligned}
(4)
where
\begin{aligned} \sigma _i&\in (0,1],\quad i=1,\ldots ,k, \end{aligned}
(5)
\begin{aligned} \pi&: {\mathbb {R}}^k_+ \rightarrow (0,1], \end{aligned}
(6)
\begin{aligned} \beta&: [0,\infty ) \rightarrow [0,\infty ) \text { is continuous, strictly increasing and bounded.} \end{aligned}
(7)
Here the $$\sigma _i$$ are survival probabilities under “standard” conditions and the (uniform, i.e., i-independent) factor $$\pi (\varvec{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 $$\sigma _k \pi (\varvec{N})N_k$$. The function $$\beta$$ is a composite model ingredient. It describes the number of offspring at the census moment per emerging adult. The monotonicity of $$\beta$$ 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).

## 3 The main results

So far we did not specify any property of $$\pi$$ 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 $$\pi$$.

### Hypothesis 1

(dissipativity (Hale 1988)) There exists $$R>0$$ such that
\begin{aligned} B^+_R:=\{ \varvec{N}\in {\mathbb {R}}^k_+:\ \sum _{i=1}^k N_i \le R\} \end{aligned}
(8)
is forward invariant and attracts all orbits.

By “attracts all orbits” we mean that for every $$\varvec{N}_0 \in {\mathbb {R}}^k_+$$ there exists $$m = m(\varvec{N}_0)$$ such that the solution of (1) with initial condition $$\varvec{N}(0) =\varvec{N}_0$$ satisfies $$\varvec{N}(t) \in B^+_R$$ for $$t\ge m(\varvec{N}_0)$$. We adopt Hypothesis 1 for the rest of the paper. In Sect. 4 below we provide easily verifiable conditions on $$\pi$$ and $$\beta$$ that guarantee that Hypothesis 1 is satisfied.

At the end of Sect. 1 we already observed that (due to symmetry, see (Diekmann and van Gils 2003; 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 $$\sigma$$’s) before reaching the highest age.

### Theorem 1

The set
\begin{aligned} \varOmega _k := \{\varvec{N} \in {\mathbb {R}}^k_+ : N_k > \sigma _i\cdots \sigma _{k-1}N_i \text { for } i=1,\ldots ,k-1\} \end{aligned}
(9)
is forward invariant under the full life cycle map, and if $$\varvec{N}(0) \in \varOmega _k$$ and for given $$i\in \{1,\ldots ,k-1\}$$, $$N_i(0) > 0$$, then
\begin{aligned} \frac{N_i(mk)}{N_k(mk)} \end{aligned}
is a strictly decreasing function of m. This function has limit zero for $$m\rightarrow \infty$$ in case $$\limsup _{m\rightarrow \infty } N_k(mk)>0$$.

The set $$\varOmega _3$$ is depicted in Fig. 1. Note that all year classes go extinct if $$N_k(mk) \rightarrow 0$$ for $$m\rightarrow \infty$$. Hypothesis 1 rules out the possibility that both $$N_i(mk)$$ and $$N_k(mk)$$ grow beyond any bound for $$m\rightarrow \infty$$ while their ratio goes to zero. So we conclude that for initial conditions belonging to $$\varOmega _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
\begin{aligned} \varvec{N}(3) = L(\varvec{h}(\varvec{N}(2))) L(\varvec{h}(\varvec{N}(1))) L(\varvec{h}(\varvec{N}(0))) \varvec{N}(0). \end{aligned}
(10)
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
\begin{aligned} \theta _1 := \sigma _1\sigma _2\sigma _3 \pi (\varvec{N}(0))\pi (\varvec{N}(1)) \pi (\varvec{N}(2)) \end{aligned}
(11)
in common and are in fact of the form $$\beta (\cdot )\theta _1$$ with the argument of $$\beta$$ being $$\theta _1 N_1(0)$$ for the first element, $$\theta _2 N_2(0)$$ with $$\theta _2:=\sigma _2\sigma _3 \pi (\varvec{N}(0)) \pi (\varvec{N}(1))$$ for the second element, and $$\theta _3 N_3(0)$$ with $$\theta _3 := \sigma _3 \pi (\varvec{N}(0))$$ for the third element.
To verify the claim, we elaborate (10) from right to left, using $$\pi _1$$ to denote $$\pi (\varvec{N}(0))$$, $$\pi _2$$ to denote $$\pi (\varvec{N}(1))$$ and $$\pi _3$$ to denote $$\pi (\varvec{N}(2))$$. All factors in the products in the vectors below follow the order of events over time.
\begin{aligned} \begin{pmatrix} N_1\\ N_2\\ N_3 \end{pmatrix}&\mapsto \begin{pmatrix} \beta (\sigma _3 \pi _1N_3) \sigma _3\pi _1 N_3\\ \sigma _1\pi _1 N_1\\ \sigma _2\pi _1 N_2 \end{pmatrix} \mapsto \begin{pmatrix} \beta (\sigma _3 \pi _2 \sigma _2\pi _1N_2) \sigma _3 \pi _2 \sigma _2\pi _1N_2\\ \sigma _1\pi _2\beta (\sigma _3\pi _1 N_3)\sigma _3\pi _1 N_3\\ \sigma _2\pi _2\sigma _1\pi _1N_1 \end{pmatrix}\nonumber \\&\quad \mapsto \begin{pmatrix} \beta (\sigma _3 \pi _3 \sigma _2 \pi _2 \sigma _1 \pi _1N_1) \sigma _3 \pi _3 \sigma _2 \pi _2 \sigma _1 \pi _1N_1\\ \sigma _1\pi _3\beta (\sigma _3 \pi _2 \sigma _2\pi _1N_2)\sigma _3 \pi _2 \sigma _2\pi _1N_2 \\ \sigma _2\pi _3\sigma _1\pi _2\beta (\sigma _3\pi _1N_3)\sigma _3\pi _1 N_3 \end{pmatrix}. \end{aligned}
(12)
Noting that $$\theta _1=\sigma _1\sigma _2\sigma _3\pi _1\pi _2\pi _3$$, inspection reveals that the claim is correct. It follows that
\begin{aligned} \frac{N_1(3)}{N_3(3)} = \frac{\beta ( \theta _1 N_1(0))}{\beta ( \theta _3 N_3(0))}\frac{N_1(0)}{N_3(0)}. \end{aligned}
(13)
Now for $$\varvec{N}(0) \in \varOmega _3$$, $$\sigma _1\sigma _2 N_1(0) < N_3(0)$$, so that
\begin{aligned} \theta _1 N_1(0)&= \sigma _3 \pi (\varvec{N}(0))\Big ( \sigma _2\sigma _1 \pi (\varvec{N}(1)) \pi (\varvec{N}(2))\Big )N_1(0)\\&\le \sigma _3 \pi (\varvec{N}(0)) \sigma _2\sigma _1N_1(0)\\&< \sigma _3 \pi (\varvec{N}(0)) N_3(0),\\&= \theta _3 N_3(0), \end{aligned}
and we conclude that, since $$\beta$$ is strictly increasing,
\begin{aligned} \beta ( \theta _1 N_1(0)) < \beta ( \theta _3 N_3(0)). \end{aligned}
Hence, when $$N_1(0)>0$$,
\begin{aligned} \frac{N_1(3)}{N_3(3)}< \frac{N_1(0)}{N_3(0)} < \frac{1}{\sigma _1\sigma _2}. \end{aligned}
Essentially the same argumentation yields that, when $$\varvec{N}(0)\in \varOmega _3$$ and $$N_2(0)>0$$, then
\begin{aligned} \quad \frac{N_2(3)}{N_3(3)}< \frac{N_2(0)}{N_3(0)} < \frac{1}{\sigma _2}. \end{aligned}
It follows that $$\varOmega _3$$ is forward invariant under the full life cycle map and that for $$i= 1,2$$ the sequences
\begin{aligned} \frac{N_i(3m)}{N_3(3m)} \end{aligned}
are decreasing in m, strictly when positive. It remains to be shown that either their limit is zero or all year classes go extinct.
Assume that $$N_1(3m)/N_3(3m) \rightarrow \alpha$$ for $$m\rightarrow \infty$$. Extending the definitions of $$\theta _i$$ to
\begin{aligned} \theta _1(m)&:= \sigma _1\sigma _2\sigma _3 \pi (\varvec{N}(3m)) \pi (\varvec{N}(3m+1)) \pi (\varvec{N}(3m+2)),\\ \theta _3(m)&:= \sigma _3 \pi (\varvec{N}(3m)), \end{aligned}
then by iteration (13) leads to
\begin{aligned} \frac{N_1(3(m+1))}{N_3(3(m+1))} = \frac{\beta ( \theta _1(m) N_1(3m))}{\beta ( \theta _3(m)N_3(3m))}\frac{N_1(3m)}{N_3(3m)}. \end{aligned}
(14)
Passing to the limit $$m\rightarrow \infty$$ in (14), it follows that if $$\alpha >0$$ then necessarily, since $$\beta$$ is continuous,
\begin{aligned} \frac{\beta (\theta _1(m) \alpha N_3(3m))}{\beta (\theta _3(m)N_3(3m))}\rightarrow 1 \end{aligned}
for $$m\rightarrow \infty$$. This is very well possible if $$N_3(3m)\rightarrow 0$$ for $$m\rightarrow \infty$$ (and in that case, since $$\varOmega _3$$ is forward invariant, also $$N_1(3m)$$ and $$N_2(3m)$$ converge to zero for $$m\rightarrow \infty$$, so all year classes go extinct). But if $$\limsup _{m\rightarrow \infty } N_3(3m)>0$$, this is not possible, since $$\alpha <1$$, $$\theta _1(m)\le \theta _3(m)$$ because $$\sigma _1\sigma _2 \le 1$$, $$\pi (\varvec{N}(3m+1)) \pi (\varvec{N}(3m+2))\le 1$$, and $$\beta$$ is strictly increasing. $$\square$$

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 $$\varOmega _k$$, one can apply (1) once and check whether $$\varvec{N}(1)$$ belongs to $$\varOmega _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 $$\pi$$ is stronger than the positive density dependence incorporated in $$\beta$$, in the sense that an increase of the oldest age group $$N_k$$ can lead to a decrease of $$\pi (\varvec{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, \ldots , N_{k-1}$$ the map $$N_k \mapsto \pi (\varvec{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

Prediction ToolFor an arbitrary initial condition$$\varvec{N}(0)$$, find the indexjthat maximizes the diagonal elements in the full life cycle map, i.e.,
\begin{aligned} j = \arg \max _i \frac{N_i(k)}{N_i(0)}. \end{aligned}
The prediction is that the year class having agejin 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 $${\mathbb {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 $${\mathbb {R}}$$ to $${\mathbb {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 = \cdots = 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.

The full life cycle map given by (12) is represented by a diagonal $$3\times 3$$ matrix M acting on a 3 vector,
\begin{aligned} M = \begin{pmatrix} \theta _1 \beta (\theta _1 N_1) &{} 0 &{} 0\\ 0 &{} \theta _1 \beta (\theta _2 N_2) &{} 0 \\ 0 &{} 0 &{} \theta _1 \beta (\theta _3 N_3) \end{pmatrix}, \end{aligned}
where, as before,
\begin{aligned} \theta _1&= \sigma _1\sigma _2\sigma _3\pi _1\pi _2\pi _3,\\ \theta _2&= \sigma _2\sigma _3\pi _1\pi _3,\\ \theta _3&= \sigma _3 \pi _1. \end{aligned}
(Of course, all the $$\pi _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} = \theta _1 \beta (\theta _3 N_3) = 1$$ when evaluated in the fixed point. $$N_3$$ is itself a fixed point of the SYC map
\begin{aligned} N_3 \mapsto \theta _1 \beta (\theta _3 N_3) N_3, \end{aligned}
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
\begin{aligned} L = \begin{pmatrix} \theta _1 \beta (0) &{} 0 &{} 0\\ 0 &{} \theta _1 \beta (0) &{} 0\\ L_{3,1} &{} L_{3,2} &{} L_{3,3} \end{pmatrix}. \end{aligned}
Since $$\beta$$ is increasing, $$0 \le \theta _1 \beta (0) < \theta _1\beta ( \theta _3 N_3) = M_{3,3} = 1$$.

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 $$\theta _1\beta$$. The Jacobi matrix is again lower-diagonal, and of the form
\begin{aligned} K = \begin{pmatrix} \theta _1(p-1) \cdots \theta _1(0) \beta ^p(0) &{} 0 &{} 0\\ 0 &{} \theta _1(p-1) \cdots \theta _1(0) \beta ^p(0) &{} 0\\ K_{3,1} &{} K_{3,2} &{} K_{3,3} \end{pmatrix}. \end{aligned}
By the same arguments as before, $$\theta _1(p-1) \cdots \theta _1(0) \beta ^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. $$\square$$

## 4 Model ingredients

To obtain more insight which choices of model ingredients $$\pi$$ and $$\beta$$ ensure persistence of populations and induce dissipative dynamics as defined in Hypothesis 1, we scale the variables in the following way. Let
\begin{aligned} N_2^\text {old}= \sigma _1 N_2^\text {new},\quad N_3^\text {old}= \sigma _1\sigma _2 N_3^\text {new}, \end{aligned}
and define functions $$\beta ^\text {new}$$ and $$\pi ^\text {new}$$ according to
\begin{aligned} \beta ^\text {new}(x)&= \sigma _1\sigma _2\sigma _3 \beta ^\text {old}(\sigma _1\sigma _2\sigma _3 x),\\ \pi ^\text {new}(\varvec{N}^\text {new})&= \pi ^\text {old}(\varvec{N}^\text {old}). \end{aligned}
In terms of these new variables, and dropping the superscripts for convenience, the map is seen to be simplified to
\begin{aligned} \varvec{N}(t+1) = \pi (\varvec{N}(t))\begin{pmatrix} 0 &{} 0 &{} \beta (\pi (\varvec{N}(t))N_3(t))\\ 1 &{} 0 &{}0\\ 0 &{} 1 &{} 0 \end{pmatrix}\varvec{N}(t). \end{aligned}
(15)
Let us now choose the following form for $$\pi$$:
\begin{aligned} \pi (\varvec{N}) = p(E),\quad \text {with } E:=c_1 N_1 + c_2 N_2 + c_3 N_3. \end{aligned}
(16)
So $$\pi (\varvec{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
\begin{aligned} E^\text {new}= \frac{1}{c_3}E^\text {old}= c_1^\text {new}N_1 + c_2^\text {new}N_2+ N_3, \end{aligned}
and
\begin{aligned} p^\text {new}(E^\text {new}) = p^\text {old}(E^\text {old}). \end{aligned}
Let us now consider the SYC dynamics. To facilitate the description of that case, we introduce
\begin{aligned} f(x) =\beta (p(x)x)p(x) x. \end{aligned}
Dropping again the “new” superscripts, the SYC dynamics are given by
\begin{aligned} N_3&\mapsto f(N_3)\\&\mapsto p(c_1f(N_3)) f(N_3)\\&\mapsto p(c_2 p(c_1 f(N_3)) f(N_3)) p(c_1 f(N_3)) f(N_3). \end{aligned}
To ensure boundedness of the population consisting of just 1 year class, the graph of the SYC full life cycle map,
\begin{aligned} N_3\mapsto p(c_2 p(c_1 f(N_3)) f(N_3)) p(c_1 f(N_3)) f(N_3) \end{aligned}
(17)
should be below the $$45^\circ$$ line for large values of $$N_3$$. Our model ingredients $$\beta$$ 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
\begin{aligned} \beta (0) (p(0))^3. \end{aligned}
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^\circ$$ line (see the green and orange curves in Fig. 3 for examples).
Some possible choices for p and $$\beta$$ include
\begin{aligned} p(E)&= \frac{1}{1+E},\quad p(E) = e^{-E},\\ \beta (x)&= \beta _0\left( 1-\frac{\zeta }{x} \right) _+. \end{aligned}
This choice for $$\beta$$ 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,
\begin{aligned} \frac{dx}{dt} = - \frac{x}{1+\zeta x}P, \end{aligned}
for a fixed period of time. This leads to an implicitly defined function $$\beta$$. For bifurcation diagrams of Single Year Class Maps, see Chapter V of (Davydova 2004).
Concerning dissipativity (Hypothesis 1), let for $$\varvec{N} \in {\mathbb {R}}^3_+$$
\begin{aligned} |\varvec{N}| := N_1 + N_2 + N_3. \end{aligned}
Then (15) implies that
\begin{aligned} |\varvec{N}(t+1)| \le \pi (\varvec{N}(t)) \max \{1,\beta (\infty )\} |\varvec{N}(t)|. \end{aligned}
(18)
If for some choice of $$R>0$$ and $$\epsilon > 0$$ the inequality
\begin{aligned} \pi (\varvec{N}) \beta (\infty )\le 1-\epsilon \end{aligned}
holds for all $$|\varvec{N}| \ge R$$, then the set
\begin{aligned} \left\{ \varvec{N}\in {\mathbb {R}}^3_+\ :\ |\varvec{N}| \le \max \{\beta (\infty ),1\}R\right\} \end{aligned}
is forward invariant and attracts all orbits. Indeed, any orbit starting outside the set $$\{ \varvec{N} \in {\mathbb {R}}^3_+\ :\ |\varvec{N}| \le R\}$$ has to enter this set and once inside this set we can use (18) and the fact that $$\pi (\varvec{N}) \le 1$$ to conclude that, if $$\beta (\infty ) > 1$$, it may again leave the ball with radius R but not the ball with radius $$\beta (\infty )R$$.

## 5 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. 2003, 2005; Diekmann 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 $$\pi$$ and the functions $$\beta$$ 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 well-known 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).

## Notes

### Acknowledgements

Part of the work was done during the Mathematical Biology semester at the Mittag-Leffler Institute. We thank the organisers, in particular Torbjörn Lundh and Mats Gyllenberg, for making it such a success. We also thank two anonymous reviewers for their careful reading and suggestions to improve the manuscript.

## References

1. Behncke H (2000) Periodical cicadas. J Math Biol 40:413–431
2. Blackwood JC, Machta J, Meyer AD, Noble AE, Hastings A, Liebhold AM (2018) Competition and stragglers as mediators of developmental synchrony in periodical cicadas. Am Nat 192(4):479–489
3. Bulmer MG (1977) Periodical insects. Am Nat 111:1099–1117
4. Cushing JM (2009) Three stage semelparous Leslie models. J Math Biol 59:75–104
5. Cushing JM (2015) On the fundamental bifurcation theorem for semelparous Leslie models. In: Bourguignon JP, Jeltsch R, Pintoa A, Viana M (eds) Mathematics of planet earth: dynamics, games and science, CIM mathematical sciences series, chapter 11. Springer, BerlinGoogle Scholar
6. Cushing JM, Henson SM (2012) Stable bifurcations in nonlinear semelparous Leslie models. J Biol Dyn 6:80–102
7. Danks HV (1992) Long life cycles in insects. Can Entomol 124(1):167–187
8. Davydova NV (2004) Old and young. Can they coexist? Ph.D. thesis, Utrecht University https://dspace.library.uu.nl/handle/1874/891. Accessed 26 Apr 2019
9. Davydova NV, Diekmann O, van Gils SA (2003) Year class coexistence or competitive exclusion for strict biennials? J Math Biol 46:95–131
10. Davydova NV, Diekmann O, van Gils SA (2005) On circulant populations. I. The algebra of semelparity. Linear Algebra Its Appl 398:185–243
11. de Roos AM, Persson L (2013) Population and community ecology of ontogenetic development. Monographs in population biology 51. Princeton University Press, Princeton
12. Diekmann O, van Gils SA (2003) Invariance and symmetry in a year-class model. In: Buescu J, Castro S, Dias AP, Labouriau I (eds) Bifurcation, symmetry and patterns, Birkhäuser trends in mathematics. Birkhäuser, BaselGoogle Scholar
13. Diekmann O, Davydova NV, van Gils SA (2005) On a boom and bust year class cycle. J Differ Equ Appl 11(4):327–335
14. Hale JK (1988) Asymptotic behavior of dissipative systems. American Mathematical Society, Providence
15. Hellövaara K, Väisänen R, Simon C (1994) Evolutionary ecology of periodical insects. Trends Ecol Evol 9(12):475–480
16. Hofbauer F, Hofbauer J, Raith P, Steingberger T (2004) Intermingled basins in a two species system. J Math Biol 49(3):293–309
17. Hoppensteadt FC, Keller JB (1976) Synchronization of periodical cicada emergences. Science 194(4262):335–337
18. Koenig WD, Liebhold AM (2005) Effects of periodical cicada emergences on abundance and synchrony of avian populations. Ecology 86(7):1873–1882
19. Koenig WD, Liebhold AM (2013) Avian predation pressure as a potential driver of periodical cicada cycle length. Am Nat 181(1):145–149
20. Kon R (2012) Permanence induced by life-cycle resonances: the periodical cicada problem. J Biol Dyn 6(2):855–890
21. Kon R (2017) Non-synchronous oscillations in four-dimensional nonlinear semelparous Leslie matrix models. J Differ Equ Appl 23(10):1747–1759
22. Lloyd M, Dybas HS (1966) The periodical cicada problem. Evolution 20:133–149
23. Machta Jonathan, Blackwood Julie C, Noble Andrew, Liebhold Andrew M, Hastings Alan (2019) A hybrid model for the population dynamics of periodical cicadas. Bull Math Biol 81(4):1122–1142
24. Marshall DC (2008) Periodical cicadas: Magicicada spp. (Hemiptera: Cicadidae). In: Capinera JL et al (eds) Encyclopaedia of entomology, 2nd edn. Springer, Dordrecht, pp 2785–2794Google Scholar
25. Mjølhus E, Wikan A, Solberg T (2005) On synchronization in semelparous populations. J Math Biol 50:1–21
26. Webb GF (2001) The prime number periodical cicada problem. Disc Cont Dyn Sys B 1:387–399
27. White J, Lloyd M (1979) Seventeen year cicadas emerging after 18 years: a new brood? Evolution 33:1193–1199
28. Wikan A (2017) An analysis of a semelparous population model with density-dependent fecundity and density-dependent survival probabilities. J Appl Math 2017:14 (Article ID 8934295)
29. Williams KS, Simon C (1995) The ecology, behavior, and evolution of periodical cicadas. Ann Rev Entomol 40:269–295
30. Williams KS, Smith KG, Stephen FM (1993) Emergence of 13-yr periodical cicadas (Cicadidae: Magicicada): phenology, mortality, and predators satiation. Ecology 74(4):1143–1152
{"extraction_info": {"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 0, "mathjax_display_tex": 2, "mathjax_asciimath": 0, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 0, "equation": 0, "x-ck12": 0, "texerror": 0, "math_score": 0.931999683380127, "perplexity": 2786.923350250813}, "config": {"markdown_headings": true, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "plain"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "warc_path": "s3://commoncrawl/crawl-data/CC-MAIN-2019-43/segments/1570987838289.72/warc/CC-MAIN-20191024012613-20191024040113-00153.warc.gz"}