纽约大学：《数学建模》课程教学资源（讲义）Prime Number Selection of Cycles in a Predator-Prey Model

RESEARCHARTICLEPrime Number Selection of Cycles in aPredator-Prey ModelERICGOLES.tOLIVERSCHULZ.*ANDMARIOMARKUSCenterforMathematicalModellingofComplexSystems,FCFM,UniversityofChile,Casilla170-3Santiago,Chile*Max-Planck-Institut fuir molekulare Physiologie, Postfach 500247,D-44202 Dortmund,GermanyReceivedSeptember14,2000;revised January30,2001;acceptedJanuary30,2001The factthat some species of cicadas appear every7, 13, or 17 years and that these periods are prime numbershas been regarded as a coincidence.We found a simpleevolutionarypredator-preymodel that yields primeperiodic preys having cycles predominantly around the observed values.An evolutionary game on a spatialarray leads to travelling waves reminiscent of those observed in excitable systems.Themodel marks an en-counteroftwo seemingly unrelated disciplines:biology and numbertheory.Arestrictiontothelatter,providesan evolutionarygenerator of arbitrarilylargeprime numbers.2001 John Wiley& Sons, Inc.1.INTRODUCTIONsitoids that attack eggs or adults,whichmayhavebecomehe appearance of three species of the genus Magici-extinct. Models of periodical cicadas presented so far [4,5]cada synchronously every13 or17years and asyn-show that synchronized periodical behavior is possibleforchronously every7years [1,2] has notbeenconvinc-periods longer than 1o years. An alternative mechanism toingly explained.These cicadas spend most of their lives be-the predator hypothesis is given (without model calcula-low the ground, emerging and dying within a few weeks.tions) by Yoshimura [6,7]; he argues that prime numbers areThepresent work is based on the hypothesis that the cycleselected because these cycles are the least likely to co-length is a prime number in order to optimally escapeemerge and hybridize, so that they prevent genetic break-predators. For example, a prey with a 12-year cycle willdown by breeding synchrony.Thismechanism has beenmeet-everytime itappears-properlysynchronized preda-compared [8] with thatproposed byCox and Carlton [9],tors appearing every1, 2, 3, 4, 6or12years,whereasawhich also involves advantageof primecycles dueto lessmutant with a 13-year period has the advantage of beingfrequenthybridizationsubject to fewer predators. According to R. MacArthur [1],Thepurpose of thepresent work is twofold: (i)wewantthis idea maybethe only application of number theoryinto underlinethat in spite ofmissingbiological data support-mathematical biology, a drawback, however, is that there ising the predator hypothesis, a simple,purely temporalas yet no evidencefor relevantperiodicpredators ofcicadas.model doeslead tolockingintoprime-periodicpreycycles,Nevertheless, Lloyd and Dybas [3] pointed out that thewhereas a spatiotemporal model leads to a maximum prob-predator hypothesis can be maintained by assuming para-ability of such cycles; (i) one can use these biological ideas 2001 John Wiley & Sons, Inc., Vol. 6, No. 4COMPLEXITY33

Prime Number Selection of Cycles in a Predator-Prey Model ERIC GOLES,† OLIVER SCHULZ,* AND MARIO MARKUS* †Center for Mathematical Modelling of Complex Systems, FCFM, University of Chile, Casilla 170-3, Santiago, Chile *Max-Planck-Institut fu¨r molekulare Physiologie, Postfach 500247, D-44202 Dortmund, Germany Received September 14, 2000; revised January 30, 2001; accepted January 30, 2001 The fact that some species of cicadas appear every 7, 13, or 17 years and that these periods are prime numbers has been regarded as a coincidence. We found a simple evolutionary predator-prey model that yields prime￾periodic preys having cycles predominantly around the observed values. An evolutionary game on a spatial array leads to travelling waves reminiscent of those observed in excitable systems. The model marks an en￾counter of two seemingly unrelated disciplines: biology and number theory. A restriction to the latter, provides an evolutionary generator of arbitrarily large prime numbers. q 2001 John Wiley & Sons, Inc. 1. INTRODUCTION T he appearance of three species of the genus Magici￾cada synchronously every 13 or 17 years and asyn￾chronously every 7 years [1,2] has not been convinc￾ingly explained. These cicadas spend most of their lives be￾low the ground, emerging and dying within a few weeks. The present work is based on the hypothesis that the cycle length is a prime number in order to optimally escape predators. For example, a prey with a 12-year cycle will meet—every time it appears—properly synchronized preda￾tors appearing every 1, 2, 3, 4, 6 or 12 years, whereas a mutant with a 13-year period has the advantage of being subject to fewer predators. According to R. MacArthur [1], this idea may be the only application of number theory in mathematical biology; a drawback, however, is that there is as yet no evidence for relevant periodic predators of cicadas. Nevertheless, Lloyd and Dybas [3] pointed out that the predator hypothesis can be maintained by assuming para￾sitoids that attack eggs or adults, which may have become extinct. Models of periodical cicadas presented so far [4,5] show that synchronized periodical behavior is possible for periods longer than 10 years. An alternative mechanism to the predator hypothesis is given (without model calcula￾tions) by Yoshimura [6,7]; he argues that prime numbers are selected because these cycles are the least likely to co￾emerge and hybridize, so that they prevent genetic break￾down by breeding synchrony. This mechanism has been compared [8] with that proposed by Cox and Carlton [9], which also involves advantage of prime cycles due to less frequent hybridization. The purpose of the present work is twofold: (i) we want to underline that in spite of missing biological data support￾ing the predator hypothesis, a simple, purely temporal model does lead to locking into prime-periodic prey cycles, whereas a spatiotemporal model leads to a maximum prob￾ability of such cycles; (ii) one can use these biological ideas © 2001 John Wiley & Sons, Inc., Vol. 6, No. 4 C O M P L E X I T Y 33

