DCBA: Simulating Population Dynamics P Systems with Proportional Object Distribution

DCBA: Simulating Population Dynamics P Systems with Proportional Object Distribution M.A. Mart´ınez-del-Amor1, I. P´erez-Hurtado1, M. Garc´ıa-Quismondo1, L.F. Mac´ıas-Ramos1, L. Valencia-Cabrera1, A. Romero-Jim´enez1, C. Graciani-D´ıaz1, A. Riscos-N´u˜nez1, M.A. Colomer2, M.J. P´erez-Jim´enez1 1Research Group on Natural Computing Department of Computer Science and Artificial Intelligence University of Seville Avda. Reina Mercedes s/n, 41012 Sevilla, Spain E-mail: [email protected], [email protected], [email protected], [email protected], [email protected], [email protected], [email protected], [email protected], [email protected] 2Department of Mathematics University of Lleida Avda. Alcalde Rovira Roure, 191, 25198 Lleida, Spain E-mail: [email protected] Summary. Population Dynamics P systems refer to a formal framework for ecological modelling. The semantics of the model associates probabilities to rules, but at the same time, the model is based on P systems, so the rules are applied in a maximally parallel way. Since the success of the first model using this framework [5], initially called multienvironment probabilistic P systems, several simulation algorithms have been defined in order to better reproduce the behaviour of the ecosystems with the models. BBB and DNDP are previous attempts, which define blocks of rules having the same left-hand side, but do not define a deterministic behaviour when different rules are competing for the same resources. That is, different blocks of rules present in their lefthand side common objects, being applicable at the same time. In this paper, we introduce a new simulation algorithm, called DCBA, which performs a proportional distribution of resources. Keywords: Membrane Computing, Population Dynamics, Simulation Algorithm, Probabilistic P systems, DCBA, P-Lingua, pLinguaCore 1 Introduction Membrane Computing has a far–reaching background on the modelling of biochemical phenomena, within the framework of Computational Systems Biology 28 M.A. Mart´ınez-del-Amor et al. [2, 7, 17, 19], being complementary and an alternative to more classical approaches (i.e. ODEs, Petri Nets, etc). However, in 2011 a Membrane Computing modelling framework for ecosystem dynamics was introduced [3]. Based on this framework, several ecosystem models have already been presented. Some examples are the population dynamics of Gypaetus barbatus [4] and Rupicapra p. pyrenaica [8] in the Catalan Pyrenees, as well as the population density of Dreissena polymorpha in Ribarroja reservoir [3]. Some of the assets of this framework are the ability to analyse the simultaneous evolution of a high number of species, as well as the management of a large number of auxiliary objects. These objects could represent, for instance, grass, biomass or animal bones. The results obtained from the application of the framework on different ecosystems prove its versatility and adaptability. Thus, a straightforward interpretation of the results of the simulations of its models can be easily obtained by checking the states and multisets associated to each one of the membranes. Although this framework allows a direct interpretation of the simulations of its models, the simulation itself is a complicated problem to solve from a practical point of view. Therefore, algorithms capable of capturing the semantics described by the framework are necessary. These algorithms should be able to select rules in the models according to their associated probabilities while keeping the maximal semantics of P systems. In this scenario, the concept of rule block takes form. A rule block is a set of rules whose left hand side (that is, the necessary and sufficient condition for them to be applied) is exactly the same. That is, given a P system configuration, either all or none of the rules in the block can be applied. According to the semantics associated to the modelling framework, one or more blocks are selected on each step of computation. The probability for a block to be selected is calculated out of the probabilities of its rules. Once a rule block is selected, its rules are applied a number of times in a probabilistic manner according to their associated probabilities, also known as local probabilities. Henceforth, the condition of the sum of the probabilities associated to all rules in each block being equal to 1 is imposed. The way in which the blocks and rules in the model are selected depends on the specific simulation algorithm employed. These algorithms should be able to deal with issues such as the possible overlapping of left hand sides from different blocks, which might result on the competition of blocks and rules for objects. So far, several algorithms have been developed in order to capture the semantics defined by the modelling framework. Some of these algorithms are the Binomial Block Based algorithm (BBB) and the Direct Non Deterministic algorithm with Probabilities (DNDP). A comparison on the performance of these algorithms can be found on [9]. DCBA: Simulating Population Dynamics P systems 29 The algorithms mentioned above share a common drawback. This drawback involves the distortion of the way in which blocks and rules are selected. That is, instead of blocks and rules being selected according to its probabilities in a uniform manner, this selection process is biased towards those with the highest probabilities. This paper introduces a new algorithm, known as Direct distribution based on Consistent Blocks Algorithm (DCBA). This algorithm is introduced to solve the aforementioned distortion, thus not biasing the selection process towards the most likely blocks and rules. The rest of the paper is structured as follows: Section 2 introduces preliminary concepts, such as the formal modelling framework of PDP systems and the DNDP algorithm. Section 3 describes the DCBA algorithm, together with a test example to show the differences with DNDP, and some details on the implementation in the PlinguaCore software framework. Section 4 shows the behaviour of DCBA when simulating a real ecosystem model. The simulated model has been adapted and improved from the original version. The paper ends with some conclusions and ideas for future work in Section 5. 2 Preliminaries 2.1 The P system based framework Definition 1. A Population Dynamics P system of degree (q, m)with q≥1, m≥1, taking Ttime units, T≥1, is a tuple (G, Γ, Σ, T, RE, µ, R, {fr,j :r∈R, 1≤j≤m},{Mij : 0 ≤i≤q−1,1≤j≤m}) where: •G= (V, S)is a directed graph. Let V={e1, . . . , em}whose elements are called environments; •Γis the working alphabet and Σ$Γis an alphabet representing the objects that can be present in the environments; •Tis a natural number that represents the simulation time of the system; •REis a finite set of communication rules between environments of the form (x)ej p(x,j,j1,...,jh) −−−→ (y1)ej1. . . (yh)ejh where x, y1, . . . , yh∈Σ,(ej, ejl)∈S(l= 1, . . . , h)and p(x,j,j1,...,jh)(t)∈[0,1], for each t= 1, . . . , T. If p(x,j,j1,...,jh)(t) = 1, for each t, then we omit the probabilistic function. These rules verify the following: ?For each environment ejand for each object x, the sum of functions associated with the rules from REwhose left-hand side is (x)ejcoincides with the constant function equal to 1. •µis a membrane structure consisting of qmembranes, with the membranes injectively labeled by 0, . . . , q −1. The skin membrane is labeled by 0. We also associate electrical charges from the set {0,+,−} with membranes. 30 M.A. Mart´ınez-del-Amor et al. •Ris a finite set of evolution rules of the form r:u[v]α i→u0[v0]α0 iwhere u, v, u0, v0are multisets over Γ,i∈ {0,1, . . . , q −1}, and α, α0∈ {0,+,−}. •For each r∈Rand for each j,1≤j≤m,fr,j is a computable function whose domain is {1, . . . , T}and its range is [0,1], verifying the following: ?For each u, v ∈Γ∗,i∈ {0, . . . , q −1}and α, α0∈ {0,+,−}, if r1, . . . , rz are the rules from Rwhose left-hand side is u[v]α iand the right-hand side have polarization α0, then Pz j=1 frj(t) = 1, for each t, 1≤t≤T. ?If (x)ejis the left-hand side of a rule r∈RE, then none of the rules of R has a left-hand side of the form u[v]α 0, for any u, v ∈Γ∗and α∈ {0,+,−}, having x∈u. •For each j(1 ≤j≤m),M0j,...,Mq−1,j are strings over Γ, describing the multisets of objects initially placed in the qregions of µ, within the environment ej. In other words, a system as described in the previous definition can be viewed as a set of menvironments e1, . . . , emlinked between them by the arcs from the directed graph G. Each environment ejcontains a P system, Πj= (Γ, µ, R, M0j, . . . Mq−1,j ), of degree q, such that M0j,...,Mq−1,j describe the initial multisets for this environment, and every rule r∈Rhas a computable function fr,j (specific for environment j) associated with it. The tuple of multisets of objects present at any moment in the menvironments and at each of the regions of each Πj, together with the polarizations of the membranes in each P system, constitutes a configuration of the system at that moment. At the initial configuration of the system we assume that all environments are empty and all membranes have a neutral polarization. We assume that a global clock exists, marking the time for the whole system, that is, all membranes and the application of all rules (both from REand R) are synchronized in all environments. The P system can pass from one configuration to another by using the rules from R=RE∪Sm j=1 RΠjas follows: at each transition step, the rules to be applied are selected according to the probabilities assigned to them, and all applicable rules are simultaneously applied in a maximal way. When a communication rule between environments (x)ej p(x,j,j1,...,jh) −−−→ (y1)ej1. . . (yh)ejh is applied, object xpasses from ejto ej1, . . . , ejhpossibly modified into objects y1, . . . , yh, respectively. At any moment t, 1 ≤t≤T, for each object xin environment ej, if there exist communication rules whose left-hand side is (x)ej, then one of these rules will be applied. If more than one communication rule can be applied to an object, the system selects one randomly, according to their probability which is given by p(x,j,j1,...,jh)(t). For each j(1 ≤j≤m) there is just one further restriction, concerning the consistency of charges: in order to apply several rules of RΠjsimultaneously to the same membrane, all the rules must have the same electrical charge on their right-hand side. DCBA: Simulating Population Dynamics P systems 31 2.2 DNDP simulation algorithm In this section, the Direct Non-deterministic Distribution with probabilities algorithm (DNDP) [14, 13] is briefly described (algorithm 1). The aim of this algorithm is to perform a non-deterministic object distribution, so rules having common objects in their left-hand sides (object competition) will have the same opportunities to consume objects. The input consists on a PDP system of degree (q, m), and a number Tof time units. The algorithm simulates Ttransition steps of the PDP system. Therefore, it only simulates one computation of the PDP system, by selecting and executing rules in a non-deterministic maximal consistent parallel way. Algorithm 1 DNDP MAIN PROCEDURE Require: A PDP system of degree (q, m) with q≥1, m ≥1, taking Ttime units, T≥1. 1: C0←initial configuration of the system 2: for t= 0 to T−1do 3: C0 t←Ct 4: Initialization 5: First selection phase (consistency). 6: Second selection phase (maximality). 7: Execution of selected rules. 8: Ct+1 ←C0 t 9: end for Similarly to the previous algorithms [14], the transitions of the P system are simulated in two phases, selection and execution, to synchronize the consumption and production of objects. However, selection is divided in two micro-phases: the first one calculates a multiset of mutually consistent applicable rules, and the second assure maximal application by eventually increasing the multiplicity of some rules in the previous multiset, obtaining a multiset of maximal mutually consistent applicable rules. The algorithm is described below, but for more details refer to [13]. First of all, in order to simplify the selection and execution phases, the initialization process constructs two ordered set of rules, Ajand Bj, gathering only rules from REand RΠapplicable in environment ej, in the sense of having the same charge in the left-hand side than the membranes in the configuration. In the first selection phase, a multiset of consistent applicable rules, denoted by R1 jfor each environment ej, is calculated. Moreover, a multiset of possible applicable rules, denoted by R0 j, is also created. We will say that two rules are consistent if they are associated to the same membrane, and they update it to the same charge. It is used in order to store rules having 0 as the number of applications when using the random number generator function. Hence, this multiset allows to have elements with multiplicity 0. 32 M.A. Mart´ınez-del-Amor et al. First, a random order is applied to Aj∪Bj, and stored in an ordered set Dj. Moreover, a copy of the configuration Ct, called C0 t, is created and it is updated each time that a rule is selected (removing the left-hand side). Then, a rule ris applicable if the following holds: it is consistent with the previously (according to the order in Dj) selected rules in R1 j, and the number of possible applications M in C0 tis greater than 0. If a rule ris applicable, a binomial distributed random number of applications nis calculated according to the probability. On the one hand, since C0 thas been updated by the previously selected rules, the number ncannot exceed Mto guarantee a correct object distribution. On the other hand, if the generated number nis 0, the corresponding rule is added to the multiset R0 j, giving another chance to be selected in the next phase (maximality). Note that only rules from R1 jare considered for the consistency condition, since rules from R0 jare not applied in the first selection phase. In the second selection phase, the consistent applicable rules are checked again in order to achieve maximality. Only consistent rules are considered, and they are taken from Rj=R0 j∪R1 j. If one rule r∈Rjhas a number of applications M greater than 0 in C0 t, then Mwill be added to the multiplicity of the rule. In order to fairly distribute the objects among the rules, they are iterated in order with respect to the probabilities. Moreover, one rule from the multiset R0 jcan be checked, so it is possible that another rule from R1 j, inconsistent to this one, have been previously selected. In this case, the consistent condition has to be tested again. An example of several executions of the DNDP algorithm is showed in section 3.3, together with a comparison with the new algorithm introduced in this paper. 3 Direct distribution based on Consistent Blocks Algorithm (DCBA) 3.1 Definitions for blocks and consistency The selection mechanism starts from the assumption that rules in Rcan be classified into blocks of rules having the same left-hand side, following the definitions 2, 3 and 4 below. Definition 2. The left and right-hand sides of the rules are defined as follows: (a) Given a rule r∈RΠof the form r:u[v]α h→u0[v0]α0 h: •The left-hand side of ris defined as LHS(r)=(h, α, u, v), where h∈L, α∈ {0,+,−} and u0, v0∈Γ∗. This corresponds to multiset uin the parent membrane of h, multiset vin membrane h, and membrane hwith charge α. •The right-hand side of ris defined as RHS(r)=(h, α0, u0, v0), where h∈L, α0∈ {0,+,−} and u0, v0∈Γ∗. This corresponds to multiset u0in the parent membrane of h, multiset v0in membrane h, and membrane hwith charge α0. DCBA: Simulating Population Dynamics P systems 33 (b) Given a rule r∈REof the form r: (x)ej→(y1)ej1. . . (yk)ejk: •The left-hand side of ris defined as LHS(r) = (ej, x), corresponding to the multiset with only one occurrence of object xin environment ej. •The right-hand side of ris defined as RHS(r) = (ej1, y1). . . (ejk, yk), corresponding to the kmultisets with single objects y1. . . yk, for each environment ej1. . . ejkrespectively. Definition 3. Rules from RΠcan be classified in blocks associated to (h, α, u, v) as follows: Bh,α,u,v ={r∈RΠ:LHS(r)=(h, α, u, v)}. Definition 4. Rules from REcan be classified in blocks associated to (ej, a)as follows: Bej,a ={r∈RE:∃a∈Σ, LHS(r)≡(a)ej}. Recall that, according to the semantics of the model, the sum of probabilities of all the rules belonging to the same block is always equal to 1 – in particular, rules with probability equal to 1 form individual blocks. Note that rules with overlapping (but different) left-hand sides are classified into different blocks. Definition 5. A block Bh,α,u,v is consistent if and only if ∃α0,∀r∈Bh,α,u,v, charge(RHS(r)) = α0. Definition 6. A consistent block Bh,α,α0,u,v, with h∈H,α, α0∈ {0,+,−}, u, v ∈Γ∗, is of the form Bh,α,α0,u,v ={r∈R:∃u0, v0∈Γ∗:r≡u[v]α h→u0[v0]α0 h}. Remark 1. Note that all the rules r∈Bh,α,α0,u,v are consistent, in the sense that each membrane hwith charge αgoes to the same charge α0when any rule of Bh,α,α0,u,v is applied. Definition 7. Two blocks Bh1,α1,β1,u1,v1and Bh2,α2,β2,u2,v2are mutually consistent with themselves, if and only if (h1=h2∧α1=α2)⇒(β1=β2). Definition 8. A set of blocks B={B1, B2, . . . , Bs}is self consistent (or mutually consistent) if and only if ∀i, j(i6=j⇒Biand Bjare mutually consistent). Remark 2. In such a context, a set of blocks has an associated set of tuples (h, α, α0), that is, a relationship of H×Cin C. Then, a set of blocks is mutually consistent if and only if the associated relationship H×Cin Cis functional. 3.2 DCBA pseudocode This new simulation algorithm for PDP systems has the same general scheme than its predecessor, DNDP (algorithm 1). The main loop (algorithm 2) is divided into two stages: selection and execution of rules, similarly to the DNDP algorithm. 34 M.A. Mart´ınez-del-Amor et al. Algorithm 2 DCBA MAIN PROCEDURE Require: A Population Dynamics P system of degree (q, m), T≥1 (time units), and A≥1 (Accuracy). The initial configuration is then called C0. 1: INITIALIZATION .(Algorithm 4). 2: for t←0 to T−1do 3: C0 t←Ct 4: Calculate probability functions fr,j (t) associated to the rules. 5: SELECTION of rules. .(Algorithm 3) 6: EXECUTION of rules. .(Algorithm 8) 7: Ct+1 ←C0 t 8: end for Note that the algorithm selects and executes rules, but not blocks of rules. Blocks are used by DCBA in order to select rules, and this is made in three micro-stages as seen in algorithm 3. Phase 1 calculates a proportional object distribution to the blocks. Phase 2 assures the maximality by checking the maximal applications of each block. And finally, phase 3 passes from block applications to rule applications by calculating random numbers following the multinomial distribution with the corresponding probabilities. Algorithm 3 SELECTION 1: Selection PHASE 1: distribution .(Algorithm 5) 2: Selection PHASE 2: maximality .(Algorithm 6) 3: Selection PHASE 3: probabilities .(Algorithm 7) Before starting to select and execute rules in the system, some data initialization is required (see algorithm 4). For instance, the selection stage uses a table in order to distribute the objects among the blocks. This table T, also called static table, is used in each time step, so it is initialized only once, at the beginning of the algorithm. The static table has one column per each consistent block of rules, and one row per each pair of object and compartment (i.e. each membrane and the environment in the skeleton). An expanded static table Tjis also constructed for each environment, to consider also blocks from environment communication rules. Finally, two multisets, Bjand Rj sel, are initialized for selected blocks and rules, respectively. Remark 3. The columns of the static table contain the information of their lefthand side of the blocks. The rows of the static table contain the information of the competitions for objects: each block competing for a given object will have a value different to −in the corresponding row. DCBA: Simulating Population Dynamics P systems 35 Algorithm 4 INITIALIZATION 1: Construction of the static distribution table T: •Columns: consistent blocks of rules from RΠ:Bh,α,α0,u,v •Rows: pairs (obj, membr) and (obj0, e), for all object obj ∈Γ,obj0∈Σand membrane membr ∈µ, being ea way to generically identify the environment of the skeleton of the P systems in the multienvironment system. •Values: place 1/k in the element (x, y) of the table T, if the corresponding object to the row xis in their left-hand side of the block given by column y, with multiplicity k. Otherwise, keep unmarked with −. 2: for j= 1 to mdo .(Construct the expanded table Tj) 3: Tj←T..(Initialize the table with the original T) 4: Add to table Tja column for each communication rule block from REassociated to the environment ej, and place the value 1 in the corresponding row for (obj0, e), being obj0the object appearing in the left-hand side. 5: end for 6: Initialize the multisets Bj← ∅ and Rj sel ← ∅ The distribution of objects among the blocks with overlapping left-hand sides is performed in selection phase 1 (algorithm 5). The expanded static table Tjis used for this purpose in each environment. Three filters are defined in order to adapt the Tjto the current state of the system. That is, to select which rule blocks are going to receive objects. The first filter will delete columns of the table corresponding to non applicable rule blocks due to the charges in the left-hand side. The second filter will delete the columns of the rule blocks with no applications in a configuration, because of the objects in the left-hand side. The goal of the third filter is to save space in the table, deleting rows with no correspondence with the non-filtered columns. These three filters are applied at the beginning of phase 1, and the result is a dynamic table Tt j(for the environment jand time step t). Filter functions for selection Phase 1 function Filter 1(table T, configuration C).(By columns and charges) Delete columns from table T, according to the charge of the membrane in the left-hand side of the corresponding block and in the configuration C. return T end function function Filter 2(table T, configuration C).(By columns and multiplicity) Delete columns from table T, such that for any row (obj, membr) or (obj0, e), the multiplicity of that object in Cmultiplied by 1/k (value in the table), returns a number κ, 0≤κ < 1. If all the values for that column are −, it is also filtered. return T end function function Filter 3(table T, configuration C).(By rows and multiplicity) Delete rows from Tof pairs (obj, membr),(obj0, e) according to the multisets of C, those having multiplicity 0. return T end function 42 M.A. Mart´ınez-del-Amor et al. Objects Consistent Blocks b2,0,0,∅,{a4,b4,c2}b2,0,0,∅,{a4,d}b2,0,0,∅,{b5,d2}b1,−,−,{b},{a7}b2,0,0,{a3},∅b2,0,−,{a,b},∅ < a,2>1/4 1/4 - - - - < b,2>1/4 - 1/5 - - - < c,2>1/2 - - - - - < d,2>- 1/1 1/2 - - - < a,1>- - - 1/7 1/3 1/1 < b,1>- - - - 1/1 < b,e>- - - 1/1 - - Table 2: Static table and the block, and the relation with the sum of the row. This relation somehow captures the proportion of objects to be initially assigned to each block. Then, the minimum number of each block (given by the column) is calculated. Objects Consistent Blocks Sum b2,0,0,∅,{a4,b4,c2}b2,0,0,∅,{a4,d}b2,0,0,∅,{b5,d2}b2,0,0,{a3},∅ < a,2>* 90 0.25 |11 0.25 |11 - - 0.5 < b,2>* 72 0.25 |10 - 0.2 |6 - 0.45 < c,2>* 66 0.5 |33 - - - 0.5 < d,2>* 30 - 1.0 |20 0.5 |5 - 1.5 < a,1>* 60 - - - 0.33 |20 0.33 Applications 10 11 5 20 Table 3: Selection Phase 1 - Distribution The next phase, maximality, starts from the remaining objects, selecting new applications of the blocks in a maximal way. The result of this phase is showed in table 4. This table presents the remaining objects (the ones not assigned in phase 1) and the possible blocks to be selected. The blocks are chosen in a random way, as shown in algorithm 6, and the possible applications of the block are calculated. This process guarantees a maximal set of blocks to be selected, with a maximal number of applications of each block. The last row, applications, shows that the block b2,0,0,∅,{a4,b4,c2}is applying 1 time, additional to the number of applications calculated in the distribution phase. Then the phase 3, probability, take place. For each block selected in the previous phases, its number of applications is divided among the rules being part of the DCBA: Simulating Population Dynamics P systems 43 Objects Consistent Blocks b2,0,0,∅,{a4,b4,c2}b2,0,0,∅,{a4,d}b2,0,0,∅,{b5,d2} < a,2>* 6 - - - < b,2>* 7 - - - < c,2>* 46 - - - < d,2>* 9 - - - Applications 1 - - Table 4: Selection Phase 2 - Maximality block, according to their probabilities. As a result, the number of applications of each rule is obtained, as showed in table 5. Rules Simulation 1 Simulation 2 Simulation 3 Simulation 4 Simulation 5 r1.17 10 7 6 7 r1.23 0 4 1 2 r1.31 1 5 3 1 r211 11 11 12 12 r35 5 5 6 6 r4- - - - - r520 20 20 20 20 r6- - - - - Table 5: Simulating Πtest using the DCBA algorithm It is noteworthy that the selection of rules belonging to block 1 {r1.i,1≤i≤3}, in table 5, always follows a multinomial distribution respecting the 3 probabilities. This solves the drawback we showed on table 1. Moreover, it can be seen that the maximality sometimes can give one more application to blocks 2 and 3, in spite of keeping the original 10 applications for block 1 from phase 1. In any case, the number of applications is proportionally distributed, avoiding the distortion of using a random order over the blocks (or rules), as made in the DNDP algorithm. 3.4 Implementation in pLinguaCore In [11], a Java library called pLinguaCore was presented under GPL license. It includes parsers to handle input files and built–in simulators to handle different P System based models. It is not a closed product because developers with knowledge of Java can add new components to the library. Within the scope of this paper, pLinguaCore has been upgraded to provide an implementation of the DBCA, 44 M.A. Mart´ınez-del-Amor et al. thus extending its existing probabilistic model simulation algorithms support. Along with the inclusion of other extensions, regarding to models such as Spiking Neural P Systems and Numerical P Systems, current version of the library, named pLinguaCore 3.0, and featuring an implementation of the introduced DBCA can be downloaded from [21]. In what follows, details of the implementation of the DBCA in pLinguaCore are shown. Data structures, methods, code optimization and bug fixes are reviewed. Going Top-down, Java classes involved in the implementation: •DynamicMatrix. It provides an implementation for the main operations of the DBCA. DynamicMatrix is built as a dynamic map indexed by MatrixKey class objects. MatrixKey objects are implemented as a pair of (MatrixRow, MatrixColumn) class objects. Associated to each MatrixKey object within the map, multiplicity kof the object specified by the MatrixRow row in the left hand side of the rule specified by the MatrixColumn column is stored. Note that kis stored instead of 1/k for accuracy reasons. As different filters are applied over the DynamicMatrix object, a couple of lists of MatrixRow and MatrixColumn objects respectively are associated to the matrix to keep track of its valid cells. Removal of elements from these lists is performed when filters are applied, while the DynamicMatrix object itself is reset in every step of the main loop of selection phase. Thus, DynamicMatrix object can be viewed as a hash table of multiplicities that allows a significant reduction of the required amount of memory for execution of the DBCA. Also, attributes that stores the sum of the multiplicities of the objects in the matrix by row as well as the minimum of the columns are included in the DynamicMatrix class. Inconsistent blocks are controlled by means of a list of pairs of MatrixColumn objects. DynamicMatrix class directly extends from StaticMatrix class. Methods in DynamicMatrix implements the DBCA different phases themselves, remarkably: –initData() initializes valid rows and columns lists in the DynamicMatrix object, clearing up them; also application of rules data structure is initialized. –filterColumns1() computes valid columns and associates them to the DynamicMatrix object; applies Filter 1 to these columns; –filterColumns2() applies Filter 2 over valid columns associated to the DynamicMatrix object, removing the required ones. –checkMutualConsistency() checks mutual consistency over blocks of the DynamicMatrix object; if any inconsistency is found, an exception is thrown and execution of the simulator is halted; a message listing the mutual inconsistent blocks found is shown to the user. DCBA: Simulating Population Dynamics P systems 45 –initFilterRows() computes valid rows and associates them to the DynamicMatrix object; applies Filter 3 to these rows. –filterRows() applies Filter 3 to valid rows, removing the required ones; this method is called inside the main loop of the selection phase, while the previous one is called outside, at the beginning of this phase. –normalizeRowsAndCalculateMinimums() implements the main loop of selection phase. –maximality() implements maximality phase. –executeRules() implements execution phase; remarkably, multinomial distribution is computed by computing binomial distributions, implemented through the specialized CERN Java library (cern.jet.random.Binomial). •StaticMatrix. Provides an implementation for the static matrix used by the DBCA. Similarly to DynamicMatrix class, cells within the matrix are stored as a map indexed by MatrixKey class objects, each one of them associated to a multiplicity. A couple of immutable lists of MatrixRow and MatrixColumn class objects determines the structure of the matrix. Contents of the cells are fixed once initialized. •MatrixRow. Provides an implementation for rows featured in DynamicMatrix, StaticMatrix and MatrixKey objects. Implemented by a pair of String objects representing object and membrane label respectively, it also provides a method for computing the validity of the row, i.e. to determine if the row has to be kept within the DynamicMatrix object with respect to a given environment. •MatrixColumn. Provides an implementation for columns featured in DynamicMatrix, StaticMatrix and MatrixKey objects. An abstract class, its extended and implemented by a couple of classes representing the two kinds of rule blocks: – SkeletonRulesBlock, which implements blocks of skeleton rules. – EnvironmentRulesBlock, which implements blocks of environment rules. Both classes have the same structure: a single object to store the common left hand rule side of the rule, plus a collection to store the several right hand rule side objects that conforms the block. Also, each one provides an specific method for computing the validity of the corresponding column within the dynamic matrix. To conclude, let us note that while conducting the DBCA implementation, several bugs have been fixed in pLinguaCore, notably some of them regarding to the way in which rules are parsed and stored, thus applying beyond the scope of 46 M.A. Mart´ınez-del-Amor et al. the DBCA an affecting to implementation of probabilistic models simulators as a whole: •Multisets of objects are now taken into account while checking rule blocks. In previous versions of pLinguaCore, when checking of the consistency of probabilities of a rule block was conducted (i.e. checking that sum of probabilities of the rules must equal to one), multiplicities of objects in the left hand side of the rules were ignored. •Issues with “intentional duplicate rules” solved assigning an unique identifier for every rule within the scope of probabilistic models. Issues found were: – Instantiation of parameters in syntactically different rule schemes for some models produced duplicated rules and caused the parser to throw an error and halt. As this duplicity proved intentional, the parser was modified subsequently to take it into account. – Probability was not taken into account when differencing rules. This made the parser to discard a rule syntactically identical, except for its probability, to a previous parsed one. 4 Validation 4.1 Improved model for the scavenger bird ecosystem In this section, it is presented a novel model for an ecosystem related to the Bearded Vulture in the Pyrenees (NE Spain), by using PDP systems. This model is an improved model of which is provided in [5]. The Bearded Vulture (Gypaetus barbatus) is an endangered species in Europe that feeds almost exclusively on bone remains of wild and domestic ungulates. In this model, the evolution of six species is studied: The Bearded Vulture and five subfamilies of domestic and will ungulates upon which the vulture feeds. The model consists of a PDP system of degree (2,1), Π= (G, Γ, µ, R, T, {fr:r∈R},M1,M2) where: •Gis an empty graph because RE=∅. •In the alphabet Γ, we represent the six species of the ecosystem (index iis associated with the species and index jis associated with their age, and the symbols X,Yand Zrepresent the same animal but in different states); it DCBA: Simulating Population Dynamics P systems 47 also contains the auxiliary symbol B, which represents 0.5 kg of bones, and C, which allows a change in the polarization of the membrane labeled by 2 at a specific stage. Γ={Xi,j, Yi,j, Zi,j : 1 ≤i≤7,0≤j≤ki,4}∪{B, C} The species are the following: – Bearded Vulture (i= 1) – Pyrenean Chamois (i= 2) – Red Deer Female (i= 3) – Red Deer Male (i= 4) – Fallow Deer (i= 5) – Roe Deer (i= 6) – Sheep (i= 7) •µ= [ [ ]2]1is the membrane structure, and the corresponding initial multisets are: –M1={Xqi,j i,j : 1 ≤i≤7,0≤j≤ki,4} –M2={C, Bα} where α=d 21 P j=1 q1,j ·1.10 ·682e Value αrepresents an external contribution of food which is added during the first year of study so that the Bearded Vulture survives. In the formula, q1,j represents the number of jyears of age of Bearded Vultures, the finality of constant factor 1.10 is to guarantee enough food for 10% population growth. At present, the population growth is estimated an average 4%, but this value can reach higher values. Thus, to avoid problems related with the underestimation of this value the first year we estimated the population growth (overestimated) at 10%. The constant value 682 represents the amount of food needed per year for a Bearded Vulture pair to survive. •Each year in the real ecosystem is simulated by 3 computational steps, so T= 3 ·Y ears, where Y ears is the number of years to simulate. •The rules Rto apply are: – Reproduction rules for ungulates Adult males r0,i,j ≡[Xi,j]1 1−ki,13 −−−→[Yi,j]1:ki,2≤j≤ki,4,2≤i≤7 Adult females that reproduce r1,i,j ≡[Xi,j]1 ki,5ki,13 −−−→[Yi,j, Yi,0]1:ki,2≤j < ki,3,2≤i≤7, i 6= 3 Red Deer females produce 50% of female and 50% of male springs r2,j ≡[X3,j]1 k3,5k3,130.5 −−−→ [Y3,jY3,0]1:k3,2≤j < k3,3 r3,j ≡[X3,j]1 k3,5k3,130.5 −−−→ [Y3,jY4,0]1:k3,2≤j < k3,3 48 M.A. Mart´ınez-del-Amor et al. Fertile adult females that do not reproduce r4,i,j ≡[Xi,j]1 (1−ki,5)ki,13 −−−→ [Yi,j]1:ki,2≤j < ki,3,2≤i≤7 Not fertile adult females r5,i,j ≡[Xi,j]1 ki,13 −−−→[Yi,j]1:ki,3≤j≤ki,4,2≤i≤7 Young ungulates that do not reproduce r6,i,j ≡[Xi,j]1 1 −−−→[Yi,j]1: 0 ≤j < ki,2,2≤i≤7 – Growth rules for the Bearded Vulture r7,j ≡[X1,j]1 k1,6+k1,10 −−−→ [Y1,k1,2−1Y1,j]1:k1,2≤j < k1,4 r8,j ≡[X1,j]1 1−k1,6−k1,10 −−−→ [Y1,j]1:k1,2≤j < k1,4 r9≡[X1,k1,4]1 k1,6 −−−→[Y1,k1,2−1Y1,k1,4]1 r10 ≡[X1,k1,4]1 1−k1,6 −−−→[Y1,k1,4]1 – Mortality rules for ungulates Young ungulates which survive r11,i,j ≡Yi,j[ ]2 1−ki,7−ki,8 −−−→ [Zi,j]2: 0 ≤j < ki,1,2≤i≤7 Young ungulates which die r12,i,j ≡Yi,j[ ]2 ki,8 −−−→[Bki,11 ]2: 0 ≤j < ki,1,2≤i≤7 Young ungulates which are retired from the ecosystem r13,i,j ≡Yi,j[ ]2 ki,7 −−−→[ ]2: 0 ≤j < ki,1,2≤i≤7 Adult ungulates that do not reach the average life expectancy Those which survive r14,i,j ≡Yi,j[ ]2 1−ki,10 −−−→[Zi,j]2:ki,1≤j < ki,4,2≤i≤7 Those which die r15,i,j ≡Yi,j[ ]2 ki,10 −−−→[Bki,12 ]2:ki,1≤j < ki,4,2≤i≤7 Ungulates that reach the average life expectancy Those which die in the ecosystem r16,i ≡Yi,ki,4[ ]2 ki,9+(1−ki,9)ki,10 −−−→ [Bki,12 ]2: 2 ≤i≤7 Those which die and are retired from the ecosystem r17,i ≡Yi,ki,4[ ]2 (1−ki,9)(1−ki,10) −−−→ [ ]2: 2 ≤i≤7 – Mortality rules for the Bearded Vulture r18,j ≡Y1,j[ ]2 1−k1,10 −−−→[Z1,j]2:k1,2≤j < k1,4 r19,j ≡Y1,j[ ]2 k1,10 −−−→[ ]2:k1,2≤j < k1,4 DCBA: Simulating Population Dynamics P systems 49 r20 ≡Y1,k1,4[ ]2 1 −−−→[Z1,k1,2−1]2 r21 ≡Y1,k1,2−1[ ]2 1 −−−→[Z1,k1,2−1]2 – Feeding rules r22,i,j ≡[Zi,jBki,14 ]2 1 −−−→ Xi,j+1[ ]+ 2: 0 ≤j≤ki,4,1≤i≤7 – Balance rules Elimination of remaining bones r23 ≡[B]+ 2 1 −−−→[ ]2 Adult animals that die because they have not enough food r24,i,j ≡[Zi,j]+ 2 1 −−−→[Bki,12 ]2:ki,1≤j≤ki,4,1≤i≤7 Young animals that die because the have not enough food r25,i,j ≡[Zi,j]+ 2 1 −−−→[Bki,11 ]2: 0 ≤j < ki,1,1≤i≤7 Change the polarization r26 ≡[C]+ 2 1 −−−→[C]2 •The constants associated with the rules have the following meaning: –ki,1: Age at which adult size is reached. This is the age at which the animal consumes food as an adult does, and at which, if the animal dies, the amount of biomass it leaves behind is similar to the total left by an adult. Moreover, at this age it will have surpassed the critical early phase during which the mortality rate is high. –ki,2: Age at which it begins to be fertile. –ki,3: Age at which it stops being fertile. –ki,4: Average life expectancy in the ecosystem. –ki,5: Fertility ratio (number of descendants by fertile females). –ki,6: Population growth (this quantity is expressed in terms of 1). –ki,7: Animals retired from the ecosystem in the first years, age < ki,1(this quantity is expressed in terms of 1). –ki,8: Natural mortality ratio in first years, age < ki,1(this quantity is expressed in terms of 1). –ki,9: 0 if the live animals are retired at age ki,4, in other cases, the value is 1. –ki,10: Mortality ratio in adult animals, age ≥ki,1(this quantity is expressed in terms of 1). –ki,11: Amount of bones from young animals, age < ki,1. –ki,12: Amount of bones from adult animals, age ≥ki,1. –ki,13: Proportion of females in the population (this quantity is expressed in terms of 1). 50 M.A. Mart´ınez-del-Amor et al. –ki,14: Amount of food necessary per year and breeding pair (1 unit is equal to 0.5 kg of bones). •In [5], they can be found actual values for the constants associated with the rules as well as actual values for the initial populations qi,j for each species i with age j. There are two sets of initial populations values, one beginning on year 1994 and another one beginning on year 2008. 4.2 Simulation results In [5], a simulator for the model was presented. The authors show a comparison of the results provided by the simulator and actual data obtained from the ecosystem. That simulator was written in C++ and the rules were implemented directly on the source code. So, that is a simulator implemented ad hoc for the model. The simulator does not implement any described simulation algorithm for P systems and does not implement any generic method to define P systems. We have found that ad hoc simulators like the one presented in [5] have a strong coupling design and it is a problem for debugging. So, if the simulator does not reproduce the expected behaviour of the model, what is causing the problem?. In that situation, we could think that: 1. The model is wrong. 2. The rules are not correctly written in the source code. 3. The semantics of the model is not correctly implemented in the source code. It is very difficult to find the cause of the problem with a strong coupling software design. Moreover, if we think that the cause of the problem is, for instance, 2, but it is really 1 or 3, then we can introduce new errors trying to correct it. From a software engineering point of view it is very important to decouple software components, that is the point of view of P-Lingua and pLinguaCore [21]: •The model is designed on a paper. •The rules are written on a P-Lingua file. So, the parser checks the syntactical/semantics errors. •The semantics of the model is implemented on the pLinguaCore library following a good described simulation algorithm. PLinguaCore is a simulation library that accepts the input written in P-Lingua and provides simulations of the defined P systems. For each type of P system, there are one or more simulation algorithms implemented in pLinguaCore. It is a software framework, so it can be expanded with new simulation algorithms. Thus, we have expanded the pLinguaCore library to include the DCBA simulation algorithm for PDP systems, the current version of pLinguaCore is 3.0 and it can be downloaded from [21]. In this section, we use the model of the Bearded Vulture described above to compare the simulation results produced by the pLinguaCore library using two different simulation algorithms: DNDP [14] and DCBA. We also compare the DCBA: Simulating Population Dynamics P systems 51 results of the implemented simulation algorithms with the results provided by the C++ ad hoc simulator and with the actual ecosystem data obtained from [5]. In [22] it can be found the P-Lingua file which defines the model and instructions to reproduce the comparisons. We have set the initial population values with the actual ecosystem values for year 1994. For each simulation algorithm we have made 1000 simulations of 14 years, that is, 42 computational steps. The simulation workflow have been implemented on a Java program that runs over the pLinguaCore library (this Java program can be downloaded from [22]). For each simulated year (3 computational steps), the Java program counts the number of animals for each species i, that is: Xi= ki,4 P j=0 Xi,j. After 1000 simulations, the Java program calculates average values for each year and species and writes the output to a text file. Finally, we have used the GnuPlot software [20] to produce population graphics. In figures 1, 2, 3, 4 ,5, 6 and 7 the population graphics for each species and simulation algorithm are represented. (a) Using DCBA (b) Using DNDP Fig. 1: Evolution of the Bearded Vulture birds