Synchronized Cicadae Emergence

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2006

Introduction

A fascinating feature of some insect populations is the phenomenon of synchronized emergence in which a very large population of mature adults emerges from the ground every k years, with very small numbers emerging in intervening years.  A dramatic example in the College Park, MD area is provided by the 17-year cicada.  I first experienced this phenomenon in 1987, a few years after arriving here.  Vast numbers of large, slow flying, extremely noisy insects appeared that spring, much to the delight of my 5 year-old daughter.  The noise produced by those insects in heavily wooded areas was truly deafening for about a fortnight, at the end of which there was a carpet of cicada corpses and many large empty shells attached to the trees.  My daughter had a lot of fun collecting those shells.  Although I had seen very few cicadas before, long-time residents assured me that this phenomenon occurs regularly with a 17 year cycle.  In fact, Bob Dylan's song Day of the Locusts was apparently inspired by the previous emergence in 1970, which occurred when he was awarded an honorary degree at Princeton.  (Although commonly called locusts, cicadas actually belong to the genus Magicicada and are more closely related to leafhoppers.)  Sure enough, it happened again in 2004; most students reading this notebook probably experienced the 2004 emergence of Brood X.  I hope that I will be around to see the next emergence!

To understand this phenomenon, we must know something about the life cycle of cicadae.  Adult cicadas spend only a couple of weeks above ground, during which their primary tasks are to mate, lay eggs in trees in slits cut in the ends of tree branches, and die.  The eggs soon fall to the ground and hatch nymphs.  The nymphs then burrow underground and attach themselves to a tree root.  There they spend the majority of their lives, sucking nourishment from their root.  Once mature, after 17 years for Magicicada septendecim, adults tunnel toward the surface and emerge when the ground temperature reaches 64°F, usually during mid-May in Maryland.  Cicadas are a good source of protein and are easy prey for birds, squirrels, or even dogs.  Periodical cicadas, those with synchronized emergence, rely upon a survival strategy of predator satiation — because an entire brood emerges almost simultaneously, there are simply too many for the predators to consume, which ensures that sufficient numbers survive to reproduce.  In this notebook we use a model proposed by Hoppensteadt and Keller to study how synchronized emergence can arise for long-cycle cicadae from balancing the limited carrying capacity of the environment against predator satiation.  However, this model does not explain why the life spans of periodical cicadae are relatively large prime numbers.  There is some speculation that evolution favored large prime numbers to minimize the probability that different species emerge simultaneously and hybridize.  Long maturation cycles may also reduce the likelihood that occasional environmental fluctuations in temperature or other conditions might exterminate an entire brood and lead to extinction of a species.  However, convincing explanations for these prime numbers are not yet available.

Further information and pictures can be obtained at http://www.ummz.lsa.umich.edu/magicicada/Periodical/.

Model

The simulation of this phenomenon that is presented here is based upon the model proposed by Hoppensteadt and Keller and is discussed in more detail in the textbook of Murray.  We provide no new development of the model, but rather use it to illustrate procedural programming with Mathematica.  Furthermore, the notebook format offers the opportunity to present more detailed graphics and tools for independent variations, which are not available using the traditional textbook format.

Let k represent the life span of a particular species.  The population dynamics are governed by several factors.  Each year a fraction μ of the nymphs survives, so that after k years the fraction of the nymphs that had become established on roots and that eventually emerges as adults is approximately (μ^k).  After emerging, the adults are subject to predation, primarily by birds, and have an average fertility f.  Now suppose that the predation one year is relatively light and the fertility high, so that there are a large number of nymphs seeking spots on the roots.  However, the capacity of the root system to accommodate and nourish the nymphs is limited, so that the excess nymphs fail to become established and die.  Thus, we require models of the fertility, predation, and carrying capacity.

adults = list containing the numbers of adults which emerge each year
predation = list containing the maximum number of adults that might be eaten each year by predators
offspring = list containing the numbers of offspring produced each year
capacity = list containing the numbers of offspring which can be accommodated each year
nymphs = list containing the numbers of nymphs that become established each year on the root system

Capacity

We assume that there is a maximum environmental capacity d and that the residual capacity is reduced by the number of nymphs established in previous years that survive to the present.  Thus, the formula

capacity[[t]]  =  
      Max[ 0,  
        d  -  (∑_(i = 1)^(Min[t, k] - 1))(μ^i) 
              nymphs[[t - 
                  i]] ]

sums over the previous generations of surviving nymphs and incorporates the constraint that the capacity must be nonnegative.

Emergent adults

The number of adults that emerge in any particular year is given by the number of nymphs established in the previous cycle which have survived the entire cycle, such that

adults[[t]]  =  (μ^k) 
        nymphs[[t - k]]

Offspring

The number of offspring is given by a fertility factor f times the number of adults which survive predation

offspring[[t]]  =  
      Max[ 0,  
        f (adults[[t]]  -  
              predation[[
                t]]) ]

Of course, the number of offspring must be nonnegative.

Establishment of nymphs

The number of nymphs which can be established in any given year is the smaller of the residual capacity and the total number of offspring produced.  Of course, the number of nymphs must be nonnegative.

nymphs[[t]]  =  
      Max[0,  Min[capacity[[t]], 
          offspring[[
            t]]] ]

Predation

The number of  adults which might be taken by predation depends upon the predator population, which in turn depends upon both the food supply available the previous year and other environmental constraints.  We assume that there is a predator population even without cicadas, but that the emergence of cicadas provides an increase in the next year's predator population that is proportional to this extra food supply.  However, there is also a natural attenuation of the predator population.  Let ν be the survival probability and a be the feeding factor.  Then a model of the effect of the cicada population on the predator population is