for a nonbiological,purely number-theoretical goal, namelyNow we proceed to analyze the model described in thetoconstructanevolutionaryalgorithmyieldingprimenum-previous paragraph.Let Icm(X,Y):least common multiple,gcd(X,Y): greatest common divisor of X and Y.In XY years,bersofanysize.Inaddition,onecanconsiderthepresentthe predator appears Ytimes, both predator and prey ap-work from the viewpoint of an evolutionary game [1o] in-pear XY/lcm(X,Y) times; thus, predator without prey ap-cluding spatial dimensions; such a viewpoint sets a link topears Y -- XYI lcm(X,Y) times. Considering that gcd(X,Y) -phenomena reported for spatial versions of the Prisoner'sIcm(X,Y)=XY,wethus obtain thepredatorfitnessF(X,Y) =dilemma [11-13], the hawk-dove game [14], and prebiotic2gcd(X,Y)/Y-1.Analogously,one obtains thepreyfitnessevolution[15].F(X,Y) =1 - 2gcd(X,Y)/XWewill nowshowthefollowing:if weallowrandom2.SIMULATIONOFTEMPORALPROCESSESmutations—which canbeof anysize,aslongas theylead toWe now set up a model considering first a predator of pe-mutantswithin the ranges2≤X≤L/2,L/2+2≤Y≤riod Xinteracting with a prey of period Y. We assign a mo-Lthen a sequence of suchmutations will finallylock thementaryfitness f(t) of the prey in a year t as follows: it isprey intoa stable prime period Y. We call H the domainzeroifitis notpresent, it is-1if both predator and prey aredefinedby theseallowedranges of Xand Y,Thedefinition ofpresent, and it is +1 if the prey is present but the predator isHontheX- YplaneisgiveninFigure1.not.Note that we punish prey that appears and meets aTo be more precise, wewill show that if Yis not a prime,predator (f,(t) = -1), compared with the case of nonemer-then there exists a sequence of mutations that will changegence of the prey (f,(t)=O); we do this because emergenceY, and if Y is prime, then no mutation will change it. Let ususes up metabolic resources,becauseof metamorphosis,firstassume that Y=Yis nota prime; F(X,Y)has themating,and death; these resources are lost if the prey ismaximumvalue2gcd(XY)/Y-1atthepredatorperiodeaten up by a predator, whereas they are preserved if theXm=gd(Y) (gd: greatest divisor). Note that 1 Fy,resp.F>FrThus,in the caseoffitnessDefinition of the domain H within which mutations of X (predatorequality, the resident and not the mutant is selected. Herecycle length; abscissa) and Y (prey cycle length; ordinate) areand in the rest of this work, we assume that all interactingallowedpopulations are synchronized, thus being all present at t=0.34COMPLEXITY 2001 John Wiley & Sons, Inc

