A bio-inspired computing model as a new tool for modeling ecosystems: The avian scavengers as a case study

A bio-inspired computing model as a new tool for modeling ecosystems:The avian scavengers as a case study M. Àngels Colomer a, Antoni Margalida b, Delfí Sanuy c, Mario J. Pérez-Jiménez d a Department of Mathematics, University of Lleida, Av. Alcalde Rovira Roure, 191. 25198 Lleida, Spain b Bearded Vulture Study and Protection Group, Apdo. 43, 25520 El Pont de Suert Lleida, Spain c Department of Animal Production, University of Lleida, Av. Alcalde Rovira Roure, 191. 25198 Lleida, Spain d Research Group on Natural Computing, Dpt. of Computer Science and Artificial Intelligence, University of Sevilla, Avda. Reina Mercedes s/n, 41012 Sevilla, Spain Keywords: P systems Ecosystem Avian scavengers Conservation abstract The models used for ecosystems modeling are generally based on differential equations. However, in recentyearsnewcomputationalmodelsbasedonbiologicalprocesses,orbioinspiredmodels,havearisen, amongwhichare Psystems.These areinspiredby thefunctionsof cellsandpresent importantadvantages with respect to traditional models, such as a high computational efficiency, modularity and their ability to work in parallel. They are simple, individual-based models that use biological parameters that can be obtained experimentally. In this work, we present the framework for a model based on P systems applied to the study of an ecosystem in which three avian scavengers (predators) interact with 10 wild and domestic ungulates (preys). The computation time for 100 repetitions, corresponding to 14 simulation years each, with an initial population composed of 385,422 individuals, was 30min. Our results suggest that the model presented, based on P systems, correctly simulates the population dynamics in the period of time analyzed. We discuss the usefulness of this tool in simulating complex ecosystems dynamics to aid managers, conservationists and policy-makers in making appropriate decisions for the improvement of management and conservation programs. 1. Introduction Mathematical models describing predator–prey relationships are used to study dynamics between two populations when one of them depends on the other for food and survival. The relationships among species, particularly vertebrates, in asymmetrical intraguild predator systems are complex, and generalizations remain elusive (Litvaitis and Villafuerte, 1996), but have important implications for conservation biology (Polis and Holt, 1992). Regarding the use of management and conservation measures for endangered species, it is particularly important to have an available model that allows us to reliably predict the dynamics of such populations, especially of the most endangered species (e.g., Meretsky et al., 2000; Ortega et al., 2009; Oro et al., 2008; Chapron et al., 2009; Grande et al., 2009). This tool can be useful in decision-making and in optimizing management of the given ecosystem. There are many processes in an ecosystem that run in parallel and are interrelated, hence it is essential to be able to organize these interrelations in a schematic and graphic way. Food webs are characterized by many weak interactions and a few strong interactions, which appear to promote community persistence and stability (Bascompte et al., 2005). In this sense, quantification of the strength of interactions between species is essential for understanding how ecological communities are organized and how they respond to human exploitation. Generally, the tools used for ecosystems modeling are based on differential equations. For example, the Lotka-Volterra model has been one of the most frequently used for the modeling of two species, predator and prey (Murray, 2002), and Verhulst’s model (logistic differential equations) has been utilized as a fundamental growth model in ecological studies because of its mathematical simplicity and single biological definition (Sakanoue, 2007, 2009). Fuzzy set theory has been used to estimate the parameters of the models based on differential equations studying the interaction between a prey and its predator (Da Silva et al., 2008). Despite the satisfactory results presented by these models, some authors question the deterministic and continuous approach they imply. Russell et al. (2009), for example, developed a stochastic model for three species by using ordinary logistic differential equations analyzed through numerical analyses techniques. Viability models do not consider optimal solutions but rather define all possible evolutions of a dynamic system under given constraints (Mullon et al., 2004). Viability theory considers that the evolution of a system is non-deterministic but belongs to a set of possibilities, depending on its state. Due to the mathematical complexity, the implementation requires a system of low dimensionality (<4 interacting components), which limits its effectiveness given the complexity of ecological models. Recently, some morecomplexmodelshavearisen(i.e.,takingintoaccountagreater number of parameters and variables), that are difficult to calibrate and validate. Moreover, their application requires division into smaller sizes defined arbitrarily and approximating the global model (e.g., Fulton et al., 2003, 2004). In this sense, Lawrie and Hearne (2008) propose a non-arbitrary algorithm for the division of these large models. Membrane Computing is an emergent branch of Natural Computing (P˘ aun, 1998; Ciobanu et al., 2006; P˘ aun et al., 2010) that was introduced with the purpose of defining computing devices, called P systems, which abstract from the structure and the function of the living cells. Rather than being an alternative to more classical modeling frameworks, such as ODEs (Ordinary Differential Equations), P systems constitute a complementary approach to be used when the classical modeling approaches fail. The most important property of these models is their capacity to work in parallel and to capture the randomness of the natural environmental processes using stochastic strategies based on Gillespie’s theory of stochastic kinetics (Gillespie, 1976, 1977) and the semantics defined by using probabilistic functions (Cardona et al., 2009, 2010b). In previous works we tested the utility of this new framework to managersand conservationists byapplying these modelson a community of scavengers (Cardona et al., 2009, 2010a) and the Zebra mussel Dreissena polymorpha (Cardona et al., 2010b). Here, taking a scavenger community that depends on carrion provided by wild and domestic ungulates as a model, we introduce a new model based on P systems generalizing the previous models and enabling thesimultaneousanalysisofsomeinterrelatedtrophicchains.First, interactionsamongspeciesareshown by means of Networks. P systems associate a rule to each interaction observed in the network quantifying the existent interaction. In order to evaluate its robustness, the model is checked and validated by using the experimental information from the ecosystem, which corresponds to a period of 14 years, taking three avian scavengers (as predator species) and ten ungulate species (as preys) as case studies. We discuss the results obtained from a methodological pointof-view and the usefulness of this tool in simulating complex ecosystemsdynamics to aidmanagers,conservationistsandpolicymakers in making appropriate decisions for the improvement of management and conservation programs. 2. Material and methods 2.1. Ecosystem to be modeled The study was carried out in the Pyrenean and Prepyrenean mountainsof Catalonia (NESpain,Fig. 1). Theecosystemto be modeledis composedof13species: threeavianscavengers(the Bearded vulture Gypaetus barbatus, the Egyptian vulture Neophron percnopterus and the Griffon vulture Gyps fulvus) as predator species, and six wild ungulates (the Pyrenean chamois Rupicapra pyrenaica, the Red deer Cervus elaphus, the Fallow deer Dama dama, the Roe deer Capreolus capreolus, the Wild boar Sus scrofa and the Mouflon Ovis orientalis) and four domestic ungulates that are found in an extensive or semi-extensive regime (the sheep Ovis aries, the goat Capra hircus, the cow Bos taurus and the horse Equus caballus) providing carrion for the avian scavengers and considered as prey species. Prey species are herbivores and their remains form the primaryfood resource forthe avian scavengersin the study area (>80% of the diet is based on these species, see Donázar, 1993; Margalida et al., 2009). Fig. 1. Study area and subpopulations considered in the ecosystem. Pyrenees (dark grey); Prepyrenees (light grey). The study area is characterized by the presence of two differentiated subpopulations that are interconnected. The Pyrenees region (57,244km2) is characterized by the presence of abundant wild ungulatesthroughouttheyear,astrongpresenceofdomesticungulates during the summer (as a consequence of transhumance) and low human density. In the Pyrenean region, the level of annual rainfall ranges from 800 to 1200mm and the maximum average temperatures do not exceed 25◦C in the summer and do not fall under−5◦C in thewinter. The orography presentsaltitudinal zones between 1000 and 3000m where alpine terrain, which is characterized by the presence of meadows above 2200m, dominates and subalpine ground is characterized by wooded formations of Mountain pine (Pinus uncinata). Under 1600m, montane terrain is found where wooded formations are dominated by European beech (Fagus sylvatica). The Prepyrenean region (77,372km2)is more populated by humans, with a more abundant domestic ungulate population that is regular during the winter, with low densities of wild ungulates. In the Prepyrenean region, the level of annual rainfall ranges from 800 to 1200mm and maximum average temperatures do not exceed 30◦C in the summer or fall under −2◦Cin the winter. The orography presents altitudinal zones between 600 and 1500m corresponding to submontane ground, which is dominated in its lower area by Portuguese oak (Quercus faginea) and Holm oak (Q. ilex sp. ballota) woods, and montane ground which is dominatedbyDownyoak(Q. pubescens). Occasionally, some mountain ranges are found in the montane region reaching 2000m. In regards to the distribution of the scavenger species, the bearded vulture population in the Pyrenees vs. Prepyrenees is 37.8% vs. 62.2% (n=37), the Egyptian vulture 10.2% vs. 89.8% (n= 59) and the Griffonvulture21.5% vs. 78.5%(n=822).Sinceindividuals can move from one area to another according to the resources available, the ecosystem would function as a single set (global ecosystem) composedoftwosubsets(subpopulationsseparatedbybiogeographical criteria). In this way, whenever there is a lack of trophic resources in one of the subareas, the individuals can move to the other one. In both, the ecosystem load capacity has been limited to the appropriate areas and habitats for the different species as well as the maximum density that can be reached (see Appendix A). The three scavenger species are cliff-nesting and only the Egyptian vulture is migratory (their presence in the study area is limited Elementary membrane Membranes Regions Skin Environment Environment Fig. 2. Representation of a membrane structure. to March–August). Concerning their trophic ecology, the diet of the Bearded vulture is based principally on bone remains of wild and domestic ungulates (principally sheep and Pyrenean chamois), although during the chick-rearing period small animals are important for the energetic requirements of the chicks (Margalida et al., 2005, 2009). The foraging areas are about 20km around the nest (A.M. unpubl. data). The Griffon vulture feeds mainly on wild and domestic ungulates and meat remains provided by sheep, pigs, cows and horses (Donázar, 1993). Grazing areas have a radius of approximately 25km around the hill. Finally, the Egyptian vulture is a more opportunistic species having a more heterogeneous diet which is based mainly on small corpses of mammals and birds. It is more dependent on rubbish dumps and supplementary feeding areas, and grazing areas are limited to radii under 10km from the nest (Donázar, 1993). Concerning prey species, domestic ungulates are present in the Pyrenean subpopulation during the summer (May–October) as a consequence of transhumant movements and become food resources for avian scavengers during this period (see Olea and Mateo-Tomás, 2009) whereas in the Prepyrenean region, they are present throughout the year in an extensive or semi-extensive regime. Wild ungulates are present in the area all year with the densities of Pyrenenan chamois, Red deer and Mouflon being more important in the Pyrenean region and Wild boar more important in the Prepyrenean region (see Appendix A). 2.2. A P system based model for ecosystems In this section, we will present the semantics and syntax used to model ecosystems by means of P systems and the specific model proposed for the ecosystem for the avian scavengers. 2.2.1. P systems Thestartingpointofthisnewmodelofcomputationistheobservation that the cell is the smallest living thing as well as a tiny machine with a complex structure, and the assumption that the processes taking place in the compartmental structure of a living cell can be interpreted as computations. The challenge is to take the cell itself as a support for computations, and to find those elements useful for computations in the structure and the functioning ofthe cell asa whole. Thedevicesof this modelare called Psystems, consisting of a cell-like membrane structure, in the compartments of which one places multisets of objects that evolve according to given rules. The main components of P systems are the membrane structure, multisets of objects, and evolution rules (Fig. 2). •Amembrane structure consists of several membranes arranged in a hierarchical structure inside a main membrane (the skin), anddelimitingregions(thespace in-betweenamembraneandthe immediatelyinnermembranes,ifany).Eachmembraneidentifies a region inside the system. •Regions defined by a membrane structure contain objects corresponding to chemical substances present in the compartments of a cell. The objects can be described by symbols or by strings of symbols, in such a way that multisets of objects are placed in regions of the membrane structure. •The objects can evolve according to given evolution rules, associated with the regions (and hence, with the membranes). The functioning of a P system is defined as follows: •Aconfigurationofa cell-like membrane systemconsistsofamembrane structure and a family of multisets of objects associated with each region of the structure. At the beginning, there is a configuration called the initial configuration of the system. •In each time unit, we can transform a given configuration to another configuration by applying the evolution rules to the objects placed inside the regions of the configurations, in a nondeterministic, maximally parallel way (the rules are chosen in a non-deterministic way, and in each region, all objects that can evolve must do so). In this way, we get transitions from one configuration of the system to the next. •Acomputation of the system is a sequence of configurations such that each one is obtained from the previous one by a transition, and shows how the system is evolving. The approach has a series of features that overcome several drawbacks of classical models based on differential equations: modularity (intrinsic to a membrane system), scalability/extensibility (further membranes and/or further evolution rules can be added without essentially changing the way a system works), understandability (evolution rules directly correspond to chemical reactions or interactions among species), programmability (a rewriting-based model can be easily transformed into a program, with certain programming languages, such as JAVA, C++, CLIPS), while preserving other desirable features of differential equations models, such as non-linearity of evolution. 2.2.2. The formal model In this section, we define a P system based framework where additional features, such as probabilistic functions and three electrical charges that better describe specific properties, are used. Askeleton of an extended P system with active membranes of degree q≥1, ˘=(,,R), can be viewed as a set of (polarized) membranes hierarchized by a structure of membranes (a rooted tree) labeled by 0, 1, ...,q−1. All membranes in are supposed to be (initially) neutral and they have associated with them R, a finite set of evolution rules of the form u[v]˛ i→u[v]ˇ ithat can modify their polarization but not their label. is an alphabet that represents the objects (i.e., Bearded vulture, Pyrenean chamois, etc., see Fig. 3). A probabilistic functional extended P system with active membranes of degree q≥1 taking Ttime units, ˘=(,,R,T,{fr:r∈R}, M0,...,Mq−1), can be viewed as a skeleton (,,R) with the membranes hierarchized by the structure labeled by 0, 1, ...,q−1. Tis anaturalnumberthat represents the simulation timeofthesystem. For each rule r∈Rand a,1≤a≤T,fr(a) is a whole number between 0 and 1, which represents a probabilistic constant associated with rule rat moment a. In a generic way, we denote r:u[v]˛ i fr(a) −→ u[v]˛ i. The tuple of multisets of objects present at any moment in the qregions of the system constitutes the configuration of the system at that moment. The tuple (M0,...,Mq−1) is the initial configuration of ˘. 0 1 2 3 45 6 a c c b a e b0 0 0 0 00 0Charge Membrane label Objects 0 213 456 Rules: ][[] ][][ ][ ][ 0 1 0 1 2 6 0 6 2 0 3 0 3 1 aar aar afder k→≡ →≡ →≡ + ... d Membranes hierarchized Fig. 3. A skeleton of an extend P system ˘(,,R) of degree 7 where =[[]1[[]4]2[[]5[]6]3]0and ={a,b,c,d,e,f}. The P system can pass from one configuration to another by using the rules from Ras follows: •A rule u[v]˛ i fr(a) −→ u[v]˛ iis applicable to a membrane labeled by i, and with ˛as electrical charge if multiset uis contained in the membrane immediately outside of membrane i,itistosay membrane father of membrane i, and multiset vis contained in the membrane labeled by ihaving ˛as electrical charge. When that rule is applied, multiset u(respectively v) in the father of membrane i(respectively in membrane i) is removed from that membrane, and multiset u(respectively v) is produced in that membrane, changing its electrical charge to ˛. •M() is the set formed by the multisets of .If u, v∈M(),i∈{0,...,q−1},˛∈{0,+,−} and r1,...,rz are the rules applicable whose left-hand side is u[v]˛ iat given moment a, then it should be verified that fr1(a)+...+frz(a)=1, and the rules will be applied according to the corresponding probabilities fr1(a),...,f rz(a). Amultienvironment probabilistic functional extended P system with active membranes of degree (m,q) taking Ttime units (˙, G, RE,,,R,T,{frj :r∈R˘,1≤j≤m},M ij : 0≤i≤q−1,1≤j≤m) can be viewed as a set of menvironments e1,...,emlinked by the arcs from the directed graph G. Each environment ejcontains a probabilistic functional extended P system with active membranes of degree q,˘j=(,,R,T,{frj :r∈R˘,1≤j≤m},Mij :0≤i≤q−1, 1≤j≤m) each of them with the same skeleton, ˘=(,,R), and such that M0,j,...,Mq−1,jdescribes their initial multisets. is an alphabet that represents the objects of that can be present in the different environments (Fig. 4). The communication rule between environments in REare of the form re:(x)ej px,j,k −→ (y)ek, and for each x∈,1≤j≤m,1≤a≤T, it verifies m  k=1 px,j,k(a)=1. When a rule of this type is applied the object xmoves from environment ejto environment ekconverted into y, according to the probability pj,k. We assume that a global clock exists, marking the time for the whole system (for its compartments), that is, all membranes and the application of all rules are synchronized. In the P systems, a configuration consists of multisets of objects present in the menvironments and at each of the regions of the P systems located in the environment. The P system can pass from one configuration to another by using the rules from R=RE∪∪ m j=1R˘jas follows: at each transition step, the rules to be applied are selected according to the probabilities assigned to them, and all applicable rules are simultaneously applied and all occurrences of the left-hand side of the rules are consumed, as usual. 2.2.3. The model Firstly, we graphically present the problem to be modeled by means of networks that are descriptors of ecological systems that can show the composition of numerous elements and the interactions among them (Bascompte, 2009). The network approach provides a powerful representation of the ecological interace1 e2 e3 e4 e1e2 e3 e4 Associate graph, G Fig. 4. Multienvironment probabilistic functional extended P system with active membranes of degree (4,7) (four environments and seven membranes). Fig. 5. Networks of the energetic needs and contributions: (a) Pyrenean subpopulation and (b) Prepyrenean subpopulation. Nodes represent scavengers’ energetic needs in the form of meat and the energetic contributions of each ungulate species. The nodes’ numeric labels correspond to ungulate species: (1) Pyrenean chamois, (2) Red deer, (3) Fallow deer, (4) Roe deer, (5) Mouflon, (6) Wild boar, (7) sheep, (8) cow, (9) goat and (10) horse. A link between nodes means that the ungulate forms part of the scavenger feeding. The thickness of the link denotes the energetic contributions and needs of all the species as a whole and is expressed in percentages. tions among species and highlights their global interdependence (Ulanowicz, 2004; Bascompte, 2009; Miehls et al., 2009). In this sense, the strength of the interaction among the three species of avian scavengers (predators) on the ungulate community (preys) is quantified for each subpopulation, and is measured as the biomass provided by prey species and the energetic requirements of these three avian scavenger species that are expressed as megacalories per year according to the population size. We consider the vegetable biomass on which wild and domestic ungulates depend as not being a limiting factor but exceeding the needs (energetic requirements) of these species (García, 2008). Subsequently, the model presented associates a rule quantifying the interaction to each of the relations shown in the networks (Fig. 5). In the networks, the energetic needs in the form of meat of the three species of avian scavengers studied inhabiting the Pyrenees (Fig. 5a) and the Prepyrenees (Fig. 5b), are shown by means of nodes as is the proportion of biomass made up by wild and domestic ungulates. Themodelproposedmustconsider:(a)thepopulationdynamics of the nine wild species (three avian scavengers and six ungulates) and four domestic species, (b) the interactions among the 13 species,(c)the presence of two zones inthestudyarea,(d) the communication protocol between the two areas and (e) the ecosystem maximum load capacity for each of the areas. In order to model the ecosystem, we use a multienvironment probabilistic functional extended P system with active membranes of degree (2,2) (two membranes and two environments) taking Ttime units (simulation years). (˙, G, RE,,,R,T{frj :r∈R˘,1≤j≤2},M ij : 0≤i≤1,1≤j≤2) The skeleton consists of the working alphabet, , formed by all objects that belong to the initial configuration and the objects that appear in the evolution of the P system (all of them appear in Appendix A). The membrane structure is formed by the skin membrane labeled 0 and an in membrane labeled 1, both having neutral charge. The set of rules are shown in the appendix; there are 49 types of rules. The probabilistic functional extended P system ˘=(,,R,T, {fr:r∈R˘},M0,M1) is defined as follows: for each r,fris a constant function, the initial configuration is M0=Xqi,j ,d iand M1={R0, F0}. The objects Xi,j,1 are associated with one animal belonging to the species i,jyears old in the instance 1and qi,jis the amount of objects Xi,j,1, we use the object difor a control of the maximum load of the animals of species i.F0is used for generating external contributions of different kinds of food, and finally, object R0allows us to synchronize the P system. The alphabet that can be present in the environment is ˙= {Zi,j,s,Z i,j,s}and the relationships in graph G are of each environment to itself and to other environments. The algorithmic scheme of the model is structured following a series of modules which are run sequentially corresponding to the passing of 1 year in the ecosystem (Fig. 6). Except for reproduction, all the processes are continuous and annual. In the model, processes are discretized in Cardona et al. (2009), which verified that the order in which modules are applied does not affect the final result obtained from the model. Some ungulates (Red deer and all domestic ungulates) of the ecosystem have been classified into two different groups according to their management. Specifically, in the case of the Red deer, sexes have been separated because mortality in males is higher than in females due to hunting activities. Thus, the total number of groups of animals considered is 18. The spatio-temporal distribution of domestic animals has been determined bearing in mind that some animals spend the entire year in the study area whereas others spend variable periods of time there (e.g., summer transhumance). Regarding transhumance, sheep, cows and horses are in mountain passes for 6 months of the year, making use of summer grazing (from mid May to mid October depending on the climate; Roigé, 1995). Data about transhumance were taken from the information provided by the Departament d’Agricultura i Ramaderia of the autonomous government of Catalonia. The model used is based on objects that evolve by means of a series of rules among which those associated with the interactions shown in the Networks are found (see Fig. 5). Objects Xi,j,y,Y i,j,y,Z i,j,y,Z i,j,y and Wi,j,yrepresent the same animal at its different stages throughout the module series. Within the objects, the first index codifies the group to which the animal belongs, the second one codifies its age and the third one the simulation year. The number of years (T) to be simulated is an input of the model. 2.2.4. Reproduction module At the initial instance, an object of type Xis associated with each animal. When rules from the reproduction module are applied to objects X, they evolve to objects of type Y. Objects Xi,j,y associated with females that reproduce when they reach fertility also generate objects Yi,0,yassociated to newborn animals. In this module, objects associated with the amount of food produced by the ecosystem itself (grazing or external contributions of No Yes (no change environment) REPRODUCTION MORTALITY FEEDING + DENSITY REGULATION (1) CHANGE ENVIRONMENT UPDATING FEEDING + DENSITY REGULATION (2) 1 step 8 steps (synchronizing) 5 steps 1 step 1 step 1 step 2 steps Fig. 6. Modules that form the model. meat and bones by man) are also generated. This module takes one simulation step. Breeding parameters of each species were obtained from the literature and unpublished data (see Donázar, 1993; Grande, 2006; Margalida et al., 2003; Oro et al., 2008; Le Gouar et al., 2008, see Table 1,Appendix A). 2.2.5. Mortality module Hunting, in the case of wild animals, and mortality due to both natural causes and human management, in the case of domestic animals, provide the animal biomass on which avian scavengers feed.Theirsurvivalhasbeenestimated by using bibliographical references (Montserrat and Villar, 2007; Blasco et al., 1992; Casasús et al., 1999; Margalida et al., 2009) as well as by means of personal surveys(authorsunpubl.data).Moreover,apercentageofdeadanimals that may be accessible to scavengers was estimated. The input of this module is formed by objects of the type Yi,j,y, which do not evolve but rather move to another membrane and generate a new object D for each animal that does not die. Objects associated with animals that die evolve to objects associated with meat (C,M) and bones (H,B). This module takes one simulation step. 2.2.6. Feeding and density regulation module (1) Annual energetic requirements (expressed as calories or megacalories) as well as the maximum load capacity in the area under study have been estimated for all the species. In this sense, there is no overexploitation of mountain grazing and cattle raising loads have decreased in recent years such that the vegetable biomass available is not a limiting factor (García, 2008). Regarding the avian scavenger community, its population has increased in the last 20 years and maximum values have been estimated according to the unoccupied space as well as to the maximum density that each of these species can achieve (see Table 1,Appendix A). Whether or not the maximum load capacity of the ecosystem for eachspecieshasbeenreachedisdeterminedbyusingobjectsDpreviously generated for each surviving animal. Furthermore, objects Yevolve to objects Zto begin the feeding process. In the second step of this module, objects Zevolve to objects Wif there is enough physical space and food. 2.2.7. Change environment module When one of the subareas reaches its maximum load capacity, it has been considered that any of the species modeled can move to another subarea according to its ecological requirements. This module will apply if there is some object Zthat has not evolved in the previous step; that is, if the resources have been insufficient for all of the animals. In this case, the model simulates animals’ movements to find the necessary resources for their survival. Objects Zgo out to the environment by taking two simulation steps and subsequently move to another environment by evolving to objects Z. Next, in two more steps, they enter the inner membrane of the P system where the objects associated with the possible unused resources are found. This module takes five simulation steps. 2.2.8. Feeding and density regulation module (2) This module will apply if there are enough resources for animals coming from other areas. In such a case, objects Zevolve to objects of type Win one simulation step. 2.2.9. Updating module This module will apply after eight simulation steps; that is, after the running of the third module. After running one cycle within the moduleseries, the initial configuration must be re-established such that a new year (cycle) begins. Objects associated with the surviving animals, Wi,j,y, evolve to objects Xi,j+1,y+1. The remaining food is removed and the objects associated with those animals which have not survived evolve to objects representing food. This module takes one simulation step. In summary, the running of one cycle in the module succession takes 11 simulation steps and represents the passing of a 1-year periodin the ecosystem.Afterthe running ofthe cycle, the Psystem returns the number of living animals of each species as well as the resources in the form of meat that each has contributed. 3. Results For execution of the model, MeCoSim software (free software under licence), developed by members of the Natural Computation Group at the University of Sevilla (GNU GPL; http://www.plingua.org), has been used. The population trend of the three scavenger species and the six wildungulatesobtainedbythemodelwithrespecttodataobtained by direct censuses from 1994 to 2008 is shown in Figs. 7 and 8. In order to check and validate the model, only initial (1994) and final (2008) data were available for the six wild ungulates whereas inter-annual censuses were available for avian scavengers. The population trend of the species present in the ecosystem throughout the period under study has been obtained by running the simulator 100 times for 14 years with the same input data. The simulator executions have allowed us to estimate the population confidence intervals of the different species. The ecosystem mod- M.À. Colomer et al. / Ecological Modelling 222 (2011) 33–47 39 Table 1 Values of parameters used in the model for each species (F=female, M= male, A= spend the entire year in the mountain, P= spend part of the year in the mountain). g1g2g3g4g5g6g7k1k2k3m1m2m3m4f1f2f3f4f5f6f7f8f9 Gypaetus barbatus 1 1 1 6 20 21 0 0.65 0.35 1 0.06 0.08 0 1 00001350450 0 Neophron percnopterus 1 0.5 1 5 24 25 1 0.80 0.57 1 0.28 0.08 0 1 0000001000 0 Gyps fulvus 1 1 1 5 24 25 0 0.75 0.56 1 0.06 0.07 0 1 0000002300 0 Rupicapra pyrenaica 1 1 1 2 18 18 0 0.55 0.75 1 0.6 0.06 0 1 346240000.50.5 Cervus elaphus (Female) 1 1 1 2 17 17 0 1 0.75 1 0.34 0.06 0 1 7 13 15 60 0 0 0 0.6 0.6 Cervus elaphus (Male) 1 1 1 2 20 20 0 0 0 0 0.34 0.06 0 1 12 15 24 96 0 0 0 0.6 0.6 Dama dama 1 1 1 2 12 12 0 0.75 0.55 1 0.5 0.06 0 1 1 14 2 37 0 0 0 0.25 0.25 Capreolus capreolus 1 1 1 1 10 10 0 0.67 1 1 0.58 0.06 0 1 141190000.25 0.25 Ovis orientalis 1 1 1 2 12 12 0 0.5 0.9 2 0.6 0.06 0 1 346220000.60.6 Sus scrofa 11114600.50.55 4 0.14 0.1 0 1 4 6 12 60 0 0 0 0.25 0.25 Ovis aries (Adult) 01128800.96 0.75 1 0.15 0.03 0 0 347380000.70.7 Ovis aries (Young) 00.5128800.96 0.75 1 0.15 0.03 0 0 347380000.70.7 Bos taurus (Adult) 012291400.90.910.057 0.045 0 0 10 60 6 518 0 0 0 0.6 0.6 Bos taurus (Young) 00.42291400.90.910.057 0.045 0 0 10 60 6 518 0 0 0 0.6 0.6 Capra hircus (Adult) 01128800.97 0.9 1 0.12 0.015 0 0 349370000.60.6 Capra hircus (Young) 00.5128800.97 0.9 1 0.12 0.015 0 0 349370000.60.6 Equus caballus (Adult) 013392000.97 0.9 1 0.034 0.0142 0 0 10 60 9 891 0 0 0 0.8 0.8 Equus caballus (Young) 0 0.55 3392000.97 0.9 1 0.034 0.0142 0 0 10 60 9 891 0 0 0 0.8 0.8 g1: 1 wild animal and 0 domestic animals. g2: proportion of time that animals remain in the mountains during the year. g3: age at which adult size is reached. This is the age at which the animal consumes an adult diet, and at which if the animal dies, the amount of biomass it leaves is similar to the total left by an adult. Moreover, at this age it will have surpassed the critical early phase during which the mortality rate is high. g4: age at which fertility begins. g5: age at which fertility ends. g6: average life expectancy in the ecosystem. g7: 1 if an important proportion of the diet of the species can be based on other small species and 0 for the remainder. k1: proportion of females in the population (per one). k2: fertility ratio (proportion of fertile females that reproduce). k3: number of descendants for fertile females that reproduce. m1: natural mortality ratio in first years, age <g4(per one). m2: mortality ratio in adult animals, age ≥g4(per one). m3: percentage of domestic animals removed from non-stabilized populations at early ages. m4: is equal to 1 if the animal dies at the age of g6and is not removed, and is equal to 0 if the animal does not die at the age of g6but is removed from the ecosystem. f1: amount of bones from young animals, age <g4. f2: amount of meat from young animals, age <g4. f3: amount of bones from adult animals, age <g4. f4: amount of meat from adult animals, age <g4. f5: amount of bones necessary per year and animal (kg). f6: amount of grass necessary per year and animal (kg). f7: amount of meat necessary per year and animal (kg). f8: Percentage of useful bones. f9: Percentage of useful meat. eled in this work is made up of 13 species totaling 18 animal types. Themodelshowsthepopulationtrendobtainedby equallydividing the evolutions of each animal. The initial population is composed of 385.422 individuals and the computation time (on a personal computer) for 100 repetitions, corresponding to 14 simulation years each, was 30min. The comparison between the real population tendency, estimated by means of censuses carried out, and that obtained by the simulatorisshown in Fig. 7. Theadjustmentofthepopulation trend shownforthethreeavianscavengerspecieswithrespecttothedata obtainedbythesimulatorshowsthatthemodelfunctionsproperly. With respect to the results obtained for wild ungulates (Fig. 8), it is observed that in simulation year 10, Roe deer reaches its maximum load capacity in the Prepyrenean area (zone 2) and some of the animals move to the Pyrenean area, causing in the latter an important population increase in simulation year 11. Thus, the model is able to show the dispersive capacity of some species. A similar situation is observed for Red deer in simulation year 10 and Pyrenean chamois in simulation year 9. Note that as simulations continue beyond the initial year of simulation the confidence interval increases. It is a random model such that results are spread horizontally. 4. Discussion The model presented based on P Systems correctly simulates the population dynamics in the period of time analyzed. Our model considers the population dynamics and the simultaneous interaction among the 18 animal types. In Cardona et al. (2009) we described a model based on P systems (the ecosystem modeled is composed by one avian scavenger species and five prey species) that considers neither the ecosystem maximum load capacity nor the appearance of density-dependent regulatory phenomena in the species. In Cardona et al. (2010a), a new model was presented that overcame some limitations of the previous model by widening the number of species (three avian scavenger species and ten prey species) including some specific characteristics of the speciesformingtheecosystem. Nonetheless, this model considered the ecosystem to be a closed entity, that is, when the necessary resources for a species to survive (space, feeding,...) are insufficient, the animal dies without considering the possibility that the animals move to some other ecosystem. In this study, we improve the model of Cardona et al. (2010a) considering the heterogeneity of the landscape and the possibility of spatial movements of the species when food resources are enough to cover the energetic requirements of the species. Models based on individuals such as that presented in this study are generally more flexible and enable the consideration of the heterogeneity of the population and the environment. Our model is composed of modules that are applied sequentially. In our model the advantages with respect to other models are (1) it is not necessary to divide the problem to be analyzed (Fulton et al., 2003, 2004; LawrieandHearne,2008);(2)thenumberofinter-andintraspecific interactions among individuals is not limited (Mullon et al., 2008); (3) it is able to capture the randomness inherent to the processes; and (4) in the function of ecosystem dynamics, it is able to update the parameters in a functional way, or in other words, to readjust (Cardona et al., 2010a). Another advantage in the application of our model is that it is easily programmable and is fast compared to the time of computationfor a deterministicmodelwithsimilarcharacteristics involving 9h in one simulation year with 240,000 individuals (Morales et al., 2006). Amongthespecializedprogramsthatenablethestudyofspecies viability (e.g., Gapps, Inmat, Ramas, Vortex) the most widely used Gypaetus barbatus 0 5 10 15 20 25 30 35 40 14131211109876543210 Year Territories Pyrenees Prepyrenees Total Experimental data Neophron percnopterus 0 10 20 30 40 50 60 70 14131211109876543210 Year Territories Pyrenees Prepyrenees Total Experimental data Gyps fulvus 0 200 400 600 800 1000 1200 1400 1514131211109876543210 Year Pairs (n) Pyrenees Prepyrenees Total Experimental data Fig. 7. Comparison of the results obtained by the simulator (lines) with the data obtained experimentally (dots). Unbroken line represents the whole population whereas broken lines represent the results obtained in both subpopulations, Pyrenees: circles; Prepyrenees: squares. is Vortex, which considers a higher number of factors. A comparative study showed that the results obtained from the diverse programs were different if the models were not normalized (Brook et al., 1999). Like the model presented in this paper, Vortex enables the modeling of different populations and considers fertility ratios, male-female relationships, number of descendants and interactions among different populations. Nonetheless, Vortex considers the possibility that natural disasters may take place whereas the model presented here can only consider such a possibility by including a new module. Finally, Vortex can simultaneously model different populations that interact with each other although they must correspond to the same species. The model used in this work enables the study of the dynamics of different populations whether or not they are composed of the same species, even allowing interaction and competition among them. In addition, the availability of energetic resources that are essential for the dynamics of populations that are interrelated as well as spatio-temporal regulation are considered. The flexibility of our model enables the increase in the number of species without having to modify the model but simply by adding the new information on the biological parameters of the new species to be included to the simulator or modifying the infor- Fig. 8. Results provided by the simulator. The starting point is the population of every species in the year 1994. Zone 1: Pyrenees, Zone 2: Prepyrenees.