A Review of Membrane Computing Models for Complex Ecosystems and a Case Study on a Complex Giant Panda System

Research Article A Review of Membrane Computing Models for Complex Ecosystems and a Case Study on a Complex Giant Panda System Yingying Duan, 1 Haina Rong , 1 Dunwu Qi, 2 Luis Valencia-Cabrera , 3 Gexiang Zhang , 1 , 4 and Mario J. P´ erez-Jim´ enez 3 1 School of Electrical Engineering, Southwest Jiaotong University, Chengdu 611756, Sichuan, China 2 Chengdu Research Base of Giant Panda Breeding, Chengdu 6110081, Sichuan, China 3 Research Group on Natural Computing, Department of Computer Science and Artificial Intelligence, Universidad de Sevilla, Avda. Reina Mercedes s/n, 41012 Sevilla, Spain 4 Research Center for Artificial Intelligence, Chengdu University of Technology, Chengdu 610059, Sichuan, China Correspondence should be addressed to Haina Rong; [email protected] and Gexiang Zhang; [email protected] Received 1 March 2020; Revised 8 August 2020; Accepted 18 August 2020; Published 24 September 2020 Academic Editor: Carlos Gershenson Copyright ©2020 Yingying Duan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Ecosystem modelling based on membrane computing is emerging as a powerful way to study the dynamics of (real) ecological populations. These models, providing distributed parallel devices, have shown a great potential to imitate the rich features observed in the behaviour of species and their interactions and key elements to understand and model ecosystems. Compared with differential equations, membrane computing models, also known as P systems, can model more complex biological phenomena due to their modularity and their ability to enclose the evolution of different environments and simulate, in parallel, different interrelated processes. In this paper, a comprehensive survey of membrane computing models for ecosystems is given, taking a giant panda ecosystem as an example to assess the model performance. This work aims at modelling a number of species using P systems with different membrane structure types to predict the number of individuals depending on parameters such as reproductive rate, mortality rate, and involving processes as rescue or release. Firstly, the computing models are introduced conceptually, describing the main elements constituting the syntax of these systems and explaining the semantics of the rules involved. Next, various modelled species (including endangered animals, plants, and bacteria) are summarized, and some computer tools are presented. Then, a discussion follows on the use of P systems for ecosystem modelling. Finally, a case study on giant pandas in Chengdu Base is analysed, concluding that the study in this field by using PDP systems can provide a valuable tool to deepen into the knowledge about the evolution of the population. This could ultimately help in the decision-making processes of the managers of the ecosystem to increase the species diversity and modify the adaptability. Besides, the impacts of natural disasters on the population dynamics of the species should also be considered. The analysis performed throughout the paper has taken into consideration this fact in order to increase the reliability of the prospects making use of the models designed. 1. Introduction Membrane computing is a fast-growing branch of natural computing [1]. The computational devices within this paradigm are called membrane systems or P systems, and they have attracted major interest since their appearance. Nowadays, membrane computing community is actively combining deep theoretical research studies with practical applications. On the one hand, theoretical studies focus on how to design membrane computing models according to the compartmentalized structure and the functioning of biological membranes within living cells and how to evaluate their computational power and computational complexity [2], addressing efficiency aspects. Different types of membrane structures abstracted from biological cells can be distinguished in membrane systems. The most widely studied are cell-like P systems [3–6] (inspired in the Hindawi Complexity Volume 2020, Article ID 1312824, 26 pages https://doi.org/10.1155/2020/1312824 compartmentalized hierarchical structure inside a cell), tissue-like P systems [7, 8] (with the focus on the interconnection among cells, not entering into details of the internal compartments inside each cell), and spiking neural P systems [9–12] (inspired from the transmission of electrical pulses, also known as spikes, among neurons). Along the last twenty years, plenty of research results have proved that a number of variants of these computational models are equivalent in computing power to Turing machines (computationally complete), and many of them obtained efficient solutions (polynomial solutions, even linear in some cases) to a variety of computationally hard, NP-complete [13–15], or PSPACE problems [16–18]. On the other hand, applied research of membrane computing models aims to effectively apply the introduced computational models to handle several practical problems, from basic ones to the modelling of complex systems. For automatic design of membrane computing models (ADMCM), regarded as an automatic computation device, the models can adaptively achieve basic arithmetic operations [19]. Thus, in [20], Huang et al. applied P systems with the Q-bit representation to compute power n 2 , something also achieved in [21] by Ou et al., who applied P systems to calculate power n 2 with a different approach using an elitist genetic algorithm (GA). Within this same research line, in [22, 23], five types of automatic P systems were used to compute addition, subtraction, product, division, and power. In some studies, spiking-like P systems are used as a computing device to solve several arithmetic problems; for instance, the addition of nnatural numbers and the product of two arbitrary natural numbers with a given length of binary bits [24, 25]. Of course, apart from providing automatic design and arithmetic operations, P systems are applied to hard problems such as vertex cover [26, 27], quadratic assignment [28], graph coloring problems [29–31], and non-semilinear sets [32, 33]. Besides, they have been applied to solve image processing problems [34–39], complex optimization problems [40–42], complex market interactions [43], and intelligent control problems for robots [44–46]. Many applications have ultimately been implemented by using simulation tools. For the state of the art of such simulators and environments for the virtual experimentation with P system-based models, refer to [47] and the associated website (https://www.gcn.us.es/SimulationMC). The applications summarized above plus some others show that membrane computing models are useful tools to solve many practical problems of very different nature. It seems clear that each type of computationally hard problem addressed by membrane computing models somehow implies an extension of application fields, providing theoretical and practical foundations, opening paths that can be worth exploring, and widening the application space. Regarding the contributions of membrane computing to model complex systems, relevant achievements have been made along the last decade, with a special attention to the study of real ecosystems (focusing on endangered or invasive species, among others) and population dynamics in general [48]. Thus, endangered species in ecosystems present reproduction rates usually low, along with mortality rates abnormally high, due to reasons derived from the biology of the species, the increased threats by other species, and the effects of human activities or natural disasters. It is often the case that these endangered species are not considered to be free of danger even when the threat is vanishing because of their scarcity in the undisturbed fragments so that isolated populations sometimes cannot survive after destruction and might become extinct. Hence, the qualitative and quantitative understanding of the inherent laws or processes underlying the disturbances in the population size and distribution (i.e., the population dynamics) has become critical for the successful management and conservation of endangered species [49]. As briefly mentioned, apart from their natural mortality, most species have suffered the effects of environmental disasters. As certain studies point out, some types of such natural disasters have caused or may potentially cause the risk of species biodiversity collapse, or even some species extinction in certain regions [50–53]. For example, in [54–58], the authors studied the impacts of climate change on parasite biodiversity [54], the divergence of species responses [55], information fusion of numerous natural or environmental factors by using a combination of mathematical models [56], the parameters related with the population dynamics of songbirds [57], and the impact for agricultural welfare [58]. In order to accurately assess the impacts of these factors, several mathematical models are used to analyse the effect of these elements on population dynamics, e.g., differential equations [59], generalized linear models (GLMs) [60], generalized additive models (GAMs) [61], ecological niche factor analysis (ENFA) [62], or machine-learning methods such as Bayesian approaches [63], population viability analysis [64, 65], agent-based models [66], maximum entropy method (Maxent) [67], and neural networks [68]. According to the predicted results of these models, the parameters related can affect the population change in varying degrees. It is highlighted by the authors the fact that it would be worth providing projections for endangered species population dynamics under the influence of potential natural disasters in order to protect them. Besides, it would be advisable to address the research on the population dynamics of the species following different approaches; there is no single Swiss knife resulting in the most adequate tool to apply in all situations. This work focuses on membrane computing (MC) models of species in certain ecosystems, thus involving mostly two main fields: membrane computing and population dynamics. Concerning the former one, MC has among its strengths the availability of rigorous and complete, strongly founded theoretical-practical developments; in addition, it provides parallel distributed devices in a framework with flexible evolution rules. With respect to the latter one, the population dynamics of species obviously involves several biological processes, including some common ones (such as feeding, reproduction, and mortality), and some more specific of certain studies, such as rescue, release, and biochemical reactions of bacteria; besides, the dynamics of the species could also be affected by the potential impacts of natural disasters on the populations under study. 2Complexity Based on the characteristics of membrane systems and the features of endangered species, the evolutionary behaviour of such species can be expressed by the rules of membrane systems. Hence, along this paper, we recapitulate ecosystem models using different types of P systems. Firstly, we analyse between the elements related with the species and their interactions in the ecosystem and those related with the definition of P systems. So far, there are several references to study the applications of P systems on different real ecosystems. The common characteristics can be summarized as follows: (a) each individual is represented by an object in the P system; (b) different behaviours or processes affecting the individuals of the species are abstracted as the rules in the P systems; and (c) a certain living environment in the ecosystem is abstracted as a membrane structure. Regarding the dynamics of these systems, at every moment, all individuals will evolve synchronously. As these species inhabit different geographical regions but are subject to the same processes (with the values of certain parameters possibly changing), this scenario can be represented by using multienvironment P systems, where communication of individuals is possible among different environments. For some species, such as plants and bacteria, the evolutionary processes are also modelled according to their characteristics and possibly distinguishing environments with different parameters. Through the analysis above, we have depicted some of the most relevant facts taken into account when modelling ecosystems based on P systems in order to accurately predict data about the biological evolution of the species, aiming to capture in these models (i.e., to mimic) the relevant elements of the real biological phenomena under study. The rest of this paper is arranged as follows. Section 2 introduces membrane computing models for ecosystems. After that, Section 3 summarizes the applications of such models to the population dynamics of certain ecosystems. Then, Section 4 lists several simulation software tools to perform virtual experiments for P system models of ecosystems. Later on, in Section 5, a case study on the population dynamics of giant pandas is analysed. Finally, some conclusions and possible further developments are discussed in Section 6. 2. Membrane Computing Models for Ecosystem Modelling As outlined in Section 1, different types of mathematical models have been applied to ecosystems. These models are representations imitating the real systems under study, using a certain formalism. In particular, some of these approximations are computational models, which means they follow the rules of some computing paradigms that regulate their behaviour and can be computed directly in their computational devices or be simulated following the same exact rules; on the contrary, noncomputational models (e.g., differential equations) require the use of approximated methods in order to be computed by some computing devices. When membrane computing is used to create a representation of an ecosystem, incorporating their main parameters, individuals, processes, etc., involved in their dynamics, this is considered a computational model because it is a model of the ecosystem that is based on a computational paradigm (in this case, membrane computing), whose computation follows the exact rules of the formal model, not requiring any approximate method to be computed. These models of ecosystems are based on membrane computing, so they are usually called membrane computingbased models. Such models are represented by computational devices called membrane systems (commonly referred to as P systems) (commonly referred to as P systems). As computational devices, membrane systems or P systems are abstracted from the structure and functioning of living cells. There are three main classes of P systems: (1) celllike P systems, inspired from living cells; (2) tissue-like P systems, inspired from the interactions of cells in tissues; and (3) SN P systems (spiking neural P systems), inspired from the communication of electrical pulses (also known as spikes) in biological neural networks. Concerning the population dynamics of ecosystems, both cell-like and tissue-like P systems have been useful to inspire the appearance of a new class of P system, originally called the multienvironment P system. Regarding the membrane structures employed in the models built with these devices, they define a structure combining the hierarchical structure of cell-like P systems with an upper layer of the so-called environments, each one containing a single cell-like system and introducing the possibility of communication among environments, similar to that in tissuelike P systems. This structure, as we can observe, is richer than the two previous models, thanks to the combination of the two constituent parts: the internal hierarchical structure present in cell-like P systems (which is not present in tissuelike P systems, where only cells without the internal structure are present); the existence of different regions with cells inside and allowing communication among environments as a graph (which is not possible in cell-like systems). Of course, with this new class given by the multienvironment P systems, we can also model ecosystems with a single environment involved, defining such environment as the only node in the graph, with cell-like P systems inside, and no other nodes to communicate with through edges of the graph. In contrast, there are more interesting problems that require additional environments. In such cases, a more complex graph must be present to allow the movement of certain elements among environments. In these scenarios, such communication among environments will play a crucial role. In addition, each environment will contain its own hierarchical structure given by the cell-like P system it holds. With respect to the dynamics of the systems involved in these computational models, mostly two main paths have been followed when dealing with multienvironment systems: a stochastic approach followed by the so-called multicompartmental systems (there are zero or several celllike P systems inside each environment, which are able to communicate with entire P systems from one environment to another), with rules subject to chemical laws and kinetic constants; a probabilistic approach (there is one and only one cell-like P system in each environment and only single object communicates among environments), with rules Complexity 3 subject with probabilities, whose values are typically based on evidence; for further details, refer to [19]. Nowadays, the stochastic approach is mostly associated with computational models at a microlevel (e.g., involving molecular interactions), not being widely used to model ecosystems at a macrolevel. Consequently, in the following, we only consider the computational models of P systems following the probabilistic approach, termed as population dynamics P systems (PDP systems) [69]. As just mentioned, PDP systems are a variant of P systems introducing probability mechanisms into membrane systems. The framework constituted by PDP systems was conceived for multiple environments, and therefore belongs to the class of multienvironment P systems, but can obviously support the particular case of having the number of environments which equals to 1 and consequently not having communication among environments. Now, while the framework is uniform, from a practical point of view, it is easy to distinguish systems depending on the number of environments min the system; thus, if m�1, we may call them as single-environment PDP systems (with a cell-like P system inside a single environment, see Figure 1(a)); for m>1, we would talk strictly about multienvironment PDP systems (with several environments, each one containing a single P system inside, see Figure 1(b)). These two types of PDP systems are introduced as follows, including syntactic and semantic aspects. Definition 1 (see [69]). A single-environment P system of degree q, with q≥1, is a tuple: π�Γ,μ, M1,. . . , Mq, R, pr |r∈R 􏼈 􏼉􏼐 􏼑,(1) where Γis a finite alphabet constructed by all the objects in the PDP systems and µis a membrane structure (MS), consisting of qmembranes, labelled as 1, 2, ...,q. The skin membrane is marked as 1. We associate electrical charges with membranes from the set {−, 0, +}, negative, neutral, and positive. M i , 1 ≤i≤q, are finite multisets over Γ, representing multisets of objects initially placed in the nregions delimited by the membranes of the hierarchical structure µ.R i (R i ∈R), 1≤i≤n, are a finite set of rules of the following form: u[v]α i⟶pru′v′ 􏼂 􏼃β i,(2) where u,v,u′, and v′are multisets over Γ(possibly empty, except uand/or v) and pr is a real number between 0 and 1 associated with the rule, and αand βare electric charges, with α,β∈{−, 0, +}. Besides, in each computation step, the same left-hand side of the rule can produce different evolutionary states (e.g., surviving or not), and the sum of these rules sharing their left-hand side (including the electrical charge) must always be equal to 1. 2.1.Rule Analysis. In order to model real-life ecosystems, the rules are abstracted from the behaviour of the species considered, food distributions, natural disasters, bacterium reactions, etc. From the point of view of the certainty of the application of the rules, two kinds of operational rules might be somehow distinguished. The first type would be rules without explicit probability written; this is equivalent to a probability of 1; that is, these rules will be executed whenever selected, which will happen whenever they are applicable, and no rules are competing for the same objects involved (if more than one rule is competing, they would be chosen nondeterministically). The other type of rules, in this sense, would be the rules with explicit probability lower than 1; that is, rules which are once selected will be executed depending on their probability. Let us consider an example taking Figure 1(a) as an example, involving two rules following the general schema presented for any rule in R i : r1≡u[v]α 2⟶ pr u′v′ 􏼂 􏼃β 2, r2≡u[v]α 2⟶ 1−pr u′[λ]β 2. (3) The pattern transformation of the system above would be as follows: let us suppose an object vhas moved from region 1 into region delimited by membrane 2 and then starts evolving by using some of the two rules in 2. Thus, if rule r 1 is selected according to its probability pr, then object vis rewritten into object v′in region 2 (of course, this could be any multiset), with usimultaneously evolving to u′in region 1; if r 2 is selected instead, then object vis removed from the skin membrane. The system halts when reaching a given condition, typically a number of iterations or cycles of the evolution of the system, because, usually, when modelling complex systems, there is no beginning or end (different from P systems generating numbers, computing functions, or solving computationally hard problems); instead, in this case, the result of the computation is indeed the observation of the system itself, including whichever elements (individuals and other possible variables involved) subject to study. Definition 2 (see [70]). A multienvironment P system of degree (m,q) with m≥1 and q≥1, taking T≥1 time units, is a tuple: Π�G, Γ,􏽘, T, Ej􏼌􏼌􏼌􏼌􏼌1≤j≤m 􏼚 􏼛, RE, 􏼒 Πk�Γ,μ, R, Mi,j 􏼌􏼌􏼌􏼌􏼌1≤i≤q, 1≤j≤m 􏼚 􏼛, 􏼚 pr,j 􏼌􏼌􏼌􏼌􏼌r∈R^ 1≤j≤m 􏼚 􏼛,1≤k≤m􏼛􏼓, (4) where G�(V,S) is a directed graph such that (x,x)∈S, for each x∈V. Let V�{e 1 ,e 2 ,. . .,e m }, whose elements are called environments. Γis the working alphabet, and Σ⊊Γ is an alphabet describing the objects that can be presented in different environments. R E is a finite set of communication rules between two environments, which is of the form rej,ejl ≡(x)p⟶ (x,ji,j2,...,jh)x1 ′  􏼁ej1,. . . , xh ′  􏼁ejh,(5) where x,x1 ′,. . .,xh ′∈Σ, (e j ,ejl)∈S(l �1, ...,h), and p(x,j1,j2,&,jh,)(t)∈[0, 1]. For the same left-hand side xej, the 4Complexity sum of functions associated with the rules from R E is equal to 1. (i) Π k � (Γ,μ, R, Mi|, j1≤i≤q, 1<j≤m 􏼈 􏼉)is a P system with a skeleton (Γ,μ,R) included in each one of the m P systems placed inside the menvironments (with the same alphabet, membrane structure, and rules). Thus, each environment e j contains exactly one P system with the same skeleton given by (Γ,μ,R). The only difference among them will be derived from different parameter values inside each environment, with these different values potentially implying different initial multisets Μ i,j inside the P system of each environment and possibly different probabilities p t,j affecting the skeleton rules Rinside each ecosystem: u[v]α i⟶ pr,j u′v′ 􏼂 􏼃β i.(6) (i) Mi,j, 1 ≤i≤q, 1≤j≤m, are the multisets of objects initially present inside each of the qmembranes of the menvironments (ii) Ej, 1 ≤j≤m, are the multisets of objects initially present in the menvironments Taking Figure 1(b) as an example, in the system depicted, there is a P system with a structure similar to the one in Figure 1(a) inside each environment. There is an object x from an environment ej, 1 ≤j≤4, which can move to another environment ek(maybe to more than one at the same time) using rule rej,ejl; during its transmission, object xfrom environment ejcan be rewritten as xi ′, 1 ≤i≤h, in environments ej1to ejh. When studying real-world ecosystems, the first type (single environment) can be appropriate to study the population dynamics of species in a region (it might cover more regions with an enriched structure in μ, but would not be too intuitive nor adequate). On the contrary, the second type (pure multienvironment P systems) is used to model various species distributed in more than one environment, subject to the same rules but possibly different conditions given by parameters (possibly different values for environmental factors, soil conditions, etc.). Thus, for the latter, each environment ejcontains one P system with the P system skeleton (Γ,μ, R)in the structure of Πk(with k aiming to identify different multisets and probability values inside the P system of each environment). Besides, it is worth emphasizing the important role played by the information communication among environments (by sending one object to one or more neighbouring environments, possibly transforming this object into a different one inside each target environment). At the same time, the mP systems placed in different regions are executed synchronously (let us also note that, inside each one of these mP systems, a second level of parallelism is present through the parallel execution of rules in the ninternal membranes of each system). As it can be shown, a single-environment PDP system is a special restricted case of a multienvironment PDP system. In the following, the main uses of PDP systems to model ecosystems are analysed. Thus, a synthesis of several papers published since 2013 illustrates that different types of P systems have been used for predicting population dynamics of ecosystems. Well-known membrane systems can be used for modelling ecosystems in order to assess the projected number of individuals of certain species and their distribution (in terms of ages and locations). So far, a number of endangered species (listed in the literature in Section 1) have been studied. They focus on different species and study different processes and phenomena, and there are some differences in the definition of rules such as counts or types and subtly different membrane structures (e.g., single vs. multienvironment or different number of membranes in the P system skeleton). However, the general structure of the systems and their dynamics can be extracted for a general protocol, as explained in [71]. A simplified version of such a protocol is outlined in the following steps, with a brief explanation about the application of this protocol to the particular scenario studied in this paper. Skin membrane e1 U, V α 21 0 r1, r2 Objects Charge Rules Inner membrane Label Environment (a) e1r1 r2 4 2 p (u, e 3 ) (u)e 2 (u 3 )e 3 e3 r1 3 e4 e2 r4 3 (b) Figure 1: A portion of classified PDP systems used for modelling ecosystems. (a) A single-environment PDP system with two membranes. (b) A multienvironment PDP system with four environments with the same P system skeleton placed inside each environment (their internal structure being omitted). The system shows population activity of the four environments e1to e4. Complexity 5 Step 1: obtain biological data of the species studied. In the present study, some pedigree data are available. The information includes the number of female (male) individuals, age, and birthday or death date. Depending on the biology of the species (usually animals), we need to further know other information which is not recorded in datasets, typically related with processes of interest, for example, their living habits or feeding needs. Step 2: define a conceptual model. According to the evolutionary behaviour of the species, i.e., feeding, reproduction, mortality, and so on, a preliminary general model (conceptual model) is abstracted from these basic processes, and then each module of this model is given a certain priority and sequencing. Step 3: define the computational model. Starting from the conceptual model, the computational model is built based on a mathematical framework, in our case, PDP systems. Natural evolutionary behaviours from the conceptual model are symbolized, representing the underlying processes with the elements of the computational model. The necessary mapping for this model involves (a) designing the membrane structure of the skeleton P systems (cell-like structure) to place inside the environments, including the initial multisets representing individual objects (e.g., an animal ⟶an object) and symbolic food; (b) designing the evolutionary rule sets capturing the main processes affecting the biology of the species under study, according to their living behaviours. Eventually, a complete multienvironment PDP system (or some other equivalent complete model for the ecosystem) is established. This model should be ready to analyse in terms of its functioning under different scenarios. Step 4: choose simulation software. The previous model designed might be analysed with manual traces to validate against real data and later use to formulate hypothesis and check the behaviour of the system under potential scenarios of interest. However, the manual analysis of big complex systems is not only tedious or error-prone but also impractical and directly intractable in certain case studies. Thus, we need simulation tools, where we can debug the models, experimentally validate them, and finally use them for intensive virtual experiments under different scenarios of major interest for the ecosystems under study. In the case of PDP systems and similar types of membrane systems, the most widely used software has been the framework provided by P-Lingua and MeCoSim. For the introduction about this software, refer to Section 4. Step 5: output predicted datasets. Taken the statistical dataset of a certain year as the input, along with all the parameters related with the biology of the species and the conditions of the ecosystem, the system predicts a set of experimental results by using the MeCoSim environment and then running the model loaded for a certain number of cycles (usually years) to get the output. The protocol briefly described above provides an organized generic sequence of steps to design a model based on the PDP system and use it in a practical way to get new insights from the study of the phenomena under study. This not only provides a theoretical but also practical framework for the use of size-based indicators to monitor the ecosystem changes of species. From a conservation and management viewpoint, a key advantage followed with these models based on PDP systems and the tools available (where many potential scenarios can be analysed by simply changing the input data) is that predictions can be obtained according to the evolutionary behaviour of species rather than relying solely on historical baselines that may not be relevant under the current or future environmental conditions. The applications of models in the context studied can also include the analysis of how several parameters, including growth, reproduction, and mortality ratios, or climatic factors, among others, affect predicted changes in the number of species. 3. Application of P Systems to the Study of the Population Dynamics of Ecosystems Concerns about endangered species arise because most of these known species have smaller range, lower reproduction rate, and higher extinction rate, thus causing the sharp decline of population size or extinction. After some disaster or difficult situation, even if the situation gets better later, endangered species are often not considered to be free of threat because the total population might recover but some populations might be scarce even in the undisturbed fragments, thus potentially causing that such populations may not remain after destruction. The scenarios to consider might be very complex to assess in order to incorporate many parameters and processes affecting the species. In this context, it is worth searching for new types of assessing approaches aiming to capture crucial aspects summarizing the change law of the population dynamics. Particularly, based on the fact that P systems can incorporate and quantify several biological behaviours of species (including reproduction, mortality, and the direct and indirect interactions over migration from one place to another one, among others), such devices are applied to model the population dynamics of complex ecosystems. The use of P systems in ecosystem modelling concerns mainly the prediction of population size of numerous endangered species using P system devices such that a prospect of the species population in the coming years can be obtained based on the known processes and parameters and the absence of data about such years. In most cases, multienvironment P systems have been designed and simulated. Numerous systems have been studied, with different processes and factors involved. In every ecosystem analysed, the elements to consider had distinct features based on input parameters such as (1) membrane structures to capture physical and abstract compartments and (2) rules about the natural biological behaviour of the species under study, related processes, human intervention, and so on. In Table 1, we have listed relevant parameters of the P systems (rules, 6Complexity Table 1: Summary of studies that have used a new frontier approach, termed PDP systems with different constraints, to assess the number of endangered species under conditions of different types. Case study Region/condition Comments Colomer et al. [71] Bearded vulture Region: the cliff-nesting and territorial mountains in the Catalan Pyrenees (Northeastern Spain) Condition: single-environment (2, 1) with two electrical charges (0 or +), where the skin region is used to fix reproduction and mortality and the inner one to fix feeding; five wild and domestic ungulates are included as carrion (prey) species. Cardona et al. [72] Bearded vulture Region: Catalan Pyrenees (NE) Condition: single-environment The structure of this system is the same as that of [69]. The only difference is this system is a dynamic P system with the probabilistic approach, while the former used stochastic constants (a rule can be used when the reaction condition reaches a given constant). Cardona et al. [73] Scavenger birds Region: Catalan Pyrenees (NE) Condition: single-environment (2, 1) with two charges. This system considers notnomadic species (also called invasion alien species—see part (b) in Section 3) and density regulation in order to coexist. Subsequently, this model contains 13 species including two new scavenger birds in competition. Colomer et al. [70] Pyrenean chamois Region: Catalan Pyrenees (NE) Condition: multienvironment (11, 4, 1) with three electrical charges (−, 0, +). The model mainly considers four influencing factors: introduced disease such as pestivirus infection, climate change (refer to part (a) in Section 3), hunting, and migrations between areas. Colomer et al. [74] Bearded vulture Region: the cliff-nesting and territorial mountains in the Catalan Pyrenees (Northeast, Spain). Condition: multienvironment The computational model of the probabilistic P system is the same as that of [70] (refer to the third case in this table for the detailed introduction about the model of a P system). Cardona et al. [75] Scavengers/zebra mussel Region: Catalan Pyrenees (NE Spain)/a fluvial reservoir (Riba-roja-Ebro river, NE Spain) Condition: multienvironment For the scavengers, the structure is the same as [69]; hence, many details have been skipped. For mussels, the structure is (5, 17, 1) with tree electrical charges. This model mainly focuses on factors such as water temperature and its effect on reproduction (see part (a) in Section 3 for impacts of the factor), fixation of the mussel to the substrate, movement of larvae, and density regulations. Colomer et al. [76] Scavenger birds Region: Catalan Pyrenees (Spain)/Pyrenean and pre-Pyrenean mountains Condition: multienvironment (2, 2) with the environment change module, where any of species will move to another area when the capacity reaches a threshold. The model studied: (a) 13 species, including three avian scavengers (three types of vultures) as predator species plus six wild ungulates and four domestic ungulates as prey species; (b) the interactions between species; (c) the communication between two areas; and (d) load capacity regulation. Colomer et al. [77] Plant communities Region: (sub) Alpine (NE Spain) Condition: multienvironment (5, 5) with climatic variability (part (a) in Section 3) and orographic factors (part (c)). More importantly, the model first emphasizes on the impact of the plant community module on population dynamics. The remaining modules are similar to those in the previous models. Colomer et al.[78] A carnivore that predates on ungulates and five ungulates Region: Catalan Pyrenees (NE) Condition: single-environment (11, 2) with three electrical charges. This model mainly considers the impacts of environment factors such as weather, orography, and soil conditions on carnivore size. Margalida et al. [79] Scavenger birds Region: Catalan Pyrenees (NE) Condition: multienvironment The model only considers wild ungulates due to the limitation of domestic carcasses. Undoubtedly, this causes an impact on the biomass. The model of the (2, 2) structure verified that when considering only wild ungulates, the ecosystem cannot offer enough food for predators. Complexity 7 membranes, initial configuration, etc.) used in a number of models of ecosystems for different species modelled by using P systems. The listed works are sorted chronologically, i.e., according to the time sequence of different papers studying species (fourteen years from 2005 to 2018, where the first paper in this list about modelling ecosystems using P systems was published in 2005). In the following section, we will mainly introduce the modelling process of each species studied following the basic sequence summarized in the table mentioned. We will mostly focus on three types of animals (bearded vulture, zebra mussel, and Pyrenean chamois) and two additional not-animal species (Arabidopsis thaliana and Vibrio fischeri). Two main aspects will be described for each case study: (a) the information about the geographical environment of the ecosystems and their processes involved; (b) the biological factors affecting the population size in the ecosystem, analysing how to model the biology of the species under study depending on their features. 3.1. Endangered Species 3.1.1. Bearded Vulture. Bearded vulture, Gypaetus barbatus, is a near-threatened species at a worldwide level but a most severely endangered species at a local level in southern Europe. In the ecosystem under study, bearded vultures are generally living in the southern slope of the central Pyrenees (Aragon region, Spain), a mountainous area belonging to the Eurosiberian biogeographic region, which encompasses the three geomorphological regions of the Pyrenees: Axial Pyrenees, Internal Sierras, and External Sierras. Bearded vultures are distributed in different regions within these areas [86]. Bearded vulture is a cliff-nesting and territorial large scavenger. This species is the only vertebrate that feeds almost exclusively on bone remains of herbivores living in the three habitats mentioned, i.e., animals such as red deer, fallow deer, roe deer, and sheep. The remains of these animals were predicted to be the major limiting factor for the survival of avian scavengers during winter and summer Table 1: Continued. Case study Region/condition Comments Margalida and Colomer [80] European vultures (i) Bearded vulture (ii) Egyptian vulture (iii) Cinereous vulture Regions: 10 municipalities in Catalonia, Northern Spain Food source: the four scenarios of food availability considered Condition: multienvironment Taking 10 areas and 4 avian scavengers as the research object, the model considers the impact of climate variations, such as seasons (summer and winter) (part (a) in Section 3), food shortage, density regulation, and changes in species habitats (insufficient resources), on population dynamics. Colomer et al. [81] Zebra mussel Region: reservoir of Ribarroja Condition: multienvironment (40, 17), where the first 20 membranes are used for 20 weeks of reproductive cycle, 16 for the weeks of the second reproductive cycle, and the last two membranes are used to handle regulation and mortality. Huang et al. [82] Captive giant panda Two regions: Chengdu Research Base of Giant Panda Breeding (GPBB)/China Conservation and Research Centre for Giant Panda (CCRCGP) (Wolong) Condition: single-environment (2, 1), where two membranes are used to evolve and store object information; the evolution process of the species: RMF + rescue module, where RMF is also modified as RFM, FMR, or other forms, showing the robustness of the system independently on the order of the modules. Tian et al. [99] Captive giant panda Two regions: GPBB/CCRCGP Condition: single-environment The membrane structure is the same as in [82], and the only difference is that the release module is added to the previous module, that is, RMF + rescue module + release module. Bernardini and Gheorghe [7] The quorum-sensing regulatory networks of the bacterium Vibrio fischeri Region: marine Condition: single-environment Evolutionary rule choices: in the stochastic way (9, 1), where multisets of objects are used to model bags or soups of chemicals, whereas rules are used to model generic biochemical processes. Romero-Campero and P´erezJim´ enez [83] Quorum sensing in Vibrio fischeri Region: marine Condition: multienvironment Rule choices: stochastic approach (N, 25), multicompartmental P system, where N bacteria are randomly placed inside a multienvironment with 25 different regions, that is, there is an uncertain number of bacteria in each region. Valencia-Cabrera et al. [84] Gene regulatory networks Condition: single-environment The first membrane computing model applied to reconstruct the behaviour of logic networks of species with PDP systems. Valencia-Cabrera et al. [85] Gene regulatory networks Case study: Arabidopsis thaliana Condition: single-environment Based on [84], P systems are used to reproduce a logic gene network of (real) Arabidopsis thaliana in order to regulate the flowering processes. 8Complexity [80]. Bearded vulture has a mean lifespan, in wild birds, of 21.4 years [87]. In general, the mean age of the successful reproduction is 11.4 years [88]. With every spawning, usually, only one chick survives due to the aggression, although the species can produce (as frequently does) two eggs. Recently, field technicians from the Conservation of the Bearded Vulture have carried out annual breeding surveys, indicating that the fertility ratio of the species in the Pyrenees is around 30%, which makes this species become one of the rarest raptors. Taking into consideration all the evolutionary characteristics and the core parameters affecting the changes in the population size of bearded vulture, different types of P systems were used to model ecosystems related to bearded vulture. In the initial phase, a bearded vulture model was presented by a single-environment P system, whereas, also, different rule selection methods started to be explored. Thus, in [80], based on the principle of biochemistry reactions, a P system with stochastic constants is used to model bearded vulture, that is, a rule will be executed if the condition of intrinsic reactivity given a certain threshold is met. However, this technique cannot exploit, in general, the full range of rules of the system, involving a number of different processes subject to different natural laws. In order to cope with the problems of limiting the use of rules, a probability-based population dynamic P system is proposed. In such systems, the rules are applied in a probabilistic way, except for these rules without probabilities (implicit probability 1). Experimental results show that, in comparison with the previous P system based on stochastic constants, this system can simulate more precisely the trends of population dynamics of bearded vultures, with a higher accuracy with respect to the validation dataset. According to the analysis for bearded vultures in a region, the reproduction of the vultures in a small area can easily lead to the loss of genetic diversity of these vultures; this phenomenon may happen not only during the founding event but also during subsequent generations when the population remains small and the exchange of individuals with other populations is minimal. In order to modify the breeding rate by increasing genetic diversity and enhance the survival rate of individuals, there should exist different communications among bearded vultures living in different regions. Based on this idea, Colomer et al. [70] applied a multienvironment P system to model bearded vultures of the Catalan Pyrenees. In this system, each environment contained 17 different types of animals corresponding to 13 species. Besides, there was communication between bearded vultures of different regions, that is, individuals can migrate from one region to another one, thus increasing genetic diversity and modifying the survival rate of bearded vultures. In order to properly predict the change in the population size of the vultures, several nature disasters were added to a multienvironment population dynamics P system [80]. Through experimental validation based on simulations and contrast with real data, it was concluded that this model presents good prediction results, regarding experimental data. However, under certain conditions, a significant difference between model and real data was observed. For bearded vultures, it was due to the fact that the initial models did not consider any process to capture the regulation of the populations, which was later incorporated to adequately represent the carrying capacity of a region of a certain ecosystem. 3.1.2. Zebra Mussel. Zebra mussel, Dreissena polymorpha, is a freshwater mussel living in several of the major river basins, including Ribarroja reservoir in the north of Spain. This species is an invasive species, implying dramatic changes in the ecosystems where they settle, in terms of species distribution, water and soil conditions, etc. Its appearance in Spain and several European countries resulted in adverse impacts on industry, economy, and ecology [89, 90]. Zebra mussel is a dioecious species with an r-selected reproductive strategy, consisting in external fertilization and planktonic larval stages. Its success colonizing new environments may be attributed to high fecundity, efficient larval dispersal, few natural controls, and its ability to adhere to hard substrates [91]. Zebra mussel has become a dangerous threat by feeding competition and alternation of river sediments to native mussels. As these native mussels are threatened or endangered, current control strategies in Spain water bodies are therefore limited to avoid spreading of zebra mussel by regulating boating and fishing activities. For these reasons, different biochemical and histological biomarkers have been undertaken to study the impacts on the population dynamics of zebra mussel, thus aiming to control the dispersal of the species over others. In general, traditional approaches applied logistic regression [92], classification and regression tree model [93], rule-based genetic algorithms [94], and maximum entropy method (Maxent) [95] to analyse the alteration of population dynamics of zebra mussel, obtaining a series of good results. However, the use of such equations in the case of the zebra mussel ecosystem imposed some restrictions on its ecological analysis. Hence, P systems were applied to predict the changes in the population of larvae and adult individuals of zebra mussel (i.e., the population dynamics of the species). The most relevant advantage derived from modelling zebra mussel using P systems is that they make it possible not only to mimic the evolutionary features of the population as a whole, exploring the birth or mortality trend of such population, but also add the traceability (at the level of animals instead of populations) of each adult or larvae individual during the evolution of the system. In [75], according to the categories of species and distribution of their regions, a multienvironment population dynamics P system with 5 cells and 17 areas is used to model zebra mussel of the Ribarroja reservoir. Since zebra mussel must breed at a strict temperature, this parameter is also considered in this P system. Comparing with statistical data, the deviation rate is controlled within reasonable 10% of error. Subsequently, Colomer et al. [81] used a PDP system with 40 membranes and 17 areas to model zebra mussel in the reservoir or Ribarroja ecosystem with a more in-depth analysis. This model included a total of 18 environments Complexity 9 r6≡Xi,j[ ]0 2⟶Yi,j 􏽨 􏽩− 2,1≤i≤2, ki,13 ≤j≤ki,5. (13) (vi) Neonatal individuals (gender determination: female or male): r7≡[Y]− 2⟶Yi,0 􏽨 􏽩+ 2,1≤i≤2.(14) (3) Rescue rules: (i) Number of giant pandas rescued from the wild field: r8≡[A]− 2⟶ pccACc[ ]+ 2, cmin ≤c≤cmax.(15) (ii) Sex for rescued giant pandas: r9≡C⟶ pgiCi 􏼂 􏼃0 1,1≤i≤2.(16) (iii) Age for rescued giant pandas: r10 ≡Ci⟶ pgiCi,j+1+j 3􏼕0 1,1≤i≤2,0≤j≤cmax age. 􏼢 (17) (4) Mortality rules: (i) Survival individuals of infancy giant pandas: r11 ≡Yi,j ⟶ 1−ki,6Zi,j􏼕+ 2,1≤i≤2,0≤j<ki,1. 􏼔 (18) (ii) Mortality individuals of infancy giant pandas: r12 ≡Yi,j ⟶ ki,6λ􏼣+ 2 ,1≤i≤2,0≤j<ki,1. 􏼢 (19) (iii) Survival individuals of young giant pandas: r13 ≡Yi,j ⟶ 1−ki,7Zi,j􏼕+ 2,1≤i≤2, ki,1≤j<ki,2. 􏼔 (20) (iv) Mortality individuals of young giant pandas: r14 ≡Yi,j ⟶ ki,7λ􏼣+ 2 ,1≤i≤2, ki,1≤j<ki,2. 􏼢 (21) (v) Survival individuals of adult giant pandas: r15 ≡Yi,j ⟶ 1−ki,8Zi,j􏼕+ 2,1≤i≤2, ki,2≤j<ki,3. 􏼔 (22) (vi) Mortality individuals of adult giant pandas: r16 ≡Yi,j ⟶ ki,8λ􏼣+ 2 ,1≤i≤2, ki,2≤j<ki,3. 􏼢 (23) (vii) Survival individuals of middle-age giant pandas: r17 ≡Yi,j ⟶ 1−ki,9Zi,j􏼕+ 2,1≤i≤2, ki,3≤j<ki,4,1. 􏼔 (24) (viii) Mortality individuals of middle-age giant pandas: r18 ≡Yi,j ⟶ ki,9λ] + 2,1≤i≤2, ki,3≤j<ki,4,1. 􏼢 (25) (ix) Survival individuals of middle-aged and old giant pandas: r19 ≡Yi,j ⟶ 1−ki,10Zi,j 􏽨 􏽩+ 2,1≤i≤2, ki,4.1≤j<ki,4,2. (26) (x) Mortality individuals of middle-aged and old giant pandas: r20 ≡Yi,j ⟶ ki,10 λ􏼣+ 2 ,1≤i≤2, ki,4.1≤j<ki,4,2. 􏼢 (27) (xi) Survival individuals of old giant pandas: r21 ≡Yi,j ⟶ 1−ki,11Zi,j 􏽨 􏽩+ 2,1≤i≤2, ki,4.2≤j<ki,5. (28) (xii) Mortality individuals of old giant pandas: r22 ≡Yi,j ⟶ ki,11 λ􏼣+ 2 ,1≤i≤2, ki,4.2≤j<ki,5. 􏼢 (29) (xiii) Longevity giant pandas: r23 ≡Yi,ki,5⟶λ 􏼔 􏼕+ 2,1≤i≤2. (30) (5) Feeding rules: giant pandas in different periods need to acquire different quantities of food; therefore, we divide feeding rules into three periods such as infancy, young, and other periods. (i) Feeding rules for infancy giant pandas: r24 ≡Zi,jSfi,1Bfi,2Ofi,3 􏽨 􏽩+ 2⟶Wi,j 􏽨 􏽩− 2,1≤i≤2,0≤j<ki,1. (31) 16 Complexity (ii) Feeding rules for young giant pandas: r25 ≡Zi,jSfi,4Bfi,5Ofi,6 􏽨 􏽩+ 2⟶Wi,j 􏽨 􏽩− 2,1≤i≤2, ki,11 ≤j<ki,2. (32) (iii) Feeding rules for giant pandas during other periods: r26 ≡Zi,jSfi,7Bfi,8Ofi,9 􏽨 􏽩+ 2⟶Wi,j 􏽨 􏽩− 2,1≤i≤2, ki,2≤j<ki,5. (33) (6) Update rules: (i) Food removal rules: r27 ≡[S]− 2⟶[λ]0 2, r28 ≡[B]− 2⟶[λ]0 2, r29 ≡[O]− 2⟶[λ]0 2. (34) (ii) Cycle update rules: r30 ≡Wi,j 􏽨 􏽩− 2⟶Xi,j+1[]0 2,1≤i≤2, ki,2≤j<ki,5, r31 ≡F[]− 2⟶[F]0 2, r32 ≡A[]− 2⟶[A]0 2. (35) (1) Parameters’ Description. In the following, some relevant parameters related with the life stages of giant pandas in our model are explained. Thus, symbol ki,1indicates that captive giant pandas are in subadulthood (i.e., this parameter defines the age at which subadult condition is reached); symbol ki,2sets the boundary age when these pandas reach adulthood; symbol ki,3 represents the age when pandas become middle-aged; ki,4,1 represents the initial age when the giant pandas can be considered that some pandas are a bit old (passing from middleaged to the initial phase of old stage); and finally, symbol ki,4,2 stands for clearly the old age stage. Additionally, symbol ki,5 denotes the longevity age (lifespan) of giant pandas. Symbols ki,6–ki,11 denote the mortality of giant pandas for their different age groups, where ki,6stands for the infancy stage, ki,7for subadulthood, ki,8for adulthood, ki,9middle-aged, ki,10 early old stage, and ki,11 for old-age stage. Symbols ki,12 and ki,13 set the range of fertile age of giant pandas for breeding purposes (i.e., they indicate the beginning and the end of the reproductive age, respectively). These ranges are usually coupled with data related with the previous parameters determining the age groups, but in this refined model, they have been considered separately in order to allow different boundaries in ranges used for mortality and the ages affecting the reproduction module, which makes the model more flexible in terms of the possible scenarios to define. Regarding the parameters related with fertility, symbol g1denotes the probability of a fertile female giant panda to give birth to one single giant panda in the reproductive period of a cycle (a year), while g2denotes the probability for this individual to give birth to twins. With respect to feeding parameters, g3is the number of supplied bamboos within a year, while g4and g5represent the corresponding amount of bamboo shoots and other sources of food (fruits, etc.) per year, respectively. This information is related with the total provision of the centres. However, other parameters regulate the amount of actual food required per individual a year. Thus, symbol fi,1is the amount of bamboo shoots consumed by an infant giant panda individual a year (similarly, fi,2and fi,3refer to the needs of bamboo and others, respectively). The same applies to fi,4,fi,5, and fi,6for subadults, and the corresponding symbols fi,7,fi,8, and fi,9refer to the needs of adults. Concerning the rescue module, symbols cmin and cmax define, respectively, the minimum and maximum number of rescued wild giant pandas per year, with cmax age standing for the maximum age of the rescued wild individuals; additionally, pccis the probability of rescuing cwild individuals in a year, while pgiis the probability of such rescued one to have gender i and pajthe probability of such giant panda to have age j. The values of all these constants have been obtained experimentally after cleaning the data and extracting information through statistical measures from the raw data. However, we must be cautious about the scope of the study and the later applicability to other scenarios, given the variability observed and the relatively reduced size of the samples, in terms of the number of years and specificity of the population. For each probabilistic parameter, a significantly large fluctuation was observed in big values, whereas in small ones, there was no obvious change observed in size. In both cases, it became impossible to obtain truly reliable estimates that can be extrapolated for scenarios apart from the one considered or if significantly different conditions or population distributions appear. That being said, some parameters included in the model show the severe effects derived from the natural stochasticity present in the ecosystem. More specifically, there are very fluctuating factors, such as the number of giant pandas rescued in a year, which can influence the evolution of the population and ultimately depends, among other things, on the presence of natural disasters. For sure, these parameters considered have considerable significance in population and conservation biology [108]. (2) Method Recap. This section ended with the definition of the PDP system providing a computational model of the ecosystem subject to study followed by a detailed explanation about the parameters involved. This would be the last step before starting the process to experimentally validate the model against the real data available and according to the judgment of the experts and managers of the giant panda base. However, before getting to this point, several steps were needed as we summarized in the following: Step 1: obtain the data set. Data sets about captive giant pandas come mainly from GPBB. Step 2: initialization: designing a successful membrane system can require plenty of parameters such as mortality, reproduction rate, and rescue rate. In order to provide a proper setup for the model design, these parameters should be initialized first. Complexity 17 Step 3: design of a basic conceptual model: some behavioural packages such as mortality module are abstracted from the observed daily behaviour of the species. According to its evolutionary cycle, certain sequencing is performed in order for the system to evolve successfully. The model described by the whole picture including different building blocks is called a conceptual model. Step 4: design of a computational model: it consists in applying a formal framework through a mathematical model capturing the details of the conceptual model from the previous given step. More specifically, our model must be computational so that it can be directly computed by an abstract machine, not approximated. Step 5: output: simulate and obtain a predicted number of giant pandas. In summary, for our given example—a giant panda population prediction method based on membrane systems—, the detailed introduction of this method is as follows: we first need to count the basic information of giant pandas in the researched region (i.e., counts, age, sex, and so on); then, we design a conceptual model with the execution sequence according to the fragmented habits (reproduction, feeding, death, and rescue) of all researched pandas; next, we can also design a computational model containing the elements set by the formal model, including the proper structure, alphabets, initial multisets, and a set of rules abstracted from the evolutionary habits of the species according to the conceptual model given and the detailed observation and study from the expert on the problem domain; and finally, we will obtain a series of computational results by running simulation software, performing plenty of virtual experiments, and processing the data obtained. The theoretical analysis indicates that, in the absence of real data as a reference, this method can effectively help in analysing potential variation trends in the population size and distribution (in age and gender) so that the evolution of the population of giant pandas can be projected under many possible scenarios given different plausible conditions. 5.3. Experiments. This section is devoted to the detailed description of the experiments conducted on the model designed. Along this work, the framework provided by P-Lingua and MeCoSim [106] was used to debug the model and run our virtual experiments. P-Lingua provides a standard specification language to define P systems of different types, including the probabilistic framework mentioned in this paper. Subsequently, we translate the model presented 0 20 40 60 80 100 120 140 160 180 200 Number of giant pandas 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 20172005 Year (a) 0 50 100 150 2008 2011 2014 20172005 (b) 0 50 100 2008 2011 2014 20172005 (c) Figure 4: Population size of endangered giant pandas published in the last 13 years, where all the used data are shown through three histograms. (a) Ecological data: the total number of real giant pandas in each year. (b) Female data: the number of female individuals in a selected year. (c) Male data: the number of male individuals in a selected year. (That is, NA�NB+NC;Npresents the number of giant pandas.). Such datasets are called statistical datasets which are generally used as the input for prediction. 18 Complexity Survival rate Reproduction rate Mortality rate 0.4 0.45 0.5 0.55 0.6 0.65 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0 0.01 0.02 0.03 0.04 0.05 0.06 2008 2010 2012 2014 20162006 Year 2008 2010 2012 2014 20162006 Year 2008 2010 2012 2014 20162006 Year Figure 5: The statistical survival rate, reproduction rate, and mortality rate of endangered giant pandas from GPBB for 13 years from 2005 to 2017, given three group trajectories reflecting the dominant cause of changes in population size: natural factors; human behaviour; or enigmatic factors. The purpose of three graphs is to offer the basis reference for setting parameters in the process of running PDP systems such that the data in Table 3 are tuned on the basis of Figure 5. (a) Survival rate ranges. (b) Reproduction rate ranges. (c) Mortality rate ranges. Prediction data Prediction Statistic Minimal deviation ratio (MIDR): 0 N1 = N2 Maximal deviation ratio (MADR): 5.99% 60 70 80 90 100 110 2007 2008 2009 20102006 Year (a) Prediction data Prediction Statistic MIDR: 7.29% MADR: 14.41% 80 100 120 140 160 2010 2011 2012 20132009 Year (b) Prediction data Prediction Statistic MIDR: 6.98% MADR: 20.42% 2013 2014 2015 20162012 Year 100 120 140 160 180 200 (c) Prediction data Prediction Statistic MIDR: 2.69% MADR: 3.08% 2018 and 2019: no statistical data 0 100 200 300 2016 2017 2018 20192015 Year (d) Prediction data Prediction Statistic Blue line: Prediction data No red line: no statistical data from 2018 to 2022: 2019 2020 2021 20222018 Year 200 250 300 350 (e) Deviation ratio 2005 2008 2011 2014 –0.1 0 0.1 0.2 0.3 23451 Year (f) Figure 6: Prediction data changes and statistical data changes based on different initial years taken as input data. There is a significant and rough difference in the average per-year prediction data and statistical data between the size-rise and the size-decline phases across the 5 years analysed. Each pair of values corresponds to the predicted or statistic result of the same year. (a) Taking 2005 as an input, a PDP system is used to predict five-year population size from 2006 to 2010, respectively; (b) taking 2008 as an input, prediction years from 2009 to 2013; (c) taking 2011 as an input, prediction years from 2012 to 2016; (d) taking 2014 as an input, prediction years from 2015 to 2019; (e) taking 2017 as an input, prediction years from 2018 to 2022; (f) this graph describes the comparison of the deviation ratio of five groups of datasets (also, see Figure 7). In the five figures, we list the minimal deviation ratio and maximal ratio of each input year (shown by the arrows), where ratio � (N1−N2)/N2, in which N1represents prediction data and N2represents statistical data in reality. In (d) and (e), because the pedigree data only counted the data before 2017, there are no data for these years from 2018 to 2022 (lack of red line). Complexity 19 Table 3: Values of the ecological parameters used in this model for giant pandas (GP for short) (F�female, M �male,and A�±1.35 ×107); for other explanations, refer to part 5.2. Species i ki,1ki,2ki,3ki,4,1ki,4,2ki,5ki,6ki,7ki,8ki,9ki,10 ki,11 ki,12 ki,13 g1g2g3g4g5fi,1fi,2fi,3fi,4fi,5fi,6fi,7fi,8fi,9 GP (F) 1 1 4 8 17 27 34 0.09 0.001 0.007 0.008 0.1 0.15 6 20 0.191 0.098 A A A 0 0 182 2920 2920 292 11680 10950 1276 GP (M) 2 1 4 6 17 27 36 0.05 0.001 0.005 0.0058 0.034 0.091 5 20 0.191 0.098 A A A 0 0 182 2920 2920 292 11680 10950 1276 20 Complexity into P-Lingua language and prepare the setup for the software application based on MeCoSim; then, the model and the data are loaded, the virtual experiments are performed by running the simulations, and finally, the results obtained are processed and analysed. In the following, we highlight the main facts related with a summary of these processes, where we (a) first introduce experimental design and simulation and (b) then analyse the experimental results. 5.3.1. Experimental Design and Simulation. As introduced at the beginning of this paper, PDP systems scale individuallevel processes up to ecosystem structure and dynamics. Here, we present in three steps how pedigree data about giant panda can be used to inform the model, where population size of giant pandas in each year is shown in Figure 4. First, the individual pedigree data of captive giant pandas that can be used to parameterize the model include food consumption; number of rescued individuals; gender and age; and division of individuals in age groups ranging from offspring and subadults to elderly. Once parameterized, the PDP model can be used to simulate the ecosystem under certain scenarios, aiming to predict the number of individuals along the years. This process requires setting of certain biological parameters. Thus, first, we need to calculate the number of births observed and the mortality of the individuals, among others. These data are used to calculate the fixed size-specific survival, reproductive, and mortality rates (see Figure 5). Naturally, reproductive and mortality rates definitely influence the degree of change (increasing or decreasing) on the population size and distribution, in and out of the studied species. As each individual goes through each module of the model changing its status depending on the application of probabilistic rules, at each time step, the number of individuals varies, and this evolution is subject to natural variability among repetitions of the experiment. The predicted changes in the population size through time are uncertain numerically, but by performing a number of experiments, they are controlled within a given confidence interval. The numerical density of species is summed across all individuals at different age groups to get the final results (but of course, also the details per age and gender remain available for possible later studies). These summarized population sizes are outputted at each time step along with predicted changes. The variation in the population size can be described by fitting, at each time step, a straight line (see Figure 6), with the parameters being used in the PDP-based model which are derived from the data shown in Table 3. The predicted evolution of the population and the corresponding resulting changes in its number of individuals can then be confronted with empirical data for comparison or repeating the above process in conjunction with a statistical procedure to formally estimate parameters and their uncertainty (see Figure 7). 5.3.2. Analysis of the Experimental Results. In Figure 6, by comparing varying trajectories of the number of predicted individuals in the species with the trends in the real data, statistically obtained, we identified the changing regulation of population size across the GPBB. In the five subcharts given by this figure, we observe strong evidence that the Deviation rate 2005 (2006–2010) 2008 (2009–2013) 2011 (2012–2016) 2014 (2015–2017) 6.98% 9.6% 4.46% –1.59% 2.69% 20.42% 17.19% 7.29% 0.99% No value No value 5.61% 14.41% 3.08% 0% 4.28% 9.48% 10.16% 9.9% 10.11% –5% 0 5% 10% 15% 20% 25% 30% 23451 Year Figure 7: How the difference occurs in the process of prediction by PDP models (dotted line to solid line). The changes of five-year time interval (2006–2010) deviation rates (2005 as the input); the changes of five-year (2009–2013) deviation rates (2008 as the input); the changes of five-year (2012–2016) deviation rates (2011 as the input); the changes of five-year (2015–2017) deviation rates (2014 as the input); the deviation trajectory generated by the prediction errors. Overlapping prediction parts, i.e., 2009–2010 {input 2005 (5.61% and 0.99%), 2008 (10.11% and 9.9%)}, 2012–2013 {input 2008 (14.41% and 6.98%), 2011 (10.11% and 9.48%)}, and 2015–2016{input 2011 (17.19% and 6.98%), 2014 (4.46% and 2.69%)}, indicate the input data of different years can predict different results for the same year. In this graph, “No value” means no deviation rate due to the lack of statistical data in reality. Complexity 21 average rate of changes in giant panda populations is 1.86% (data from 2005 as the input) per year, as shown in Figure 6(a), almost consistent with the real statistical data, and similar rates of changes occur across other four groups of experiments we provided, i.e., 10.26% (2008 as the input), 12.84% (2011 as the input), and 3.41% (2014 as the input) (see four other figures in Figure 6). These rates obtained suggest that, on average, the number of giant pandas will be predicted with a small error rate within several years. According to these results, we find that the deviation in the number of individuals per year is controlled within 10% except for the special point (Figure 6) although different inputs can also lead to different prediction results in terms of the number of population sizes per year (Figure 7). These deviation results from the model show that there is no single, fully integrated model that can simulate with the same precision all possible scenarios, and the variation in the numbers emerges from the combination of uncertainty in parameters including growth rate, birth rate, and mortality rate rather than caused by a single fact and taking into account the difficulty of parameterizing interactions. 6. Conclusions and Future Work In this article, we provided an overview of membrane computing models for complex ecosystems and a case study on a complex giant panda ecosystem. Membrane computing models used for modelling ecosystems are very promising, yielding truly distributed and parallel implementations. Distribution is mainly manifested in that the living space of each species is closed, and parallelism is mainly manifested in the simultaneous evolution of different species in different regions at the same time, which are in line with the development trend of ecosystems [109–111]. Various P system models have been used to model a large number of species. The differences between these P systems are mainly reflected in the structures and rules of the systems. For the four species most intensively studied, the differences are mainly reflected in the types of species and the types of environment in which species live; in terms of types of rules present in the systems, the distinction gets mainly reflected in the evolutionary behaviours of species and the types of natural disasters they suffer. Most of the ecosystems, described in more detail within the overview, were referred to a variety of models simulated within the framework provided by P-Lingua and MeCoSim. The experimental results show that P systems can approximately predict the trend of population size by mimicking the evolution state of species. As a case study, a single-environment PDP system is used to model giant pandas. It can be seen from the experiment results that the deviation rates in many years between predicted data and statistical data are controlled within 10% expect for those in several years. As P systems can be used to try to predict the number of species in the next few decades based on the current evolutionary behaviours of species, the datasets obtained through the virtual experiments based on the PDP system models provided can assist the decision makers with further prospects enriching the information available in order to make more informed decisions in the future. Finally, as a future research direction, we may consider the following points: (a) designing a multienvironment P system to model captive or wild giant pandas in different regions, significantly increasing the complexity of the model with respect to the previous model presented in this work. (b) Including potential impacts due to natural or not so natural factors (e.g., climate change influence, introduction of invasive species, and habit destruction produced by environmental disasters such as earthquakes); these elements might be considered in the designed models, possibly having a great influence on the number of individuals of the species, especially on wild species, given a certain setup of conditions for different parameters involved and given a certain population inside each area, including the detailed information or estimation of the distribution of gender and age. Data Availability The data used in this paper come from Chengdu Research Base of Giant Panda Breeding. They can be made available upon request to the corresponding author. Conflicts of Interest The authors declare that there are no conflicts of interest regarding the publication of this paper. Authors’ Contributions The research structure was conceived and designed by Y. D. and G. Z.; D. Q. provided the experimental data; L. V. and M. J. wrote the program and designed the experiment method; Y. D. and H. R. wrote the paper and analysed the experimental data; and G. Z. made revisions to the final manuscript. The final manuscript was read and corrected by all authors. Acknowledgments This work was partially supported by the National Natural Science Foundation of China (Grant nos. 61672437, 61972324, and 61702428), New Generation Artificial Intelligence Science and Technology Major Project of Sichuan Province (Grant no. 2018GZDZX0043), and Artificial Intelligence Key Laboratory of Sichuan Province (Grant no. 2019RYJ06). The authors also acknowledge the support of the research project TIN2017-89842-P (MABICAP), cofinanced by Ministerio de Econom´ ıa, Industria y Competitividad (MINECO) of Spain, through the Agencia Estatal de Investigaci´ on (AEI), and by Fondo Europeo de Desarrollo Regional (FEDER) of the European Union. References [1] G. Pǎun, G. Rozenberg, and A. Salomaa, “Membrane computing with external output,” Fundamenta Informaticae, vol. 41, no. 3, pp. 313–340, 1998. 22 Complexity [2] M. J. P´ erez-Jim´ enez, “A computational complexity theory in membrane computing,” Lecture Notes in Computer Science, pp. 125–148, Springer, Berlin, Germany, 2010. [3] A. Pǎun and G. Pǎun, “The power of communication: P systems with symport/antiport,” New Generation Computing, vol. 20, no. 3, pp. 295–305, 2002. [4] G. Pǎun, “P systems with active membranes: attacking NPcomplete problems,” Journal of Automata, Languages and Combinatorics, vol. 6, no. 1, pp. 75–90, 2001. [5] A. J. Tanentzap, S. Walker, R. T. Theo Stephens, and W. G. Lee, “A framework for predicting species extinction by linking population dynamics with habitat loss,” Conservation Letters, vol. 5, no. 2, pp. 149–156, 2012. [6] D. Orellana-Mart´ ın, L. Valencia-Cabrera, A. Riscos-N´uñez, and M. J. P´ erez-Jim´ enez, “Minimal cooperation as a way to achieve the efficiency in cell-like membrane systems,” Journal of Membrane Computing, vol. 1, no. 2, pp. 85–92, 2019. [7] F. Bernardini and M. Gheorghe, “Population P systems,” Journal of Universal Computerence, vol. 10, no. 5, pp. 509– 539, 2004. [8] C. Mart´ ın-Vide, G. P˘ aun, J. Pazos, and A. Rodr´ ıguez-Pat´ on, “Tissue P systems,” Theoretical Computer Science, vol. 296, no. 2, pp. 295–326, 2003. [9] L. Pan and G. P˘aun, “Spiking neural P systems with antispikes,” International Journal of Computers Communications & Control, vol. 4, no. 3, pp. 273–328, 2009. [10] A. P˘ aun and G. Pǎun, “Small universal spiking neural P systems,” BioSystems, vol. 90, no. 1, pp. 48–60, 2007. [11] T. Song, L. Pan, and G. P˘ aun, “Asynchronous spiking neural P systems with local synchronization,” Information Sciences, vol. 219, pp. 197–207, 2013. [12] R. T. A. P˘ aun, F. G. Cabarle, and H. N. Adorna, “Generating context-free languages using spiking neural P systems with structural plasticity,” Journal of Membrane Computing, vol.1, no. 3, pp. 161–177, 2019. [13] M. J. P´erez-Jim´enez, “The P versus NP problem from the membrane computing view,” European Review, vol. 22, no. 1, pp. 18–33, 2014. [14] H. Sharaf, A. Badr, and I. Farag, “Using P system with innate immunity to solve NP complete problems,” Computing & Information Systems, vol. 14, no. 2, 2010. [15] P. Sos´ ık, “P systems attacking hard problems beyond NP: a survey,” Journal of Membrane Computing, vol. 1, no. 3, pp. 198–208, 2019. [16] P. Sos´ ık and A. Rodr´ ıguez-Pat´ on, “Membrane computing and complexity theory: a characterization of PSPACE,” Journal of Computer and System Sciences, vol. 73, no. 1, pp. 137–152, 2007. [17] A. Alhazov, C. Mart´ ın-Vide, and L. Pan, “Solving graph problems by p systems with restricted elementary active membranes,” in Aspects of Molecular Computing, pp. 1–22, Springer, Berlin, Germany, 2003. [18] A. Leporati, L. Manzoni, G. Mauri, A. E. Porreca, and C. Zandron, “Characterizing PSPACE with shallow nonconfluent P systems,” Journal of Membrane Computing, vol. 1, no. 2, pp. 75–84, 2019. [19] G. Zhang, F. Zhou, X. Huang et al., “A novel membrane algorithm based on particle swarm optimization for solving broadcasting problems,” Journal of Universal Computer Science, vol. 18, no. 13, pp. 1821–1841, 2012. [20] X. Huang, G. Zhang, H. Rong, and F. Ipate, “Evolutionary design of a simple membrane system,” in International Conference on Membrane Computing, pp. 203–214, Springer, Berlin, Germany, 2011. [21] Z. Ou, G. Zhang, T. Wang, and X. Huang, “Automatic design of cell-like P systems through tuning membrane structures, initial objects and evolution rules,” International Journal of Unconventional Computing, vol. 9, no. 5-6, pp. 425–443, 2013. [22] Y. Chen, G. Zhang, T. Wang, and X. Huang, “Automatic design of a P system for basic arithmetic operations,” Chinese Journal of Electronics, vol. 23, no. 2, pp. 302–304, 2014. [23] G. Zhang, H. Rong, Z. Ou, M. J. P´ erez-Jim´ enez, and M. Gheorghe, “Automatic design of deterministic and nonhalting membrane systems by tuning syntactical ingredients,” IEEE Transactions on NanoBioscience, vol. 13, no. 3, pp. 363–371, 2014. [24] X. Peng, X. Fan, J. Liu, and H. Wen, “Spiking neural p systems for performing signed integer arithmetic operations,” Journal of Chinese Computer Systems, vol. 34, no. 2, pp. 360–364, 2013. [25] X. Zhang, X. Zeng, and L. Pan, “A spiking neural p system for performing multiplication of two arbitrary natural numbers,” Chinese Journal of Computers, vol. 32, no. 12, pp. 2362–2372, 2009. [26] C. Lu and X. Zhang, “Solving vertex cover problem by means of tissue P systems with cell separation,” International Journal of Computers Communications & Control, vol. 5, no. 4, pp. 540–550, 2010. [27] T. Song, H. Zheng, and J. He, “Solving vertex cover problem by tissue P systems with cell division,” Applied Mathematics & Information Sciences, vol. 8, no. 1, p. 333, 2014. [28] Y. Niu, K. G. Subramanian, I. Venkat, and R. Abdullah, “A tissue P system based solution to quadratic assignment problem,” International Journal of Foundations of Computer Science, vol. 23, no. 7, pp. 1511–1522, 2012. [29] D. D´ ıaz-Pernil, M. A. Guti´ errez-Naranjo, M. J. P´ erezJim´enez, and A. Riscos-N´uñez, “A linear-time tissue P system based solution for the 3-coloring problem,” Electronic Notes in Theoretical Computer Science, vol. 171, no. 2, pp. 81–93, 2007. [30] Y. Niu, Y. Jiang, and J. Xiao, “Time-free solution to 3-coloring problem using tissue P systems,” Chinese Journal of Electronics, vol. 25, no. 3, pp. 407–412, 2016. [31] J. Cooper and R. Nicolescu, “Alternative representations of P systems solutions to the graph colouring problem,” Journal of Membrane Computing, vol. 1, no. 2, pp. 112–126, 2019. [32] A. Alhazov and S. Cojocaru, “Small asynchronous P systems with inhibitors defining non-semilinear sets,” Theoretical Computer Science, vol. 701, pp. 12–19, 2017. [33] P. Sos´ ık, “A catalytic P system with two catalysts generating a non-semilinear set,” Romanian Journal of Information Science and Technology, vol. 16, no. 1, pp. 3–9, 2013. [34] H. A. Christinal, D. Diaz-Pernil, M. A. Gutierrez-Naranjo, and M. J. P´ erez-Jim´ enez, “Thresholding 2d images with celllike p systems,” Romanian Journal of Information Science and Technology, vol. 13, no. 2, pp. 131–140, 2010. [35] H. Peng, J. Wang, M. J. P´erez-Jim´enez, and P. Shi, “A novel image thresholding method based on membrane computing and fuzzy entropy,” Journal of Intelligent & Fuzzy Systems, vol. 24, no. 2, pp. 229–237, 2013. [36] T. Song, S. Pang, S. Hao, A. Rodr´ ıguez-Pat´ on, and P. Zheng, “A parallel image skeletonizing method using spiking neural p systems with weights,” Neural Processing Letters, vol. 50, no. 2, pp. 1485–1502, 2019. Complexity 23 [37] B. Wang, L. L. Chen, and J. Cheng, “New result on maximum entropy threshold image segmentation based on P system,” Optik, vol. 163, pp. 81–85, 2018. [38] R. Yahya, S. Hasan, L. E. George, and B. Alsalibi, “Membrane computing for 2D image segmentation,” International Journal of Advances in Soft Computing and Its Applications, vol. 7, no. 1, pp. 35–50, 2015. [39] G. Zhang, M. Gheorghe, and Y. Li, “A membrane algorithm with quantum-inspired subalgorithms and its application to image processing,” Natural Computing, vol. 11, no. 4, pp. 701–717, 2012. [40] D. Bahuguna, S. Abbas, and J. Dabas, “Partial functional differential equation with an integral condition and applications to population dynamics,” Nonlinear Analysis: Theory, Methods & Applications, vol. 69, no. 8, pp. 2623–2635, 2008. [41] S. Strohm and R. C. Tyson, “The effect of habitat fragmentation on cyclic population dynamics: a reduction to ordinary differential equations,” Theoretical Ecology, vol. 5, no. 4, pp. 495–516, 2012. [42] G. Zhang, H. Rong, F. Neri, and M. J. P´ erez-jim´ enez, “An optimization spiking neural P system for approximately solving combinatorial optimization problems,” International Journal of Neural Systems, vol. 24, no. 5, Article ID 1440006, 2014. [43] E. S´anchez-Karhunen and L. Valencia-Cabrera, “Modelling complex market interactions using PDP systems,” Journal of Membrane Computing, vol. 1, no. 1, pp. 40–51, 2019. [44] C. Buiu, C. Vasile, and O. Arsene, “Development of membrane controllers for mobile robots,” Information Sciences, vol. 187, pp. 33–51, 2012. [45] A. B. Pavel and C. Buiu, “Using enzymatic numerical P systems for modeling mobile robot controllers,” Natural Computing, vol. 11, no. 3, pp. 387–393, 2012. [46] X. Wang, G. Zhang, F. Neri et al., “Design and implementation of membrane controllers for trajectory tracking of nonholonomic wheeled mobile robots,” Integrated Computer-Aided Engineering, vol. 23, no. 1, pp. 15–30, 2016. [47] L. Valencia-Cabrera, D. Orellana-Mart´ ın, M. ´ A. Mart´ ınezdel-Amor, and M. J. P´erez-Jim´enez, “An interactive timeline of simulators in membrane computing,” Journal of Membrane Computing, vol. 1, no. 3, pp. 209–222, 2019. [48] L. Valencia-Cabrera, C. Graciani, I. P´erez-Hurtado-deMendoza, and M. J. P´ erez-Jim´ enez, “A decade of ecological membrane computing applications,” Bulletin of the International Membrane Computing Society, pp. 39–50, 2018. [49] J. N. Welch and C. Leppanen, “The threat of invasive species to bats: a review,” Mammal Review, vol. 47, no. 2, pp. 277–290, 2017. [50] J. H. R. Lambers, “Extinction risks from climate change,” Science, vol. 348, no. 6234, pp. 501-502, 2015. [51] A. Soultan, M. Wikelski, and K. Safi, “Risk of biodiversity collapse under climate change in the Afro-Arabian region,” Scientific Reports, vol. 9, no. 1, p. 955, 2019. [52] T. J. Stohlgren and J. L. Schnase, “Risk analysis for biological hazards: what we need to know about invasive species,” Risk Analysis, vol. 26, no. 1, pp. 163–173, 2006. [53] D. Tilman, R. M. May, C. L. Lehman, and M. A. Nowak, “Habitat destruction and the extinction debt,” Nature, vol. 371, no. 6492, p. 65, 1994. [54] C. J. Carlson, K. R. Burgio, E. R. Dougherty et al., “Parasite biodiversity faces extinction and redistribution in a changing climate,” Science Advances, vol. 3, no. 9, Article ID e1602422, 2017. [55] S. Fei, J. M. Desprez, K. M. Potter et al., “Divergence of species responses to climate change,” Science Advances, vol. 3, no. 5, Article ID e1603055, 2017. [56] C. Proistosescu and P. J. Huybers, “Slow climate mode reconciles historical and model-based estimates of climate sensitivity,” Science Advances, vol. 3, no. 7, Article ID e1602821, 2017. [57] B. Saether, J. Tufto, S. Engen, K Jerstad, O. W Rostad, and J. E Skˆ atan, “Population dynamical consequences of climate change for a small temperate songbird,” Science, vol. 287, no. 5454, pp. 854–856, 2000. [58] M. Stevanovic, A. Popp, H. Lotze-Campen et al., “The impact of high-end climate change on agricultural welfare,” Science Advances, vol. 2, no. 8, Article ID e1501452, 2016. [59] R. A. Fisher, “The wave of advance of advantageous genes,” Annals of Eugenics, vol. 7, no. 4, pp. 355–369, 1937. [60] J. A. Nelder and R. W. M. Wedderburn, “Generalized linear models,” Journal of the Royal Statistical Society. Series A (General), vol. 135, no. 3, pp. 370–384, 1972. [61] T. W. Yee and N. D. Mitchell, “Generalized additive models in plant ecology,” Journal of Vegetation Science, vol. 2, no. 5, pp. 587–602, 1991. [62] A. H. Hirzel, J. Hausser, D. Chessel, and N. Perrin, “Ecological-niche factor Analysis: how to compute habitat-suitability maps without absence data?” Ecology, vol. 83, no. 7, pp. 2027–2036, 2002. [63] F. Gamboa and E. Gassiat, “Bayesian methods and maximum entropy for ill-posed inverse problems,” The Annals of Statistics, vol. 25, no. 1, pp. 328–350, 1997. [64] C. Gonz´ alez-Crespo, E. Serrano, S. Cahill et al., “Stochastic assessment of management strategies for a Mediterranean peri-urban wild boar population,” PLoS One, vol. 13, no. 8, Article ID e0202289, 2018. [65] E. Serrano, A. Colom-Cadena, E. Gilot-Fromont et al., “Border disease virus: an exceptional driver of chamois populations among other threats,” Frontiers in Microbiology, vol. 6, 2015. [66] A. J. McLane, C. Semeniuk, G. J. McDermid, and D. J. Marceau, “The role of agent-based models in wildlife ecology and management,” Ecological Modelling, vol. 222, no. 8, pp. 1544–1556, 2011. [67] S. J. Phillips, R. P. Anderson, and R. E. Schapire, “Maximum entropy modeling of species geographic distributions,” Ecological Modelling, vol. 190, no. 3-4, pp. 231–259, 2006. [68] M. J. V. Zanden, J. D. Olden, J. H. Thorne, and N. E. Mandrak, “Predicting occurrences and impacts of smallmouth bass introductions in north temperate lakes,” Ecological Applications, vol. 14, no. 1, pp. 132–148, 2004. [69] M. Cardona, M. A. Colomer, M. J. P´ erez-Jim´ enez, D. Sanuy, and A. Margalida, “Modeling ecosystems using p systems: the bearded vulture, a case study,” in International Workshop on Membrane Computing, pp. 137–156, Springer, Berlin, Germany, 2008. [70] M. A. Colomer, S. Lav´ ın, I. Marco et al., “Modeling population growth of Pyrenean chamois (Rupicapra p. pyrenaica) by using P-systems,” in International Conference on Membrane Computing, pp. 144–159, Springer, Berlin, Germany, 2010. [71] M. ` A. Colomer, A. Margalida, and M. J. P´ erez-Jim´ enez, “Population dynamics P system (PDP) models: a standardized protocol for describing and applying novel bioinspired computing tools,” PLoS One, vol. 8, no. 4, Article ID e60698, 2013. [72] M. Cardona, M. A. Colomer, M. J. P´erez-Jim´enez, D. Sanuy, and A. Margalida, “A P System modeling an ecosystem related 24 Complexity to the bearded vulture,” in Proceedings of the Sixth Brainstorming Week on Membrane Computing, pp. 51–66, 2008. [73] M. Cardona, M. A. Colomer, A. Margalida et al., “AP system based model of an ecosystem of some scavenger birds,” in International Workshop on Membrane Computing, pp. 182–195, Springer, Berlin, Germany, 2009. [74] M. A. Colomer, M. A. Mart´ ınez-del-Amor, I. P´ erez-Hurtado, M. J. P´ erez-Jim´ enez, and A. Riscos-N´uñez, “A uniform framework for modeling based on P systems,” in Proceedings of the 2010 IEEE Fifth International Conference on Bio-Inspired Computing: Theories and Applications (BIC-TA), pp. 616–621, Changsha, China, 2010. [75] M. Cardona, M. A. Colomer, A. Margalida et al., “A computational modeling for real ecosystems based on P systems,” Natural Computing, vol. 10, no. 1, pp. 39–53, 2011. [76] M. ` A. Colomer, A. Margalida, D. Sanuy, and M. J. P´ erezJim´ enez, “A bio-inspired computing model as a new tool for modeling ecosystems: the avian scavengers as a case study,” Ecological Modelling, vol. 222, no. 1, pp. 33–47, 2011. [77] M. A. Colomer, C. Fondevilla, and L. Valencia Cabrera, “A new P system to model the subalpine and alpine plant communities,” in Proceedings of the Ninth Brainstorming Week on Membrane Computing, pp. 91–112, Seville, Spain, 2011. [78] M. A. Colomer, I. P´ erez-Hurtado, M. J. P´ erez-Jim´ enez, and A. Riscos-N´uñez, “Comparing simulation algorithms for multienvironment probabilistic P systems over a standard virtual ecosystem,” Natural Computing, vol. 11, no. 3, pp. 369–379, 2012. [79] A. Margalida, M. A. Colomer, and D. Sanuy, “Can wild ungulate carcasses provide enough biomass to maintain avian scavenger populations? An empirical assessment using a bio-inspired computational model,” PLoS One, vol. 6, no. 5, 2011. [80] A. Margalida and M. A. Colomer, “Modelling the effects of sanitary policies on European vulture conservation,” Scientific Reports, vol. 2, p. 753, 2012. [81] M. Colomer, A. Margalida, L. Valencia, and A. Palau, “Application of a computational model for complex fluvial ecosystems: the population dynamics of zebra mussel Dreissena polymorpha as a case study,” Ecological Complexity, vol. 20, pp. 116–126, 2014. [82] Z. Huang, G. Zhang, and D. Qi, “Application of probabilistic membrane systems to model giant panda population data,” Computer Systems & Applications, vol. 26, no. 8, pp. 252–256, 2017. [83] F. J. Romero-Campero and M. J. P´erez-Jim´enez, “A model of the quorum sensing system in vibrio fischeri using P systems,” Artificial Life, vol. 14, no. 1, pp. 95–109, 2008. [84] L. Valencia-Cabrera, M. Garcia-Quismondo, M. J. P´ erezJim´enez et al., “Modeling logic gene networks by means of probabilistic dynamic P systems,” International Journal of Unconventional Computing, vol. 9, no. 5-6, pp. 445–464, 2013. [85] L. Valencia-Cabrera, M. Garc´ ıa-Quismondo, M. J. P´ erezJim´enez et al., “Analysing gene networks with PDP systems. Arabidopsis thailiana, a case study,” in Proceedings of the Eleventh Brainstorming Week on Membrane Computing, pp. 257–272, 2013. [86] J. A. Gil, G. Ch´eliz, ´ I. Zuberogoitia, and P. L´opez-L´opez, “First cases of polygyny for the bearded vulture gypaetus barbatus in the central Pyrenees,” Bird Study, vol. 64, no. 4, pp. 565–568, 2017. [87] C. J. Brown, “Population dynamics of the bearded vulture Gypaetus barbatus in southern Africa,” African Journal of Ecology, vol. 35, no. 1, pp. 53–63, 1997. [88] R. J. Antor, A. Margalida, H. Frey, R. Heredia, L. Lorente, and J. A. Ses´ e, “First breeding age in captive and wild bearded vultures Gypaetus barbatus,” Acta Ornithologica, vol. 42, no. 1, pp. 114–118, 2007. [89] G. L. Mackie and D. W. Schloesser, “Comparative biology of zebra mussels in Europe and North America: an overview,” American Zoologist, vol. 36, no. 3, pp. 244–258, 1996. [90] A. Ardura, A. Zaiko, Y. J. Borrell, A. Samuiloviene, and E. Garcia-Vazquez, “Novel tools for early detection of a global aquatic invasive, the zebra mussel Dreissena polymorpha,” Aquatic Conservation: Marine and Freshwater Ecosystems, vol. 27, no. 1, pp. 165–176, 2017. [91] J. D. Ackerman, B. Sim, S. J. Nichols, and R. Claudi, “A review of the early life history of zebra mussels (Dreissena polymorpha): comparisons with marine bivalves,” Canadian Journal of Zoology, vol. 72, no. 7, pp. 1169–1179, 1994. [92] S. Hallstan, U. Grandin, and W. Goedkoop, “Current and modeled potential distribution of the zebra mussel (Dreissena polymorpha) in Sweden,” Biological Invasions, vol. 12, no. 1, pp. 285–296, 2010. [93] M. Li, W. Ju, and S. Kumar, “Modeling potential habitats for alien species Dreissena polymorpha in Continental USA,” Acta Ecologica Sinica, vol. 28, no. 9, pp. 4253–4258, 2008. [94] Q. Chen and A. Mynett, “Applications of soft computing to environmental hydroinformatics with emphasis on ecohydraulics modelling,” in Practical Hydroinformatics, pp. 405–420, 2009. [95] J. M. Drake and J. M. Bossenbroek, “Profiling ecosystem vulnerability to invasion by zebra mussels with support vector machines,” Theoretical Ecology, vol. 2, no. 4, pp. 189–198, 2009. [96] F. Masini and S. Lovari, “Systematics, phylogenetic relationships, and dispersal of the chamois (Rupicapra spp.),” Quaternary Research, vol. 30, no. 3, pp. 339–349, 1988. [97] C. Luzzago, E. Ebranati, O. Cabez´ on et al., “Spatial and temporal phylogeny of border disease virus in Pyrenean chamois (Rupicaprap. pyrenaica),” PLoS One, vol. 11, no. 12, Article ID e0168232, 2016. [98] J. Guo, “Wildlife conservation: giant panda numbers are surging—or are they?” Science, vol. 316, no. 5827, pp. 974-975, 2007. [99] H. Tian, G. Zhang, H. Rong et al., “Population model of giant panda ecosystem based on population dynamics P system,” Journal of Computer Applications, vol. 38, no. 5, pp. 1488– 1493, 2018. [100] M. C. Andersen, H. Adams, B. Hope, and M. Powell, “Risk assessment for invasive species,” Risk Analysis, vol. 24, no. 4, pp. 787–793, 2004. [101] M. Garc´ ıa-Quismondo, R. Guti´ errez-Escudero, M. A. Mart´ ınez-del-Amor, E. Orejuela-Pinedo, and I. P´erezHurtado, “P-Lingua 2.0: a software framework for cell-like P systems,” International Journal of Computers Communications & Control, vol. 4, no. 3, pp. 234–243, 2009. [102] M. A. Mart´ ınez-del-Amor, I. P´ erez-Hurtado, M. J. P´ erezJim´ enez, and A. Riscos-N´uñez, “A P-Lingua based simulator for tissue P systems,” The Journal of Logic and Algebraic Programming, vol. 79, no. 6, pp. 374–382, 2010. [103] L. F. Mac´ ıas-Ramos, I. P´ erez-Hurtado, M. Garc´ ıa-Quismondo et al., “AP-Lingua based simulator for spiking neural P systems,” In International Conference on Membrane Computing, pp. 257–281, 2011. Complexity 25