for a nonbiological, purely number-theoretical goal, namely to construct an evolutionary algorithm yielding prime num￾bers of any size. In addition, one can consider the present work from the viewpoint of an evolutionary game [10] in￾cluding spatial dimensions; such a viewpoint sets a link to phenomena reported for spatial versions of the Prisoner’s dilemma [11–13], the hawk-dove game [14], and prebiotic evolution [15]. 2. SIMULATION OF TEMPORAL PROCESSES We now set up a model considering first a predator of pe￾riod X interacting with a prey of period Y. We assign a mo￾mentary fitness fy(t) of the prey in a year t as follows: it is zero if it is not present, it is 11 if both predator and prey are present, and it is +1 if the prey is present but the predator is not. Note that we punish prey that appears and meets a predator (fy(t) = 11), compared with the case of nonemer￾gence of the prey (fy(t) = 0); we do this because emergence uses up metabolic resources, because of metamorphosis, mating, and death; these resources are lost if the prey is eaten up by a predator, whereas they are preserved if the prey stays as larvae below the ground. The momentary predator fitness fx(t) is defined analogously as for the prey, but with opposite signs. The fitness Fx, resp. Fy, is defined by the sum over the fx(t), resp. fy(t), t = 0, . . . ,XY, divided by the number of predator, resp. prey, generations. (Note that this yields an average valid for t→`, because the process is pe￾riodic with period XY). We divided by the number of gen￾erations because we found that no such division favors short generation times, instead of properly timed genera￾tions (prime cycle lengths). This is important if one consid￾ers that each generation uses considerable metabolic en￾ergy, as mentioned above, and these expenses should be minimized in the long run. Assuming that prey larvae can survive for a long time below the ground, our model favors seldom emergences, as long as the prey is safe when it ap￾pears. A similar reasoning applies to the predator, e.g., fungi attacking eggs of prey; such a predator can survive for a long time as spores, and our model favors seldom appearances, as long as they get nourishment when they appear. Never￾theless, this model does not cause the cycles to become longer and longer, because prey cycles eventually get locked into a prime number, as we will show below, bringing evo￾lution to a stop. We compare now a prey mutating to a cycle Y8 with the resident prey (cycle length Y ) at constant X. Analogously, we compare mutant cycles X8 with resident cycles X at constant Y. (Note that we do not restrict cycle length changes to 51.) A mutant prey (resp. predator) substitutes the resident if and only if Fy8 > Fy, resp. Fx8 > Fx. Thus, in the case of fitness equality, the resident and not the mutant is selected. Here and in the rest of this work, we assume that all interacting populations are synchronized, thus being all present at t = 0. Now we proceed to analyze the model described in the previous paragraph. Let lcm(X,Y ): least common multiple, gcd(X,Y ): greatest common divisor of X and Y. In XY years, the predator appears Y times, both predator and prey ap￾pear XY/lcm(X,Y ) times; thus, predator without prey ap￾pears Y 1 XY/lcm(X,Y ) times. Considering that gcd(X,Y ) ? lcm(X,Y ) = XY, we thus obtain the predator fitness Fx(X,Y ) = 2gcd(X,Y )/Y 1 1. Analogously, one obtains the prey fitness Fy(X,Y )=1 1 2gcd(X,Y )/X. We will now show the following: if we allow random mutations—which can be of any size, as long as they lead to mutants within the ranges 2 # X # L/2, L/2 + 2 # Y # L—then a sequence of such mutations will finally lock the prey into a stable prime period Y. We call H the domain defined by these allowed ranges of X and Y. The definition of H on the X 1 Y plane is given in Figure 1. To be more precise, we will show that if Y is not a prime, then there exists a sequence of mutations that will change Y, and if Y is prime, then no mutation will change it. Let us first assume that Y = YN is not a prime; Fx(X,YN) has the maximum value 2gcd(XM,YN)/YN 1 1 at the predator period XM = gd(YN) (gd: greatest divisor). Note that 1 < XM < L/2 + 1 → (XM,YN) e H. A sequence of random mutations keeping Y = YN constant will eventually lead to XM. However, (XM,YN) is abandoned if mutations lead to (XM,YN 5 1). In fact, gcd(XM,YN) = XM, implying that Fy(XM,YN) = 11; gcd(XM,YN 5 1) cannot be equal to XM (the reason is: (YN 5 1)/XM = YN/XM 5 1/XM, the first term being an integer, but the sec￾ond not, so that XM is not a divisor of YN 5 1) and gcd(XM,YN Definition of the domain H within which mutations of X (predator cycle length; abscissa) and Y (prey cycle length; ordinate) are allowed. FIGURE 1 34 C O M P L E X I T Y © 2001 John Wiley & Sons, Inc

unstable with respect topreymutations.This can be shownFIGURE2easily:gcd(jY,Y,) =Ypwhereas gcd(jYpY,-k) with k =1,2,3,cannotbelargerthanYp-k,thusF,(jYpYp)-1 =F,(XyYN). Sothepreyofeachcell arereplacedbythefittestamongthefar, we have shown that there exists a sequence of muta-neighbors. The neighbourhood is defined by the cell itself,tionssuchthat apreywitha nonprime cycleY=Yisand the eight cells around it.We call C, (i= 1, ...,9) the9extinguished.Assume now that Y=Y,is a prime; anyX-XneighborsofCThemomentaryfitnessF(t)ofapredatorinsuch that (XYp) eH satisfies the condition 1≤X≤Yp-thetimesteptiscomputedhereasfollows:f(t)=oifthe1 and is thus relatively prime to Y,: therefore gcd(XY,) = 1,predator does not appear in C, in that time step;if theso that starting from (Xr,Yp) eH, there exist no predatorpredatorappears in C,and thenumbervof cells in the3×mutants that are fitter than a resident. On the other hand,3neighbourhood of C,ocuppied bypreyis notzero(1≤gcd(XY)≥1,whereY isapreymutant, ascompared to≤9), then f(t) =v; f(t) = -p if the predator appears in C,gcd(XYp)=1,sothatF(XY)≤F(XYp), i.e., no preyand v=0.p is a natural numberdescribing a"punishment"mutant is fitter than a resident.In conclusion,if H containsfor a predator that appears but finds no prey at all Theaprimepreycyclelength,any initial random choiceof (X,Y)momentary fitness f,(t)of the prey is computed analo-e H will lead and lock to a prime Y aftera sufficiently largegously,but withopposite signs.Thefitness F,resp.Fy,ofanumberofmutations.Thislocking intoprimes is illustratedpredator,resp.prey, in C, are given by the sum of the f(t),in Figure 2.Figure 2a is intended to illustrate a biologicalresp.f(t), over all t, tranging from 1 to theproduct of all 9process, so we chose a low value of L (L =22) and obtainedcycle lengths interacting in the neighborhood of C;thislockingatY=17.Incontrast,Figure2bisintended todem-onstrate a case of number-theoretical relevance, so wesum is then divided bythe number of generations of thechose L=2.2X109and obtained locking of Yinto the Eulerpredator,resp.theprey.In order to updateXand Yfor a cellprime E.C, we perform the evaluation just described in all 9 neigh-At this point we would like to comment on our restric-boring cells C,We then replace X, resp.Y, in C by the valuetion to the domain H. We imposed this restriction becauseof X, resp.Y, of the cell C, that yielded the largest fitness Fthepoints(jYp,Yp),whereYpis primeand j=1,2,3,*,areresp.Fy2001 John Wiley & Sons, Inc.35COMPLEXITY

5 1) can certainly not be larger than XM; thus gcd(XM,YN 5 1) 11 = Fy(XMYN). So far, we have shown that there exists a sequence of muta￾tions such that a prey with a nonprime cycle Y = YN is extinguished. Assume now that Y = YP is a prime; any X = XR such that (XR,YP) e H satisfies the condition 1 # XR # YP 1 1 and is thus relatively prime to YP; therefore gcd(XR,YP) = 1, so that starting from (XR,YP) e H, there exist no predator mutants that are fitter than a resident. On the other hand, gcd(XR,Y8) $ 1, where Y8 is a prey mutant, as compared to gcd(XR,YP) = 1, so that Fy(XR,Y8) # Fy(XR,YP), i.e., no prey mutant is fitter than a resident. In conclusion, if H contains a prime prey cycle length, any initial random choice of (X,Y ) e H will lead and lock to a prime Y after a sufficiently large number of mutations. This locking into primes is illustrated in Figure 2. Figure 2a is intended to illustrate a biological process, so we chose a low value of L (L = 22) and obtained locking at Y = 17. In contrast, Figure 2b is intended to dem￾onstrate a case of number-theoretical relevance, so we chose L = 2.2 2 109 and obtained locking of Y into the Euler prime E. At this point we would like to comment on our restric￾tion to the domain H. We imposed this restriction because the points (jYP,YP), where YP is prime and j = 1, 2, 3, . . . , are unstable with respect to prey mutations. This can be shown easily: gcd(jYP,YP) = YP, whereas gcd(jYP,YP 1 k) with k = 1, 2, 3, . . . cannot be larger than YP 1 k; thus Fy(jYP,YP) < Fy(jYP,YP 1 k). This means that convergence to a prey with period YP is not possible if mutations to the points (jYP,YP) are permitted. In order to discard these points, one could restrict the system to points above the diagonal X = Y. This, however, is not plausible if we want the model to be appli￾cable to biological systems; in fact, this would mean that the limits for prey cycle lengths depend on predator cycles lengths, and vice versa. Such dependences are avoided by re￾stricting mutations to the rectangle H. In addition, H ful￾fills the requirement that it contains XM = gd(YN) (see above). Note that, as a purely number-theoretical game, i.e., without biological considerations, one may loosen the re￾striction to the domain H and allow any mutations that lead to (X,Y ) pairs above the main diagonal X = Y. The restriction then would read 2 # X < Y. 3. SPATIO-TEMPORAL SIMULATIONS We will now modify the model so that prime numbers are selected by competition between neighboring residents in a spatially extended system, instead of competition between mutants and residents. In this spatio-temporal model there are no random mutations; instead we start the process with cycles X and Y that are randomly distributed in space. Fig￾ure 3 shows results obtained with a cellular automaton (CA) evolving from such a random configurations in a two￾dimensional habitat with cyclic boundary conditions. The CA is updated after a time Dt, which is equal to the lcm of all X and Y at time t in the CA. At time t + Dt, the predator and the prey of each cell are replaced by the fittest among the neighbors. The neighbourhood is defined by the cell itself, and the eight cells around it. We call Ci (i = 1, . . . , 9) the 9 neighbors of C. The momentary fitness Fx(t) of a predator in the time step t is computed here as follows: fx(t) = 0 if the predator does not appear in Ci in that time step; if the predator appears in Ci and the number v of cells in the 3 2 3 neighbourhood of Ci ocuppied by prey is not zero (1 # v # 9), then fx(t) = v ; fx(t) = 1p if the predator appears in Ci and v = 0. p is a natural number describing a “punishment” for a predator that appears but finds no prey at all. The momentary fitness fy(t) of the prey is computed analo￾gously, but with opposite signs. The fitness Fx, resp. Fy, of a predator, resp. prey, in Ci are given by the sum of the fx(t), resp. fy(t), over all t, t ranging from 1 to the product of all 9 cycle lengths interacting in the neighborhood of Ci; this sum is then divided by the number of generations of the predator, resp. the prey. In order to update X and Y for a cell C, we perform the evaluation just described in all 9 neigh￾boring cells Ci . We then replace X, resp. Y, in C by the value of X, resp. Y, of the cell Ci that yielded the largest fitness Fx, resp. Fy. Evolution to the prey perod Y = 17 (a) and to Y = E = 2147483647 (b; prime number discovered by Euler) through mutation￾selection sequences (abscissa: number of time steps t; left ordi￾nate: prey period Y; right ordinate: predator period X). FIGURE 2 © 2001 John Wiley & Sons, Inc. C O M P L E X I T Y 35

cycles in thebackground X, Y and those in the wave XwFIGURE3YwThis is graphically explained in Figure 4.A one-dimensional waveis composed offivezonesdefinedonthecells denoted by i= 1, 2,...,n,... n+m,..., where m>2isthe width of the wave: (i) Xg,Ygat i= 1,2,*..,n, (i) X,Yw(aati=n+l; (i)XwYwati=n+2,...,n+m, (iv)XwYati= n+ m+ l; (u) Xg, Ygat i= n + m+ 2, ...Note that thet=ot=172t=187predator wave (i= n+ 2, ..., n+ m+ 1) is displaced one cellto the right relatively to the prey wave (i= n+ l,., n+ m);this determines the moving direction to the right, as we willexplainnow.Consideringthatpredator-preyinteractions(b)occurhereonlywiththetwoimmediateneighborsofeachcell, the predator at i= n+ m+1 (resp. at i= n+ 1) can feedt=0t=489t=493on two types of prey and thus has a larger fitness than thepredator at i= n+ m+ 2 (resp. at i= n + 2), which can onlySpatiotemporal dynamics of the cellular automaton starting fromfeed on onetypeofpreytherefore,thepredator wavewilla random distribution and evolving to gliders (a) and to a periodicattractor consisting of a rotating spiral surrounded by pulsatingmove one cell to the right in the next time step.The prey atsmaller domains (b). At t= 0, one predator (with cycle X) and onei=n+m+1 (resp.theprey at i=n+1) can be eaten bytwoprey (with cycle Y) are placed in each cell by random choicetypes of predators and thus has a lowerfitness than the preywithin VE<X<E-2, E-8<Y<E+ 8, where Eis the Euler-primeat i= n+m (resp.at i= n),which can onlybe eaten by one(see Fiqure 2b).Number of cells:64 x 64.Only Yis shown heretypeofpredator;therefore,also theprey wave will move oneDifferent grey shadings correspond to different Y, black corre-sponding to Y = E.celltotheright.Intwodimensionsthemechanismismorecomplicated (especiallyforgliders moving diagonally,as inFigure3a),buttheycanbeunderstoodbythesametypeofThe CA just described yields a large diversity of coexist-reasoning steps.It is remarkable that one-dimensionaling attractors, depending solely on the choice of the initial,waves and spiral waveshere showabehavior similartorandomly distributed cycle lengths. As in the"Prisoner'swaves in excitable media, such as chemical reactions, heartDilemma" CA[11], thefate of each cell depends here on 25muscle,and epidemics (seee.g-,Refs.18-22 and referencesneighbors; in fact, the fitness of each cell is evaluated by thetherein). This functions as follows: the spatial domain,encounters among the inhabitants of this cell and of its8which is in the excited state, becomes refractory; the do-neighbors,butaftereachtimesteptheinhabitantsofacellmain in the refractory state becomes excitable;thedomainare replaced by those having the largest fitness, consideringin the excitable state becomes excited; and then the cyclethe fitness of that cell and that of each of its 8 neighbors.starts again. Also in the spirals found here, there are threeThis contrasts with Conway's"Game of Life"[16], wherespatial domains characterized by cycle pairs (X,Y), (X2,Y2),only 9 cells specify a cell's fate.After a sufficient number ofand (Xg,Yg)that undergo a similar cycle:by virtue of ourCAupdatings,weobtainhomogeneity,travelingwaves(in-fitnesscriterion,thefirstreplacesthesecond,thesecondcluding"gliders" as those shown in Figure 3a), spiral wavesreplaces the third, and the third replaces the first; this se-(Figure 3b)and periodical,complicatelypulsating areas,quential replacements cause the revolving of the wave.suchasthosesurroundingthespiral inFigure3b.Because(Note that grey levels only display Y in Figure 3). Spiralall solutions displayperiodicity (traveling waves do so be-causeofthe cyclicboundaryconditions),theexistence ofanattractor is numerically well defined.Wefound a predomi-FIGURE4nant appearance of prime-periodic preys, either in thewaves or in the background. In the examples given in FigureXa.YgXg.YwXw.YwXw,Y,Xg.Y3 it is the Euler prime E that appears (black) in the back-ground of the figure.Our particular choice of initial, ran-n.1#+2namdomly distributed Ywhich is given in the caption-aimedPrey wavetoward the locking into the large prime E.Alarger interval ofinitial Yvalues causes locking into other primes.As in thePredator waveworkofHasseletal.[17],theemergentpatchinessallows,asinnature,thepersistence ofotherwiseunstablepopula-Spatial distribution of cycle lengths for a one-dimensional travel-tions, coexisting in space.ing wave (Xw. Yw) and for the backround (Xg, Ys). i, cell index.Thetraveling waves can be better understoodby consid-Thewavemovestotheright.ering a CA with only one spatial dimension.We call the36COMPLEXITY 2001 John Wiley & Sons, Inc

The CA just described yields a large diversity of coexist￾ing attractors, depending solely on the choice of the initial, randomly distributed cycle lengths. As in the “Prisoner’s Dilemma” CA [11], the fate of each cell depends here on 25 neighbors; in fact, the fitness of each cell is evaluated by the encounters among the inhabitants of this cell and of its 8 neighbors, but after each time step the inhabitants of a cell are replaced by those having the largest fitness, considering the fitness of that cell and that of each of its 8 neighbors. This contrasts with Conway’s “Game of Life” [16], where only 9 cells specify a cell’s fate. After a sufficient number of CA updatings, we obtain homogeneity, traveling waves (in￾cluding “gliders” as those shown in Figure 3a), spiral waves (Figure 3b) and periodical, complicately pulsating areas, such as those surrounding the spiral in Figure 3b. Because all solutions display periodicity (traveling waves do so be￾cause of the cyclic boundary conditions), the existence of an attractor is numerically well defined. We found a predomi￾nant appearance of prime-periodic preys, either in the waves or in the background. In the examples given in Figure 3 it is the Euler prime E that appears (black) in the back￾ground of the figure. Our particular choice of initial, ran￾domly distributed Y—which is given in the caption—aimed toward the locking into the large prime E. A larger interval of initial Y values causes locking into other primes. As in the work of Hassel et al. [17], the emergent patchiness allows, as in nature, the persistence of otherwise unstable popula￾tions, coexisting in space. The traveling waves can be better understood by consid￾ering a CA with only one spatial dimension. We call the cycles in the background XB, YB and those in the wave XW, YW. This is graphically explained in Figure 4. A one￾dimensional wave is composed of five zones defined on the cells denoted by i = 1, 2, . . . ,n, . n + m,., where m > 2 is the width of the wave: (i) XB, YB at i = 1, 2, . . . ,n; (ii) XB, YW at i = n + 1; (iii) XW, YW at i = n + 2, . . . ,n + m; (iv) XW, YB at i = n + m + 1; (v) XB, YB at i = n + m + 2, . . . . Note that the predator wave (i = n + 2, . . . , n + m + 1) is displaced one cell to the right relatively to the prey wave (i = n + 1, . . . , n + m); this determines the moving direction to the right, as we will explain now. Considering that predator-prey interactions occur here only with the two immediate neighbors of each cell, the predator at i = n + m + 1 (resp. at i = n + 1) can feed on two types of prey and thus has a larger fitness than the predator at i = n + m + 2 (resp. at i = n + 2), which can only feed on one type of prey; therefore, the predator wave will move one cell to the right in the next time step. The prey at i = n + m + 1 (resp. the prey at i = n + 1) can be eaten by two types of predators and thus has a lower fitness than the prey at i = n + m (resp. at i = n), which can only be eaten by one type of predator; therefore, also the prey wave will move one cell to the right. In two dimensions the mechanism is more complicated (especially for gliders moving diagonally, as in Figure 3a), but they can be understood by the same type of reasoning steps. It is remarkable that one-dimensional waves and spiral waves here show a behavior similar to waves in excitable media, such as chemical reactions, heart muscle, and epidemics (see e.g., Refs. 18–22 and references therein). This functions as follows: the spatial domain, which is in the excited state, becomes refractory; the do￾main in the refractory state becomes excitable; the domain in the excitable state becomes excited; and then the cycle starts again. Also in the spirals found here, there are three spatial domains characterized by cycle pairs (X1,Y1), (X2,Y2), and (X3,Y3) that undergo a similar cycle: by virtue of our fitness criterion, the first replaces the second, the second replaces the third, and the third replaces the first; this se￾quential replacements cause the revolving of the wave. (Note that grey levels only display Y in Figure 3). Spiral Spatiotemporal dynamics of the cellular automaton starting from a random distribution and evolving to gliders (a) and to a periodic attractor consisting of a rotating spiral surrounded by pulsating smaller domains (b). At t = 0, one predator (with cycle X) and one prey (with cycle Y ) are placed in each cell by random choice within √E < X < E − 2, E −8< Y < E + 8, where E is the Euler-prime (see Figure 2b). Number of cells: 64 × 64. Only Y is shown here. Different grey shadings correspond to different Y, black corre￾sponding to Y = E. FIGURE 3 Spatial distribution of cycle lengths for a one-dimensional travel￾ing wave (XW, YW) and for the backround (XB, YB). i, cell index. The wave moves to the right. FIGURE 4 36 C O M P L E X I T Y © 2001 John Wiley & Sons, Inc

results, a5x5grid hasmorethanonemaximum:localFIGURE5maximumagainatY=17andaglobal maximumatY=29,the latter being a prime cycle length that has not been ob-PY)served for cicadas; note that for this small grid, nonprimes3(Y=25 and Y=49) are selected with higher probabilities0.2than for larger grids. Moreover, we found that such non-0.5prime selection becomes increasingly more pronounced as0.09the grid isfurtherreducedThus,a 10 Xl0grid (as in0.06Figures 5 and 6a) or a 20 x 20 grid (as in Figure 6b) are largerthan the smallest square grid for which our results are valid.0.040.024.DISCUSSION2720S50Qualitatively,thepeak shapes illustrated inFigure5and6ProbabilitiesP()thatpreyperiodsYareselectedafterthecellularareexplained as follows.Prey with sufficiently small Yareautomaton has converged to an attractor. At t= O, cycle lengthsexposed to predators having periods jY≤U, j=1, 2,3, are randomly chosen within 2 ≤ X, Y ≤ 50 and randomly dis-tributed among 10 x 10 cells. 10,000 different initial spatial con-figurations are evaluated. Symbols corresponding to prime, resp.FIGURE6nonprime, Y values are shown in black, resp., gray. Squares:p= 5, i.e., predators emerging but getting no prey lose the aver-P(Y)123age of what they would gain (between 1 and 9) if they found preyin the neighborhood; analogously: prey emerging but meeting no0.00predators gain the average of what they would lose (between 1and 9) if predators emerge in the neighborhood. Triangles: p = 0,0.0i.e., predators emerging but finding no prey are undisturbed,whereas prey emerging but meeting no predators are not re-(a)0.04warded. Stars: modified model, so that f(t) = +1 if the predatoremerges and finds prey, independently of the number v of prey-0.02populated neighboring cells (satiation effect) and f(t) = -1 ifemerging prey meets predators, no matter how many.p(Y)waves are also obtained in host-parasitoid dynamics [17]0.2and in prebiotic evolutionary automata [15].We investigate now the capability of our CAto produceprime numbers.For this, we determine the probabilitiesP(Y) that different prime Yvalues appear in the attractors of(b)0.1theCA.Figure5showsresultsofsuchanalyses.Oneclearlysees thereamuchmorefrequentappearanceof primepreycycles (dark symbols), especially with periods around 17,compared with nonprime prey cycles (grey symbols, all be-ing very close to the abscissa). This result is robust to drasticvariations of the model, these variations being displayed byPY!0.2the symbol shapes in Figure 5 (squares, triangles, and stars)and by Figure 6. Note that in the spatiotemporal modelleading to Figures5 and 6,we only restrict populations byupper bounds (2 ≤ X,Y≤ U; U= 50 in the case of Figures 5,(c)0.16b, and 6c; U=100 in the case of Figure 6a)We investigated the influenceof theupper bound Uandof the grid size of the CA for p=5 (ranges: 25≤ U≤200; gridsize between 10 × 10 and 30 X 30).Wefound that none ofAYthese parameters change the peak-like shape of P(Y) vs.Y.143 0Moreover, thevalueof YforwhichP(Y)is maximum isaResults fromthe same model as that correspondingtotheprime number,independently of U and the grid size.Insquares in Figure 5 (10 × 10 cells, U= 50), but with different gridFigure 6 we show examples for the influence of these pa-sizes and upper bounds of cycles. (a) 10 × 10 cells, U= 100; (b)rameters.Note that the maximum is attained at Y=17 in20 × 20 cells, U = 50; (c) 5 × 5 cells, U = 50.Figure 6a and at Y=13 in Figure 6b.In contrast to these@ 2001 John Wiley & Sons, Inc.37COMPLEXITY

waves are also obtained in host-parasitoid dynamics [17] and in prebiotic evolutionary automata [15]. We investigate now the capability of our CA to produce prime numbers. For this, we determine the probabilities P(Y ) that different prime Y values appear in the attractors of the CA. Figure 5 shows results of such analyses. One clearly sees there a much more frequent appearance of prime prey cycles (dark symbols), especially with periods around 17, compared with nonprime prey cycles (grey symbols, all be￾ing very close to the abscissa). This result is robust to drastic variations of the model, these variations being displayed by the symbol shapes in Figure 5 (squares, triangles, and stars) and by Figure 6. Note that in the spatiotemporal model leading to Figures 5 and 6, we only restrict populations by upper bounds (2 # X,Y # U; U = 50 in the case of Figures 5, 6b, and 6c; U = 100 in the case of Figure 6a). We investigated the influence of the upper bound U and of the grid size of the CA for p = 5 (ranges: 25 # U # 200; grid size between 10 2 10 and 30 2 30). We found that none of these parameters change the peak-like shape of P(Y ) vs. Y. Moreover, the value of Y for which P(Y ) is maximum is a prime number, independently of U and the grid size. In Figure 6 we show examples for the influence of these pa￾rameters. Note that the maximum is attained at Y = 17 in Figure 6a and at Y = 13 in Figure 6b. In contrast to these results, a 5 2 5 grid has more than one maximum: local maximum again at Y = 17 and a global maximum at Y = 29, the latter being a prime cycle length that has not been ob￾served for cicadas; note that for this small grid, nonprimes (Y = 25 and Y = 49) are selected with higher probabilities than for larger grids. Moreover, we found that such non￾prime selection becomes increasingly more pronounced as the grid is further reduced. Thus, a 10 2 10 grid (as in Figures 5 and 6a) or a 20 x 20 grid (as in Figure 6b) are larger than the smallest square grid for which our results are valid. 4. DISCUSSION Qualitatively, the peak shapes illustrated in Figure 5 and 6 are explained as follows. Prey with sufficiently small Y are exposed to predators having periods jY # U, j = 1, 2, 3, . . . ; Probabilities P(Y) that prey periods Y are selected after the cellular automaton has converged to an attractor. At t = 0, cycle lengths are randomly chosen within 2 # X, Y # 50 and randomly dis￾tributed among 10 × 10 cells. 10,000 different initial spatial con￾figurations are evaluated. Symbols corresponding to prime, resp. nonprime, Y values are shown in black, resp., gray. Squares: p = 5, i.e., predators emerging but getting no prey lose the aver￾age of what they would gain (between 1 and 9) if they found prey in the neighborhood; analogously: prey emerging but meeting no predators gain the average of what they would lose (between 1 and 9) if predators emerge in the neighborhood. Triangles: p = 0, i.e., predators emerging but finding no prey are undisturbed, whereas prey emerging but meeting no predators are not re￾warded. Stars: modified model, so that fx(t) = +1 if the predator emerges and finds prey, independently of the number v of prey￾populated neighboring cells (satiation effect) and fy(t) = −1 if emerging prey meets predators, no matter how many. FIGURE 5 Results from the same model as that corresponding to the squares in Figure 5 (10 × 10 cells, U = 50), but with different grid sizes and upper bounds of cycles. (a) 10 × 10 cells, U = 100; (b) 20 × 20 cells, U = 50; (c) 5 × 5 cells, U = 50. FIGURE 6 © 2001 John Wiley & Sons, Inc. C O M P L E X I T Y 37

the number of such predators decreases on increasing Y,REFERENCESroughly explaining the left branch of the P(Y) peak. Note in1.R.M.May.Periodical cicadas.Nature277,1979,pp.347349this contextthat we showed above that the pairs (jYpYp),2.J.D.Murray.Mathematical biology,Springer, Berlin,1989where Y, is prime, are unstable. By virtue of this, the prey3.M.Lloyd and H.S.Dybas.The periodical cicada problem. Il.with Y=17 is subject to the predators with X= 17 and withEvolution.Evolution 20,1966,pp.4665054.F.C.HoppensteadtandJ.B.Keller.SynchronizationofperiodicalX = 34; however, the number of (initially equally probable)cicada emergences.Science 194,1976, pp.335-337.favorablepredators is sufficientlylarge to allowfora maxi-5.M.G. Bulmer. Periodical insects. Am Nat 111, 1977,pp.1099mumofP(Y)at Y=17inFigures5and6a.Inorderto1117explain qualitatively the right branch of the peak, let us6. J.Yoshimura.The evolutionary origins of periodical cicadascompare a largeprimeY (sayY=47in Figure6b)with aduringiceages.AmNat149,1997,pp.112-1247.A.Barnett.Odd timing keeps cicadas in the prime.New Sci153,medium prime (say Y= 13). As the CA proceeds in time, theMarch 1997.p.18.prey with Y= 47 will favor the survival of one predator in its8.R.T.CoxandC.E.Carlton.Acommentaryonprimenumbersandspatial neighborhood, namely that having X = 47, becauselivecyclesof periodical cicadas.Am Nat152,1998,pp.162164.only this predator can feed well on this prey; however, this9.R.T.Cox and C.E.Carton.Paleoclimatic influences in the evolutionpredator will meet the prey in each prey generation and willof periodical cicadas. Am Midland Nat 120, 1988. pp. 183-193.10.J.Maynard Smith.Evolution and the theory of games,Cam-thus stronglydecimatetheprey, ifthis prey appears,thenF,bridge Univ. Press, 1982.= v. The prey with the smaller Y= 13 willfavor three types11.M.ANowak and R.M.May.Evolutionarygamesand spatialof predators in its spatial neighborhood,namelythosewithchaos.Nature 359,1992.pp.826-829X=13,26,and39.Ifwemakethesimplificationthatthat12.M.A.Nowak andK.Sigmund.Astrategyofwin-staylose-shiftthese three are equally distributed, i.e, each has an averagethat out-performs tit-for-tat in the Prisoner's dilemma game.Nature364,1993,pp.5658.number v/3 in the prey's neighborhood and if we consider13.A.V.M.Herz.Collective phenomena in spatially extended evolu-that X= 26 only hits the prey every second prey generationtionary games.JTheor Biol 169,1994,pp.6587.and X=39 only everythird preygeneration,then thefitness14.M.A.Nowak and R.M.May.The spatial dilemmas of evolution.of the prey (which is evaluated per preygeneration) is F,=Int JBifurcation Chaos 3,1993,pp.35-78.15.M.C.Boerlijst and P.Hogeweg.Spiralwave structure in prebi--(1+1/2+1/3)/3=-u(11/18).Thus,thepreywithY=13otic evolution:hypercycles stable against parasites.Physica Disfitter thanthat with Y=47,exemplifyingthe decreaseof48,1991,pp.1728P(Y) at large Y(right side of the peak). Weleave it as an open16.E.Berlekamp.J.Conway.and R.Guy.Winningways, Vol.2task toformalizethese explanationsofthepeaks,which areAcademicPress,NewYork,1982based so far only on estimates and on empirical observa-17.M.PHassell,HN.CominsandR.M.May.Spatialstructureandtions in the course of our CA simulations.chaos in insect population dynamics.Nature 353, 1991, pp.255-258Ourresults,bothfromthepurelytemporal and fromthe18V.S.Zykov.Simulationsof wavephenomena inexcitablemediaspatiotemporal model, suggest that there are generic prop-ManchesterUniv.Press,1988.erties of this type of dynamics that favor prime numbers.19.A.V.Holden,M.Markus,and H.G.Othmer,Eds.Nonlinear waveAlthough there are traditional methods for prime numberprocesses in excitable media, Plenum, New York, 1991detection (see e.g., Ref.23)that arefaster than the methods20.M.Markus,and B.Hess.Isotropic cellular automaton formod-ellingexcitablemedia.Nature347,1990,pp.56-58presented here, it isremarkablethatthe generation ofprime21.M.Markus,Zs.Nagy-Ungvarai,andB.Hess.Phototaxisofspiralnumbers canbeperformed usingabiological modelwaves.Science257,1992,pp.225-22722.M.Markus,G.Kloss,andI.Kusch.Disorderedwaves in a homogeneous, motionless excitable medium.Nature 371,1994,ACKNOWLEDGEMENTSpp.402404.We thank the Deutsche Forschungsgemeinschaft and23. P.Ribenboim. The new book of prime number records,FONDAP (Chile)forfinancial support.Also,we thank MalteSpringer,NewYork,1991Schmick and David Bressoud for fruitful interactions38COMPLEXITY 2001 John Wiley & Sons, Inc

the number of such predators decreases on increasing Y, roughly explaining the left branch of the P(Y ) peak. Note in this context that we showed above that the pairs (jYP, YP), where YP is prime, are unstable. By virtue of this, the prey with Y = 17 is subject to the predators with X = 17 and with X = 34; however, the number of (initially equally probable) favorable predators is sufficiently large to allow for a maxi￾mum of P(Y ) at Y = 17 in Figures 5 and 6a. In order to explain qualitatively the right branch of the peak, let us compare a large prime Y (say Y = 47 in Figure 6b) with a medium prime (say Y = 13). As the CA proceeds in time, the prey with Y = 47 will favor the survival of one predator in its spatial neighborhood, namely that having X = 47, because only this predator can feed well on this prey; however, this predator will meet the prey in each prey generation and will thus strongly decimate the prey; if this prey appears, then Fy ≈ 1v. The prey with the smaller Y = 13 will favor three types of predators in its spatial neighborhood, namely those with X = 13, 26, and 39. If we make the simplification that that these three are equally distributed, i.e., each has an average number v/3 in the prey’s neighborhood and if we consider that X = 26 only hits the prey every second prey generation and X = 39 only every third prey generation, then the fitness of the prey (which is evaluated per prey generation) is Fy ≈ 1v(1 + 1/2 + 1/3)/3 = 1v(11/18). Thus, the prey with Y = 13 is fitter than that with Y = 47, exemplifying the decrease of P(Y ) at large Y (right side of the peak). We leave it as an open task to formalize these explanations of the peaks, which are based so far only on estimates and on empirical observa￾tions in the course of our CA simulations. Our results, both from the purely temporal and from the spatiotemporal model, suggest that there are generic prop￾erties of this type of dynamics that favor prime numbers. Although there are traditional methods for prime number detection (see e.g., Ref. 23) that are faster than the methods presented here, it is remarkable that the generation of prime numbers can be performed using a biological model. ACKNOWLEDGEMENTS We thank the Deutsche Forschungsgemeinschaft and FONDAP (Chile) for financial support. Also, we thank Malte Schmick and David Bressoud for fruitful interactions. REFERENCES 1. R.M. May. Periodical cicadas. Nature 277, 1979, pp. 347–349. 2. J.D. Murray. Mathematical biology, Springer, Berlin, 1989. 3. M. Lloyd and H.S. Dybas. The periodical cicada problem. II. Evolution. Evolution 20, 1966, pp. 466–505. 4. F.C. Hoppensteadt and J.B. Keller. Synchronization of periodical cicada emergences. Science 194, 1976, pp. 335–337. 5. M.G. Bulmer. Periodical insects. Am Nat 111, 1977, pp. 1099– 1117. 6. J. Yoshimura. The evolutionary origins of periodical cicadas during ice ages. Am Nat 149, 1997, pp. 112–124. 7. A. Barnett. Odd timing keeps cicadas in the prime. New Sci 153, March 1997, p. 18. 8. R.T. Cox and C.E. Carlton. A commentary on prime numbers and live cycles of periodical cicadas. Am Nat 152, 1998, pp. 162–164. 9. R.T. Cox and C.E. Carlton. Paleoclimatic influences in the evolution of periodical cicadas. Am Midland Nat 120, 1988, pp. 183–193. 10. J. Maynard Smith. Evolution and the theory of games, Cam￾bridge Univ. Press, 1982. 11. M.A. Nowak and R.M. May. Evolutionary games and spatial chaos. Nature 359, 1992, pp. 826–829. 12. M.A. Nowak and K. Sigmund. A strategy of win-stay, lose-shift that out-performs tit-for-tat in the Prisoner’s dilemma game. Nature 364, 1993, pp. 56–58. 13. A.V.M. Herz. Collective phenomena in spatially extended evolu￾tionary games. J Theor Biol 169, 1994, pp. 65–87. 14. M.A. Nowak and R.M. May. The spatial dilemmas of evolution. Int J Bifurcation Chaos 3, 1993, pp. 35–78. 15. M.C. Boerlijst and P. Hogeweg. Spiral wave structure in prebi￾otic evolution: hypercycles stable against parasites. Physica D 48, 1991, pp. 17–28. 16. E. Berlekamp, J. Conway, and R. Guy. Winning ways, Vol. 2, Academic Press, New York, 1982. 17. M.P. Hassell, H.N. Comins, and R.M. May. Spatial structure and chaos in insect population dynamics. Nature 353, 1991, pp. 255–258. 18. V.S. Zykov. Simulations of wave phenomena in excitable media, Manchester Univ. Press, 1988. 19. A.V. Holden, M. Markus, and H.G. Othmer, Eds. Nonlinear wave processes in excitable media, Plenum, New York, 1991. 20. M. Markus, and B. Hess. Isotropic cellular automaton for mod￾elling excitable media. Nature 347, 1990, pp. 56–58. 21. M. Markus, Zs. Nagy-Ungvarai, and B. Hess. Phototaxis of spiral waves. Science 257, 1992, pp. 225–227. 22. M. Markus, G. Kloss, and I. Kusch. Disordered waves in a ho￾mogeneous, motionless excitable medium. Nature 371, 1994, pp. 402–404. 23. P. Ribenboim. The new book of prime number records, Springer, New York, 1991. 38 C O M P L E X I T Y © 2001 John Wiley & Sons, Inc