A new simulation algorithm for multienvironment probabilistic P systems

A new simulation algorithm for multienvironment probabilistic P systems M.A. Mart´ ınez-del-Amor #1, I. P´ erez-Hurtado #2, M.J. P´ erez-Jim´ enez #3, A. Riscos-N´ u˜ nez #4, M. Angels Colomer ∗5 #Dpt. of Computer Science and Artificial Intelligence, University of Seville Avda. Reina Mercedes s/n, 41012 Sevilla, Spain 1[email protected] 2[email protected] 3[email protected] 4[email protected] ∗Department of Mathematics, University of Lleida Avda. Alcalde Rovira Roure, 191, 25198 Lleida, Spain 5[email protected] Abstract—Multienvironment P systems are the base of a general framework for modeling ecosystems dynamics. On one hand, this modeling framework represents the structural and dynamical aspects of real ecosystems in a discrete, modular and compressive way. On the other hand, the inherent randomness and uncertainty of biological systems are captured by using probabilistic strategies. Nowadays, the simulation of these P systems based models is fundamental for experimentation and validation. In this paper, we introduce a new simulation algorithm, called DNDP, which performs object distribution and maximal consistency in the application of rules, that are crucial aspects of these systems. The paper also depicts a parallel implementation of the algorithm, and a comparison with the existing algorithm in PLinguaCore is provided. In order to test the performance of the presented algorithm, several experiments (simulations) have been carried out over four simple P systems with the same skeleton and different number of environments. I. INTRODUCTION Since Gh. P˘ aun introduced Membrane Computing in 1998 [9], this bio-inspired computing paradigm has proved to be an active research area within Natural Computing. The aim of this area is to define computational models which abstract from the functioning (chemical reactions) and structure (chemical substances, membranes and compartments) of living cells. Devices of this model are called P systems, and they consist of a cell-like membrane structure, in the compartments of which are placed multisets of objects that evolve according to given rules in a synchronous (assuming a global clock) non-deterministic maximally parallel way. P systems have several syntactic ingredients: a membrane structure consisting of a hierarchical arrangement of membranes embedded in a skin membrane, and delimiting regions or compartments where multisets of objects (corresponding to chemical substances) and sets of evolution rules (corresponding to reaction rules) are placed. The skin membrane delimits the internal region of the P system with the environment. Every membrane has an associated label, and depending on the P system model, also a polarization that can be modified in the computation. P systems also have two main semantic ingredients: their inherent parallelism and non-determinism. Recent research works have been focused on using P systems as a modeling tool for biological phenomena, within the framework of computational Systems Biology and Population Dynamics [10], [1], [2], [6], [11], [12], being complementary and an alternative to more classical approaches (i.e. ODEs, Petri Nets, etc). They are used as a formalism for describing, and simulating, the behavior of biological systems, with the advantage of providing a discrete and modular formal model. In [3], a P systems based general framework for modeling ecosystems dynamics is presented. This computational modeling has been used for real ecosystems, such as the scavenger bird in the Catalan Pyrenees [4] and the zebra mussel in Ribarroja reservoir [3]. These P system based models are used to study the ecosystem dynamics, analyzing the simultaneous evolution of a high number of species. The modularity of these models also enables to easily adding or removing new ingredients and characteristics. The aim of this P system based modeling tool is to help the ecologists to adopt a priori management strategies for the real system by executing virtual experiments. However, since no in vivo nor in vitro implementations of P systems are yet available, the computation and analysis of these models is currently performed by simulators. Thus, the design of simulators and other related software tools becomes a critical point in the process of model validation, as well as for virtual experimentation. In this sense, a software tool was developed and presented in [3]. It provides, among others, the following features: a graphical user interface (users be ecologists or model designers), definition of model and ecosystem’s parameters, execution of simulations, creation of statistical data in form of tables, graphs, etc. The core of this application is PLinguaCore [7], a software library for Membrane Computing. The models are defined by plain-text files using P-Lingua specification language. The application loads that file, configures the corresponding parameters, calls PLinguaCore to execute, and collects the results of the simulation. Inside the PLinguaCore library, many simulation algorithms are defined for the different P system models. Specifically, for probabilistic P systems, the implemented simulation algorithm is based on binomial distribution and blocks of rules. This simulation algorithm is efficient enough for small and medium instance sizes, but it lacks maximal consistency application of rules and calculation of probability functions. In this paper, we propose a new simulation algorithm, inspired on the algorithm presented by V. Nguyen et al. in [8], that removes the concept of block of rules, and solves the above mentioned restrictions. We also analyze and validate the algorithm and study the new features and constraints. The rest of the paper is structured as follows. Section II describes the modeling framework based on probabilistic P systems. Section III depicts the details of the proposed simulation algorithms. In section IV we show some results by using a toy P system, validating the algorithm. The paper ends with some conclusions and ideas for future work. II. THE PSYSTEM BASED FRAMEWORK Definition 1: A multienvironment functional probabilistic P system with active membranes of degree (q, m)with q≥ 1, m ≥1, taking Ttime units, T≥1, is a tuple (G, Γ,Σ, RE,Π,{fr,j :r∈RΠ,1≤j≤m}, {Mij : 0 ≤i≤q−1,1≤j≤m}) where: •G= (V, S)is a directed graph such that (x, x)∈S, for each x∈V. 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; •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 ejand for each x, the sum of functions associated with the rules from REwhose left–hand side is (x)ejcoincide with the constant function equal to 1. •Π = (Γ, µ, RΠ)where –µis a membrane structure consisting of qmembranes, with the membranes injectively labeled with 0, . . . , q −1. The skin membrane is labeled with 0. We also associate electrical charges from the set {0,+,−} with membranes; and –RΠis a finite set of evolution rules of the form r: u[v]α i→u0[v0]α0 iwhere u, v, u0, v0∈M(Γ), i ∈ {0,1, . . . , q −1}, and α, α0∈ {0,+,−}; •For each r∈RΠand for each j,1≤j≤m,fr,j is a computable function whose domain is {1,2, . . . , T}and its range is contained in [0,1], verifying the following: ?For each u, v ∈M(Γ),i∈ {0, . . . , q −1}and α∈ {0,+,−}, if r1, . . . , rzare the rules from RΠwhose left–hand side is u[v]α i, then Pz j=1 frj(t)=1, for each t, 1≤t≤T. •For each j(1 ≤j≤m),M0j,...,Mq−1,j are strings over Γ, describing the multisets of objects initially placed in the qregions of µ. A multienvironment probabilistic functional extended P system with active membranes of degree (q, m)taking Ttime units (G, Γ,Σ, RE,Π,{fr,j :r∈RΠ,1≤j≤m}, {Mij : 0 ≤i≤q−1,1≤j≤m}) can be viewed as a set of menvironments e1, . . . , emlinked between them by the arcs from the directed graph G. Each environment ejcontains a functional probabilistic P system with active membranes of degree q, each of them with the same skeleton, Π, and such that M0j, . . . Mq−1,j describe their initial multisets. When a communication rule between environments (x)ej p(x,j,j1,...,jh) −−−→ (y1)ej1. . . (yh)ejh is applied, object xpass from ejto ej1, . . . , ejhpossibly modified into objects y1, . . . , yh, respectively. In any moment t,1≤t≤T, at which an object xis in environment ej, one and only one rule will be applied according to its probability which is given by p(x,j,j1,...,jh)(t). We assume that a global clock exists, marking the time for the whole system (for its compartments), that is, all membranes and the application of all rules are synchronized. The tuple of multisets of objects present at any moment in the menvironments and at each of the regions of the P systems located within them, and 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. 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 and all occurrences of the left–hand side of the rules are consumed, as usual. III. BINOMIAL BLOCK BASED SIMULATION ALGORITHM In this section we describe the first simulation algorithm developed for the model presented in [3], multienvironment probabilistic P systems. It is available in the current release of PLinguaCore library [7], and it is implemented following a strategy based on the binomial distribution and blocks of rules. Let us consider a multienvironment probabilistic functional extended P system with active membranes of degree (q, m) with q≥1, m ≥1, taking Ttime units, T≥1,Π= (Γ,Σ, G, RE,Π,{fr,j :r∈RΠ,1≤j≤m},{Mij : 0 ≤ i≤q−1,1≤j≤m}), as defined in section II, where {fr,j :r∈RΠ,1≤j≤m}are constant functions with range in [0,1]. The computation of the system is a sequence of configurations Ct,0≤t≤Tconstructed by the application of rules from R=RE∪Sm j=1 RΠj. The algorithm presented below shows the pseudocode of the binomial block based simulation algorithm. Roughly speaking, it is divided into two main phases: selection and execution. In the first one, the rules are selected, determining the number of times that they will be applied in the simulated step, according to their left-hand sides and the available objects in the configuration. In the second phase, the chosen rules are applied the elected number of times by adding the multisets of the right-hand side to the configuration, and possibly changing the polarization of membranes. Let us now describe the algorithm in detail. The input data for selection phase consists of the configuration in time t,Ct, and the set of defined rules R. The output data of this stage is a multiset of the form Rsel ={hr, nri}, where r∈Rand nr∈Nis the number of times to be executed. Note that only in the execution phase, the righthand side of rules are added to the configuration. Adding new objects before finishing selection phase leads to inconsistent states, since the simulation would select and execute rules from different transitions of the system. 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 and 3 below. Recall that, according to the semantics of the model, the sum of probabilities of all the rules from a 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 in different blocks. Definition 2: Given a rule r∈RΠof the form r:u[v]α h→u0[v0]α0 h: •The left-hand side of the rule ris defined as LHS(r) = u[v]α h, corresponding to the multiset uin the parent membrane of h, the multiset vin the membrane h, and the membrane hwith charge α. •The right-hand side of the rule ris defined as RHS(r) = u0[v0]α0 h, corresponding to the multiset u0in the parent membrane of h, the multiset v0in the membrane h, and the membrane hwith charge α0. Definition 3: Two rules r1and r2have the same left-hand side if LHS(r1) = u1[v1]α1 h1,LHS(r2) = u2[v2]α2 h2and the following holds: •h1=h2and α1=α2 •the strings v1and v2(resp. u1and u2) represent the same multiset of objects placed in the region h1(resp. in the parent region of h1). The rules selection mechanism iterates through the randomly ordered set of blocks of rules, and within each block, rules are selected in a maximal way (i.e. they will consume as many objects from the configuration as possible). More precisely, given a block, the number of times that a rule r is applied is determined according to a binomial distribution B(N, pr)(see definition 4 below) where Nis the number of copies of the multisets in LHS(r)contained in the configuration, and pris the probability associated with r. Definition 4: In probability theory, the binomial distribution B(n, p)is the discrete probability distribution of the number of successes in a sequence of nindependent Bernoulli experiments (success/failure experiment), each of which yields success with probability p. In order to guarantee that we do not consume more resources than what it is available, after determining the number of applications for the first rule of the block (Nr=B(N, pr)), the first parameter of the binomial distribution is reduced so that for the next rule the maximum number of applications matches the available objects (N←N−Nr). This is done for all rules in the block (randomly ordered) except for the last one, which skips the binomial distribution and takes directly all the remaining applications. It is worth noting that this process is equivalent to calculating a multinomial distribution (definition 5). Definition 5: In probability theory, the multinomial distribution M(n, p1, p2, . . . , pk),Pk i=1 pi= 1 is a generalization of the binomial distribution, where each independent experiment can have k > 0different results. The output of this function is composed by nivariables, with 1≤i≤k, and Pk i=1 ni=n. When the selection phase finishes, it returns the multiset Rsel containing the applicable rules and the number of times they will be executed in the simulated step. Execution phase iterates through the selected rules, and adds the corresponding right-hand sides to the configuration (taking into account the number of applications). Next we show the pseudocode for the binomial block based algorithm. Input: A multienvironment functional P system with active membranes of degree (q, m)with q≥1, m ≥1, taking T time units, T≥1. 1: for t←0to T−1do 2: Rsel ←Selection-phase (Ct, R) 3: Ct+1 ←Execution-phase (Ct, Rsel) 4: print Ct+1 5: end for Selection phase 1: Rules from R=RE∪RΠj,1≤j≤mare classified into sets (blocks) Blso that all the rules belonging to a block have the same left-hand side. 2: Let Fb(N, p)be a function that returns a random natural number using the binomial distribution B(N, p). 3: A random order on the family of all blocks of rules is considered. 4: for all blocks of rules Bl={r1, . . . , rs}, according to the selected random order do 5: A random order to the rules {r1, . . . , rs}is selected 6: Let us suppose that the common left-hand side is u[v]α hand their respective probabilistic constants are cr1, . . . , crs 7: The highest number Nis computed so that uNappears in the parent membrane of hand vNappears in membrane h(only if the charge is α) in the configuration Ct. 8: if N > 0then 9: d←1 10: for k←1to s−1, according to the selected order do 11: crk←crk d 12: nrk←F(N, crk) 13: N←N−nrk 14: q←1−crk 15: d←d∗q 16: end for 17: nrs←N 18: for k←1to sdo 19: Ct←Ct−nrk∗LHS(rk) 20: Rsel ←Rsel ∪ hrk, nrki 21: end for 22: end if 23: end for 24: return Rsel Execution phase 1: for all hr, ni ∈ Rsel do 2: if n > 0then 3: Ct←Ct+n∗RHS(r) 4: Update charges of membranes from Ctusing RHS(r) 5: end if 6: end for 7: Ct+1 ←Ct 8: return Ct+1 This simulation algorithm is useful for the majority of models, such as [4], [3]. However, it has several disadvantages that restricts the P systems models to be correctly simulated: •It needs to classify the rules by their left-hand side. The creation of blocks of rules, with the restriction that probabilities have to sum 1, leads to simulate only models with rules distributed in this special way. •It does not handle rules with intersections on their left-hand side (object competition). Rules with common objects but not the same left-hand side are classified in different blocks, so the common objects will be not distributed since these blocks are maximally executed. •It does not check the consistency of charges in the selection of rules. As seen in section II, rules are executed in a maximally consistent way. If there are rules changing the charge of a membrane, and others changing to a different one, they cannot be executed in the same step. Fixing one value to the charge, all the rules consistent to it must, if possible, be executed. •It does not evaluate probabilistic functions related to rules. Only constant probabilities are considered, what is not the case of the new P systems based models These constraints lead to develop new simulation algorithms, looking for flexibility on the semantics to simulate (i.e. do not restrict blocks of rules) and providing supplementary features for the new requirements of models (maximal consistency and probabilistic functions). IV. DND-P SIMULATION ALGORITHM A. Inspiration: Direct Non-deterministic Distribution algorithm (DND) The algorithm to develop has to solve the restrictions mentioned above: classification of rules, object competition, maximal consistency and calculation of probability functions, and possibly improve the performance. In order to solve the first two points, the algorithm should work with the rules individually without classifications and blocks, and also to assure the object distribution to the rules (according to the probabilities) and maximal execution. In [8], V. Nguyen et al. introduce an algorithm performing non-deterministic, maximally parallel object distribution for transition P systems and its hardware implementation. Assuming that it is possible to have more than one solution to the object distribution problem, several approaches are also analyzed: •Indirect approaches: These approaches consider both solutions and non-solutions to the problem in the searching process. For example, the Incremental Approach constructs a solution starting from a non-solution. •Direct approaches: These approaches consider only solutions in the exploration process. For example, the Direct Non-deterministic Distribution algorithm (DND) constructs an object distribution in one step, with a possible second step to fix the maximal parallelism. In order to do that, the algorithm executes two phases: –Forward phase: it is a loop iterating the rules of a region in a random order, choosing a random value for each and storing it for the next rules that can have intersections. –Backward phase: it is a loop iterating the rules of such region in the previous random order, checking if there exists more applications to the rules, and assigning that number to each. It can be considered that the forward phase implements a non-deterministic object distribution, and the backward phase makes the maximal parallelism, both selecting the rules individually (without using blocks). Therefore, the idea of the DND algorithm is suitable for solving the restriction of blocks of rules and object competition in probabilistic P systems. B. Direct Non-deterministic Distribution with probabilities algorithm (DND-P) In this section we describe the pseudocode of our proposed simulation algorithm for multienvironment functional P systems with active membranes. It is called DNDP, which comes from the inspiration on the algorithm DND [8], and the extension for probabilities. Similarly to the previous algorithm (section III), the transitions of the P system are simulated in two phases, selection and execution, in order to synchronize the consumption and creation of objects. However, selection is divided in two micro-phases. The first one calculates a multiset of consistent applicable rules. The second changes the multiplicity of the rules in the previous multiset to assure maximal application, obtaining a multiset of maximally consistent applicable rules. Before selecting and executing rules, some previous calculations are required. For each environment j, a set of rules (Dj), one per each environment, is constructed by using a random order. It is of the form: rules from REapplied to the environment, and rules from RΠwhose probabilities are non zero, and that are applied according to the charges of membranes in the corresponding configuration and environment. Finally, a new feature is provided in this initialization procedure, which is that the probabilities of rules are recalculated for each moment t. In the first selection phase, a multiset of consistent applicable rules, R1 sel, is calculated. For each environment, rules from Djare processed in a iterative process where the rules are selected a number of times, updating the configuration with the left-hand side. In this loop, two copies of the configuration are needed: one to be updated each time that a rule is selected (removing the left-hand side) and the original one before starting the process. They are called C0 tand Ctrespectively. Suppose that the rules in R1 sel ={< r1, n1>, < r2, n2> , . . . , < rx, nx>}have been selected before processing the rule r. Then, ris selected if: •The application of the rule ris consistent with the application of the previously selected rules. Assuming that rule ris of the form ur[vr]αr hr→u0 r[v0 r]α0 r hr, if there exists a selected rule < ri, ni>∈R1 sel, ri=ui[vi]αi hi→ u0 i[v0 i]α0 i hisuch that hr=hiand α0 r6=α0 i, then the rule is considered non consistent, and so, discarded to be applied. •The number of applications of rin the auxiliary configuration C0 tis greater than 0. If one rule < ri, ni>∈R1 sel has common objects with rin the left-hand side, the total amount of applications in C0 thas been (possibly) decreased nitimes, so the rule rwill use the rest of objects in the configuration. In this way, the object competition is solved by serializing the selection of rules. Note that the rules are randomly chosen, so that for each simulation, the first rule trying to get objects will be different. Once a rule is selected (because it is consistent with the rest of selected rules and it is applicable), a random number of applications is calculated. In this version, this number is generated by using a binomial distribution. However, the function can be changed in order to use other probability distributions without any repercussion to the code. It is considered that the number of times to apply the rule depends on the associated probability and on the initial number of applications in Ct, according to the multinomial distribution. For example, consider a rule rxwith probability px, and rule rywith py, both with the same left-hand side, and px+py= 1. Also consider that the total amount of instances of the left-hand side is N. For the rule rx, the number of applications will be calculated by B(N, px) = Nx, so that there are N0=N−Nx instances left. If we calculate the number for ryas B(N0, py), then it is not a real multinomial distribution. The correct calculation would be B(N, py) = Ny. When Nyis greater than the amount of objects left (Ny> N0), the maximum amount will be directly assigned (Ny=N0). When the binomial random number generator returns 0, the corresponding rule is added to a multiset called R2 sel. This multiset is used to give a new chance to the rules to be selected when applying maximality in the second phase. However, it has to be handled in a special way at this phase. The rule cannot reserve the change of a membrane charge since it has not be actually selected. Hence, R1 sel is the multiset of consistent rules selected a number of times greater than 0, and R2 sel is the multiset of rules selected 0times. 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 Rsel =R1 sel ∪R2 sel. If one rule r∈Rsel has a number of applications (N) different to 0in C0 t, this number will be assigned directly to the rule r. 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 R2 sel can be checked, so it is possible that other rule from R1 sel, inconsistent to this one, has been previously selected. In this case, the consistent condition has to be tested again. Finally, execution phase is similar to the binomial algorithm. It will iterate all the rules in R1 sel (maximal applicable consistent rules), and it will add the right-hand side of them to the configuration C0 t. At the end of the process, C0 tis actually the next configuration: the left-hand side of rules has been removed in the first and second selection phases, and the righthand side of rules is added in the execution stage. Next, the pseudocode for the DNDP algorithm is depicted. Input: A multienvironment functional P system with active membranes of degree (q, m, n)with q≥1, m ≥1, n ≥1, taking Ttime units, T≥1, and a natural number K≥1. 1: for t←0to T−1do 2: Ct←configuration of the system at the moment t 3: C0 t←Ct 4: Initialization 5: First selection phase. It generates a multiset of consistent applicable rules. 6: Second selection phase. It generates a multiset of maximally consistent applicable rules. 7: Execution of selected rules. 8: Ct+1 ←C0 t 9: end for Initialization 1: RΠ←ordered set of rules of Π 2: for j←1to mdo 3: RE,j ←ordered set of rules from RErelated to the environment j 4: Aj←ordered set of rules from RE,j whose probability at the moment tis >0 5: Mj←ordered set of pairs hlabel, chargeifor all the membranes from Ctcontained in the environment j 6: Bj← ∅ 7: for each hh, αi ∈ Mj(following the considered order) do 8: Bj←Bj∪ordered set of rules u[v]α h←u0[v0]β h from RΠwhose probability at the moment tis greater than 0for the environment j 9: end for 10: end for First selection phase (consistency) 1: for j←1to mdo 2: R1 sel,j ←the empty multiset 3: R2 sel,j ←the empty multiset 4: Dj←Aj∪Bjwith a random order 5: for each r∈Dj(following the considered order) do 6: if ris consistent with the rules in R1 sel,j then 7: N0←max{number of times that ris applicable to C0 t} 8: if N0>0then 9: if pr,j (t)=1then 10: n←N0 11: else 12: N←max{number of times that ris applicable to Ct} 13: n←Fb(N, pr,j (t)) 14: if n>N0then 15: n←N0 16: end if 17: end if 18: if n > 0then 19: C0 t←C0 t−n·LHS(r) 20: R1 sel,j ←R1 sel,j ∪ {< r, n >} 21: else 22: R2 sel,j ←R2 sel,j ∪ {< r, n >} 23: end if 24: end if 25: end if 26: end for 27: end for Second phase of rules selection (maximality) 1: for j←1to mdo 2: Rsel,j ←R1 sel,j +R2 sel,j with an order by the rule probabilities, from highest to lowest 3: for each < r, n >∈Rsel,j (following the selected order) do 4: if n > 0∨(ris consistent with the rules in R1 sel,j) then 5: N0←max{number of times that r is applicable to C0 t} 6: if N0>0then 7: R1 sel,j ←R1 sel,j ∪ {< r, N0>} 8: C0 t←C0 t−N0·LHS(r) 9: end if 10: end if 11: end for 12: end for Execution of selected rules 1: for each < r, n >∈R1 sel,j do 2: C0 t←C0 t+n·RHS(r) 3: Update the electrical charges of C0 taccording to RHS(r) 4: end for C. A parallel implementation In this section we provide an implementation to the pseudocode following the imperative programming paradigm. We have included this implementation in the pLinguaCore library in order to evaluate it. The key aspect of this implementation is the provided parallel solution. Considering that in the left-hand side of rules, only one environment is affected, the selection of rules can be executed in parallel for each environment, without compromising the concurrency. However, the selection of rules have to be synchronized, so that the execution phase cannot start before finishing the selection in every environment. First, the set of rules is processed. For each transition and for each environment, a list of applicable rules to the configuration is constructed, according only to the membrane charge in the left-hand side. When these lists are completed, they are used to fill an array (one per environment), called Dj, of triples < rule, probability, number−applications >. At this point, the probabilities of rules are calculated using the associated function, and stored in the triple with the corresponding rule. The construction of this array, assigning a random order, corresponds to the set of rules Djdefined in the pseudocode. Furthermore, a boolean array B, with one entry per membrane, is used to check the consistency of rules in a efficient way. In the pseudocode, the consistency condition is assured by testing each rule with the rest of the previously selected. In this implementation, each time that a rule r=u[v]α h→ u0[v0]α0 his selected a number of times greater than 0, the value in Bfor the active membrane his set to true (B[h]←true). Moreover, the charge is changed in C0 t. Then, if an inconsistent rule is checked, it can be seen by Bthat a rule has reserved the change of the membrane charge, that is stored in C0 t. The first selection phase works directly with the array Dj, so it is a loop using a shuffle iterator (an index that takes random values between 0and maxk). Like in the pseudocode, after checking the consistency condition, the number of applications is calculated. If it is non-zero, the binomial distribution is applied, getting a number n. When the rule is processed, nis stored in the corresponding triple, which is changed by the last one (maxk) in order to do not select again. In the case that the rule is inconsistent or not applicable, it has to be discarded. That is, the rule is changed by the last one (maxk) as before, but also changed with the last selected rule (maxDj), avoiding to be processed in the second phase. At the end of the first selection phase, Djcorresponds to the consistent applicable rules, both with zero and nonzero selected numbers (Rsel). The second selection phase is similar to the one depicted in the pseudocode. It iterates the rules in Dj, according to the order given by probabilities, to fix the number of applications of each one by checking the consistency and maximality. The algorithm ends with the execution phase, which adds the right-hand side of rules to C0 t. Next we show the algorithm that simulates a computation of a functional multienvironment probabilistic P system with active membranes of grade (q, m), with q≥1,m≥1, using Ttime units. (Γ,Σ, G, RE,Π,{fr,j :r∈RΠ,1≤j≤ m},{Mij : 0 ≤i≤q−1,1≤j≤m}) The algorithm must receive the following input parameters: •The initial configuration of the system, C0. •The sets of rules REyRΠ. •The values q,myT. •The set of functions {fr,j :r∈RΠ,1≤j≤m}. The following computable functions are considered for random number generation: •Fb(N, p)is a function that returns a random natural number n∈ {0, . . . , N}, according to the binomial distribution B(N, p), where N∈ {0,...,264}and p∈ {R∩[0,1]}. Real numbers are encoded in floating point (float) with a precision of 32 bits. •Fu(N)is a function that returns a random natural number n∈ {0, . . . , N −1}, according to the uniform distribution U(N), where N∈ {1,...,264}. Procedure DNDP (C0,RE,RΠ,q,m,T,K,{fr,j :r∈ RΠ,1≤j≤m}) 1: {LR,E,j, LRΠ,j,h,α : 1 ≤j≤m, 0≤h≤q−1, α ∈ {−,0,+}} ← Initialization(RE,RΠ,q,m) 2: for t←0to T−1do 3: C0 t←Ct 4: for j←1to mdo 5: Throw new thread Selection-execution(j,t,C0 t,Ct, LR,E,j,{LRΠ,j,h,α 0≤h≤q−1, α ∈ {−,0,+}}, q,m,K,{fr,j :r∈RΠ}) 6: end for 7: Barrier-synchronization {All threads wait until everyone reaches this point} 8: Ct+1 ←C0 t 9: print Ct+1 10: end for Thread Selection-execution (j,t,C0 t,Ct,LR,E,j, {LRΠ,j,h,α 0≤h≤q−1, α ∈ {−,0,+}},q,m,K, {fr,j :r∈RΠ}) 1: Dj, maxDj, B ←Initialize-selection-phase(j,t,C0 t,Ct, LR,E,j,{LRΠ,j,h,α 0≤h≤q−1, α ∈ {−,0,+}},q,m, K,{fr,j :r∈RΠ}) 2: Dj, maxDj, B ←Selection-first-phase(j,Dj,maxDj, B,K,C0 t,Ct) 3: Dj, maxDj, B ←Selection-second-phase(j,Dj,maxDj, B,C0 t) 4: Barrier-synchronization { All threads wait until everyone reaches this point} 5: Execution(j,Dj,maxDj,C0 t) Function Selection-first-phase (j,Dj,maxDj,B,K,C0 t, Ct) 1: maxk←maxDj 2: while maxk>0do 3: i←Fu(maxk) 4: < r, p, n > ←Dj[i] 5: if Is-consistent(r,j,Bj,C0 t)then 6: N0←Count-applications(r,j,C0 t) 7: if N0>0then 8: if p= 1 then 9: ni←N0 10: else 11: N←Count-applications(r,j,Ct) 12: ni←Fb(N, p) 13: if ni> N0then 14: ni←N0 15: end if 16: end if 17: if ni>0then 18: Remove-left-hand-rule-objects(r,j,ni,C0 t) 19: Update-charge(r,j,B,C0 t) 20: n←n+ni 21: end if 22: Swap(Dj,i,maxk−1) 23: maxk←maxk−1 24: else 25: Swap(Dj,i,maxk−1) 26: Swap(Dj,maxk−1,maxdj−1) 27: maxk←maxk−1 28: maxDj←maxDj−1 29: end if 30: else 31: Swap(Dj,i,maxk−1) 32: Swap(Dj,maxk−1,maxdj−1) 33: maxk←maxk−1 34: maxDj←maxDj−1 35: end if 36: end while 37: return Dj,maxDj,B Function Selection-second-phase (j,Dj,maxDj,B,C0 t) 1: Sort triples of the form < r, p, n > from and to Dj, from 0to maxDj, in decrement order by p. 2: for i←0to maxDjdo 3: < r, p, n >←Dj[i] 4: if n > 0∨Is-consistent(r,j,Bj,C0 t)then 5: N0←Count-applications(r,j,C0 t) 6: if N0>0then 7: n←n+N0 8: Remove-left-hand-rule-objects(r,j,N0,C0 t) 9: if n= 0 then 10: Update-charge(r,j,B,C0 t) 11: end if 12: end if 13: end if 14: end for 15: return Dj,maxDj,B Procedure Execution (j,Dj,maxDj,C0 t) 1: for i←0to maxDjdo 2: < r, p, n >←Dj[i] 3: if n > 0then 4: Add-left-hand-rule-objects(r,j,n,C0 t) 5: end if 6: end for Function Initialization (RE,RΠ,q,m) 1: for j←1to mdo 2: The set of rules of the form (x)ej p(x,j,j1,...,js) −−−→ (x1)ej1,...,(xs)ejs⊆REis ordered lexicographically with respect of x, x1, . . . , xs. 3: Be LRE,j the list of rules with the previous order. 4: end for 5: for h←0to q−1do 6: for all α∈ {−,0,+}do 7: The set of rules of the form u[v]α h→u0[v0]α0 h⊆ RΠjis ordered lexicographically with respect of u, v, u0, v0. 8: Be LRΠ,j,h,α the list of rules with the previous order. 9: end for 10: end for 11: return {LR,E,j, LRΠ,j,h,α 0≤j≤m−1,0≤h≤ q−1, α ∈ {−,0,+}} Function Initialize-selection-phase (j,t,C0 t,Ct,LR,E,j, {LRΠ,j,h,α 0≤h≤q−1, α ∈ {−,0,+}},q,m,K, {fr,j :r∈RΠ}) 1: Dj←New array of triples < rule − id, float, integer >, representing < rule, probability, selection−number > respectively, of size: Pq−1 h←0maximumα∈{−,0,+}{Length(LRΠ,j,h,α )}+ Length(LRE,j ) 2: maxDj←0 3: for h←0to q−1do 4: α←charge of membrane hin environment j, from Ct 5: for all r∈LRΠ,j,h,α do 6: p←fr,j (t) 7: if p > 0then 8: Dj[maxDj]←< r, p, 0> 9: maxDj←maxDj+ 1 10: end if 11: end for 12: Bj[h]←false 13: end for 14: for all r∈LRE,j do 15: p←p(x,j,j1,...,js)(t) 16: if p > 0then 17: Dj[maxDj]←< r, p, 0> 18: maxDj←maxDj+ 1 19: end if 20: end for 21: return maxDj,Dj,Bj Function Is-consistent (r,j,B,C) 1: check ←true 2: if ris of the form u[v]α h→u0[v0]α0 hthen 3: β←charge of membrane hin Πj∈C 4: if B[h] = true ∧α06=βthen 5: check ←false 6: end if 7: end if 8: return check Function Count-applications (r,j,C) 1: if ris of the form (x)ej p(x,j,j1,...,js) −−−→ (x1)ej1,...,(xs)ejs then 2: n←multiplicity of object xin environment ej∈C 3: else if ris of the form u[v]α h→u0[v0]α0 hthen 4: n←minimum value such that unappears in the multiset of the upper compartment of hin Πj∈C and vnappears in the multiset of hin Πj∈C 5: end if 6: return n Procedure Remove-left-hand-rule-objects (r,j,n,C) 1: if ris of the form (x)ej p(x,j,j1,...,js) −−−→ (x1)ej1,...,(xs)ejs then 2: Remove the multiset xnfrom environment ej∈C 3: else if ris of the form u[v]α h→u0[v0]α0 hthen 4: Remove the multiset unfrom the multiset of the parent membrane of hin Πj∈C 5: Remove the multiset vnfrom the multiset of membrane hin Πj∈C 6: end if Procedure Update-charge (r,j,C) 1: if ris of the form u[v]α h→u0[v0]α0 hthen 2: Assign α0to the electrical charge of membrane hin Πj∈C 3: end if Procedure Swap (D,i,j) 1: aux ←D[i] 2: D[i]←D[j] 3: D[j]←aux Procedure Add-left-hand-rule-objects (r,j,n,C) 1: if ris of the form (x)ej p(x,j,j1,...,js) −−−→ (x1)ej1,...,(xs)ejs then 2: Add the multisets xn 1, . . . , xn sto the environments ej1, . . . , ejs∈C, respectively. 3: else if ris of the form u[v]α h→u0[v0]α0 hthen 4: Add the multiset u0nto the multiset of the parent membrane of hin Πj∈C 5: Add the multiset v0nfrom the multiset of membrane h in Πj∈C 6: end if V. PERFORMANCE EVALUATION In order to test the performance of our parallel implementation, we have created four simple P systems with the same skeleton, that has no biological meaning. The four test P systems are multienvironment functionalprobabilistic P systems with active membranes of degrees (2,2),(4,2),(8,2) and (16,2) (2 membranes per P system, and only varying the number of environments per each) with electrical charges (0,+,−), of the following form: Π=(G, Γ, µ, R, T, {fr:r∈R},M1,M2) , where: •Gis a empty graph because RE=φ. •Γ = {a, b, x, y, z, r, fi, gi: 1 ≤i≤1000} •µ= [[]2]1is the membrane structure. –M1={a232 }∪{fi, gi: 1 ≤i≤1000} –M2={b232 } •The rules Rto apply are: r1≡a, fi[b]0 2 0.8 −−−→ a, fi[b, x]0 2,1≤i≤1000 r2≡a, fi[b]0 2 0.2 −−−→ a, fi[b, y]0 2,1≤i≤1000 r3≡a, gi[b]0 2 0.9 −−−→ a, gi[b, z]0 2,1≤i≤1000 r4≡a, gi[b]0 2 0.1 −−−→ a, gi[b, r]0 2,1≤i≤1000 These four P systems have been simulated using PLinguaCore with the three available implementations written in Java: binomial block based (binomial in short), DNDP sequential and DNDP parallel algorithms. The simulations have been executed in a computer with an Intel core2 Quad Q9550 system running at 2.83GHz with 4GB of main memory, and using Ubuntu Linux Server 8.04 with the Java Virtual Machine 1.6.0. Figure 1 shows the simulation time for the three implementations simulating the example depicted above. It can be seen that the DNDP sequential code is a little bit faster than the binomial one. As the complexity of the P system increments (in this example, the number of environments), the DNDP reaches better performance. In this experiment, we report up to 1,07x of speedup for the P system with 16 environments. Furthermore, the parallel implementation of the DNDP algorithm has a better overall performance than the other two. Using a 4 cores based computer, the Java Virtual Machine distributes the threads (assigned to each environment of the system) among them. In this way, for the P system with 16 environments, we report a 1,72x of speedup with respect to the DNDP sequential, and 1,84x with the binomial. Binomial DNDP seq DNDP par 0 1 2 3 4 5 6 7 8 9 10 2 env 4 env 8 env 16 env Simulation time (sec) Fig. 1. Simulation time of the three different implementations in a 4 cores based computer VI. CONCLUSIONS AND FUTURE WORK In this paper we have introduced a novel algorithm for multienvironment functional probabilistic P systems with active membranes. The new applications of this model to the field of Population (Ecosystems) Dynamics leads to search for new and efficient algorithms. We have presented the previous algorithm, called binomial block based algorithm, and studied the deficiencies in the simulation. Then, we have based on the idea presented by V. Nguyen et al. [8] in order to improve it. Finally, we have developed an algorithm, called DNDP, performing object dis-