Evolutionary ecology in-silico Does mathematical modelling help in understanding the generi

a r X i v :q -b i o /0504020v 1 [q -b i o .P E ] 16 A p r 2005Evolutionary ecology in-silico:Does mathematical modelling help in understanding the “generic”trends?Debashish Chowdhury 1,∗and Dietrich Stauffer 2,†1Physics Department,Indian Institute of Technology,Kanpur 208016,India.2Institute for Theoretical Physics,University of Cologne,50923K¨o ln,Germany.Abstract:Motivated by the results of recent laboratory experiments (Yoshida et al.Nature,424,303-306(2003))as well as many earlier field observations that evolutionary changes can take place in ecosystems over relatively short ecological time scales,several “unified”mathematical models of evolutionary ecology have been developed over the last few years with the aim of describ-ing the statistical properties of data related to the evolution of ecosystems.Moreover,because of the availability of sufficiently fast computers,it has become possible to carry out detailed computer simulations of these models.For the sake of completeness and to put these recent developments in the proper perspective,we begin with a brief summary of some older models of ecological phenomena and evolutionary processes.However,the main aim of this article is to review critically these “unified”models,particularly those published in the physics literature,in simple language that makes the new theories accessible to wider audience.∗E-mail:debch@iitk.ac.in†E-mail:stauffer@thp.uni-koeln.de11IntroductionEnormous progress has been made in the twentieth century in the domain of sub-cellular and cell biology,particularly in area of molecular genetics and genomics.One of the challenges of the twenty-first century will be to link the insight gained from the molecular level research on uni-cellular as well as multi-cellular organisms to biological research at higher levels of organization, namely,those at the levels of colonies,communities and,finally,eco-systems [1,2].Admittedly,at present,we are far from that goal.In traditional paleobiology,analysis of the fossil data has always been the most popular way of understanding the causes and consequences of extinction of species as well as those of biotic recoveries from mass extinctions[3,4,5, 6,7,8].Unfortunately,the available record of the history of life,written on stone in the form of fossils,is incomplete and ambiguous[5].Laboratory experiments have also played equally important role so far in ecology and evolutionary biology.However,an alternative enterprise seeks to recreate the evolution on a computer by simulating theoretical models;this is often referred to as in-silico experiments.Models are normally useful in understanding the real world.In princi-ple,models can be verbal or symbolic,graphical or abstract,qualitative or quantitative.However,throughout this paper,by the term model we shall always mean mathematical models that not only indicate qualitative features of various quantities of interest but can also make quantitative predictions. Mathematical modelling often helps in getting insight into ecological phe-nomena and evolutionary processes[9,10].Most of the ecological and evolu-tionary models are too complicated to be solved analytically;for such models computer simulation is one of the most powerful tools of analysis.It has been realized in recent years that ecosystems are examples of com-plex adaptive systems[11,12,13,14].Over the last ten years statistical physi-cists have used the conceptual toolbox of their profession for understanding some aspects of the dynamical evolution of eco-systems which include,for example,the“generic”trends in the statistics of the data on speciation and extinction[15,16,17].Significant progress has been made over the last three years in developing detailed models that incorporate not only ecological phe-momena on short periods of time but also evolutionary processes on longer time scales.In this article we present a critical overview of the current status of the“unified generic theories”of evolutionary ecology.Models intended to describe the spatio-temporal patterns in ecological2and evolutionary processes must set the spatio-temporal scale unambiguously [18,19].For example,let us consider a single aerial photograph of a land-scape.While the boundaries of the photograph determine the spatial extent of the observation,the size of the pixels(grain size)in the photograph imposes the limit on the spatial resolution.Similarly,if a sequence of photographs of the same landscape is taken at regular intervals,the time difference between thefirst and the last photograph is a measure of the temporal extent of the observations while the time difference between the successive photographs determines the corresponding temporal resolution[20].For example,in ecol-ogy the temporal resolution can be days while the temporal extent can be up to decades,whereas the temporal resolution in evolutionary biology,par-ticularly empirical observations from fossil data varies,usually,from tens of thousands to millions of years.In this review we shall consider different classes of models with widely different scales of spatio-temporal resolution.The“ecological”models,that describe population dynamics in detail us-ing,for example,the Lotka-Volterra equations(discussed in section3)usually ignore the slow macro-evolutionary changes in the eco-system;hardly any ef-fects of these would be observable before the computer simulations would run out of computer time.On the other hand,in order to simulate the billion-year old history of life on earth with a computer,the elementary time steps in“evolutionary”models have to correspond to thousands of years,if not millions;consequently,thefiner details of the ecological processes over shorter periods of time cannot be accounted for by these models in any ex-plicit manner.However,despite the practical difficulties,it is desirable,at least in principle,to develop one single theoretical model which would be able to describe the entire dynamics of an eco-system since thefirst appearance of life in it up till now and in as much detail as possible.This dream has now come closer to reality,mainly because of the availability of fast com-puters.It has now become feasible to carry out computer simulations of eco-system models where,each time step(i.e.,temporal resolution)would correspond to typical times for“micro”-evolution while the total duration (i.e.,temporal duration)of each of the simulations is long enough to capture “macro”-evolution.The mathematical models in evolutionary ecology can be broadly classi-fied into different classes with different levels of detailed description,as shown in Fig.1.In the earliest mathematical models of population dynamics,only one predator species and one prey species were considered.However,for mod-3elling the population dynamics of more than two species,one needs to know the food web which is a graphical way of describing the prey-predator rela-tions,i.e.,which species eats which one and which compete among theselves for the same food resources[53,54,55,56].More precisely,a food web is a directed graph where each node is labelled by a species’name and each di-rected link indicates the direction offlow of nutrient(i.e.,from a prey to one of its predators).In the early works of this type,the food web was assumed to be static,i.e.,independent of time.An altogether different class of mod-els were developed to study macro-evolution;in such models,because of the Darwinian evolution,the food web is a dynamic network.In recent times, models of evolutionary ecology have been developed by a synthesis of eco-logical models of population dynamics and macroevolutionary models with evolving food webs.However,a more detailed theoretical description has also been attempted by incorporating individual organisms explicitly in the model where the birth,ageing and death of each individual occurs naturally.tems usually reach a stationary state after sufficiently long time.In contrast to these models,the models closer to biological reality might never reach a stationary state.After all,life forms in nature have evolved over billions of years from simple bacteria and archaea to complex structures like the bodies of dinosaurs and the brains of our readers.So far as the results of theoretical modeling are concerned,most of the recent works in physics literature have been concerned with the possibility of the existence of“power laws”in the statistics of extinction data.For example,suppose it is claimed that the relative frequency P(s)of extinctions of size s follows the power-law P(s)∝s−τwithτ≃2.If irrefutable evidences in favour of such power laws can be gathered,either from fossil data or from mathematical modeling,it would imply that the self-organizing dynamics of eco-systems exhibitfluctuations that are statistically self-similar because a change of scale(s→s′=bs)leaves the form of the power law unchanged. 2Ageing and age-structured population in single speciesBiological ageing of adults is best measured through the mortality q(a)= [S(a)−S(a+1)]/S(a)where S(a)is the number of survivors to the age a in suitable time units(like years for humans or days forflies and worms).More accurate is the mortality functionµ(a)=−d ln S(a)/da which is defined in terms of a derivative instead of a difference.The Gompertz law states that the mortality function for adults increases exponentially with age;this is also valid for many animals,but does not hold for the youngest and perhaps the oldest ages[21].Many theories exist to explain ageing[22,23,24]:accumu-lation of hereditary detrimental mutations[25],reliability theory as known from engineering in inanimate machines[26],loss of telomeres in cell division [27,28],damage caused by oxygen radicals,trade-offbetween longevity and fecundity(disposable soma),wear and tear,etc.Only for thefirst three, several quantitative computer simulations or mathematical solutions,giving the Gompertz law,are known to us,as listed in the cited literature.These models do not explicitly incorporate inter-species interactions like, for example,prey-predator interactions.These models cannot capture macro-evolutionary phenomena like,for example,extinctions which depend crucially on the prey-predator interactions.Such interactions are an essential part of5the ecological models which we mention briefly in the next section.3Ecological models of population dynamics: prey-predator interactionsTraditionally,the population dynamics of prey-predator systems have been described quantitatively in terms of the Lotka-Volterra equations[29,30,31,32].The nonlinearity of these deterministic differential equations leads toa rich variety of dynamical behaviour of the system.Although originally only two interacting species(predator and prey)were considered,later the approach was extended to more than two(but only a few)interacting species. Pioneering mathematical work of May[33]raised the question of stability of these dynamical equations when the number of interacting species increases. This challenged the earlier common belief that diversity of species ensured enhanced stability of the eco-system[34,35,36,37,38].Unlike the models of ageing discussed in the preceeding section,these Lotka-Volterra-type models of population dynamics monitor only the increase (or decrease)of populations caused by the birth(or death)of individual organisms but,usually,do not keep track of their ageing with time.The original formulation of the Lotka-Volterra equations assume that the population of the prey as well as that of predators are uniformly distributed in space.The absence of the spatial degrees of freedom in these equations is usually interpreted in the statistical physics literature as a mean-field-like approximation.This situation is similar to a well-stirred chemical reaction where the spatialfluctuations in the concentrations of the reactants and the products is negligibly small.On the other hand,spatial inhomogeneities in the eco-systems and migration of organisms from one eco-system to another are known to play crucial roles in evolutionary ecology[39,40,41].In recent years,the spatial inhomogeneity of the populations,i.e.,varia-tion in the population of the same spacies from one spatial patch to another [42],that have been observed in real ecosystems has been captured by ex-tending the Lotka-Volterra systems on discrete lattices where each of the lattice sites represents different spatial patches or habitats of the ecosystem [43,44,45,46,47,48,49,50,51,52].For modelling the population dynamics of more than two species,one needs to know the food web which is a graphical way of describing the prey-6predator relations,i.e.,which species eats which one and which compete among theselves for the same food resources[53,54,55,56].The structure of foodwebs and their statistical properties have been investigated both using field data on real eco-systems as well as abstract mathematical models[20, 57,58,59,60,61,62,63,64,65,66,67,68],particularly,in the recent years in the light of scale-free and small world networks[69,70].A static(time-independent)food web may be a good approximation over a short period of time.But,a more realistic description,valid over longer period of time,must take into account not only the adaptations of the species and their changing food habits,but also their extinction and creation of new species through speciation or migration of alien species into a new habitat. These processes make the food web a slowly evolving graph.Such slow time evolution of the food webs are naturally incorporated in macroevolutionary models which we summarize in the next section.4Modelling macroevolution and extinction: evolving food websThese models are intended to throw light on the mechanisms of origination of species through speciation as well as their extinction arising from biotic and abiotic causes.Several models have been developed just to account for the different routes to speciation[71,72,73,74,75,76,77,78,79].However, in this section we shall focus mainly on those works that have been inspired by close similarity with concepts or phenomena in statistical physics.4.1Self-organized critical models of eco-systems Inspired by the work of Per Bak and Kim Sneppen[80],a large number of evolutionary models have been developed over the last ten years[80,81,82, 83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98]Most of these works,including that of Bak and Sneppen,claimed the discovery of self-organized criticality[99]in the statistics of the numerical data on extinction. They also drew attention to the close relation of these observations with the concept of“punctuated equilibrium”introduced many years ago,by Gould and Eldredge[100,101]in the context of extinction of species[102].74.2Modelling evolution as a walk on afitness land-scapeAn alternative approach views macroevolution as random walk in a rugged “fitness landscape”[103].In recent years this approach has been extendedby allowing slow evolution of the landscape itself to incorporate the effect ofco-evolution of species[104,105,106].The notion of“fitness”has been used often loosely to mean different things[107].4.3Modelling eco-system as network of interacting speciesA network model of ecosystems was developed by Sole and Manrubia[88,89]. The system consists of N species,each labelled by an index i(i=1,2,...N). The state of the i-th species is represented by a two-state variable S i;S i=0or1depending on whether it is extinct or alive,respectively.The inter-species interactions are captured by the interaction matrix J;the elementJ ij denotes the influence of the species j on the species i.If J ij>0while, simultaneously,J ji<0then i is the predator and j is the prey.On the other hand,if both J ij and J ji are positive(negative)the two species cooperate (compete).The dynamics of the system consists in updating the states of the system (i.e.,to determine the state at the time step t+1from a complete knowledgeof the state at the time t)in the following three steps:Step(i):One of the input connections J ij for each species i is picked up randomly and assigned a new value drawn from the uniform distribution inthe interval[−1,1],irrespective of its previous magnitude and sign(this,we believe,is not a very realistic description of the inter-species interaction).Step(ii):The new state of each of the species is decided by the equationS i(t+1)=Θ N j=1J ij S j(t)−θi (1) whereθi is a threshold parameter for the species i andΘ(x)is the standard step function,i.e.,Θ(x)=1if x>0but zero otherwise.If S(t+1)becomes zero for m species,then an extinction of size m is said to have taken place.8Step(iii):All the niches left vacant by the extinct species are refilled by copies of one of the randomly selected non-extinct species.Sole and Manrubia[88,89]recorded extinctions of sizes as large as500 and the distributions of the sizes of these extinctions could befitted to a power law of the form N(m)∼m−αwith an exponentα≃2.3.Moreover, the periods of stasis t s were also found to obey a power law N(t s)∼t−γs with the exponentγ≃3.0.However,surprisingly,in none of their papers [88,89],did Sole and his collaborators report the distributions of the lifetimes of species which,according to some claims(see,for example,refs.[15,16] for references to the experimental literature and data analysis),also follows a power law.Amaral and Meyer[108]considered a hierarchical food web which was assumed to be organized into trophic levels;a species in levelℓfeeds on some species at the levelℓ−1(except for those atℓ=1which are au-totrophic).Origination of species through speciation was assumed to take place as follows:an empty niche is occupied by a non-extinct species at the same trophic level and the prey of the new species are selected randomly from among those at the level immediately below the trophic level of the new species.However,Amaral and Meyer did not treat the population dynamics of the species explicitly.Instead,a fraction p of the species at the lowest level is randomly selected and made extinct.Then,any species in the next higher level for which all prey species became extinct are also made extinct;this procedure is repeated for all the levels upto the highest one.Although this may be a more realistic description of inter-species interactions than that in the Sole-Manrubia model,the dynamics of the model is oversimplified.From the point of view of mathematical analysis,the advantage of this model is that its properties can be obtained not only numerically[108,109]but also analytically[110].The main limitation of these models i(see also[111])is that the individual organisms do not appear explicitly.On the other hand,it is the individual organisms,rather than species,which are the primary objects of selection [112,113,114,115].Moreover,the extinction of a species is nothing but eventual demise of all the individual organisms.Furthermore,direct experi-mental evidences[116,117,118,119,120,121,122]have established that sig-nificant evolutionary changes can occur over ecologically relevant time scales. In other words,the dynamics of Ecology and Evolution are inseparable.95”Unified”models of evolutionary ecology The need for“unification”of the various ecological subdisciplines,e.g.,popu-lation ecology,community ecology and evolutionary ecology,has been felt for quite some time[123].In the recent years attempts have been made to model evolutionary ecology in terms of“unified”models that describe both micro-and macro-evolution.Some of these models describe population dynamics in terms of one single dynamical variable and,therefore,fail to account for the age-structured populations of each species.However,only in the last two years it has been possible to develop detailed models that describe the birth,ageing and death of individual organisms.This is partly because of the availability of relatively fast computers.Abramson[124]considered a simple evolving ecosystem where each site of a one-dimensional lattice offinite length L represent a species such that the species i feeds on the species i−1and its eaten by the species i+1.The species1,which feeds at a constant rate on the environment,represents the species at the lowest level of the hierarchy whereas the species L,occupying the top of the chain,is not eaten by any other species.Such a linear food web is not realistic,but the importance of the model lies in its its simplicity.A species is considered extinct when its population falls below a preas-signed threshold.Because of the one-dimensional nature of the food web, the system would break into disjoint parts if any site is allowed to remain vacant following extinction of the corresponding species.In order to avoid such a situation,each niche that is left vacant by the extinction of a species is re-filled by another new species which interacts with the two neighbouring species with interactions whose strengths are drawn from a uniformly dis-tributed random fraction.However,Abramson did not monitor the ageing of each individual organisms.Instead,he monitored only the time evolution of the total populations of each species in the eco-system.McKane,Higgs and collaborators[125,126,127,128]also modelled the population dynamics in terms of one single dynamical variable.But,unlike, Abramson[124],they took into account the hierarchical organization of the species in food webs.In their webworld model,each species is represented by a set of L features(or,phenotypic characters)chosen from a set of K possible features.Evolution of the webworld model of eco-system is implemented by speciation events during which a new species is created from a randomly chosen existing species;the new species differs from the parent species by just one randomly chosen feature.There are some close similarities between10this model and some models developed in recent years incorporating the individual organisms explicitly in the model.Over the last two years,a few“unified”models of evolutionary ecology have been developed incorporating the individual organisms explicitly[129, 130,131,132,133,134,135,136,137].Each individual organism in the Rikvold-Zia model[135]has a genome of L genes,each of which can take one of two possible values,namely,0and 1.Thus the total number of different genotypes is2L.Rikvold and Zia[135] assumed that each of the different genotypes represent a separate species.A plausible justification,suggested by Rikvold and Zia,is that each binary “gene”actually represents a group of real genes in a coarse-grained sense.The spirit in which such“coarse-grained genes”are used in this model is somewhat similar to that of using the“phenotypic”characters in the webworld model [125,126,127,128].The number of individuals of genotype I in generation t is n I(t);the total population is N tot(t)= I n I(t).In each generation,the genomes of the individual organisms are subjected to random mutation with probabilityµ/L per gene per individual where L is the total size of the genome.Thus,by working with the genomes,Rikvold and Zia account for the genetic mutations explicitly and use the genotypes to label the different species.In order to keep the model as simple as possible,Rikvold and Zia assumed the successive generations to be nonoverlapping.More precisely,the organ-isms are incapable of living through successive reproduction cycles;an indi-vidual organism produces a litter of F offspring and immediately thereafter it dies.Consequently,this model does not describe age-structured populations of any species.Rikvold and Zia[135]considered a random food web where the effects of species j on the population of the species i is modelled by the element J ij of the interaction matrix J,exactly as in the Sole-Manrubia model[88, 89].In this model J is taken to be a time-independent random matrix, with vanishing diagonal elements,whose off-diagonal elements are selected randomly from a uniform distribution over the interval[−1,1].In each generation,the probability that an individual of genotype I pro-duces a litter of F offspring before it dies is P I({n J(t)})whereas the proba-bility that it dies without giving birth to offspring is1−P I.The reproduction11probability P I is assumed to be given by[135]P I({n J(t)})= 1+exp − J J IJ n J N0 −1(2) Here the second term in the exponential,commonly called the Verhulst factor, represents an environmental carrying capacity where N0is determined by the limited shared resorces like,for example,space,water,light,etc.The tangled nature model[136,137]is slightly more general than the model studied by Rikvold and Zia because overlapping generations are al-lowed.An individual organism is removed from the system with a constant probability p kill per time step.The algorithm used for reproduction proba-bility in the tangled nature model is also slightly different from that used by Rikvold and Zia.However,to our knowledge,there is no apriori justification at present for preferring either of these.In our works[129,130,131,132,134],we have provided the most detailed description of the ecological as well as the evolutionary processes.We have incorporated not only the hierarchical architecture of the natural food webs in a simplified manner but also the emergence of this architecture through self-organization as well as the possibility of migration of populations from one“patch”to another of the same eco-system for predation or merely for occupying a habitat.We have modelled the eco-sytem as a dynamic hierarchical network(see fig.2).Each node of the network represents a niche,rather than a species. Each niche can be occupied by at most one species at a time.The“micro”-evolution,i.e.,the birth,growth(ageing)and natural death of the individ-ual organisms,in our model is captured by the intra-node dynamics.The “macro“-evolution,e.g.,adaptive co-evolution of the species,is incorporated in the same model through a slower evolution of the network itself over longer time scales.Moreover,as the model eco-system evolves with time,extinction of species is indicated by vanishing of the corresponding population;thus, the number of species and the trophic levels in the model eco-system can fluctuate with time.Furthermore,the natural process of speciation is imple-mented by allowing re-occupation of the vacant nodes by mutated versions of non-extinct species.The prey-predator interaction between two species that occupy the nodes i and k at two adjacent trophic levels is represented by J ik;the three possible values of J ik are±1and0.The sign of J ik indicates the direction of trophic flow,i.e.from the lower to the higher level.J ik is+1if i is the predator and12n j(3)2where the superscript±on J ij indicates that the sum is restricted to only the positive(negative)elements J ij.The sum S+i is a measure of the total food13currently available to the i-th species whereas−S−i is a measure of the total population of the i-th species that would be,at the same time,consumed as food by its predators.If the food available is less than the requirement,then some organisms of the species i will die of starvation,even if none of them is killed by any predator.This way the matrix J can account for the shortfall in the food supply and the consequent competition among the organisms of the species i.The intra-species competition among the organisms of the same species for limited availability of resources,other than food,imposes an upper limit n max of the allowed population of each species;n max is a time-independent parameter in the model.Our model captures the starvation deaths and killing by the predators,in addition to the natural death due to ageing.self-organizing dynamics of the eco-system.1101001000100001000001e+061e+071e+081e+091101001000n u m b e rlifetime Figure 4:Distributions of the lifetimes of the species in the latest version of our “unified”model of evolutionary ecology.The size of the eco-system is 5×5in an appropriate dimensionless unit where each unit corresponds to a large patch.The lifetimes of the species are indicated along the X-axis whereas the number of times species with a given lifetime are encountered in our simulation,is plotted along the corresponding Y-axis.The symbols +,×,∗,open and filled squares correspond,respectively,to t =103,104,105,106and 107,where t is the time in dimensionless units (but,can be interpreted,for example,as one year)for which the ecosystem evolves following the dynamics of our model.In the original version of our “unified”model,we assumed that the pop-ulation of the prey as well as that of predators are uniformly distributed in space.This situation is similar to a well-stirred chemical reaction where the spatial fluctuations in the concentrations of the reactants and the products is negligibly small.On the other hand,spatial inhomogeneities in the eco-systems and migration of organisms from one spatial “patch”to another are known to play important roles in evolutionary ecology [39,40,41].We have15。