Dynamical Probabilistic P Systems: Definitions and Applications

Dynamical Probabilistic P Systems: Definitions and Applications Dario Pescini1, Daniela Besozzi2, Claudio Zandron1, Giancarlo Mauri1 1Universit`a degli Studi di Milano-Bicocca, DISCo Dipartimento di Informatica, Sistemistica e Comunicazione Via Bicocca degli Arcimboldi 8, 20126 Milano, Italy E-mail: pescini/zandron/[email protected] 2Universit`a degli Studi di Milano, DICo Dipartimento di Informatica e Comunicazione Via Comelico 39, 20135 Milano, Italy E-mail: [email protected] Summary. We introduce dynamical probabilistic P systems, a variant where probabilities associated to the rules change during the evolution of the system, as a new approach to the analysis and simulation of the behavior of complex systems. We define the notions for the analysis of the dynamics and we show some applications for the investigation of the properties of the Brusselator (a simple scheme for the Belousov-Zabothinskii reaction), the Lotka-Volterra system and the decay process. 1 Introduction P systems [4] are a class of distributed parallel computing devices, inspired by the structure and the functioning of cells. The basic model consists of a cell-like membrane structure, composed by several compartments where multisets of objects evolve according to given rules, in a nondeterministic and maximally parallel manner. A computation device is obtained starting from an initial configuration and letting the system evolve. In the following, we assume that the reader is familiar with the basic notions and the terminology underlying P systems. We refer, for details, to [5]. Updated information about P systems can be found at http://psystems.disco.unimib.it/. Many research studies around P systems concentrates on computational power aspects. In this paper, we propose a new approach for the investigation and the application of P systems, which consists in interpreting them as tools for the description and the analysis of the dynamical behavior of complex systems. A similar approach is considered also in [1, 6]. As said, membrane systems are inspired from the functioning of the cell, hence it is natural to consider them for modelling different cellular processes and natural living systems, with the final goal of producing 276 D. Pescini, D. Besozzi, C. Zandron, G. Mauri new tools and acquiring useful information for the scientists (mainly, biologists) working on the modelled system. Some first steps in this direction have already been made, see [3] for various applications. Since we are interested in describing the evolution of a complex system, and since changes of many different conditions can have direct influence on the reaction parameters and behavior, the basic model of P systems is not suitable to describe these kind of processes. For this reason we introduce, in Section 2, dynamical probabilistic P systems, where a probability is associated to each rule and it changes during the whole process (we will talk about evolution instead of computation). In Section 3 we introduce some notions which will then be used to analyze the behavior of such systems. In Section 4 we show some applications to the Brusselator, a simplified theoretical scheme which describes the Belousov-Zhabotinskii reaction (BZ, in short), the Lotka-Volterra and decay processes. In particular, we show we can simulate the behavior of chemical oscillator reactions. Indeed, the interaction of two or more oscillating systems is of interest for many biological processes and systems, as it constitutes an important factor to keep alive an organism or a complex system constituted by several sub-components of different types. Finally, in Section 5 we present the conclusion and give some perspective for future work. 2 Dynamical Probabilistic P Systems In this section we give the definition of a probabilistic P system where the probabilities associated to the rules vary during the evolution of the system. The method for evaluating probabilities and the way the system works are explained in details. Then, we extend the definition to consider families of P systems of this type, whose members differ among each other for the choice of some parameters, but not for the main structure. We assume the reader to be familiar with the basic notions and notations of P systems [5]. Some prerequisites about multisets are here recalled. Let Vbe an alphabet, we denote by V∗the set of all strings over V, by λ the empty string and by V+=V∗\ {λ}the set of non-empty strings. A multiset over Vis a map M:V→N, where M(a) is the multiplicity of any symbol a∈V,Nis the set of natural numbers. A multiset Mover V={a1, . . . , al}can be explicitly represented by the string x=aM(a1) 1aM(a2) 2. . . aM(al) l, for all aisuch that M(ai)6= 0, and by all its possible permutations. By interpreting a multiset in the corresponding form of a string x, we can denote by |x|its length and by |x|athe number of occurrences of a symbol ain x. The set of symbols from V occurring in xis denoted by alph(x). Moreover, to every string x∈V∗we can associate the Parikh vector ΨV(x)=(|x|a1,|x|a2,...,|x|al) (the order of symbols occurring in xmatters). Definition 1. Adynamical probabilistic P system (DPP, in short) of degree nis a construct Π= (V, O, µ, M0, . . . , Mn−1, R0, . . . , Rn−1, I;E),where: Dynamical Probabilistic P Systems: Definitions and Applications 277 •Vis the alphabet of the system, O⊆Vis the set of analyzed symbols; •µis a membrane structure consisting of nmembranes labelled with the numbers 0, . . . , n −1. The skin membrane is labelled with 0; •Mi,i= 0, . . . , n−1, is the multiset over Vinitially present inside membrane i; •Ri,i= 0, . . . , n −1, is a finite set of evolution rules associated with membrane i. An evolution rule is of the form r:uk −→ v, where uis a multiset over V, vis a string over V×({here, out} ∪ {inj|1≤j≤n−1}) and k∈Ris a constant associated to the rule; •I⊆ {0, . . . , n −1} ∪ Eis the set of labels of the analyzed regions; •E={VE, ME, RE}is called the environment, it consists of an alphabet VE⊆ V, a feeding multiset MEover VEand a finite set of feeding rules REof the type r:u→(v, in0), for u, v multisets over VE. The alphabet Oand the set Ispecify which symbols and regions are of peculiar importance in Π, namely those elements whose evolution will be analyzed and simulated. Definition 2. Let Πbe a DPP. We call the parameters of Πthe set Pconsisting of: (1) the multisets M0, . . . , Mn−1, MEinitially present in µand in E, (2) the constants of all rules in R0, . . . , Rn−1. Note that the alphabets V, O, VE, the membrane structure µ, the form of the rules in R0, . . . , Rn−1, REand the set Iof analyzed regions do not belong to the set of parameters of Π. We call these components the main structure of Π. We can now extend Definition 1 and consider a family of DPPs, where the main structure is equal for all members of the family, while the parameters can change from member to member. For instance, one can choose to analyze the same DPP with some different settings of initial conditions, such as different initial multisets and/or different rule constants (this can be useful when not all of them are previously known) and/or different feeding multisets. Definition 3. Afamily of DPPs is defined as F={(Π, P)|Πis a DPP and Pis the set of parameters of Π}. Given two elements (Π1,P1),(Π2,P2)∈ F, it holds Π1=Π2for the main structure and P16=P2for the choice of (all or some) elements in P1and P2. In the following, we will talk about the evolution, not computation, of a DPP, since we are not interested in generating languages but in simulating biological or chemical systems. The family Fdescribes a general model of the biological or chemical system of interest and, for any choice of the parameters, we can investigate the evolution of the corresponding fixed DPP. A fixed initial configuration of Πdepends on the choice of P, hence it consists of the multisets initially present inside the membrane structure, the chosen rule constants and the feeding multiset, which is given as an input to the skin membrane from the environment at each step of the evolution by applying the feeding rules. Different strategies in the feeding process can be used: for instance, one can 278 D. Pescini, D. Besozzi, C. Zandron, G. Mauri use it to keep at a constant value the concentrations of chemicals involved in a certain reaction (see Section 4.1 for an application to the BZ), or to increase the concentrations of substances mimicking the biological transport from the extracellular space. We assume that, as long as the system evolves, the environment contains as many symbols as they are needed to continuously feed the system. At each step of the evolution, all applicable rules are simultaneously applied and all occurrences of the left-hand sides of the rules are consumed, hence the parallelism is maximal at both levels of objects and of rules. For simplicity, in this paper we assume that the system evolves according to a universal clock, that is all membranes and the application of all rules are synchronized. In the future, this condition will be extended to considering also non-synchronized evolutions. The applied rules are chosen according to the probability values dynamically assigned to them; the rules with the highest normalized probability value will be more frequently tossed. In simulations, the tossing process is obtained by means of a random number generator, as described below. If some rules compete for objects and have the same probability values, then objects are nondeterministically assigned to them. The probability associated to each rule in any set Ri,i= 0, . . . , n −1, is a function of its constant and of the current multiset occurring in membrane i, and it is evaluated as follows. Let V={a1, . . . , al},Mibe the multiset inside membrane i,r:uk −→ va rule in Ri; let alph(u) = {a1, . . . , as}and u=aα1 1. . . aαs s. To obtain the actual normalized probability piof applying rwith respect to all other rules that are applicable in membrane iat the same step, we need to evaluate the nonnormalized probability epi(r) of r, which depends on the constant associated to r and on the left-hand side of r, namely: epi(r) = k· s Y h=1 Mi(ah)! αh!(Mi(ah)−αh)! = =k· s Y h=1 Mi(ah)(Mi(ah)−1) . . . (Mi(ah)−αh+ 1) αh!(1) that is epi(r) is dynamically defined, according to the current multiset occurring inside membrane i, since we choose αhcopies of symbols ahamong all its Mi(ah) copies currently available in the membrane itself. If Ri={r1, . . . , rm}, the normalized probability of any rule rjis pi(rj) = epi(rj) Pm j=1 epi(rj).(2) In the simulations, the parallel application of the rules is done by splitting one parallel step into several sequential sub-steps. It is possible to separate each single parallel step into two stages, exploiting the fact that the probability distribution and the applicability of the rules are functions only of the left-hand side of the rules and their constants. In the first stage objects are assigned to rules by means of a Dynamical Probabilistic P Systems: Definitions and Applications 279 random number generator, while in the second one the multiset is updated using a stored trace of the rules previously tossed. It should be pointed out that, during the first stage, the probability distribution of the rules has to be kept constant, otherwise the application of the rules would become sequential. Remark 1. A different probability distribution over rules could be obtained by using the classical rate law of Chemistry, though the approach used in equation (1) is more accurate from the combinatorial point of view [2]. It is well known from Chemistry that the rate of a reaction ρat any time is governed by the concentration of the chemicals involved, namely ρ=k·Y j∈J [Aj]σj, where the index jvaries over all chemicals involved in the reaction, [Aj]represents the concentration of each chemical Ajand kis called the rate constant. The value σjis always experimentally determined but in the elementary reactions, where it is assumed to be the stoichiometric coefficient. Indeed, at high concentrations (multiplicities) the two approaches are undistinguishable, but at lower ones our choice is preferable since it accounts for the exact number of all possible tuples of evolving objects. 3 Analysis of the Dynamics in DPP In this section we introduce some notions that will be used for the analysis of the behavior of a DPP. The final goal is to introduce an appropriate definition of the phase space. Usually, the evolution of a physical system is completely determined by means of the motion equations, a set of differential equations inferred by the system properties. In the case of P systems this role should be accomplished by the evolution rules, which create a one-to-one mapping between the application of each rule and the relative displacement of the system in the phase space. First of all, to keep trace of the system evolution we extend the definition of the alphabet V={a1, . . . , al}of Πby introducing the parameter time, that is we define the space e V:= V×N=V× {time}. Definition 4. Let M={aα1 1, . . . , aαl l}be a multiset over V, where αi≥0 for all h= 1, . . . , l. We call a t-multiset the structure M={aα1 1, . . . , aαl l, t}over the space e V. By abuse of notation, we will denote both the multiset over Vand the t-multiset in e Vwith the same symbol M, being it clear when one considers also the time component or not. To represent a t-multiset in the space e Vwe define its position relatively to the t-multiset O={0, . . . , 0}of e V(the first lcomponents of Oare the null multiplicities of the symbols from V). We need also to extend the notion 280 D. Pescini, D. Besozzi, C. Zandron, G. Mauri of Parikh vector to the space e Vas Ψe V(M) = (α1, . . . , αl, t). This is necessary if we want to distinguish among two multisets having the same total numbers of symbols but different multiplicities for (at least) one symbol from V. Definition 5. The position of a t-multiset Mover e Vis the vector −→ M=Ψe V(M). The vector −→ O=Ψe V(O) is called the origin of e V. From Definition 5 it follows that the positions of t-multisets −→ Oand −→ Mare vectors in the space Nl×N. The next step is to introduce a scalar product in Nl, to naturally define the notion of distance between t-multisets, thus giving the structure of an euclidian space to Nl. Definition 6. Let −→ Mi,−→ Mjbe two positions in Nl×N. The distance between −→ Mi,−→ Mjis a function d:Nl+1 ×Nl+1 −→ R+defined as d2(−→ Mi,−→ Mj) = Pm k=1(αi,k − αj,k)2. Note that the two positions −→ Mi,−→ Mjin Definition 6 need not to be necessarily one the evolution of the other (that is, the multiset inside the same membrane taken into different time steps). In fact, given a family Fof DPP and two positions −→ Mi,−→ Mj, the following cases may hold: 1. −→ Mi,−→ Mjoccur in distinct time steps, in the same membrane of the same DPP with equal setting P; 2. −→ Mi,−→ Mjoccur in distinct or equal time steps, in different membranes of the same DPP with equal setting P; 3. −→ Mi,−→ Mjoccur in distinct or equal time steps, in the same membrane of the same DPP with different settings P1,P2; 4. −→ Mi,−→ Mjoccur in distinct or equal time steps, in different membranes of the same DPP with different settings P1,P2. That is, we might be interested in looking at the multiset occurring inside a membrane during its evolution, or comparing two multisets of different membranes of the same DPP (in equal or different time steps), or else two multisets inside the same (or even a different) membrane but analyzed in two different evolutions of the family of the DPP. In each of the four cases, the distance gives information about “how far” the states in the two trajectories are (that is, the t-multisets in the two evolutions). In particular, given any couple of positions −→ Mi,−→ Mjof the same DPP (for the same or different set of fixed parameters P), we can say that they are simultaneous if they exist at the same time step. This concept can be useful mainly when one considers a membrane structure with degree n≥1, where many multisets are co-evolving. Definition 7. Let −→ Mi,−→ Mjbe two positions in Nl+1. The displacement between −→ Mi,−→ Mjis a function −→ u:Nl+1 ×Nl+1 −→ Zldefined as −→ u(−→ Mi,−→ Mj)=(αi,1− αj,1, . . . , αi,l −αj,l). Dynamical Probabilistic P Systems: Definitions and Applications 281 Note that the displacement can be either a positive or negative value, and it tells how the system “moves”; in details, it tells how the multiplicities in the positions −→ Mjdiffer from those in −→ Mi. Hence, it gives more information than the distance, since it also considers the direction of the variation. Indeed, it is also possible to construct the versor bu:Nl+1 ×Nl+1 −→ Rlof the displacement which only gives the information about the direction of −→ u: bu(−→ Mi,−→ Mj) = Ãαi,1−αj,1 d(−→ Mi,−→ Mj), . . . , αi,l −αj,l d(−→ Mi,−→ Mj)!. Note that −→ u=bu·d, for construction. The last step before arriving to the definition of the phase space consists in defining the velocity, which carries on the information about the time the displacement between two t-multisets (in the same DPP, with equal initial settings) needs to take place. That is, it tells how fast the evolution from one state of the DPP to the other is. Definition 8. Let −→ Mi,−→ Mjbe positions occurring inside the same membrane of a DPP (for a fixed choice of the parameters) in distinct time steps. The velocity is a function −→ v:Nl+1 ×Nl+1 −→ Rldefined as −→ v(−→ Mi,−→ Mj) = µαi,1−αj,1 ti−tj , . . . , αi,l −αj,l ti−tj¶. It should be pointed out here that, actually, this is the definition of the average velocity, which becomes the “instantaneous” velocity when ti−tj= 1, which is the minimal time increment allowed in a discrete-time system, as a DPP is. Note that if −→ Mjis the position evolved from −→ Miin the same membrane, then the instantaneous velocity gives the variation of that multiset in a single time step. We are now ready to define the phase space for a DPP, which is constructed as the cartesian product of the phase spaces of all membranes in the DPP. Let −→ Mi= (α1, . . . , αl, t) be the position of the t-multiset inside membrane iat time t, and let −→ v(−→ Mi,−→ Mi)=(v1, . . . , vl) be its instantaneous velocity (that is, the variation of −→ Miat time twith respect to time t−1, for any t≥1). Definition 9. We call a phase point of −→ Mithe vector −→ ϕti= (α1, . . . , αl, v1, . . . , vl)∈Nl×Rl, for any fixed t∈N. The phase point represents the state of membrane iat any given time t. The evolution of the multiset in membrane ican be described by the phase curve, which is a function −→ ϕi:N−→ Nl×Rlsuch that −→ ϕi(t) = −→ ϕti. The space Nl×Rlis the set of all the points −→ ϕticorresponding to an evolution of the multiset inside any membrane. Definition 10. Let Πbe a DPP of degree n, for some n≥1. The space Φi=Nl× Rlis called the phase space of the membrane i. The space ΦΠ=Φ0×· · ·×Φn−1×ΦE is called the phase space of the DPP. 282 D. Pescini, D. Besozzi, C. Zandron, G. Mauri Hence, the phase space of a DPP describes the evolution of the whole system, with respect to both the change of all multisets and the passing of time. Actually, in analyzing the behavior of a given DPP, we will be interested in considering only the phase space restricted to the regions specified in the set I(see Definition 1). Similarly, only the evolution of symbols from Owill be analyzed for the multisets present in the regions in I. 4 Applications In this section we present some applications of the DPP model to known problems and the relative results obtained from the corresponding simulations. 4.1 The Belousov-Zhabotinskii Reaction The BZ chemical reaction is considered the prototype oscillator and exhibits an extraordinary variety of temporal and spatial phenomena. Its oscillating behavior is one of the most widely studied, both theoretically and experimentally, thus making this reaction a suitable workbench for the capabilities of DPP. Its basic mechanism consists in the oxidation of malonic acid, in acid medium, by bromate ions and catalyzed by cerium, which has two states. The sustained periodic oscillations are observed in the cerium ions. The Brusselator is a simplified theoretical scheme introduced in [7] to explain the nonlinear oscillating behavior, and after that was carefully studied in, e.g., [8]. Despite the fact that it is physically unrealistic, as it involves a trimolecular state, it is recognized to be the skeleton for the explanation of the oscillating behavior in chemical reactions. Moreover, it has a very simple description: Ak1 −→ X, B +Xk2 −→ Y+D, 2X+Yk3 −→ 3X, X k4 −→ E. In this section we describe the Brusselator in terms of DPP and we show the analysis and some results obtained from the simulations. Indeed, in order to describe a chemical or a biological system evolving over time, a kind of rule able to react to the variation of occurrences of symbols (that is, concentrations of substances) is needed. For this purpose, we believe that the dynamical probabilistic rules are really suitable, so we consider the family of DPP that, according to Definition 1, are given by ΠBZ = ({A, B, X, Y },{X, Y },[0]0, M0, R0,0; EBZ), where M0={Am1, Bm2, Xm3, Y m4}and R0={r1:Ak1 −→ X, r2:BX k2 −→ Y, r3:XXY k3 −→ XXX, r4:Xk4 −→ λ}, for some k1, . . . , k4∈R. Note that, with respect to the original equations in the Brusselator, we choose not to consider the chemicals Dand Esince they are not relevant for the system evolution. The environment EBZ is given by the alphabet {A,B}, the multiset MEBZ ={An1, Bn2}, for some n1, n2∈Nand the feeding rules REBZ ={r5:A−→ (A, in0), r6: B−→ (B, in0)}. According to Definition 2, the set of parameters of ΠBZ is Dynamical Probabilistic P Systems: Definitions and Applications 283 PBZ ={m1, . . . , m4, k1, . . . , k4, n1, n2}. A family FBZ can be given by considering different values for the elements in PBZ. The simulations based on the DPP approach have shown all the dynamical behaviors which characterize the continuously stirred BZ (see for example [1, 2, 8]), but here due to space limits we only present the quasi periodic oscillations (in Figure 1, for Pqp BZ ={100,100,1000,2000,50,0.5,5·10−5,5,100,100}) and the attractor (in Figure 2, for Patt BZ ={100,100,1000,2000,1,1,1,1,100,100}). A fading transition from one to the other is possible by tuning the parameters in PBZ . Since in the literature on the Brusselator the phase plane has been widely identified with the X-Yplane, our attention is focused on the dynamic of these symbols. A first characterization of the system dynamic can be obtained by looking directly to the temporal evolution of the two variables: Fig.1.(a) and Fig.2.(a) allow to discriminate the quasi periodic oscillation of the first case from the attracted dynamic of the second one. Fig.1.(b) and Fig.2.(b) show the phase space of membrane 0: in the first case we obtain a limit cycle, in the second case only the initial multiset (point at right-up corner) and the attractor (point at left-bottom corner) can be displayed. Fig.1.(c) and Fig.2.(c) show the evolution of multiplicities of Xand Y; the projection on X−Yplane of these pictures obviously correspond to Fig.1.(b), Fig.2.(b), respectively. Finally, Fig.1.(d) and Fig.2.(d) show the spectra: in the first case, the spectrum shows the highest peak, corresponding to the principal oscillation frequency, and some other harmonics, plus the stochastic contribute which is spread all over the other frequencies; in the second case (where the Yaxis is in logarithmic scale), the spectrum corresponds to a δof Dirac centered in the 0 frequency (the height of δis equal to the mean value of the multiplicities of X and Y), since this is the Fourier transform of a constant (in time) signal. Remark 2. To make clear the definitions of Section 3, we give some examples by extracting three t-multisets from the simulated evolution of (ΠBZ ,Pqp BZ ). Chosen the t-multisets M39 ={100,100,1921,1029,39},M40 = {100,100,2701,262,40},M53 ={100,100,109,1055,53}, their positions are M39 = (100,100,1921,1029,39),M40 = (100,100,2701,262,40),M53 = (100,100,109,1055,53). The distance between M53 and M39 is d(M53,M39) = (0 + 0 + (−1812)2+ 262)1/2≈1812.19, while the displacement is u(M53,M39) = (0,0,−1812,26). The versor associated to this displacement is ˆu(M53,M39) = (0,0,−1812/1812.19,26/1812.19) ≈(0,0,−0.99,0.0014), which says that the predominant direction of the motion is along the Xaxes (that is, the highest variation occurs for the multiplicities of the symbol X). The average velocity v(M53,M39) = (0,0,−1812/14,26/14) ≈(0,0,−129.43,1.86) is quite different from the instantaneous one, which is v(M40,M39) = (0,0,780,−767). 4.2 Decay The radioactive decay process is such that the population of a radioactive isotope varies in time with a fixed rate and usually it is described by the ordinary differential equation dX dt =cX. This process is very simple and can be described by the