predation[[
        t]]  =  ν predation[[
            t - 1]]  +  
        a adults[[
            t - 1]]

Remember that the predation variable actually gives the number of cicadas which can be eaten, but that is obviously proportional to the predator population.

References:

F.C. Hoppensteadt and J. B. Keller, Synchronization of periodical cicada emergence, Science 194, 335 (1976)

J.D. Murray, Mathematical Biology, (Springer-Verlag, Berlin, 1989), chapter 4.

G. F. Webb, The prime number periodical cicada problem, Discrete and continuous dynamical systems B1, 387 (2001)

Implementation

The following population module must be provided k initial conditions, carried by the list nymphs0, representing the initial populations of nymphs on tree roots for the first cycle.  In the absence of data to the contrary, it is simplest to begin with a uniform distribution that is a small fraction of the environmental capacity.  We must also provide numerical values for the various fertility, survival, and predation parameters described in the model.  We assume that these parameters depend upon species but are independent of time and population.

The module begins by initializing the various lists.  Note that by using the Take function, we can allow the nymphs0 list to be longer than the actual life cycle of a particular species.  For each year up to the maximum, given by t, the number of emergent adults, eaten adults, offspring, nymphs, and residual capacity are evaluated and appended to the appropriate lists.  These lists grow as the model is propagated forward in time.  Upon exit, these lists contain t entries each.

ClearAll["Global`*"];
Needs["Graphics`Graphics`"]

\(\(population[nymphs0_, t_, k_, μ_, f_, d_, ν_, a_]\  := \ 
        Module[{i, j}, 

 (*\ 
            Initialize\ first\ cycle\ *) 
nymphs\  = \ 
            Take[nymphs0, k]; 
capacity\  = 
            Table[Max[0, 
                d - ∑\+\(i = 1\)\%\(j - 1\)\(μ\^i\) 
                      nymphs[[j - 
                          i]]], {j, 2, k}]; 
          capacity\  = 
            Prepend[capacity, 
              d]; 
predation\  = \(adults = \(offspring = \ 
                Table[0, {k}]\)\); 
		
 (*\ Loop\ over\ years\ *) \ 

          Do[
adults\  = \ 
              Append[adults, \(μ\^k\) 
                  nymphs[[j - 
                      k]]]; 
predation\  = Append[predation, ν\ predation[[
                      j - 1]]\  + \ 
                  a\ adults[[j - 1]]]; 
            offspring\  = \ 
              Append[offspring, \ 
                Max[0, f\ \((adults[[
                          j]] - 
                        predation[[
                          j]])\)]\ ]; 
capacity\  = 
              Append[capacity, 
                Max[0, d - ∑\+\(i = 1\)\%\(k - 1\)\(μ\^i\) 
                        nymphs[[
                          j - i]]]\ ]; 
            nymphs\  = 
              Append[nymphs, 
                Min[\ capacity[[j]], \ 
                  offspring[[
                    j]]\ ]\ ], 
\ {j, k + 1, 
              t}]];\)\)

Results

In this section we show a few representative results based upon the parameters suggested by Murray.  For simplicity, we assume a uniform initial population of 100 nymphs per year and a maximum capacity of 10000.  (In other words, the initial population per year is only 1% of the maximum capacity).

nymphs0  =  Table[100, {20}];

7-year cicada

Consider first a species with a 7-year life cycle.

population[nymphs0, 130, 7, 0.95, 10, 10000, 0.95, 0.042]

BarChart[adults, Edges -> False];

Although the population exhibits some strong peaks over its first few cycles, eventually an equilibrium condition is reached in which the emerging cicada population is uniform.

17-year cicada

Now consider a species with a 17-year life cycle, but otherwise the same fertility and predation parameters.

population[nymphs0, 200, 17, 0.95, 10, 10000, 0.95, 0.042]

BarChart[adults, PlotRange -> {Automatic, All}];

We now find that the initial transients become entrenched, such that eventually the emergence of adult cicadas is synchronized with a 17-year period.  By examining the competition between various aspects of the model, we can begin to appreciate how this remarkable phenomenon occurs.  First we note that although predation appears to be relatively mild, it is crucial and the fertility factor amplifies its influence.

BarChart[predation, PlotRange -> {Automatic, All}];

Next we notice that the number of offspring can be much larger than the capacity of the environment to sustain nymphs.  As this capacity becomes exhausted, most of the next year's offspring die off so that the next cycle for that group will have few adults and even fewer surviving predation to produce offspring.

BarChart[offspring, PlotRange -> {Automatic, All}];

BarChart[nymphs, PlotRange -> {Automatic, All}];

BarChart[capacity, PlotRange -> {Automatic, All}];

The Hoppensteadt-Keller model illustrates how competition between several interacting species for limited resources with delayed effects can produce a wide variety of phenomena.  With realistic choices of the parameters it can explain the synchronized emergence of insect populations.  However, like any model it neglects many other effects which may also be important and it does not explain the periods themselves.  Cicada species with life cycles of 3, 4, 7, 13, and 17 years are common and there probably are other values as well.  With the present parameters one can show that periods of 3, 4, or 7 years lead to uniform steady states, while periods of 13 or 17 years lead to synchronized emergence.  For a given period k, one can derive analytic relationships between the other parameters that determine which type of outcome is realized, but is there any fundamental significance to these parameters?  Such an analysis is beyond the scope of the present discussion and we defer to Murray for further development of these aspects of the model.  However, why are there no species with life spans of 15, 16, or 18 years?  And how do these insects count years anyway?

Demonstrate that without predation the 17-year cicada population reaches a uniform steady state.

Determine the minimum period which leads to synchronized emergence given the fertility and predation parameters suggested above.