Table Of Content

Spatiotemporal complexityofa ratio-dependent predator-prey system Weiming Wang,1,2,∗ Quan-Xing Liu,1,† and Zhen Jin1,‡ 1Department of Mathematics, NorthUniversityof China, Taiyuan, Shan’xi 030051, P.R.China 2School ofMathematicsandInformationScience, Wenzhou University, Wenzhou,Zhejiang, 325035 (Dated:February9,2008) Inthispaper,weinvestigatetheemergenceofaratio-dependentpredator-preysystemwithMichaelis-Menten- typefunctionalresponseandreaction-diffusion.WederivetheconditionsforHopf,TuringandWavebifurcation on a spatial domain. Furthermore, we present a theoretical analysis of evolutionary processes that involves organismsdistributionandtheirinteractionofspatiallydistributedpopulationwithlocaldiffusion. Theresults of numerical simulations reveal that thetypical dynamics of population density variation isthe formation of 7 isolatedgroups, i.e.,stripelikeorspottedorcoexistenceofboth. Ourstudyshowsthatthespatiallyextended 0 modelhasnotonlymorecomplexdynamicpatternsinthespace,butalsochaosandspiralwaves. Itmayhelp 0 usbetterunderstandthedynamicsofanaquaticcommunityinarealmarineenvironment. 2 n PACSnumbers: 87.23.Cc,89.75.Kd,89.75.Fb,47.54.-r a J INTRODUCTION prey attached per unit time per predator as the prey density 7 changes[10].Ingeneral,functionalresponsecanbeclassified ] Ecologicalsystemsarecharacterizedbytheinteractionbe- as:(i)preydependent,whenpreydensityalonedeterminesthe E tween species and their naturalenvironment[1]. Such inter- response,i.e.,f(x,y) = p(x); (ii)predatordependent,when P action may occur over a wide range of spatial and temporal bothpredatorandpreypopulationsaffecttheresponse.Partic- . o scales [2, 3]. The study of complex population dynamics ularly,whenf(x,y) = p(x),wecallmodel(1)strictlyratio- y bi is nearly as old as population ecology. In the 1920s, Lotka dependent; and (iii) multi-species dependent, when species - and Volterra independentlydeveloped a simple model of in- otherthanthefocalpredatoranditspreyspeciesinfluencethe q teracting species that still bears their joint names. This was functionalresponse [11]. Differingfrom the prey-dependent [ anearlylinearmodel,butthepredator-preyversiondisplayed predator-preymodels, the ratio-dependentpredator-preysys- 1 neutrallystablecycles[4, 5]. Fromthenon, thedynamicre- tems have two principal predictions: (a) equilibrium abun- v lationshipbetweenpredatorsandtheirpreyhaslongbeenand dances are positively correlated along a gradient of enrich- 2 1 will continue to be one of dominantthemes in both ecology ment and (b) the “paradoxof enrichment”either completely 0 and mathematical ecology due to its universalexistence and disappearsorenrichmentislinkedtostabilityinamorecom- 1 importance[6,7,8]. plexway [12, 13]. The ratio-dependentpredator-preymodel 0 Predator-preymodelsfollowtwogeneralprinciples: oneis hasbeenstudiedbyseveralresearchersrecentlyandveryrich 7 that population dynamics can be decomposed into birth and dynamicshavebeenobserved[7,8,12,14,15]. 0 deathprocesses;theotheristheconservationofmassprinci- / o Ontheotherhand,patternformationinnonlinearcomplex ple,statingthatpredatorscangrowonlyasafunctionofwhat i systems is one of the centralproblemsof the natural, social, b they have eaten [9]. With these two principles we can write and technological sciences [16]. In particular, starting with - thecanonicalformofapredator-preysystemas q the pioneering work of Segel and Jackson [17], spatial pat- : v x˙(t)=xg(x)−f(x,y)y−µ (x)x, ternsandaggregatedpopulationdistributionsarecommonin x i (1) nature and in a variety of spatio-temporalmodels with local X ( y˙(t)=γf(x,y)y−µy(y)y. ecologicalinteractions[1,18].Promulgatedbythetheoretical r a paper of Turing [19], the field of research on pattern forma- where g(x) is the per capita prey growth rate in the absence tion modeled by reaction-diffusion systems, which provides ofthepredator,µ andµ arenaturalmortalitiesofpreyand x y ageneraltheoreticalframeworkfordescribingpatternforma- predatorrespectively,f(x,y)isthefunctionalresponse. And tionin systems frommanydiversedisciplinesincluding(but γf(x,y)isthepercapitaproductionofpredatorduetopreda- not limited to) biology [16, 20, 21, 22, 23, 24, 25], chem- tion,whichisoftencalledthenumericalresponse. Thefunc- istry[26,27,28,29,30,31],physics[32,33,34,35],andso tionalresponseplaysamainroleinsystem(1):theknowledge on, seems to be a new increasingly interesting area, particu- ofthisfunctiondeterminesthedynamicsofthewholesystem larlyduringthelastdecade. andthe transferof the biomassin the predationbecauseit is proportionalto the numericalresponse. Usually one consid- InRef.[2],Neuhausersurveyssomecurrentworkonspatial ersconsumptionto bethemajordeathcausefortheprey. In mathematicalmodelsinecology. Muchofthisworkconsists thiscase µ (x) can be neglectedand set to 0 (aslongas the ofbuildingspatialdimensionsintoexistingclassicalmodels, x predatorexists)[9]. such as the Lotka-Volterra model that describes competition Inpopulationdynamics,afunctionalresponseofthepreda- between species. But the research on the spatial pattern of tor to the prey density refers to the change in the density of ratio-dependentpredator-preymodelsseemstoberare. 2 STABILITYANDBIFURCATIONANALYSIS To perform a linear stability analysis, we linearize the dynamic system (4) around the spatially homogenous fixed In this paper, we mainly focus on the ratio-dependent point(5)forsmallspace-andtime-dependentfluctuationsand predator-prey system with Michaelis-Menten-type(or expandtheminFourierspace Michaelis-Menten-Holling)functionalresponse: N(~x,t)∼n∗eλtei~k·~x, ∂N =r(1− N)N − αN/P P +D ∇2N (6) ∂t K 1+αhN/P 1 P(~x,t)∼p∗eλtei~k·~x.  =r(1− N)N − αN P +D ∇2N,  ∂∂Pt =γP+ααNhKNP −µPP++αhDN2∇2P.1 (2) andobtainthecharacteristicequation |A−k2D−λI|=0, (7) whereN∀,(NP,sPta)nd∈f[o0r,p∞re]y2\a(n0d,0p)r.edatordensity,respectively. where D1,D2aretheirrespectivediffusioncoefficients,∇2 = ∂∂x2+ D1 0 ∂ . Allparametersarepositiveconstants,rstandsformaxi- D = , (8) m∂ya2lgrowthrateoftheprey,γ conversionefficiency,µpreda- 0 D2 ! tordeathrate,K carryingcapacity,αcapturerateandhhan- andAisgivenby dlingtime. Note that αN/P is strictly correctonly for P > 0. In ∂ f ∂ f f f 1+αhN/P N P N P A= = , (9) the case of P = 0 and N > 0 we can define f(N,0) := 1 (thelimitoff(x)forx→∞). h ∂Ng ∂Pg !(n∗,p∗) gN gP ! Let where the elementsare the partialderivativesof the reaction Nˆ = αhN, Pˆ = αhN, R= rh, kineticsevaluatedatthestationarystate(n∗,p∗). NowEq.(7) γK γ2K γ (3) canbesolved,yieldingthesocalledcharacteristicpolynomial Q= hµ, S = αh, tˆ= γt. oftheoriginalproblem(Eq.4) γ γ h Forsimplicitywe willnotwritethehat(ˆ)inthe restofthis λ2−tr λ+∆ =0, (10) k k paper. Andinthesenewvariables,from(2)-(3), we arriveat thefollowingequationscontainingdimensionlessquantities: where ∂∂Nt =R(1− NS)N − PS+NSNP +D1∇2N, trk =fN+gP−k2(D1+D2)=tr0−k2(D1+D2), (11) (4)  ∂P = SN P −QP +D ∇2P.  ∂t P+SN 2 ∆ =f g −f g −k2(f D +g D )+k4D D k N P P N N 2 P 1 1 2 Moredetails about the choice of dimensionless variables in =∆ −k2(f D +g D )+k4D D , the system (2) as well as possible implicationscan be found 0 N 2 P 1 1 2 (12) in [14]. TherootsofEq.10yieldthedispersionrelation Thedimensionlessmodel(Eq.4)hasonlythreeparameters: R,whichcontrolsthegrowthrateofprey;Q,whichcontrols 1 thedeathrateofthepredator;andS,whichmeasurescaptur- λ1,2(k)= 2 trk± trk2−4∆k . (13) ingrate. (cid:16) q (cid:17) The first step in analyzing the model is to determine the It’s well known that reaction-diffusion systems have led behavior of the non-spatial model obtained by setting space to the characterization of three basic types of symmetry- derivativesequalto zero. Thenon-spatialmodelhasatmost breaking bifurcations responsible for the emergence of spa- three equilibria (stationary states), which correspond to spa- tiotemporalpatterns.Thespace-independentHopfbifurcation tially homogeneous equilibria of the full model (Eq. 4), in breaksthetemporalsymmetryofasystemandgivesrisetoos- the positive quadrant: (0,0)(totalextinct), (S,0) (extinctof cillationsthatareuniforminspaceandperiodicintime. The thepredator)anda nontrivialstationarystate (n∗,p∗) (coex- (stationary)Turingbifurcationbreaksspatialsymmetry,lead- istencetopreyandpredator),where ingtotheformationofpatternsthatarestationaryintimeand oscillatory in space. The wave (oscillatory Turing or finite- n∗ = S(R+(Q−1)S), wavelengthHopf)bifurcationbreaksbothspatialandtempo- R (5) ralsymmetry,generatingpatternsthatareoscillatoryinspace p∗ = S(1−Q)n∗ = S2(R−S+QS)(1−Q) andtime[28]. Q RQ TheHopfbifurcationoccurswhen Easy to know that n∗ is positive for all S < R , which 1−Q impliesQ<1andthereforeensuresthepositivityofp∗ [14]. Im(λ(k))6=0, Re(λ(k))=0 at k=0. (14) 3 thenwecangetthecriticalvalueofHopfbifurcationparame- terS equals R+Q−Q2 S = . (15) H 1−Q2 AttheHopfbifurcationthreshold,thetemporalsymmetryof thesystemisbrokenandgivesrisetouniformoscillationsin spaceandperiodicoscillationsintimewiththefrequency ω =Im(λ(k))= ∆ = Q(Q−1)(R−S+QS), H 0 thecorrespondingwavpelengthips 2π 2π λ = = . (16) H ωH Q(Q−1)(R−S+QS) FIG.1: (Coloronline)Bifurcationdiagramforthesystem(2)with TheTuringbifurcatiopnoccurswhen R = 0.5, Q = 0.6, D2 = 0.2. Hopfbifurcationline: SH = 3327; Turingbifurcationline:ST = 5556+5258D1withkT2 =5.48571;Wave Im(λ(k))=0, Re(λ(k))=0 at k =kT 6=0. (17) bifurcationline: SW = (1265D1+ 156)kW2 + 3327 withkW = 0.334. Turing-Hopfbifurcationpoint:(0.08864,1.15625)andTuring-Wave thecriticalvalueofbifurcationparameterS equals bifurcationpoint:(0.11655,1.21110). D D k4 +(D R+D Q(1−Q))k2 +RQ(1−Q) S = 1 2 T 2 1 T , T Q3−(kT2D2+2)Q2+Q+kT2D2 The Hopf bifurcation line, the Wave bifurcation line and (18) theTuringbifurcationlineintersectattwocodimension-2bi- where furcation points, the Turing-Hopf bifurcation point and the Turing-Wave bifurcation point. The bifurcation lines sepa- ∆ kT2 = D D0 , ratetheparametricspaceintosixdistinctdomains.Indomain r 1 2 I,locatedbelowallthreebifurcationlines, thesteadystateis andatthe Turingthreshold,the spatialsymmetryof the sys- theonlystablesolutionofthesystem. DomainIIisregionof tem is brokenand the patternsare stationary in time and os- pureTuringinstabilities,anddomainIIIispureHopfinstabil- cillatoryinspacewiththewavelength ities. IndomainIV,bothHopfandTuringinstabilitiesoccur, andindomainV,theWaveandHopfmodesarise. Whenthe 2π λ = . (19) parameterscorrespondtodomainVI,whichislocatedabove T k T allthree bifurcationlines, all threeinstabilities occur. Fig. 2 AndtheWavebifurcationoccurswhen showsthedispersionrelationsofunstableHopfmode,transition of Turingand Wave modesfrom stable to unstable. It’s Im(λ(k))6=0, Re(λ(k))=0 at k =kw 6=0. (20) easytoknowthatallthreebifurcationsaresupercritical. thecriticalvalueofWavebifurcationparameterS equals k2(D +D )+R+Q−Q2 SPATIOTEMPORALPATTERNANALYSIS S = w 1 2 , (21) W 1−Q2 Wehaveperformedextensivenumericalsimulationsofthe where spatially extended model (4) in two-dimensional space, and thequalitativeresultsareshownhere. Allournumericalsim- k2 = Q (D −D )2Q4−2(D2+D2)Q3 w 2D22(Q+1) 1 2 1 2 ulations employ the periodic Neumann (zero-flux) boundary (cid:16) conditions with a system size of 200×200 space units and 1/2 +(D12+D22+6D1D2−4D22R)Q2−4D1D2Q+4D22R . R = 0.5, Q = 0.6, D1 = 0.02, D2 = 0.2. The spatiotem- (cid:17) poraldynamicsofadiffusion-reactionsystemdependsonthe Easytoknowthat,attheWavethreshold,bothspatialandtem- choiceofinitialconditions,whichsomeauthorshaveconsid- poralsymmetriesarebrokenandthepatternsareoscillatoryin eredinconnectionwiththeproblemofbiologicalinvasionin spaceandtimewiththewavelength afewpapers[16,36,37,38],wheretheinitialconditionsare naturallydescribedbyfinitefunctionsandthedynamicsofthe 2π λW = k . (22) community mainly consists of a variety of diffusive popula- w tionalfronts. Ingeneral,therearetwoinitialconditionsused Linearstabilityanalysisyieldsthebifurcationdiagramwith for analysis of the spatial extendedsystems. One is random R=0.5, Q=0.6, D =0.2showninFig.1. spatial distribution of the species, which seems to be more 2 4 FIG.3:(Coloronline)Snapshotsofcontourpicturesofthetimeevo- lution of the prey at different instants with S = 0.6 < ST. (A) 0 iteration; (B) 5000 iterations; (C) 45000 iterations. [Additional movieformatavailablefromtheauthor] cansee thatthe stripelikespatialpatternsarise fromtheran- dominitialconditions. Afterthestripelikepatternsform(cf. FIG.2: (Coloronline) Dispersionrelationsshowing unstableHopf mode, transition of Turing and Wave modes from stable to unsta- Fig. 4(B)) they grow steadily with time until they reach cer- ble, e.g., asD1 decreased. Parameters: S = 1.2, R = 0.5, Q = tain width—armlength, and the spatial patterns become dis- 0.6, D2 = 0.2 and (1) D1=0.15; (2) D1=0.12; (3) D1=0.10; (4) tinct. Finally,thestripelikespatialpatternsprevailthewhole D1=0.07;(5)D1=0.04;(6)D1=0.02. domain (cf. Fig. 4(C) and Fig.5). Comparing the Fig. 4(C) withFig.5,wefindthattheparameterS isclosertoS ,and T thestripelikespatialpatternsaremoredistinct. Hereweomit generalfrom the biologicalpointof view (cf. Fig. 3(A) and thepre-imageofFig.5astheyaresimilartotheFig.4(A)and 4(A), 7(A)); The other is a special choice, i.e., taking the (B).Itiseasytoseethatthestationarypatternsareessentially speciescommunityinahorizontallayerasdecreasinggradu- differentfromthepreviouscase(cf. Fig.3). ally and the verticaldistribution of species homogeneous(cf. Fig. 9(A)).Inthissectionwe choosethe former,andthe lat- ter in the next section (IV. Discussion). The equations (4) are solved numerically in two-dimensional space using a finite difference approximation for the spatial derivatives and anexplicitEulermethodforthe time integrationwith a time stepsizeof∆t=0.01andspacestepsize∆h=0.25. FromtheanalysisofsectionIIandphase-transitionbifurca- tiondiagram(cf. Fig.1),theresultsofcomputersimulations show that the type of the system dynamicsis determinedby FIG.4:(Coloronline)Snapshotsofcontourpicturesofthetimeevo- the values of S and D . We run the simulations until they 1 lution of the prey at different instants with S = 0.9 < ST. (A) reachastationarystateoruntiltheyshowabehaviorthatdoes 0 iteration; (B) 5000 iterations; (C) 45000 iterations. [Additional notseemtochangeitscharacteristicsanymore. Fordifferent movieformatavailablefromtheauthor] setsofparameters,thefeaturesofthespatialpatternsbecome essentiallydifferentifS exceedsacriticalvalueS , S and T T S respectively,theydependonD . W 1 Fig.3showstheevolutionofthespatialpatternofpreyat0, 5000and45000iterations,withrandomsmallperturbationof thestationarysolutionn∗andp∗ofthespatiallyhomogeneous systemswhenS isless thantheTuringbifurcationthreshold S . Inthiscase,onecanseethatforthesystem(4),theran- T dominitialdistributionleadtotheformationofastronglyir- regular transient pattern in the domain. After the irregular patternform(cf. Fig.3(B))itgrowsslightlyand“jumps”al- FIG.5:(Coloronline)Snapshotsofcontourpicturesofthetimeevo- ternately for a certain time, and finally the chaos spiral pat- lutionof thepreyatdifferent instantswithST < S = 1.1 < SH ternsprevailoverthewholedomain,andthedynamicsofthe (45000iterations). [Additionalmovieformatavailablefromtheau- system does not undergo any further changes (cf. Fig. 3(C) thor] andtheadditionmovieforFig.3). Figures 4 and 5 show spontaneous formation of short When S < S = 1.16 < S , we find that the spotted H W stripelikeandspottedspatialpatternsemergeandcoexiststa- patternsandthestripelikepatternscoexistinthespatiallyex- bly when the bifurcation parameter S < S (Fig. 4) and tendedmodel(cf. Fig.6). T S < S < S (Fig.5). Fromthe snapshotsormovies,one Fig.7showssnapshotsofpreyspatialpatternat0,5000and T H 5 theoriginalimages. InFig.8(mid-handcolumn),shortwave- length, represented by a large circle, corresponds to Turing structures; longer wavelength, represented by a small white circle, correspondsto traveling and/or standing waves. This technique can be particularly appropriate for characterizing quasi-orderedarraysforFig.7(C). DISCUSSIONANDCONCLUSION FIG.6:(Coloronline)Coexistenceofstationaryspottedpatternsand The numericalresults correspondperfectlyto our theoret- stripelike patterns of the prey for long time run with SH < S = ical findingsthat there are a range of parameters in S −D 1 1.16<SW.[Additionalmovieformatavailablefromtheauthor] planewherethedifferentspatialpatternsemerge(cf. Fig.1). Fig.1andtheresultsofsimulationspresentthatthechaospat- ternswillpersistinthespatiallyextendedmodel(Eq.4)when 45000iterationsfortheparameterS = 1.2 > S . Although W theparametersareinthedomainI.Theboundaryofthisdo- thedynamicsofthesystemstartsfromthesameinitialcondi- maincan be computednumericallyandis shown asthe blue tionaspreviouscases, thereisanessentialdifferenceforthe line“Turing”inFig.1. Thestationarystateofstripelikepat- spatially extended model (Eq. 4). Form Fig. 7, one can see ternsexistswhen theparametersarein the domainII,where that the regular spotted patterns prevails over the whole do- itsboundarycanalsobecomputednumericallyandisshown mainatlast,andthedynamicsofthesystemdoesnotundergo asthe black line “Hopf”in Fig. 1. Theperiodicspottedpat- anyfurtherchanges(cf. theadditionalmovieforFig.7). terns appear in the domain VI, where the boundary can be computednumericallyby the wave bifurcationand is shown astheredline“Wave”inFig.1. Moreover,thereistransverse domain IV (cf. Fig. 1) in the system between the stripelike patternsand spottedpatterns, where the spottedpatternsand thestripelikepatternscoexist(cf.Fig.6). Dothestationarypatternsarisedependentontheinitialcon- ditions?Wetestthedifferentinitialconditionsforthespatially extended system, but the final spatial patterns are the same in qualitative. In those figures we find that the spatial chaos FIG.7:(Coloronline)Snapshotsofcontourpicturesofthetimeevo- lution of the prey at different instants with S = 1.2 > SW. (A) patterns come from the destruction of the spirals, when we 0 iteration; (B) 5000 iterations; (C) 45000 iterations. [Additional choosethespecialinitialcondition(cf. Fig.9)inthedomain movieformatavailablefromtheauthor] I. This phenomenon coheres with the results of the study in Refs.[16,43]. On the other hand, discrete Fourier transform is a basic We have presented a theoretical analysis of evolutionary mathematicaltool used to decomposea signal or image into processesthatinvolvesorganismsdistributionandtheirinter- different periodic components. It has been widely used for actionofspatiallydistributedpopulationwithlocaldiffusion. thespatialpatterns[39,40,41]. Wehavealsoperformednu- Ouranalysisandnumericalsimulationsrevealthatthetypical mericalinvestigationsintotwo-dimensionalspacebyFourier dynamicsof populationdensity variation is the formationof spectra. ThenumericalcomputationoftheFouriertransform isolatedgroups(stripelikeorspottedorcoexistenceofboth). isdonebythewell-establishedtwo-dimensionalFastFourier Thisprocessdependsonseveralparameters,includingS,D 1 Transform(FFT2) algorithm[42]. Spatial Fourier transform andD . Thefieldmeaningofourresultsmaybefoundinthe 2 of the stripelike and spotted patterns in Figures 4(c), 5 and dynamics of an aquatic community which is affected by the 7(c)areshownasFig.8. Anddigitalordigitizedtransmission existenceofrelativelystablemesoscaleinhomogeneityinthe electron micrographs (TEMs) are analyzed by using Matlab fieldofecologicallysignificantfactorssuchaswatertemper- (Ver.7.0). ature,salinityandbiogenconcentration. FromFig.8,wefindthatFig.4(C)andFig.5havethesame InRef.[16],theauthorsexplainedthefieldmeaningbyus- spatialfrequencyinthelengthofthespaceunitandpresence ing aquatic community in the ocean (cf. Fig. 9). Our study of one mode with differentwave numbers. On the contrary, showsthatthespatiallyextendedmodel(Eq.4) hasnotonly Fig. 7(C) has two modes with different wave numbers. The more complexdynamicpatterns in the space, but also chaos spatialfrequencyanddirectionofanycomponentinthepower patternsandspiralwaves,soitmayhelpusbetterunderstand spectrum are given in the length and direction, respectively, thedynamicsofanaquaticcommunityinarealmarineenvi- of a vector from the origin to the point on the circle. The ronment. It is also importantto distinguish between “intrin- magnitudeis depicted by a gray scale or color scale, but the sic” patterns, i.e., patternsarising due to trophic interactions unitsaredimensionlessvaluesrelatedtothetotaldarknessof likethoseconsideredabove,and“forced”patternsinducedby 6 x 107 0.45 3.5 Spatial Frequency (1/s.u.)1.05 12..12355 Percent total power0000000...123....1234555 1.5 0.5 0.05 0 1.5 0 1.5 0 0.005 0.01 0.015 0.02 0.025 Spatial Frequency (1/s.u.) Frequency (1/s.u.) x 107 0.4 Spatial Frequency (1/s.u.)0.05 12..12355 Percent total power000000...123...123555 0.5 0.5 0.05 0 0.5 0 0.5 0 0.005 0.01 0.015 0.02 0.025 Spatial Frequency (1/s.u.) Spatial Frequency (1/s.u.) x 106 0.45 Spatial Frequency (1/s.u.)00..055 01..1255 Percent total power0000000...123....1234555 0.05 0 0.5 0 0.5 0 0.005 0.01 0.015 0.02 0.025 Spatial Frequency (1/s.u.) Spatial Frequency (1/s.u.) FIG. 8: (Color online) Patterns (left-hand column), it spatial Fourier transformation (mid-hand column), and radial average of the power spectrum(right-handcolumn).Up-rowS =0.9;Mid-rowS =1.1;Low-rowS =1.2. ratio-dependent predator-prey model (Eq. 4) also represents richspatialdynamics,suchaschaosspiralpatterns,stripelike patterns, spotted patterns, coexistence of both stripelike and spotted patterns, etc. It will be useful for studying the dy- namiccomplexityofecosystems. Thiswork was supportedby the NationalNaturalScience Foundation of China under Grant No. 10471040 and the Natural Science Foundation of Shan’xi Province Grant No. FIG.9:(Coloronline)Snapshotsofcontourpicturesofthetimeevo- 2006011009. lutionofthepreyatdifferentinstantswiththespecialinitialcondition andS =0.6<ST. (A)0iteration;(B)1000iterations;(C)400000 iterations.[Additionalmovieformatavailablefromtheauthor] ∗ Electronicaddress:[email protected] the inhomogeneity of the environment. The physical nature † Electronicaddress:[email protected] oftheenvironmentalheterogeneity,andthusthevalueofthe ‡ Electronicaddress:[email protected] dispersionofvaryingquantitiesandtypicaltimesandlengths, [1] M.Baurmann,T.Gross,andU.Feudel,J.Theor.Biol.InPress canbeessentiallydifferentindifferentcases. Neuhauserand (2006). Pacala [44] formulatedthe Lotka-Volterramodelas a spatial [2] C.Neuhauser, NoticesoftheAmericanMathematical Society 47,1304(2001). model. Theyfoundthestrikingresultthatthecoexistenceof [3] S.A.Levin,B.Grenfell,A.Hastings,andA.S.Perelson,Sci- patternsis actually harderto get in the spatialmodelthan in ence275,334(1997). thenon-spatialone.Onereasoncanbetracedtohowlocalin- [4] R.M.May,Nature261,459(1976). teractionsbetweenindividualmembersofthespeciesarerep- [5] B.E.Kendall,EncyclopediaofLifeSciences(NaturePublish- resentedinthemodel. Inthisthesis,ourresultsshowthatthe ing Group, ondon, 2001), vol. 13, chap. Nonlinear dynamics 7 andchaos,pp.255–262. Epstein,Nature406(2000). [6] A.A.Berryman,Ecology73,1530(1992). [28] L. Yang, M. Dolnik, A. M. Zhabotinsky, and I.R. Epstein, J. [7] Y.KuangandE.Beretta,J.Math.Biol.36,389(1998). Chem.Phys.117,7259(2002). [8] J. D. Murray, Mathematical Biology II: Spatial Models and [29] L.Yang,M.Dolnik,A.M.Zhabotinsky,andI.R.Epstein,Phys. BiomedicalApplications,vol.18ofBiomathematics(Springer, Rev.Lett.88,208303(2002). NewYork,3edition,2003). [30] L.YangandI.R.Epstein,Phys.Rev.Lett.90,178303(pages4) [9] C. Jost, Phd-thesis, Institute National Agronomique, Paris- (2003). Grignon(1998). [31] L.YangandI.R.Epstein,Phys.Rev.E69, 026211 (pages6) [10] S.RuanandD.Xiao, SIAMJournal onAppliedMathematics (2004). 61,1445(2001). [32] W.Just,M.Bose,S.Bose,H.Engel,andE.Scho¨ll,Phys.Rev. [11] P.A.AbramsandL.R.Ginzburg,TrendsinEcologyandEvo- E64,026219(2001). lution15,337(2000). [33] M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. 65, 851 [12] R.ArditiandL.R.Ginzburg,J.Theor.Biol.139,311(1989). (1993). [13] M.L.Rosenzweig,Science171,385(1971). [34] A. Yochelis, A. Hagberg, E. Meron, A. L. Lin, and H. L. [14] C.Jost,O.Arino,andR.Arditi,Bull.Math.Biol.61,19(1999). Swinney,SIAMJournalonAppliedDynamicalSystems1,236 [15] D.M.XiaoandL.S.Jennings, SIAMJ.Appl.Math.65, 737 (2002). (2005). [35] T.Leppa¨nen,M.Karttunen,R.A.Barrio,andK.Kaski,Phys. [16] A. B. Medvinsky, S. V. Petrovskii, I. A. Tikhonova, H. Mal- Rev.E70,066202(2004). chow,andB.-L.Li,SIAMReview44,311(2002). [36] J.Sherratt,M.Lewis,andA.Fowler,PNAS92,2524(1995). [17] L.A.SegelandJ.L.Jackson,J.Theor.Biol.37,545(1972). [37] N. Shigesada and K.Kawasaki, Biological Invasions: Theory [18] M.Pascual,Proc.Roy.Soc.Lond.B251,1(1993). andPractice(OxfordUniversityPress,Oxford,1997). [19] A.M.Turing,Philos.Trans.R.Soc.LondonB237,7(1952). [38] S.V.Petrovskii andH.Malchow, Theor.Popul.Biol.59, 157 [20] C.A.Klausmeier,Science284,1826(1999). (2001). [21] R.T.Liu,S.S.Liaw,andP.K.Maini,Phys.Rev.E74,011914 [39] R. O. Prum and R. H. Torres, Integr. Comp. Biol. 43, 591 (2006). (2003). [22] J.vonHardenberg,E.Meron,M.Shachak,andY.Zarmi,Phys. [40] R.O.Prum,R.Torres,S.Williamson,andJ.Dyck,Proc.Roy. Rev.Lett.87,198101(2001). Soc.Lond.B1414,13(1999). [23] O.Lejeune,M.Tlidi,andP.Couteron,Phys.Rev.E66,010901 [41] R.O.PrumandR.Torres,J.Exp.Biol.206,2409(2003). (2002). [42] W.L.BriggsandV.E.Henson,TheDFT(SocietyforIndustrial [24] E.Sharon,M.Marder,andH.L.Swinney,AmericanScientist andAppliedMathematics,Philadelphia,Pennsylvania,1995). 92,254(2004). [43] W. S. C. Gurney, A. R. Veitch, I. Cruickshank, and [25] M. Pascual, M. Roy, and A. Franc, Ecology Letters 5, 412 G.McGeachin,Ecology79,2516(1998). (2002). [44] C. Neuhauser and S. W. Pacala, Ann. Appl. Probab 9, 1226 [26] Q.OuyangandH.L.Swinney,Nature352(1991). (1999). [27] V.K.Vanag,L.Yang,M.Dolnik,A.M.Zhabotinsky,andI.R.