Abstract
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 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 n2, something also achieved in [21] by Ou et al., who applied P systems to calculate power n2 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 n natural 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.
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 computing-based 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) cell-like 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 tissue-like 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 tissue-like 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 cell-like 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 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 m in 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.
(a)
(b)
Definition 1. (see [69]). A single-environment P system of degree q, with q ≥ 1, is a tuple:where Γ is a finite alphabet constructed by all the objects in the PDP systems and µ is a membrane structure (MS), consisting of q membranes, 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. Mi, 1 ≤ i ≤ q, are finite multisets over Γ, representing multisets of objects initially placed in the n regions delimited by the membranes of the hierarchical structure µ. Ri (Ri ∈ R), 1 ≤ i ≤ n, are a finite set of rules of the following form:where u, , u′, and are multisets over Γ (possibly empty, except u and/or ) 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 Ri:
The pattern transformation of the system above would be as follows: let us suppose an object has 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 r1 is selected according to its probability pr, then object is rewritten into object in region 2 (of course, this could be any multiset), with u simultaneously evolving to u′ in region 1; if r2 is selected instead, then object is 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:where G = (V, S) is a directed graph such that (x, x) ∈ S, for each x ∈ V. Let V = {e1, e2, …, em}, whose elements are called environments. Γ is the working alphabet, and Σ ⊊ Γ is an alphabet describing the objects that can be presented in different environments. RE is a finite set of communication rules between two environments, which is of the formwhere x, , …, ∈ Σ, (ej, ) ∈ S (l = 1, …, h), and ∈ [0, 1]. For the same left-hand side , the sum of functions associated with the rules from RE is equal to 1.(i)Πk = is a P system with a skeleton (Γ, μ, R) included in each one of the m P systems placed inside the m environments (with the same alphabet, membrane structure, and rules). Thus, each environment ej 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 pt,j affecting the skeleton rules inside each ecosystem:(i), , are the multisets of objects initially present inside each of the membranes of the environments(ii), , are the multisets of objects initially present in the environmentsTaking 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 from an environment , , which can move to another environment (maybe to more than one at the same time) using rule ; during its transmission, object from environment can be rewritten as , , in environments to .
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 contains one P system with the P system skeleton in the structure of (with 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 P systems placed in different regions are executed synchronously (let us also note that, inside each one of these P systems, a second level of parallelism is present through the parallel execution of rules in the internal 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. 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, 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 [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 distinguishing areas (plus an additional control environment). Besides, within the environments, each corresponding inner P system contained 40 cells, devoted to handle different steps and processes related with the biology of the species and different reproductive cycles (with more regions to control the evolution of the first yearly reproductive cycle, a smaller amount for the second cycle, and different ones to control regulation processes).
The simulation results obtained showed that PDP models provide very useful tools to model complex, partially desynchronized, processes that work in a parallel way. According to the analysis above, P system-based models could predict better results when handling characteristics such as the ones present in this species (with respect to the previous approaches), thus increasing the confidence in the effectiveness of the methodology proposed.
3.1.3. Pyrenean Chamois
Pyrenean chamois, Rupicapra pyrenaica, belonging to genus Rupicapra, is a mountain ungulate distributed over most of the medium- to high-altitude mountain ranges of Andorra, France, and Spain [96]. There are about 53,000 individuals of Pyrenean chamois living in these places, and the status of the species has not always been favorable. For example, in the late 60s, the population decreased down to the edge of extinction due to indiscriminate hunting [74].
In order to mimic the evolutionary behaviour of this species and estimate the effects of introducing pestivirus affecting the species, Colomer et al. [74] modelled Pyrenean chamois using a multienvironment PDP system with thresholds such as maximum density restrictions (i.e., carrying capacity was considered). In such a PDP system, there were four environments (given the availability of data from 4 different areas in the Catalan Pyrenees). Each environment included the same skeleton with 10 regions. According to the requirements captured from ecologists studying Pyrenean chamois, the processes mainly considered were those of feeding, reproduction, and mortality, being these evolutionary behaviours properly abstracted as some rules in the PDP system designed. In the following, we will introduce some key factors that affect the population dynamics of this species.
In this P system, the model designed considered the impacts of diseases caused by border disease virus (BDV, dangerous pestivirus) on the population dynamics of chamois under study, given that the evolution of this species was highly influenced by an infection of BDV [97] in a significant part of certain regions of the Catalan Pyrenees. In addition to this pestivirus, a number of relevant environmental conditions (e.g., weather, soil, and vegetation) were added to this model because these factors can affect the biology of the species, specifically in processes such as reproduction. Besides, human activities (such as hunting and changes in land use) could be added to this model because of the influence they may have on the population dynamics in the ecosystem considered. In the design of the PDP system, these conditions would be incorporated through certain parameters and rules making use of them. Again, this model involved several processes simultaneously taking place during the evolution of the population.
According to the simulation results obtained with the model designed, belonging to ecological data of 22 years, some differences were observed with respect to the real data available to validate the model. It was concluded after a deeper study that these differences between the values obtained with the model and the statistical data could be further reduced by introducing more nature conditions into P systems. The modular, flexible, and extensible nature of membrane systems made it possible to introduce new elements without major changes in the existing model. Differently from other modelling techniques, the modularity of P system enables the reusability of existing structures and rules in the construction of increasingly complex models where new processes are incorporated iteratively.
Due to the fact that monitoring data including weather change varies continuously and due to the existing relations among different natural conditions (e.g., climate change can increase or decrease the spread of diseases and invasions of alien species may also bring new diseases), it is likely that, for assessing population dynamics, we should focus on multi-index fusion, not just single-index study. A few potential examples can enlighten this thought.
3.1.4. Giant Pandas
Giant pandas, Ailuropoda melanoleuca, are listed as endangered species by the IUCN Red Data Book (recently updated to class VU—vulnerable—in the last assessment, 2016). According to the fourth giant panda survey by the National Forestry Administration (NFA), only 1864 wildlife giant pandas and 375 captive giant pandas existed by the end of 2013. It is a sharply debated question whether panda populations are just beginning to regain lost ground or are already healthier than they have been for many years [98]. Recent studies have clarified that giant pandas in more than 200 countries and regions around the world are extinct, or almost extinct, given that present survival rates are exceptionally low.
Facing the survival status of giant pandas, approaches on assessing population size are built. As P systems can incorporate and quantify several behaviours of species (i.e., reproduction, feeding, mortality, etc.), along with direct or indirect species interactions implying communication, P systems are being used to model giant pandas. In [82], based on the ecological data of giant pandas in Chengdu Research Base of Giant Panda Breeding, a probabilistic membrane system was designed. This PDP system included a single environment, with an internal membrane structure consisting of two nested regions. The model included a series of objects (representing giant panda individuals with indexes representing ages and genders plus sources of food and control objects) and evolution rules to represent the main giant panda biological processes. The experimental results obtained showed that the model can approximately simulate the trends of population dynamics of giant pandas. In order to properly mimic the natural behaviours and lifestyles of giant pandas in the real natural environment, based on the ecological data of giant pandas in the China Giant Panda Conservation Research Centre, in [99], a release module was added to the previous population dynamics P system. The simulation results indicated that the maximum deviation rate between the prediction results and the actual data was basically controlled within 4.13%, which improved the accuracy of solutions using the P system compared with [82]. Although these systems have already been dedicated to further improve the quality of the simulation results mentioned, they are still limited to single-environment conditions, with no clear-cut studying for multienvironment P systems so far. Besides, natural disasters are also not considered in the studied systems to this day, leaving room for an important open research line being explored.
3.2. Other Biological Systems
In this section, we survey ecosystem modelling for other biological systems including the quorum-sensing regulatory networks of the bacterium Vibrio fischeri and those of Arabidopsis thaliana. We will shortly outline the contributions of P systems to studies about the behaviour of populations of other ecosystems. For such ecosystems, P systems were used to model the biological phenomena rather than being interested in the particular behaviour of each individual. In the following, model designs and experiment assays applied are briefly described, along with some conclusions about their study.
3.2.1. Vibrio fischeri
Quorum sensing is a cell density-dependent gene regulation system that can manage the expression of specific sets of genes. Certain pathogenic bacteria use quorum sensing to regulate genes encoding extracellular virulence factors [100]. The cell density control of luminescence in the symbiotic marine bacterium Vibrio fischeri is the best-studied quorum-sensing system. Hence, several references used P systems to model this bacterium, concluding that these systems can really simulate the overall effect of the colony.
In [7], Bernardini and Gheorghe used a single-environment P system with membrane structure (9, 1) to model the quorum-sensing regulatory networks of Vibrio fischeri. In this model, there are 9 compartments, each one representing a bacterium. Inside each of the compartment, rules are used to regulate the reactions of the luminescence genes. Besides, objects can move between a compartment and the environment. Single-environment techniques limited the communication between different types of bacteria. In [83], Romero-Campero and Pérez-Jiménez presented a multienvironment P system with stochastic constants in order to allow the individuals of bacterial cells to communicate with each other. Through experimental results, this multienvironment computational model of quorum sensing in Vibrio fischeri was shown to predict efficiently the change of bacterial density.
3.2.2. Arabidopsis thaliana
Gene regulatory networks are useful models based on versatile frameworks for biologists to understand the interactions among genes in living organisms. In order to better reproduce the behaviour and the dynamics of gene networks, a PDP system was provided to accurately simulate the behaviour of different types of gene networks of species. Thus, in [84], Valencia-Cabrera et al. first used membrane computing models to reproduce the behaviour of a gene network constructed by the improved LAPP method. In this model, the state of each gene in the network, at every moment, is needed to be coded by a series of binary numbers in the existing environments. Then, through interactions (regulated by the rules defined), the next state of the genes will be produced. The contribution of each interaction is calculated from the previously generated objects in order to assess the global influence. Once the global clock equals to 0, the system can stop, obtaining the optimal gene networks.
Subsequently, Valencia-Cabrera et al. [85] applied the defined LN DP systems to reconstruct Boolean gene networks and gene dynamics of Arabidopsis thaliana in order to regulate the flowering processes associated with Arabidopsis thaliana on a long-day scenario. By simulating using software MeCoSim, the designed model proved to match the output data obtained by the improved LAPP method.
4. Simulation Software
P system simulators have become important computational tools in the processes of model debugging, validation, and later virtual experimentation, among other purposes. In this section, we introduce some simulation software products to design models for ecosystems by means of P systems.
The primitive simulators were mostly ad hoc simulation tools devoted to very specific problems or pedagogical aids in the understanding of computations based on the first variants of P systems. However, very relevant achievements were made with these tools, in the first decade of membrane computing, as extensively explored in [47]. This initial exciting period led to the need of simulation assistants with a broader scope, thus providing general-purpose tools. A milestone in this evolution was the P-Lingua framework, providing a standard language, along with the corresponding software framework applied to specify and simulate a number of classes of P systems [101]. Since its appearance, P-Lingua has been successfully used to cope with ecosystem modelling problems [76], formal model checking, and several computationally hard problems [101], among others. Each model displays characteristic semantic constraints that determine the way in which the rules are selected. Hence, it is necessary for simulation software to take into account different scenarios when the computational tools of P systems come to the fore. Nowadays, many types of P system models can be specified and simulated with this framework (e.g., P systems with active membranes [100], tissue P system models [102], and spiking neural P systems (SN P) systems [103], among others). As it will also be recalled later, this framework provides the core engine used by MeCoSim.
With the development of simulators, in recent years, different software applications have been applied to the simulation and validation of biological systems. Here, we will introduce several developed software tools in the following sections.
4.1. MetaPlab (Italy, 2008)
The simulator called MetaPlab (MP for short), designed by Verona University, is used for modelling biological phenomena related to metabolism. It is a computational framework for metabolic P systems. This framework consists of the following four layers: (i) MP graph, which takes MP systems as inputs and visualizes them; (ii) MP store data structure, applied to store all the elements of MP in a suitable Java object form; (iii) data processing, plugin-based module, dealing with biological data; and (iv) MP vistas, coping mainly with the graphical representation of MP structures and dynamics. For more details, refer to the website [104]. Through experiment verification, MetaPlab can deal with dynamic computation problems, flux discovery, and regulation discovery problems.
4.2. BioSimWare (Italy, 2010)
In [105], Besozzi et al. developed simulation software called BioSimWare. This product provided a user-friendly framework for complex biological systems, ranging from cellular processes to population change. BioSimWare implemented several stochastic algorithms to simulate the dynamics of single- or multienvironment models, as well as automatic tools to analyse the effect of variation of the system parameters. BioSimWare supports the SBML format and can automatically convert stochastic models into the corresponding deterministic formulation. Prediction data showed that this software can offer a better comprehension for complex biological systems.
4.3. MeCoSim (Spain, 2010)
The simulation environment MeCoSim [106] was designed by Natural Computing Research Group of the University of Seville, Spain, and it is currently available [107]. In the context of this paper, it is worth mentioning that this software is used to run experiments to estimate the biological phenomena of the species under study in the coming years under different conditions, performing the simulation from the P system and the initial data provided by the user. To sum up the process, MeCoSim contains three files: (a) a config file, where we can define the custom app including the setting of how many simulations and cycles to run, which input tables and output fields/charts we want to have in our interface, and other input/output information and parameter-related data (more info at [106, 107]); (b) data file, whose purpose is to record experimental data sets, i.e., concrete input data, output data, and all the parameters corresponding to a specific experiment; and (c) model file, where the specification of the P system itself is given, including the membrane structure, initial multisets received from the scenario, and grammatical rules capturing the rules of the P systems (both skeleton and environment rules). In summary, this software provided a multifunction application for the study, analysis, modelling, visual simulation, model checking, and optimization of all the possible variants of P systems covered by the P-Lingua framework plus some linked external simulators. In this software, some plugins have been developed to provide some analysis and model monitoring capabilities. Based on the powerful programming function of software and mostly its ability to allow the definition of custom applications with interfaces adapted to each particular problem, it has been successfully applied as an assistant tool for the iterative design of ecosystem models such as literature studies mentioned.
In summary, all software analysed mainly focus on different types of ecological systems: MetaPlab is used for modelling the internal mechanisms of biologic systems by means of metabolic P systems, and BioSimWare is used for modelling the multicompartmental complex biological systems. While they achieved great results when applied to different systems, to the best of our knowledge, it is out of the scope of these products to provide the requirements for data verification, model checking, model optimization, and mostly the definition of custom interfaces for the particular end-user applications for ecologists or managers to perform virtual experiments based on validated models. Therefore, due to the need of providing not only a research result by us as model designers but also a final tool for the experts in the giant panda ecosystem not familiar with P systems, MeCoSim has been used to simulate the evolutionary behaviour of many species in nature and has been chosen by this team. Thus, in this work, we will predict giant panda individuals in the coming years by means of MeCoSim.
5. A Case Study: Population Dynamics of Giant Pandas
This section analyses a particular case study where the methodology and tools explained in the previous sections are applied. Specifically, the species of interest is giant panda, and the following sections will describe in further detail the purpose, process, and conclusions derived from our study.
5.1. Data Source
Biological data concern mainly on populations (with distributions of age and gender), input parameters (birth rate, mortality rate, etc.), and rules governing the biology of the species, in this case, giant pandas and its evolutionary behaviours. These elements are indeed most of the crucial parts taken into account when modelling ecosystems making use of P systems (more specifically, PDP systems). In this section, we describe the study of the related data about giant pandas in the context of the specific ecosystem considered. Giant pandas studied are from Chengdu Research Base of Giant Panda Breeding (GPBB, for short). Giant panda population of GPBB consists of individuals currently living in GPBB, those in Chengdu Zoological Garden, and those who were born in GPBB but are living outside of GPBB. The reference basis for our research is giant panda pedigree data compiled by the Chinese Association of Zoological Gardens, including data belonging to 13 years (from 2005 to 2017). Adequate processing of these data sets could lead to the extraction of relevant statistical data referred to the number of female and male giant pandas per year, including the age of each individual, being these data the basis and driving force, along with the deep study of the processes involved in the biological evolution of the species, to model the ecosystem.
5.2. Model Design
Since we only study giant pandas in a region, based on the features of giant pandas, we choose as our modelling framework a single-environment population dynamics P system (PDP system, for short), as introduced in Section 2. Typically, the modelling process based on PDP systems (as well as other with frameworks) involves two stages: the design of a conceptual model followed by that of a computational model; the former one is used to build the building blocks or modules capturing different behaviours and processes related with the species in the ecosystem (such as giant pandas in GPBB) and to give several key parameters required by the model; subsequently, the latter phase is applied, based on the schema obtained by the previous model, aiming to design the specific PDP system detailing the structure and objects involved and setting the computation rules governing the processes analysed in the biology of the species in the ecosystem and its evolution.
5.2.1. Conceptual Model of PDPs
The main goal of this stage, in our case study, is the design of a novel population dynamics P system-based model for captive giant pandas in the regions described at the beginning of this section. The processes of interest in the evolution of the ecosystem include the biological aspects related with their life cycle and other possible phenomena happening in the environment; these processes should be simulated by computers according to the model designed and the input data about giant panda populations and related parameters. PDP prediction focuses on the changes in the female or male population size and distribution of ages. In the following, design of the conceptual model is described, starting with the definition of the modules composing the computational model.
In order to define the modules, we must first describe the life cycle of giant panda, the classification of age groups, and types of food required. In this simplified case study, the whole evolution behaviour of giant panda considered consists of four processes: reproduction, mortality, feeding, and rescue. Thus, the ecosystem model to design should also contain modules for these four processes involved.
In this model, according to the specialists’ understanding of giant panda life habits, age groups are classified into six ranges; a detailed introduction is given in Table 2, where the classification standards of the female and male may have a little difference. For food required, we mainly consider three kinds of food as the necessary feeding sources: bamboo, bamboo shoots, and other sources of food (possibly including milk, fruits, and others).
The whole setup (an entire year) of this conceptual model mainly consists of the four models depicted in Figure 2, where the first three models are executed sequentially, while the rescue module runs in parallel with respect to them. It involves a series of rules inside each module, abstracted from different behaviours. At every instant, each individual will evolve according to the rules in a parallel way with the other individuals. The model graph of this conceptual model is illustrated in Figure 2.
In the following, the description of the four modules is given: reproduction module, mortality module, feeding module, and rescue module.
(1) Reproduction Module. Many new individuals are born every year. Depending on the birth rate, the number of new individuals will emerge depending on the number of female individuals in the reproductive age, with certain variations due to the probabilistic nature of the rules, capturing the inherent randomness present in the nature up to a certain degree. In further extensions of this model, we could determine the number of new individuals also depending on additional or alternative direct or indirect factors.
(2) Feeding Module. During each cycle (in our case, a year), plenty of food is provided to captive giant pandas, where the main types of food are bamboo, bamboo shoots, and others (mostly, fruits or milk). In the captive environment, all required food can be satisfied. However, it is worth incorporating this process in our models in order to also track the food consumption, and, more importantly, we can also consider analysing scenarios with certain damage producing scarcity of food sources due to some natural disasters such as earthquakes or climatic anomalies.
(3) Mortality Module. Some individuals of the population can die with a certain mortality rate in this phase, also subject to probabilistic rules. In comparison with wild giant pandas, the fundamental results showed that the mortality of captive giant pandas is significantly lower, and the longevity of captive individuals is clearly higher, given the improved life conditions and medical technologies of the captive population, also subject to improvements along the years. In this model, both parameters about mortality and longevity can be adjusted on the basis of available statistical data.
(4) Rescue Module. In nature, in the ecosystem under study, the population size is affected not only by the birth and death of individuals but also by the addition of new rescue individuals to the ecosystem. The rescue module mainly captures the change in the information about yearly rescued individuals, including not only the number of giant pandas rescued from the wild each year but also their gender and age. Historically, the number of rescued individuals ranges from zero or one panda rescued in a year to at most two rescued individuals. In addition, each individual has a certain probability of being rescued. This is based on the past evidence but could be refined to include the possibility of more rescued individuals, what despite being unlikely is not completely impossible. On the contrary, it is worth noting that the maximum life span of rescued giant pandas is the same as that of wild giant pandas, until the moment of their rescue. Consequently, the probability of rescuing very old pandas (for captive giant panda standards, let us say beyond 30 years old) is almost negligible because this could contradict the previous statement. Consequently, it is difficult to find and rescue very old giant pandas because of their weak mobility, and the same generally happens with very young individuals, which are too week to survive enough if they are in danger and still lack mobility. Therefore, the rescue situation of giant pandas at these age ranges is not considered. In simulation, the number of rescued giant pandas of each gender and age is pseudo-random, always between the limits given by their lifespan and subject to the considerations just made. The rescued giant pandas, once calculated their number, gender, and age, will be added to the current population at the end of a cycle. Because the rescue module only simulates the phenomenon that rescued individuals may occur every year, the steps of performing this module are not affected by other modules, and in the following cycle, they will be regular members of the captive population.
Before designing this model, which we have just presented, it was necessary to obtain data and qualitative information about different processes and behavioural facts related with the giant panda life cycle. The model had to consider natural factors such as reproduction habits, mortality rates, specific evolutionary behaviour, and conduct patterns, determined according to the actual situation observed. Concerning the rescue module, the data about rescued individuals per age and gender were collected from the past experience, and the causes for these rescues were analysed with the experts in charge of managing the ecosystem. Some decisions were made concerning the estimation of rescued individuals per year, and some increase in the comparative age of these individuals when incorporated into captive life was applied to simulate aging produced by the past wild life in such rescued individuals depending on the amount of years spent in wild conditions.
5.2.2. Computational (Mathematical) Model Based on PDP Systems
The main goals of the research conducted on the population of giant pandas in GPBB were to assess the evolution in the population size along the years under certain given conditions and initial populations by using a model based on P systems. The target model involved a number of processes and parameters related with the biology of the species and the ecosystem that should be translated into the concepts belonging to the formal model used: P systems (that is, elements such as environments, membranes, objects, evolution rules, and so on). Besides, the proper semantic constraints and considerations inherent to P systems had to be taken into account, along with the accuracy of these conditions (like the application of certain probabilities associated with the rules) in mimicking the natural processes involved. For PDP systems, once the P system skeleton to be placed inside the environment was defined, we mainly focused on designing the evolutionary rules to capture the processes taking place in our subject ecosystem.
Due to the fact that only one species was the subject of our study and no movement among regions was considered as part of the initial design, no complex environment interactions were required. Therefore, a single environment was enough for this case study. Inside such environment, the P system skeleton was designed with two membranes: the external membrane, denoted as a skin membrane (directly contained inside the single environment) and labelled as 1; the inner membrane (used to perform most of the evolution operations), labelled as 2. Then, inside these two membranes, all processes occur, sequencing the proper operations of three of the modules and simultaneously performing the actions related with the other one: the rescue. In every module, all the individuals subject to the rules of the module evolve in parallel. A trace of a simulation of the model along a cycle (a year in this case) is illustrated in Figure 3. In the following, the model will be described in a semiformal way through its main constituent parts.
While PDP systems are defined by a more complex structure, given the fact that a single-environment system is defined, in our case, the definition of the PDP system modelling our case study can be simplified as the tuple:where(i), where(1)In the individuals in the population, with indexes and , represents female and male, while stands for the age of the individual.(2)Initially, each giant panda individual is abstracted as an object that will evolve according to different modules capturing the biological processes, transitioning among states as follows: . Object is abstracted as a -year-old male or female panda individual before reproduction modules, object is abstracted as a j-year-old male or female panda individual in mortality modules, object is abstracted as a -year-old male or female survival panda individual, and object is abstracted as a -year-old male or female (i) panda individual after feeding.(3)Symbol stands for the -year-old male or female rescued panda individuals. After getting incorporated into the population, they also transition to regular giant pandas for the next cycle: .(4)Symbols , , and denote different food objects such as bamboos, bamboo shoots, and others (fruits or milk).(5)Symbols and are auxiliary objects, being used to trigger food production and to start the rescue process.(ii) (i.e., it is a hierarchical membrane structure with membranes labelled as 1 and 2, where the skin membrane is labelled as 1).(iii)The rules in will be described later.(iv)With respect to the initial multisets, and are associated with regions 1 and 2, respectively, where(1), with being the number of all the -year-old male or female (gender given by ) pandas(2)
In this system, the rewriting rules (including all the modules) consist of rules for initialization, reproduction, rescue, mortality, feeding, food removal, and update. The rule set includes the rules introduced in the following, grouped by modules:
(1)Initialization rules: food supply:(2)Reproduction rules:(i)Panda individuals in the pre-reproductive age:(ii)Female individuals in the reproductive period who actually reproduce:(iii)Female individuals in the reproductive period who are not reproducing:(iv)Male panda individuals:(v)Aged panda individuals (postreproductive):(vi)Neonatal individuals (gender determination: female or male):(3)Rescue rules:(i)Number of giant pandas rescued from the wild field:(ii)Sex for rescued giant pandas:(iii)Age for rescued giant pandas:(4)Mortality rules:(i)Survival individuals of infancy giant pandas:(ii)Mortality individuals of infancy giant pandas:(iii)Survival individuals of young giant pandas:(iv)Mortality individuals of young giant pandas:(v)Survival individuals of adult giant pandas:(vi)Mortality individuals of adult giant pandas:(vii)Survival individuals of middle-age giant pandas:(viii)Mortality individuals of middle-age giant pandas:(ix)Survival individuals of middle-aged and old giant pandas:(x)Mortality individuals of middle-aged and old giant pandas:(xi)Survival individuals of old giant pandas:(xii)Mortality individuals of old giant pandas:(xiii)Longevity giant pandas:(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:(ii)Feeding rules for young giant pandas:(iii)Feeding rules for giant pandas during other periods:(6)Update rules:(i)Food removal rules:(ii)Cycle update rules:(1) Parameters’ Description. In the following, some relevant parameters related with the life stages of giant pandas in our model are explained. Thus, symbol indicates that captive giant pandas are in subadulthood (i.e., this parameter defines the age at which subadult condition is reached); symbol sets the boundary age when these pandas reach adulthood; symbol represents the age when pandas become middle-aged; represents the initial age when the giant pandas can be considered that some pandas are a bit old (passing from middle-aged to the initial phase of old stage); and finally, symbol stands for clearly the old age stage. Additionally, symbol denotes the longevity age (lifespan) of giant pandas. Symbols – denote the mortality of giant pandas for their different age groups, where stands for the infancy stage, for subadulthood, for adulthood, middle-aged, early old stage, and for old-age stage. Symbols and 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 denotes 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 denotes the probability for this individual to give birth to twins.
With respect to feeding parameters, is the number of supplied bamboos within a year, while and represent 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 is the amount of bamboo shoots consumed by an infant giant panda individual a year (similarly, and refer to the needs of bamboo and others, respectively). The same applies to , , and for subadults, and the corresponding symbols , , and refer to the needs of adults.
Concerning the rescue module, symbols and define, respectively, the minimum and maximum number of rescued wild giant pandas per year, with standing for the maximum age of the rescued wild individuals; additionally, is the probability of rescuing wild individuals in a year, while is the probability of such rescued one to have gender and the probability of such giant panda to have age .
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. 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 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 individual-level 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.
(a)
(b)
(c)
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).
(a)
(b)
(c)
(d)
(e)
(f)
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 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ón (AEI), and by Fondo Europeo de Desarrollo Regional (FEDER) of the European Union.