Seventh Brainstorming Week on Membrane Computing. Sevilla, February 2-February 6, 2009, Volume II : RGNC REPORT 2/2009

RGNC REPORT 2/2009 Seventh Brainstorming Week on Membrane Computing Sevilla, February 2–February 6, 2009 Volume II Miguel Angel Mart´ınez-del-Amor, Enrique Francisco Orejuela-Pinedo, Gheorghe P˘aun, Ignacio P´erez-Hurtado, Agust´ın Riscos-N´u˜nez Editors Research Group on Natural Computing REPORTS UNIVERSIDAD DE SEVILLA Seventh Brainstorming Week on Membrane Computing Sevilla, February 2–February 6, 2009 Volume II Miguel Angel Mart´ınez-del-Amor, Enrique Francisco Orejuela-Pinedo, Gheorghe P˘aun, Ignacio P´erez-Hurtado, Agust´ın Riscos-N´u˜nez Editors Seventh Brainstorming Week on Membrane Computing Sevilla, February 2–February 6, 2009 Volume II Miguel Angel Mart´ınez-del-Amor, Enrique Francisco Orejuela-Pinedo, Gheorghe P˘aun, Ignacio P´erez-Hurtado, Agust´ın Riscos-N´u˜nez Editors RGNC REPORT 2/2009 Research Group on Natural Computing Sevilla University F´enix Editora, Sevilla, 2009 c Autores I.S.B.N.(Obra completa): 978-84-613-2836-9 I.S.B.N.(Volumen I): 978-84-613-2837-6 I.S.B.N.(Volumen II): 978-84-613-2839-0 Dep´osito Legal: SE-????–09 Edita: F´enix Editora C/ Patricio S´aenz 13, 1–B 41004 Sevilla [email protected] Telf. 954 41 29 91 Printed by Publidisa Preface These proceedings, consisting of two volumes, contain the papers emerged from the Seventh Brainstorming Week on Membrane Computing (BWMC), held in Sevilla, from February 2 to February 6, 2009, in the organization of the Research Group on Natural Computing from the Department of Computer Science and Artificial Intelligence of Sevilla University. The first edition of BWMC was organized at the beginning of February 2003 in Rovira i Virgili University, Tarragona, and the next five editions took place in Sevilla at the beginning of February 2004, 2005, 2006, 2007, and 2008, respectively. In the style of previous meetings in this series, the seventh BWMC was conceived as a period of active interaction among the participants, with the emphasis on exchanging ideas and cooperation; this time, however, there were much more presentations than in the previous years, but still these presentations were “provocative”, mainly proposing new ideas, open problems, research topics, results which need further improvements. The efficiency of this type of meetings was again proved to be very high and the present volumes prove this assertion. As already usual, the number of participants was around 40, most of them computer scientists, but also a few biologists were present. It is important to note that several new names appeared in the membrane computing community, also bringing new view points and research topics to this research area. The papers included in these volumes, arranged in the alphabetic order of the authors, were collected in the form available at a short time after the brainstorming; several of them are still under elaboration. The idea is that the proceedings are a working instrument, part of the interaction started during the stay of authors in Sevilla, meant to make possible a further cooperation, this time having a written support. A selection of the papers from these volumes will be considered for publication in a special issues of International Journal of Computers, Control and Communication. After the first BWMC, a special issue of Natural Computing – volume 2, number 3, 2003, and a special issue of New Generation Computing – volume 22, number 4, 2004, were published; papers from the second BWMC have appeared in vi Preface a special issue of Journal of Universal Computer Science – volume 10, number 5, 2004, as well as in a special issue of Soft Computing – volume 9, number 5, 2005; a selection of papers written during the third BWMC have appeared in a special issue of International Journal of Foundations of Computer Science – volume 17, number 1, 2006); after the fourth BWMC a special issue of Theoretical Computer Science was edited – volume 372, numbers 2-3, 2007; after the fifth edition, a special issue of International Journal of Unconventional Computing was edited – volume 5, number 5, 2009; finally, a selection of papers elaborated during the sixth BWMC has appeared in a special issue of Fundamenta Informaticae – volume 87, number 1, 2008. Other papers elaborated during the seventh BWMC will be submitted to other journals or to suitable conferences. The reader interested in the final version of these papers is advised to check the current bibliography of membrane computing available in the domain website http://ppage.psystems.eu. *** The list of participants as well as their email addresses are given below, with the aim of facilitating the further communication and interaction: 1. Alhazov Artiom, Hiroshima University, Japan, [email protected] 2. Ardelean Ioan, Institute of Biology of the Romanian Academy, Bucharest, Romania, [email protected] 3. Bogdan Aman, A.I.Cuza University, Iasi, Romania, [email protected] 4. Caravagna Giulio, University of Pisa, Italy, [email protected] 5. Ceterchi Rodica, University of Bucharest, Romania, [email protected] 6. Colomer-Cugat M. Angels, University of Lleida, Spain, [email protected] 7. Cord´on-Franco Andr´es, University of Sevilla, Spain, [email protected] 8. Csuhaj-Varj´u Erzs´ebet, Hungarian Academy of Sciences, Budapest, Hungary, [email protected] 9. D´ıaz-Pernil Daniel, University of Sevilla, Spain, [email protected] 10. Frisco Pierluigi, Heriot-Watt University, United Kingdom, [email protected] 11. Garc´ıa-Quismondo Manuel, University of Sevilla, Spain, [email protected] 12. Graciani-D´ıaz Carmen, University of Sevilla, Spain, [email protected] 13. Guti´errez-Naranjo Miguel ´ Angel, University of Sevilla, Spain, [email protected] Preface vii 14. Guti´errez-Escudero Rosa, University of Sevilla, Spain, [email protected] 15. Henley Beverley, University of Sevilla, Spain, [email protected] 16. Ipate Florentin, University of Pitesti, Romania, [email protected] 17. Ishdorj Tseren-Onolt, Abo Akademi, Finland, [email protected] 18. Krassovitskiy Alexander, Rovira i Virgili University, Tarragona, Spain, [email protected] 19. Leporati Alberto, University of Milano-Bicocca, Italy, [email protected] 20. Mart´ınez-del-Amor Miguel Angel, University of Sevilla, Spain, [email protected] 21. Mauri Giancarlo, University of Milano-Bicocca, Italy, [email protected] 22. Mingo Postiglioni Jack Mario, Carlos Tercero University, Madrid, Spain, [email protected] 23. Murphy Niall, NUI Maynooth, Ireland, [email protected] 24. Obtu lowicz Adam, Polish Academy of Sciences, Poland, [email protected] 25. Orejuela-Pinedo Enrique Francisco, University of Sevilla, Spain, [email protected] 26. Pagliarini Roberto, University of Verona, Italy, [email protected] 27. Pan Linqiang, Huazhong University of Science and Technology, Wuhan, Hubei, China, [email protected] 28. P˘aun Gheorghe, Institute of Mathematics of the Romanian Academy, Bucharest, Romania, and University of Sevilla, Spain, [email protected], [email protected] 29. P´erez-Hurtado-de-Mendoza Ignacio, University of Sevilla, Spain, [email protected] 30. P´erez-Jim´enez Mario de Jes´us, University of Sevilla, Spain, [email protected] 31. Porreca Antonio, University of Milano-Bicocca, Italy, [email protected] 32. Riscos-N´u˜nez Agust´ın, University of Sevilla, Spain, [email protected] 33. Rogozhin Yurii, Institute of Mathematics and Computer Science, Chisinau, Moldova, [email protected] 2 T.-O. Ishdorj et al. that they are Turing-complete when considered as number computing devices [7], when used as language generators [4, 2] and also when computing functions [13]. Investigations related to the possibility to solve computationally hard problems by using SN P systems were first proposed in [3]. The idea was to encode the instances of decision problems in a number of spikes, to be placed in an (arbitrarily large) pre-computed system at the beginning of the computation. It was shown that the resulting SN P systems are able to solve the NP-complete problem SAT (the satisfiability of propositional formulas expressed in conjunctive normal form) in a constant time. Slightly different solutions to SAT and 3-SAT by using SN P systems with pre-computed resources were considered in [8]; here the encoding of an instance of the given problem is introduced into the pre-computed resources in a polynomial number of steps, while the truth values are assigned to the Boolean variables of the formula and the satisfiability of the clauses is checked. The answer associated to the instance of the problem is thus computed in a polynomial time. Finally, very simple semi-uniform and uniform solutions to the numerical NP- complete problem Subset Sum — by using SN P systems with exponential size pre-computed resources — have been presented in [9]. All the systems constructed above work in a deterministic way. A different idea of constructing SN P systems for solving NP-complete problems was given in [11, 12], where the Subset Sum and SAT problems were considered. In these papers, the solutions are obtained in a semi-uniform or uniform way by using nondeterministic devices but without pre-computed resources. However, several ingredients are also added to SN P systems such as extended rules, the possibility to have a choice between spiking rules and forgetting rules, etc. An alternative to the constructions of [11, 12] was given in [10], where only standard SN P systems without delaying rules, and having a uniform construction, are used. However, it should be noted that the systems described in [10] either have an exponential size, or their computations last an exponential number of steps. Indeed, it has been proved in [12] that a deterministic SN P system of polynomial size cannot solve an NP-complete problem in a polynomial time unless P=NP. Hence, under the assumption that P6=NP, efficient solutions to NP-complete problems cannot be obtained without introducing features which enhance the efficiency of the system (pre-computed resources, ways to exponentially grow the workspace during the computation, nondeterminism, and so on). The present paper deals with QSAT (the satisfiability of fully quantified propositional formulas expressed in conjunctive normal form) and with Q3SAT (where the clauses that compose the propositional formulas have exactly three literals), two well known PSPACE-complete problems. For QSAT we provide a family {ΠQSAT (2n, m)}n,m∈Nof SN P systems with pre-computed resources such that for all n, m ∈Nthe system ΠQSAT (2n, m) solves all the instances of QSAT which are built using 2nBoolean variables and mclauses. Each system ΠQSAT (2n, m) is deterministic, and computes the solution in a time which is linear with respect to both nand m; however, the size of ΠQSAT (2n, m) is exponential with respect to the size of the instances of the problem. As for Q3SAT, we provide a family Solutions to QSAT and Q3SAT by SN P Systems 3 {ΠQ3SAT (2n)}n∈Nof SN P systems with pre-computed resources, such that for all n∈Nthe system ΠQ3SAT (2n) solves all possible instances of Q3SAT which can be built using 2nBoolean variables. Also in this case the systems ΠQ3SAT (2n) are deterministic and have an exponential size with respect to n. Given an instance of Q3SAT, the corresponding answer is computed in a time which is at most cubic in n. An important observation is that we will not specify how our pre-computed systems could be built. However, we require that such systems have a regular structure, and that they do not contain neither “hidden information” that simplify the solution of specific instances, nor an encoding of all possible solutions (that is, an exponential amount of information that allows to cheat while solving the instances of the problem). These requirements were inspired by open problem Q27 in [15]. Let us note in passing that the regularity of the structure of the system is related to the concept of uniformity, that in some sense measures the difficulty of constructing the systems. Usually, when considering families {C(n)}n∈Nof Boolean circuits, or other computing devices whose number of inputs depends upon an integer parameter n≥1, it is required that for each n∈Na “reasonable” description (see [1] for a discussion on the meaning of the term “reasonable” in this context) of C(n), the circuit of the family which has ninputs, can be produced in polynomial time and logarithmic space (with respect to n) by a deterministic Turing machine whose input is 1n, the unary representation of n. In this paper we will not delve further into the details concerning uniformity; we just rely on reader’s intuition, by stating that it should be possible to build the entire structure of the system in a polynomial time, using only a polynomial amount of information and a controlled replication mechanism, as it already happens in P systems with cell division. We will thus say that our solutions are exp-uniform (instead of uniform), since the systems ΠQSAT (2n, m) and ΠQ3SAT (2n) have an exponential size. The paper is organized as follows. In section 2 we recall the formal definition of SN P systems, as well as some mathematical preliminaries that will be used in the following. In section 3 we provide an exp-uniform family {ΠQSAT (2n, m)}n,m∈Nof SN P systems with pre-computed resources such that for all n, m ∈Nthe system ΠQSAT (2n, m) solves all possible instances of QSAT containing 2nBoolean variables and mclauses. In section 4 we present an exp-uniform family {ΠQ3SAT (2n)}n∈N of SN P systems with pre-computed resources such that for all n∈Nthe system ΠQ3SAT (2n) solves all the instances of Q3SAT which can be built using 2nBoolean variables. Section 5 concludes the paper and suggests some possible directions for future work. 2 Preliminaries We assume the reader to be familiar with formal language theory [16], computational complexity theory [5] as well as membrane computing [15]. We mention here only a few notions and notations which are used throughout the paper. 4 T.-O. Ishdorj et al. For an alphabet V,V∗denotes the set of all finite strings over V, with the empty string denoted by λ. The set of all nonempty strings over Vis denoted by V+. When V={a}is a singleton, then we simply write a∗and a+instead of {a}∗, {a}+. A regular expression over an alphabet Vis defined as follows: (i) λand each a∈Vis a regular expression, (ii) if E1, E2are regular expressions over V, then (E1)(E2), (E1)∪(E2), and (E1)+are regular expressions over V, and (iii) nothing else is a regular expression over V. With each regular expression Ewe associate a language L(E), defined in the following way: (i) L(λ) = {λ}and L(a) = {a}, for all a∈V, (ii) L((E1)∪(E2)) = L(E1)∪L(E2), L((E1)(E2)) = L(E1)L(E2), and L((E1)+) = (L(E1))+, for all regular expressions E1, E2over V. Non-necessary parentheses can be omitted when writing a regular expression, and also (E)+∪{λ} can be written as E∗. For a string str =y1y2. . . y2n, where y2k−1∈ {0,1},y2k∈ {0,1, x2k}, 1 ≤ k≤n, we denote by str|ithe ith symbol of the string str, 1 ≤i≤2n. For given 1≤i≤n, if str such that the 2jth symbol is x2jfor all j≤iand the 2j0th symbol equals to 1 or 0 for all j0≥i, then we denote by str|2i←xa string which is obtained by replacing the 2ith symbol of str with x2i. In particular, for a binary string bin ∈ {0,1}2n,bin|iand bin|2←x2are defined as the ith bit of bin and the string obtained by replacing the second bit of bin with x2, respectively. QSAT is a well known PSPACE-complete decision problem (see for example [5, pages 261–262], where some variants of the problem Quantified Boolean Formulas are defined). It asks whether or not a given fully quantified Boolean formula, expressed in the conjunctive normal form (CNF), evaluates to true or false. Formally, an instance of QSAT with nvariables and mclauses is a formula γn,m =Q1x1Q2x2· · · Qnxn(C1∧C2∧...∧Cm) where each Qi, 1 ≤i≤n, is either ∀or ∃, and each clause Cj, 1 ≤j≤m, is a disjunction of the form Cj=y1∨y2∨. . .∨yrj, with each literal yk, 1 ≤k≤rj, being either a propositional variable xsor its negation ¬xs, 1 ≤s≤n. For example, the propositional formula β=Q1x1Q2x2[(x1∨x2)∧(¬x1∨ ¬x2)] is true when Q1=∀and Q2=∃, whereas it is false when Q1=∃and Q2=∀. The decision problem Q3SAT is defined exactly as QSAT, the only difference being that all the clauses now contain exactly three literals. It is known that even under this restriction the problem remains PSPACE-complete (see, for example, [5, page 262]). In what follows we require that no repetitions of the same literal may occur in any clause. Without loss of generality we can also avoid the clauses in which both the literals xsand ¬xs, for any 1 ≤s≤n, occur. Further, we will focus our attention on the instances of QSAT and Q3SAT in which all the variables having an even index (that is, x2, x4, . . .) are universally quantified, and all the variables with an odd index (x1, x3, . . .) are existentially quantified. We will say that such instances are expressed in normal form. This may be done without loss of generality. In fact, for any instance γn,m =Q1x1Q2x2· · · Qnxn(C1∧C2∧. . . ∧Cm) of QSAT having nvariables and mclauses there exists an equivalent instance γ0 2n,m =∃x0 1∀x0 2. . . ∃x0 2n−1∀x0 2n(C0 1∧C0 2∧. . . ∧C0 m) with 2nvariables, where Solutions to QSAT and Q3SAT by SN P Systems 5 each clause C0 jis obtained from Cjby replacing every variable xiby x0 2i−1if Qi=∃, or by x0 2iif Qi=∀. Note that this transformation may require to introduce some “dummy” variables, that is, variables which are quantified in γ0 2n,m to guarantee the alternance of even-numbered and odd-numbered variables, but that nonetheless do not appear in any clause. For example, for the propositional formula β1=∀x1∃x2[(x1∨x2)∧(¬x1∨¬x2)] the normal form is β0 1=∃x0 1∀x0 2∃x0 3∀x0 4[(x0 2∨ x0 3)∧(¬x0 2∨¬x0 3)]; for the propositional formula β2=∃x1∀x2[(x1∨x2)∧(¬x1∨¬x2)] we have the normal form β0 2=∃x0 1∀x0 2∃x0 3∀x0 4[(x0 1∨x0 4)∧(¬x0 1∨ ¬x0 4)]. The same transformation may be applied on any instance of Q3SAT defined on nBoolean variables; in this case the result will be another instance of Q3SAT, defined on 2n variables. From now on we will denote by QSAT(2n, m) the set of all possible instances of QSAT, expressed in the above normal form, which are built using 2n Boolean variables and mclauses. Similarly, we will denote by Q3SAT (2n) the set of all possible instances of Q3SAT, expressed in normal form, which can be built using 2nBoolean variables. 2.1 Spiking neural P systems As stated in the Introduction, SN P systems have been introduced in [7], in the framework of Membrane Computing. They can be considered as an evolution of P systems, corresponding to a shift from cell-like to neural-like architectures. In SN P systems the cells (also called neurons) are placed in the nodes of a directed graph, called the synapse graph. The contents of each neuron consists of a number of copies of a single object type, called the spike. Every cell may also contain a number of firing and forgetting rules. Firing rules allow a neuron to send information to other neurons in the form of electrical impulses (also called spikes) which are accumulated at the target cell. The applicability of each rule is determined by checking the contents of the neuron against a regular set associated with the rule. In each time unit, if a neuron can use one of its rules, then one of such rules must be used. If two or more rules could be applied, then only one of them is nondeterministically chosen. Thus, the rules are used in the sequential manner in each neuron, but neurons function in parallel with each other. Note that, as usually happens in Membrane Computing, a global clock is assumed, marking the time for the whole system, and hence the functioning of the system is synchronized. When a cell sends out spikes it becomes “closed” (inactive) for a specified period of time, that reflects the refractory period of biological neurons. During this period, the neuron does not accept new inputs and cannot “fire” (that is, emit spikes). Another important feature of biological neurons is that the length of the axon may cause a time delay before a spike arrives at the target. In SN P systems this delay is modeled by associating a delay parameter to each rule which occurs in the system. If no firing rule can be applied in a neuron, then there may be the possibility to apply a forgetting rule, that removes from the neuron a predefined number of spikes. Formally, a spiking neural membrane system (SN P system, for short) of degree m≥1, as defined in [7], is a construct of the form 6 T.-O. Ishdorj et al. Π= (O, σ1, σ2, . . . , σm, syn, in, out), where: 1. O={a}is the singleton alphabet (ais called spike); 2. σ1, σ2, . . . , σmare neurons, of the form σi= (ni, Ri), 1 ≤i≤m, where: a) ni≥0 is the initial number of spikes contained in σi; b) Riis a finite set of rules of the form E/ac→ap;d, where Eis a regular expression over a, and c≥1, p≥0, d≥0, with the restriction c≥p; 3. syn ⊆ {1,2, . . . , m} × {1,2, . . . , m}, with (i, i)6∈ syn for 1 ≤i≤m, is the directed graph of synapses between neurons; 4. in, out ∈ {1,2, . . . , m}indicate the input and the output neurons of Π. A rule E/ac→ap;dwith p≥1 is an extended firing (we also say spiking) rule; a rule E/ac→apwith p= 0 is written in the form E/ac→λand is called an extended forgetting rule. Rules of the types E/ac→a;dand ac→λare said to be standard. If a rule E/ac→ap;dhas E=ac, then we will write it in the simplified form ac→ap;d; similarly, if a rule E/ac→ap;dhas d= 0, then we can simply write it as E/ac→ap; hence, if a rule E/ac→ap;dhas E=acand d= 0, then we can write ac→ap. The rules are applied as follows. If the neuron σicontains kspikes, and ak∈ L(E), k ≥c, then the rule E/ac→ap;dis enabled and can be applied. This means consuming (removing) cspikes (thus only k−cspikes remain in neuron σi); the neuron is fired, and it produces pspikes after dtime units. If d= 0, then the spikes are emitted immediately; if d= 1, then the spikes are emitted in the next step, etc. If the rule is used in step tand d≥1, then in steps t, t + 1, t + 2, . . . , t +d−1 the neuron is closed (this corresponds to the refractory period from neurobiology), so that it cannot receive new spikes (if a neuron has a synapse to a closed neuron and tries to send a spike along it, then that particular spike is lost). In the step t+d, the neuron spikes and becomes open again, so that it can receive spikes (which can be used starting with the step t+d+ 1, when the neuron can again apply rules). Once emitted from neuron σi, the pspikes reach immediately all neurons σjsuch that (i, j)∈syn and which are open, that is, the pspikes are replicated and each target neuron receives pspikes; as stated above, spikes sent to a closed neuron are “lost”, that is, they are removed from the system. In the case of the output neuron, pspikes are also sent to the environment. Of course, if neuron σi has no synapse leaving from it, then the produced spikes are lost. If the rule is a forgetting one of the form E/ac→λ, then, when it is applied, c≥1 spikes are removed. In each time unit, if a neuron σican use one of its rules, then a rule from Ri must be used. Since two firing rules E1/ac1→ap1;d1and E2/ac2→ap2;d2can have L(E1)∩L(E2)6=∅, it is possible that two or more rules can be applied in a neuron; in such a case, only one of them is chosen in a nondeterministic way. However it is assumed that if a firing rule is applicable then no forgetting rule is applicable, and vice versa. Thus, the rules are used in the sequential manner in each Solutions to QSAT and Q3SAT by SN P Systems 7 neuron (at most one in each step), but neurons work in parallel with each other. It is important to note that the applicability of a rule is established depending on the total number of spikes contained in the neuron. The initial configuration of the system is described by the numbers n1, n2, . . . , nmof spikes present in each neuron, with all neurons being open. During the computation, a configuration is described by both the number of spikes present in each neuron and the state of the neuron, that is, the number of steps to count down until it becomes open again (this number is zero if the neuron is already open). Thus, hr1/t1, . . . , rm/tmiis the configuration where neuron σicontains ri≥0 spikes and it will be open after ti≥0 steps, for i= 1,2, . . . , m; with this notation, the initial configuration of the system is C0=hn1/0, . . . , nm/0i. Using the rules as described above, one can define transitions among configurations. Any sequence of transitions starting in the initial configuration is called a computation. A computation halts if it reaches a configuration where all neurons are open and no rule can be used. Since in SN P systems the alphabet contains only one symbol (denoted by a), the input information is sometimes encoded as a sequence of “virtual” symbols, λ or ai,i≥1, where λrepresents no spike and airepresents a multiset of ispikes. The input sequence is then introduced in the input neuron one virtual symbol at one time unite, starting from the leftmost symbol of the sequence. For instance, the sequence a2λa3is composed of three virtual symbols: a2,λand a3. When providing this sequence as input to an SN P system, the virtual symbol a2(that is, two spikes) is introduced at the first computation step, followed by λ(0 spikes) at the next step, and finally by a3(three spikes) at the third step. Another useful extension to the model defined above, already considered in [10, 12, 11, 8], is to use several input neurons, so that the introduction of the encoding of an instance of the problem to be solved can be done in a faster way, introducing parts of the code in parallel in various input neurons. Formally, we can define an SN P system of degree (m, `), with m≥1 and 0 ≤`≤m, just like a standard SN P system of degree m, the only difference being that now there are ` input neurons denoted by in1, . . . , in`. A valid input for an SN P system of degree (m, `) is a set of `binary sequences (where each element of the sequence denotes the presence or the absence of a spike), that collectively encode an instance of a problem. Spiking neural P systems can be used to solve decision problems, both in a semi–uniform and in a uniform way. When solving a problem Qin the semi– uniform setting, for each specified instance Iof Qwe build in a polynomial time (with respect to the size of I) an SN P system ΠQ,I, whose structure and initial configuration depend upon I, that halts (or emits a specified number of spikes in a given interval of time) if and only if Iis a positive instance of Q. On the other hand, a uniform solution of Qconsists of a family {ΠQ(n)}n∈Nof SN P systems such that, when having an instance I ∈ Q of size n, we introduce a polynomial (in n) number of spikes in a designated (set of) input neuron(s) of ΠQ(n) and the computation halts (or, alternatively, a specified number of spikes is emitted in a 8 T.-O. Ishdorj et al. given interval of time) if and only if Iis a positive instance. The preference for uniform solutions over semi–uniform ones is given by the fact that they are more strictly related to the structure of the problem, rather than to specific instances. Indeed, in the semi–uniform setting we do not even need any input neuron, as the instance of the problem is embedded into the structure (number of spikes, graph of neurons and synapses, rules) from the very beginning. If the instances of a problem Qdepend upon two parameters (as is the case of QSAT, where nis the number of variables and mthe number of clauses in a given formula), then we will denote the family of SN P systems that solves Qby {ΠQ(n, m)}n,m∈N. Alternatively, if one does not want to make the family of SN P systems depend upon two parameters, it is possible to define it as {ΠQ(hn, mi)}n,m∈N, where hn, mi indicates the positive integer number obtained by applying an appropriate bijection (for example, Cantor’s pairing) from N2to N. In the above definitions it is assumed that the uniform (resp., semi-uniform) construction of ΠQ(n) (resp., ΠQ,I) is performed by using a deterministic Turing machine, working in a polynomial time. As stated in the Introduction, the SN P systems we will describe will solve all the instances of QSAT and Q3SAT of a given size, just like in the uniform setting. However, such systems will have an exponential size. Since a deterministic Turing machine cannot produce (the description of) an exponential size object in a polynomial time, we will say that our solutions are exp-uniform. 3 An exp-Uniform Solution to QSAT In this section we build an exp-uniform family {ΠQSAT (2n, m)}n,m∈Nof SN P systems such that for all n, m ∈Nthe system ΠQSAT (2n, m) solves all the instances of QSAT(2n, m) in a polynomial number of steps with respect to nand m, in a deterministic way. The instances of QSAT(2n, m) to be given as input to the system ΠQSAT (2n, m) are encoded as sequences of virtual symbols, as follows. For any given instance γ2n,m =∃x1∀x2. . . ∃x2n−1∀x2n(C1∧C2∧. . .∧Cm) of QSAT (2n, m), let code(γ2n,m) = α11α12 · · · α12nα21α22 · · · α22n· · · αm1αm2· · · αm2n, where each αij, for 1 ≤i≤mand 1 ≤j≤2n, is a spike variable whose value is an amount of spikes (a virtual symbol), assigned as follows: αij =     a, if xjoccurs in Ci; a2,if ¬xjoccurs in Ci; λ, otherwise. In this way the sequence αi1αi2· · · αi2nof spike variables represents the clause Ci, and the representation of γ2n,m is just the concatenation of the representations of the single clauses. As an example, the representation of γ2,2= ∃x1∀x2[(x1∨ ¬x2)∧(¬x2)] is aa2λa2. The set of all the encoding sequences of Solutions to QSAT and Q3SAT by SN P Systems 9 all possible instances of QSAT (2n, m) is denoted by code(QSAT(2n, m)). For instance, QSAT(2,1) contains the following nine formulas (the existential and the universal quantifiers are here omitted for the sake of readability): γ1 2,1= no variable appears in the clause, γ2 2,1=x2,γ3 2,1=¬x2,γ4 2,1=x1,γ5 2,1=x1∨x2, γ6 2,1=x1∨ ¬x2,γ7 2,1=¬x1,γ8 2,1=¬x1∨x2,γ9 2,1=¬x1∨ ¬x2. Therefore, code(QSAT(2,1)) = {λλ, λa, λa2, aλ, aa, aa2, a2λ, a2a, a2a2}. The structure of the pre-computed SN P system that solves all possible instances of QSAT(2n, m) is illustrated in a schematic way in Figures 1 and 2. The system is a structure of the form Π(2n,m) QSAT = (ΠQSAT (2n, m), code(QSAT(2n, m))) with: •ΠQSAT (2n, m) = (O, µ, in, out), where: 1. O={a}is the singleton alphabet; 2. µ= (H, Si∈H{mi},Sj∈HRj, syn) is the structure of the SN P system, where: –H=H0∪H1∪H2∪H3is a finite set of neuron labels, with H0={in, out, d}∪{di|0≤i≤2n}, H1={Cxi, Cxi1, Cxi0|1≤i≤2n}, H2={bin, Cbin |bin ∈ {0,1}2n}, H3={y1y2. . . y2n−1y2n|yi∈ {0,1}when iis odd and yi∈ {0,1, xi}when iis even, and there exists at least one k∈ {1,2, . . . , n}such that y2k= x2k}(we recall that even values of icorrespond to universally quantified variables). All the neurons are injectively labeled with elements from H; –md0= 2, md= 1 and mi= 0 (i∈H,i6=d0, d) are the numbers of spikes that occur in the initial configuration of the system; –Rk,k∈H, is a finite set of rules associated with neuron σk, where: Rin ={a→a, a2→a2},Rd={a→a; 2mn +n+ 6}, Rdi={a2→a2}, for i∈ {0,1, . . . , 2n−1}, and Rd2n={a2→a2, a3→λ}, RCxi={a→λ, a2→λ, a3→a3; 2n−i, a4→a4; 2n−i}, for i∈ {1,2, . . . , 2n}, RCxi1={a3→a2, a4→λ}and RCxi0={a3→λ, a4→a2}, for i∈ {1,2, . . . , 2n}, RCbin ={(a2)+/a →a}∪{a2k−1→λ|k= 1,2, . . . , 2n}, for bin ∈ {0,1}2n, Rbin ={am→a}, for bin ∈ {0,1}2n, Rstr ={a2→a}, where str ∈H3and there exists at least one i∈ {1,2, . . . , n}such that str|2i6=x2i, Rstr0={a2→a2}, where str0∈H3and str0|2k=x2k, for all 1 ≤k≤n, Rout ={(a2)+/a →a}; –syn is the set of all the synapses between the neurons. The following synapses are used in the input module (see Figure 1): (in, Cxi), (di−1, di), (di, Cxi), (Cxi, Cxi1) and (Cxi, Cxi0), for all 1 ≤i≤2n, as well as (d2n, d1) and (d, d2n). The synapses connecting the other neurons are illustrated in Figure 2: 10 T.-O. Ishdorj et al. (Cxij, Cbin), where bin ∈ {0,1}2nand bin|i=j, 1≤i≤2n,j∈ {0,1}, (Cbin, bin), where bin ∈ {0,1}2n, (bin, str), where bin ∈ {0,1}2n,str ∈H3, and str = (bin|i←x2), (strj1, strj2), where strj1, strj2∈H3and strj2= (strj1|2i←x2i), 2 ≤i≤ n, (str, out), where str ∈H3and str|2k=x2k, for all 1 ≤k≤n; 3. in, out indicate the input and output neurons, respectively; •code(QSAT(2n, m)) is the set of all the encoding sequences for all the possible instances of QSAT(2n, m), as defined above. Fig. 1. The input module of ΠQSAT (2n, m) The system is composed of four modules: input, satisfiability checking, quantifier checking, and output. To simplify the description of the system and its working, the neurons in the system are arranged in n+ 7 layers in Figures 1 and 2. The input module has three layers (the first layer contains three neurons σd0,σdand σin; the second layer contains 2nneurons σdi, 1 ≤i≤2n; the third layer contains 2nneurons σCxi, 1 ≤i≤2n). The satisfiability checking module has also three layers (the fourth layer contains 4nneurons σCxij, 1 ≤i≤2n,j= 0,1; the fifth layer contains 22nneurons σCbin,bin ∈ {0,1}2n; the sixth layer contains 22nneu- rons σbin,bin ∈ {0,1}2n). The quantifier checking module is composed of nlayers, from the 7th to the (n+ 6)th layer, where a total of 22n−2nneurons are used. Solutions to QSAT and Q3SAT by SN P Systems 11 Input module a 3ƀ a 2 a 4ƀƮ Cx 1 1 a 3ƀƮ a 4ƀ a 2 Cx 1 0 Cx 2 1 Cx 2 0 Cx 2n 1 Cx 2n 0 Ǝ a 2Ə   a ƀ a a 2k  1ƀƮ C 11111 ſ 1 k  1, ſ , 2 n C 10111 ſ 1 C 11101 ſ 1 C 10101 ſ 1 C 01 ſ 0 C 0 ſ 00 a mƀ a 11111 ſ 1 a mƀ a 10111 ſ 1 a mƀ a 11101 ſ 1 a mƀ a 10101 ſ 1 a mƀ a 01 ſ 0 a mƀ a 0 ſ 00 a 2ƀ a 1 x 2 1 ſ 1 a 2ƀ a 1 x 2 101 ſ 1 a 2ƀ a 0 x 2 0 ſ 0 a 2ƀ a 1 x 2 1 x 4 1 ſ 1 a 2ƀ a 2 1 x 2 1 x 4ſ 1 x 2n a 2ƀ a 2 1 x 2 0 x 4 1 x 6ſ 1 x 2n a 2ƀ a 2 0 x 2ſ 0 x 2n out  a 2    a  a ¹ ¹ ¹ ¸ ¸ ¸¸ ¸¸ a 3ƀ a 2 a 4ƀƮ a 3ƀƮ a 4ƀ a 2 a 3ƀƮ a 4ƀ a 2 a 3ƀ a 2 a 4ƀƮ Ǝ a 2Ə   a ƀ a a 2k  1ƀƮ k  1, ſ , 2 n Ǝ a 2Ə   a ƀ a a 2k  1ƀƮ k  1, ſ , 2 n Ǝ a 2Ə   a ƀ a a 2k  1ƀƮ k  1, ſ , 2 n Ǝ a 2Ə   a ƀ a a 2k  1ƀƮ k  1, ſ , 2 n Ǝ a 2Ə   a ƀ a a 2k  1ƀƮ k  1, ſ , 2 n ¹     Fig. 2. Structure of the SN P system ΠQSAT (2n, m) The output module only contains one neuron σout, which appears in the last layer. In what follows we provide a more detailed description of each module, as well as its working when solving a given instance γ2n,m ∈QSAT(2n, m). •Input: The input module consists of 4n+ 3 neurons, contained in the layers 1 – 3 as illustrated in Figure 1; σin is the unique input neuron. The values of the spike variables of the encoding sequence code(γ2n,m) are introduced into σin 18 T.-O. Ishdorj et al. Ź 1  0,2  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0 ź Ź 0  11,0  0,2  0,0  0,1  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0 ź Ź 0  10,0  0,0  0,2  0,2  0,3  0,1  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0 ź Ź 0  9,0  0,2  0,0  0,1  0,0  1,4  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0 ź Ź 0  8,0  0,0  0,2  0,0  0,3  0,1  0,3  0,3  0,4  0,4  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0 ź Ź 0  7,0  0,2  0,0  0,0  0,0  1,2  0,0  0,0  0,0  0,0  0,2  0,4  0,0  0,2  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0 ź Ź 0  6,0  0,0  0,2  0,0  0,2  0,0  0,3  0,3  0,0  0,0  0,1  0,3  0,0  0,1  0,1  0,1  0,0  0,1  0,0  0,0  0,0  0 ź Ź 0  5,0  0,2  0,0  0,0  0,0  0,2  0,0  0,0  0,0  0,0  0,2  0,2  0,0  0,0  0,1  0,1  0,0  0,1  0,0  0,0  0,0  0 ź Ź 0  4,0  0,0  0,2  0,0  0,2  0,0  0,0  0,0  0,0  0,0  0,1  0,1  0,0  0,0  0,2  0,2  0,0  0,1  0,0  0,0  0,0  0 ź Ź 0  3,0  0,2  0,0  0,0  0,0  0,2  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,1  0,2  0,0  0,0  0 ź Ź 0  2,0  0,0  0,2  0,0  0,2  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,1  0,0  0,0  0,2  0 ź Ź 0  1,0  0,2  0,1  0,0  0,0  0,2  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,1  0,0  0,0  0,1  0 ź Ź 0  0,0  0,0  0,3  0,0  0,2  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,1  0,0  0,0  0,1  0 ź Ź 0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,0  0,1  0,0  0,0  0,1  0 ź r d , r d0 r d1 , r in,2 r d2,1 , r in,1 , r Cx1,3 , r Cx2,1 r d1 , r in, 2 , r Cx1,s , r Cx2,4 r d2,1 , r Cx1,3 , r Cx2,1 , r Cx11,1 r Cx10,1 , r Cx21,2 , r Cx20,2 r d1 , r Cx1,s , r Cx2,2 r C11, 1 , r C10 ,1 , r C00 ,1 r d2,1 , r Cx1,2 , r Cx11,1 r Cx10,1 , r C11, 3 , r C10 , 2 , r C00 , 3 r d1 , r Cx2,2 , r C11, 1 , r C10 , 1 r d2,1 , r Cx1,2 r C11,3 , r C10 , 3 , r 11 , r 10 r d1 , r Cx2,2 , r 1x2 r d ,11s , r d1 , r Cx2,2 r d2,2 , r Cx1,2 r d2,1 , r Cx1,2 , r out Fig. 4. The transition diagram of the system illustrated in Figure 3 Solutions to QSAT and Q3SAT by SN P Systems 19 Fig. 5. Sketch of a deterministic SN P system that solves all possible instances of Q3SAT (2n) 20 T.-O. Ishdorj et al. Figure 5 outlines an SN P system which can be used to solve any instance γ2n,m of Q3SAT(2n). The input to this system is once again the instance of Q3SAT we want to solve, but this time such an instance is given by specifying — among all the possible clauses that can be built using nBoolean variables — which clauses occur in the instance. The selection is performed by putting (in parallel, in the initial configuration of the system) one spike in each of the input neurons sel1, sel2, . . . , selCl(2n)that correspond to the selected clauses. The system is a simple modification of the one used in [8] to uniformly solve all the instances of the NP-complete problem 3-SAT of a given size. To see how it works, let us consider the family {M(2n)}n∈Nof Boolean matrices, where M(2n)has 22nrows — one for each possible assignment to the variables x1, x2, . . . , x2n— and one column for each possible 3-clause that can be built using the same variables. As stated above, the number of columns is Cl(2n)∈Θ(n3), a polynomial quantity in n. In order to make the construction of the matrix M(2n)as regular as possible, we could choose to list all the 3-clauses in a predefined order; however, our result is independent of any such particular ordering, and hence we will not bother further with this detail. For every j∈ {1,2,3, . . . , 22n}and i∈ {1,2, . . . , Cl(2n)}, the x1x2x3x4··· x1∨x2∨ ¬x4· · · ¬x1∨ ¬x2∨x3··· 0 0 0 0 · · · 1· · · 1··· 0 0 0 1 ··· 0··· 1··· 0 0 1 0 ··· 1··· 1··· 0 0 1 1 ··· 0··· 1··· 0 1 0 0 ··· 1··· 1··· 0 1 0 1 ··· 1··· 1··· 0 1 1 0 ··· 1··· 1··· 0 1 1 1 ··· 1··· 1··· 1 0 0 0 ··· 1··· 1··· 1 0 0 1 ··· 1··· 1··· 1 0 1 0 ··· 1··· 1··· 1 0 1 1 ··· 1··· 1··· 1 1 0 0 ··· 1··· 0··· 1 1 0 1 ··· 1··· 0··· 1 1 1 0 ··· 1··· 1··· 1 1 1 1 ··· 1··· 1··· Assignments Clauses Fig. 6. An excerpt of matrix M(4). On the left we can see the assignments which are associated to the corresponding rows of the matrix. Only the columns corresponding to the clauses x1∨x2∨ ¬x4and ¬x1∨ ¬x2∨x3are detailed element M(2n) ji is equal to 1 if and only if the assignment associated with row j satisfies the clause associated with column i. Figure 6 shows an excerpt of matrix M(4), where each row has been labelled with the corresponding clause; only the Solutions to QSAT and Q3SAT by SN P Systems 21 Check Satisfiability(M(2n)) res ←[1 1 · · · 1] // 22nelements for all columns Cin M(2n) do if Ccorresponds to a selected clause then res ←res ∧C// bit-wise and return res Fig. 7. Pseudocode of the algorithm used to select the assignments that satisfy all the clauses of γ2n,m ∈Q3SAT (2n) columns that correspond to clauses x1∨x2∨ ¬x4and ¬x1∨ ¬x2∨x3are shown in details. Let us now consider the algorithm given in pseudocode in Figure 7. The variable res is a vector of length 22n, whose components — which are initialized to 1 — are bijectively associated with all the possible assignments of x1, x2, . . . , x2n. The components of res are treated as flags: when a component is equal to 1, it indicates that the corresponding assignment satisfies all the clauses which have been examined so far. Initially we assume that all the flags are 1, since we do not have examined any clause yet. The algorithm then considers all the columns of M(2n), one by one. If the column under consideration does not correspond to a selected clause, then it is simply ignored. If, on the other hand, it corresponds to a clause which has been selected as part of the instance, then the components of res are updated, putting to 0 those flags that correspond to the assignments which do not satisfy the clause. At the end of this operation, which can be performed in parallel on all the components, only those assignments that satisfy all the clauses previously examined, as well as the clause currently under consideration, survive the filtering process. After the last column of M(2n)has been processed, we have that the components of vector res indicate those assignments that satisfy all the clauses of the instance γ2n,m of Q3SAT given as input. Stated otherwise, res is the output column of the truth table of the unquantified propositional formula contained in γ2n,m. This algorithm can be easily transformed into an exponential size Boolean circuit, that mimics the operations performed on the matrix M(2n), described by the pseudocode given in Figure 7. Such a circuit can then be easily simulated using the SN P system that we have outlined in Figure 5 (precisely, the left side of the system, until the column of neurons that contain the final value of vector res). This part of the system is composed of three layers for each possible 3-clause that can be built using 2nBoolean variables. Two of these layers are used to store the intermediate values of vector res and the values contained in the columns of M(2n), respectively. The third layer, represented by the boxes marked with Filter in Figure 5, transforms the current value of res to the value obtained by applying 22 T.-O. Ishdorj et al. the corresponding iteration of the algorithm given in Figure 7. This layer is in turn composed by three layers of neurons, as we will see in a moment. The last part of the system is used to check the satisfiability of the universal and existential quantifiers associated with the variables x1, x2, . . . , x2n. The neurons in this part of the system compose a binary tree of depth 2n; the first layer of neurons corresponds to the bottom of the tree, and checks the satisfiability of the quantifier ∀x2n; the second layer checks the satisfiability of ∃x2n−1and so on, until the last layer, whose only neuron is σout (the output neuron), that checks the satisfiability of ∃x1. To see how the check is performed, let us consider the fully quantified Fig. 8. Example of a quantified Boolean formula formed by one clause, built using two Boolean variables. On the left, its truth table is reported with an indication of the truth assignments that make the formula true. On the right, the tree which is used to check the satisfaction of the quantifiers ∀and ∃is depicted formula ∃x1∀x2(¬x1∨x2). This formula is composed of a single 2-clause (hence it is not a valid instance of Q3SAT), built using two Boolean variables. In Figure 8 we can see the truth table of the clause, and an AND/OR Boolean circuit that can be used to check whether the quantifiers associated with x1and x2are satisfied. This circuit is a binary tree whose nodes are either AND or OR gates. Each layer of nodes is associated with a Boolean variable: precisely, the output layer is associated to x1, the next layer to x2, and so on until the input layer, which is associated to xn. If Qi=∀then the nodes in the layer associated with xiare AND gates; on the contrary, if Qi=∃then the nodes in such a layer are OR gates. The input lines of the circuit are bijectively associated to the set of all possible truth assignments. It is not difficult to see that when these input lines are fed with the values contained in the output column of clause’s truth table, the output of the circuit is 1 if and only if the fully quantified formula is true. Since in the first part of the system we have computed the output column of the truth table of the unquantified propositional formula contained in γ2n,m, we just have to feed these values to an SN P system that simulates the above AND/OR Boolean circuit to see whether γ2n,m is true or not. Implementing such a Boolean circuit by means of an SN P system is trivial: to see how this can be done, just compare the last two layers of the circuits illustrated in Figures 10 and 11. Solutions to QSAT and Q3SAT by SN P Systems 23 The overall system then works as follows. During the computation, spikes move from the leftmost to the rightmost layer. One spike is expelled to the environment by neuron σout if and only if γ2n,m is true. In the initial configuration, every neuron in the first layer (which is bijectively associated with one of the 22npossible assignments to the Boolean variables x1, x2, . . . , x2n) contains one spike, whereas neurons sel1, sel2, . . . , selCl(2n)contain one or zero spikes, depending upon whether or not the corresponding clause is part of the instance γ2n,m given as input. Stated otherwise, the user must provide one spike — in the initial configuration of the system — to every input neuron selithat corresponds to a clause that has to be selected. In order to deliver these spikes at the correct moment to all the filters that correspond to the ith iteration of the algorithm, every neuron selicontains the rule a→a; 4(i−1), whose delay is proportional to i. In order to synchronize the execution of the system, also the neurons that correspond to the ith column of M(2n)deliver their spikes simultaneously with those distributed by neurons seli, using the same rules. An alternative possibility is to provide the input to the system in a sequential way, for example as a bit string of length Cl(2n), where a 1 (resp., 0) in a given position indicates that the corresponding clause has to be selected (resp., ignored). In this case we should use a sort of delaying subsystem, that delivers — every four time steps — the received spike to all the neurons that correspond to the column of M(2n)currently under consideration. Since the execution time of our algorithm is proportional to the number Cl(2n) of all possible clauses containing 2nBoolean variables, this modification keeps the computation time of the entire system cubic with respect to n. In the first computation step, all the inputs going into the first layer of filters are ready to be processed. As the name suggests, these filters put to 0 those flags which correspond to the assignments that do not satisfy the first clause (corresponding to the first column of M(2n)). This occurs only if the clause has been selected as part of the instance γ2n,m ∈Q3SAT(2n) given as input, otherwise all the flags are kept unchanged, ready to be processed by the next layer of filters. In either case, when the resulting flags have been computed they enter into the second layer of filters together with the values of the second column of M(2n), and the input sel2 that indicates whether this column is selected or not as being part of the instance. The computation proceeds in this way until all the columns of M(2n)have been considered, and the resulting flags (corresponding to the final value of vector res in the pseudocode of Figure 7) have been computed. Before looking at how the system checks the satisfiability of the universal and existential quantifiers ∃x1,∀x2, . . . , ∀x2n, let us describe in detail how the filtering process works. This process is performed in parallel on all the flags: if the clause Cihas been selected then an and is performed between the value M(2n) ji (that indicates whether the jth assignment satisfies Ci) and the current value of the flag resj(the jth component of res); as a result, resjis 1 if and only if the jth assignment satisfies all the selected clauses which have been considered up to now. On the other hand, if the clause Cihas not been selected then the old value of resjis kept unaltered. This filtering process can be summarized by the pseudocode 24 T.-O. Ishdorj et al. Filter(seli, resj, Ci) if seli= 0 then return resj else return resj∧Ci Fig. 9. Pseudocode of the Boolean function computed by the blocks marked with Filter in Figure 5 given in Figure 9, which is equivalent to the following Boolean function: (¬seli∧resj)∨(seli∧resj∧Ci) Such a function can be computed by the Boolean circuit depicted in Figure 10, that in turn can be simulated by the SN P system illustrated in Figure 11. Note the system represented in this latter figure is a generic module which is used many times in the whole system outlined in Figure 5, hence we have not indicated the delays which are needed in neurons seliand Ci. Also neuron 1, which is used to negate the value emitted by neuron seli, must be activated together with seli, that is, after 4(i−1) steps after the beginning of the computation. The spike it contains can be reused in the namesake neuron that occurs in the next layer of filters. Fig. 10. The Boolean circuit that computes the function Filter whose pseudocode is given in Figure 9 The last part of the system illustrated in Figure 5 is devoted to check the satisfiability of the universal and existential quantifiers ∃x1,∀x2, . . . , ∀x2nassociated to the Boolean variables x1, x2, . . . , x2n. As we have seen, the final values of vector res represent the output column of the truth table of the unquantified propositional formula contained in γ2n,m. Hence, to check whether all the universal and existential quantifiers are satisfied, it suffices to feed these values as input to an SN P system that simulates a depth 2nAND/OR Boolean circuit that operates as described in Figure 8. Each gate is simply realized as a neuron that contains Solutions to QSAT and Q3SAT by SN P Systems 25 Fig. 11. An SN P system that computes the function Filter given in Figure 9, simulating the Boolean circuit of Figure 10 two rules, as depicted in the last two layers of Figure 11. At each computation step, one quantifier is checked; when the check terminates, one spike is emitted to the environment by the output neuron σout if and only if the fully quantified formula γ2n,m ∈Q3SAT(2n) given as input to the entire system is true. The total computation time of the system is proportional to the number Cl(2n) of columns of M(2n), that is, Θ(n3). As we can see, the structure of the system that uniformly solves all the instances of Q3SAT(2n) is very regular, and does not contain “hidden information”. For the sake of regularity we have also omitted some possible optimizations, that we briefly mention here. The first column of neurons in Figure 5 corresponds to the initial value of vector res in the pseudocode given in Figure 7. Since this value is fixed, we can pre-compute part of the result of the first step of computation, and remove the entire column of neurons from the system. In a similar way we can also remove the subsequent columns that correspond to the intermediate values of res, and send these values directly to the next filtering layer. A further optimization concerns the values M(2n) ji , which are contained in the neurons labelled with Ci. Since these values are given as input to and gates, when they are equal to 1 they can be removed since they do not affect the result; on the other hand, when they are equal to 0 also the result is 0, and thus we can remove the entire and gate. 5 Conclusions and Remarks In this paper we have shown that QSAT, a well known PSPACE-complete problem, can be deterministically solved in linear time with respect to the number nof variables and the number mof clauses by an exp-uniform family of SN P systems 26 T.-O. Ishdorj et al. with pre-computed resources. We have also considered Q3SAT, a restricted (but still PSPACE-complete) version of QSAT in which all the clauses of the instances have exactly three literals; we have shown that in this case the problem can be solved in a time which is at most cubic in n, independent of m. Each pre-computed SN P system of the family can be used to solve all the instances of QSAT (or Q3SAT), expressed in a normalized form, of a given size. Note that using pre-computed resources in spiking neural P systems is a powerful technique, that simplifies the solution of computationally hard problems. The pre-computed SN P systems presented in this paper have an exponential size with respect to nbut, on the other hand, have a regular structure. It still remains open whether such pre-computed resources can be constructed in a regular way by using appropriate computation devices that, for example, use a sort of controlled duplication mechanism to produce an exponential size structure in a polynomial number of steps. A related interesting problem is to consider whether alternative features can be introduced in SN P systems to uniformly solve PSPACE-complete problems. Nondeterminism is the first feature that comes to mind, but this is usually considered too powerful in the Theory of Computation. Acknowledgements Valuable comments and suggestions from professor Mario J. P´erez-Jim´enez are greatly appreciated. This work was supported by National Natural Science Foundation of China (Grant Nos. 60674106, 30870826, 60703047, and 60803113), Program for New Century Excellent Talents in University (NCET-05-0612), Ph.D. Programs Foundation of Ministry of Education of China (20060487014), Chenguang Program of Wuhan (200750731262), HUST-SRF (2007Z015A), and Natural Science Foundation of Hubei Province (2008CDB113 and 2008CDB180). The first author was partially supported by Academy of Finland, project 122426. The second author was partially supported by MIUR project “Mathematical aspects and emerging applications of automata and formal languages” (2007). References 1. J.L. Balc´azar, J. D´ıaz, J. Gabarr´o: Structural Complexity, Voll. I and II, Springer- Verlag, Berlin, 1988–1990. 2. H. Chen, R. Freund, M. Ionescu, Gh. P˘aun, M.J. P´erez-Jim´enez: On string languages generated by spiking neural P systems. Fundamenta Informaticae, 75 (2007), 141– 162. 3. H. Chen, M. Ionescu, T.-O. Ishdorj, On the efficiency of spiking neural P systems. Fourth Brainstorming Week on Membrane Computing (M.A. Guti´errez-Naranjo, Gh, P˘aun, A. Riscos-N´u˜nez, F.J. Romero-Campero, eds.), Sevilla, January 30– February 3, Sevilla University, F´enix Editora, 2006, 195–206. 4. H. Chen, M. Ionescu, T.-O. Ishdorj, A. P˘aun, Gh. P˘aun, M.J. P´erez-Jim´enez: Spiking neural P systems with extended rules: universality and languages. Natural Computing, 7, 2 (2008), 147–166. Solutions to QSAT and Q3SAT by SN P Systems 27 5. M.R. Garey, D.S. Johnson: Computers and Intractability. A Guide to the Theory of NP-Completeness. W.H. Freeman and Company, San Francisco, 1979. 6. M.A. Guti´errez-Naranjo, M.J. P´erez-Jim´enez, D. Ram´ırez-Mart´ınez: A software tool for verification of spiking neural P systems. Natural Computing, 7, 4 (2008), 485–497. 7. M. Ionescu, Gh. P˘aun, T. Yokomori: Spiking neural P systems. Fundamenta Informaticae, 71, 2–3 (2006), 279–308. 8. T.-O. Ishdorj, A. Leporati: Uniform solutions to SAT and 3-SAT by spiking neural P systems with pre-computed resources. Natural Computing, 7, 4 (2008), 519–534. 9. A. Leporati, M.A. Guti´errez-Naranjo: Solving Subset Sum by spiking neural P systems with pre-computed resources. Fundamenta Informaticae, 87, 1 (2008), 61–77. 10. A. Leporati, G. Mauri, C. Zandron, Gh. P˘aun, M.J. P´erez-Jim´enez: Uniform solutions to SAT and Subset Sum by spiking neural P systems. Natural Computing, online version (DOI: 10.1007/s11047-008-9091-y). 11. A. Leporati, C. Zandron, C. Ferretti, G. Mauri: Solving numerical NP-complete problem with spiking neural P systems. Membrane Computing, International Workshop, WMC8 (G. Eleftherakis, P. Kefalas, Gh. P˘aun, G. Rozenberg, A. Salomaa, eds.), Thessaloniki, Greece, 2007, LNCS 4860, Springer-Verlag, Berlin, 2007, 336– 352. 12. A. Leporati, C. Zandron, C. Ferretti, G. Mauri: On the computational power of spiking neural P systems. International Journal of Unconventional Computing, to appear. 13. A. P˘aun, Gh. P˘aun: Small universal spiking neural P systems. BioSystems, 90, 1 (2007), 48–60. 14. Gh. P˘aun: Computing with membranes. Journal of Computer and System Sciences, 1, 61 (2000), 108–143. See also Turku Centre for Computer Science – TUCS Report No. 208, 1998. 15. Gh. P˘aun: Membrane Computing. An Introduction. Springer-Verlag, Berlin, 2002. 16. A. Salomaa: Formal Languages. Academic Press, New York, 1973. 17. The P systems Web page: http://ppage.psystems.eu 34 A. Krassovitskiy 1. For every rule (a, b1. . . bk, c, goj)∈Ri, a, b1, . . . , bk, c ∈V∪ {ε}, k > 0 we add to P(Qi→Qj, Da,c →Da,b1Db1,b2. . . Dbk−1,bkDbk,c),where a=½a, if a∈V, t, ∀t∈V∪ {ε},if a=εc=½c, if c∈V, t, ∀t∈V∪ {ε},if c=ε 2. For every rule (a, ε, c, inj)∈Ri, a, c ∈V∪ {ε}, k > 0 we add to P (Qi→Qj, Da,c →Da,c),where aand cdefined as in the previous case. 3. Next, for every w=b1· · ·k∈Mi, i = 1, . . . , n, k > 0 we add to Pthe matrix (S→QiDε,b1Db1,b2. . . Dbk−1,bkDbk,ε). 4. In special case if ε∈Miwe add (S→QiDε,ε) to P. 5. Also, for every Da,b ∈D, a, b ∈V∪ {ε}we add (Da,b →a) to P. 6. Finally, we add (p1→ε) to P(we assume that the first cell is the output cell). The simulation of Πby the matrix grammar is straightforward. We store the label of current cell by nonterminals Q. Every nonterminal Da,c ∈D, a, c ∈V∪{ε}, represents a pair of adjacent letters, so we can use them as a context. A rule (a, b1. . . bk, c, inj)∈Ri, a, c ∈V, b1. . . bk∈Vk, can be simulated by the grammar iff the sentential form contains both Qiand Da,c.It results the label of current cell is rewritten to Qjand Da,c is rewritten to the string Da,b1Db1,b2. . . Dbk−1,bkDbk,a. Clearly, the string preserves one symbol context. In order to treat those rules which have no context we introduce productions that preserve arbitrary context (a∈V∪ {ε}and c∈V∪ {ε}). The simulation of the grammar starts with a nondeterministic choice of the axiom. Then, during the derivation any rule corresponding to the context (a, b) have to be applied (in a one to one correspondence with grammar productions). Finally, the string over Vis produced by the grammar iff Q1has been deleted from the simulated sentential form. The deletion of Q1specifies that Πreached the output cell. So, we obtain L(Π) = L(G). Hence, LStP∗(ins1,1 ∗)⊆MAT. The strictness of the inclusion follows from the fact there are languages from MAT which cannot be generated by any insertion P system from LStP∗(ins1,1 n),for any n≥1.Indeed, consider La={cakcakc|k≥1}.One may see that the matrix grammar ({Sl, Sr, S},{a, b}, S, P 0) generates La,where P0={(S→cSlcSrc), (Sl→ε, Sr→ε), (Sl→aSl, Sr→aSr)}. On the other hand, La/∈LStP∗(ins1,1 n),for any n≥1.For the contrary, assume there is such a system. We note that the system cannot delete or rewrite any letter so every insertion is terminal. And, as the language of axioms is finite, we need an insertion rule of letter a. Consider alternatives for a final insertion step in a derivation which has at most one step and derives a word cakcakc, for some k > n: (1) the last applied rule inserts the central letter c, or (2) it does not insert the On the Power of Insertion P Systems of Small Size 35 central letter c. (1) The central ccan be inserted between any two letters a. So we get a contradiction because the prefix capcmay not be equal to the suffix caqc. (2) The last applied rule can (2.1) either insert the letter c(at the end or start of the string) or, (2.2) no cis inserted by the final rule. (2.1) We get a contradiction because ccan be alternatively inserted in between two aas we assumed k > n. (2.2) The last rule cannot distinguish whether to insert abefore the central cor after the central c. So, again, we get a contradiction because the prefix capcmay not be equal to the suffix caqc. So we proved La/∈LStP∗(ins1,1 n),for any n≥1 and hence LStP∗(ins1,1 ∗)⊂MAT. Corollary 3 LSP∗(ins1,1 ∗)⊂MAT. Proof. A tree is a special case of a graph. Lemma 4 MAT ⊆hLSP∗(ins1,1 2). Proof. We prove the theorem by a direct simulation of a matrix grammar G= (N, T, S, P ). We assume that Gis in binary normal form, i.e., every matrix has the form i: (A→BC, A0→B0C0)∈P, where A, A0∈N, B, B0, C, C0∈ N∪T∪ {ε}and i= 1, . . . , n. Consider a system Π∈hLSPn+3(ins1,1 2), Π= (V, [1[2[3 Y i=1,...,n [i+3 ]i+3  ]3]2]1,{S$},∅, . . . , ∅, R1, . . . , Rn+3, h, ϕ), where V=N∪T∪ {Ci, C0 i|i= 1, . . . , n}∪{#,$}. For every matrix i: (A→BC, A0→B0C0) we add r.1.1 : (A, #Ci, α, in2),to R1; r.2.1 : (Ci, BC, α, in3), r.2.2 : (C0 i,#, α, out) to R2; r.3.1 : (Ci,#, α, ini+3), r.3.2 : (C0 i, B0C0, α, out) to R3; r.i + 3.1 : (A0,#C0 i, α, out),to Ri+3 for every α∈V\{#}.In addition we add (ε, $, ε, out) to R1. We also define the morphism hand the weak coding ϕby: h(a) = ½a, if a∈T, a# if a∈V ϕ(a) = ½a, if a∈T, εif a∈V∪ {$}. We claim that L(Π) = L(G).To do so it is enough to prove that w∈L(G) iff w0∈L(Π) and w0=ϕ(h−1(w)). 36 A. Krassovitskiy First we show that for every w∈L(G) there exists w0∈L(Π) and w0=ϕ(h−1(w)).Consider the simulation of the i-th matrix (A→BC, A0→B0C0)∈P. The simulation is controlled by letters Ciand C0 i. First, we insert #Ciin the context of a unmarked Aand send the obtained string to the second membrane. Then we use Cias a context to insert adjacently right the word BC. After that, we mark the control letter Ciand send the sentential form to the i+ 3 membrane. Here we choose nondeterministically one letter A0, mark it, write adjacently right new control letter C0 i, and, after that, send the obtained string to the third membrane. We mention that it is not possible to apply the rule r.i + 3.1 : (A0,#C0 i, α;out)ein the i+ 3 membrane and to reach the skin membrane if the sentential form does not contain the unmarked A0. So, this branch of computation cannot influence the result and may be omitted in the consideration. Next, in the third membrane, B0C0is inserted in the context of unmarked C0 iand the sentential form is sent to the second membrane. Finally, we mark C0 iand send the resulting string back to the skin membrane. We assume that at the beginning of this simulation the sentential form in the skin membrane does not contain unmarked Ci, C0 i.Hence, the insertions in the second and third membranes are deterministic. The derivation preserves the assumption, as after the sentential form is sent back to the skin membrane the introduced Ci, and C0 iare marked. At the end of computation we send the obtained sentential form out of the system by the rule (ε, $, ε, out). Let wbe a string in the skin region which contains some unmarked Aand A0. If the letter Aprecedes B, then we can write w=w1Aα1w2A0α1w3.The simulation of the matrix is the following w1Aα1w2A0α1w3 r.1.1,r.2.1,r.3.1 =⇒w1A#Ci#BCα1w2A0α1w3 r.3+i.1,r.3.2,r.2.2 =⇒ w1A#Ci#BCα1w2A0#C0 i#B0C0α1w3, where w1, w2, w3∈V∗, α1, α2∈V\{#}.We can write the derivation similarly if Bprecedes A. Hence, as a result of the simulation of i-th matrix we get both Aand A0 marked and BC,B0C0inserted in the right positions. The derivation in Πmay terminate by the rule (ε, $, ε, out) only in the first membrane. This guarantees that the simulation of each matrix has been finished. According to the definition of Π the string w0belongs to the language if w0=ϕ(h−1(w)),where wis the generated string. This is the case only if the resulting output of Πdoes not contain unmarked nonterminals. Hence we proved L(G)⊆ϕ(h−1(L(Π))). The inverse inclusion is obvious since every rule in Πhas its counterpart in G. The case when the derivation in Πis blocked corresponds to the case a matrix cannot be finished. Hence, we get MAT ⊆hLSP∗(ins1,1 2). Remark 5 One can mention that a similar result can be obtained with a smaller number of membranes at the cost of increasing the maximal length of inserted On the Power of Insertion P Systems of Small Size 37 words. I.e., for any grammar G0from MAT there is a P insertion system Π0 corresponding to hLSPn+1(ins1,1 3)such that L(G) = L(Π0),and nis the number of matrices in G0.To prove this we can use the same argument as in the previous theorem and replace rules (r. ∗.∗) by (A, #BC, α, ini+1), α ∈V\{#}to R1 (A0,#B0C0, α, in1), α ∈V\{#}to Ri+1. Corollary 6 MAT ⊆hLStP∗(ins1,1 2). Proof. Obvious, since a tree is a special case of a graph. Taking into account Lemma 4, Lemma 2, and the fact that the class of matrix grammars is closed under inverse morphisms and weak codings we get the following characterization of MAT: Theorem 7 hLS(t)P∗(ins1,1 ∗) = MAT. Next we consider computationally complete insertion P systems. In order to use concise representations of productions in 0-type grammars we need an auxiliary lemma. Lemma 8 For every 0-type grammar G0= (N0, T, S0, P 0)there exists a grammar G= (N, T, S, P)such that L(G) = L(G0)and every production in Phas the form AB →AC or AB →CB or (1) A→AC or A→CA or (2) A→δ, (3) where A, B and Care from Nand δ∈T∪N∪ {ε}. Proof. To prove the lemma it is enough to show that for any grammar in Penttonnen normal form there is an equivalent grammar having productions of the form (1)–(3). To do so it is enough to simulate the context-free productions A→BC by productions of the form (1)–(3). Let G= (N, T, S, P) be a grammar in Penttonnen normal form whose production rules in Pare of the form: AB →AC or A→BC or A→α where A, B, C and Dare from Nand α∈T∪N∪ {ε}. Let PCF ⊆Pdenotes the set of all context-free productions A→BC ∈P such that B6=C. Suppose that rules in PCF are ordered and n=card(PCF ). Consider a grammar G0= (N0, T, S, P0),where N0=N∪ {Xi, Yi, Zi|i= 1, . . . , n}, P0= (P\PCF )∪P0 CF ,and P0 CF is constructed as follows: for every i:A→BC ∈PCF add to P0 CF the following productions 38 A. Krassovitskiy r.i.1 : A→Xi, r.i.2 : Xi→XiYi, r.i.3 : XiYi→ZiYi, r.i.4 : ZiYi→ZiC r.i.5 : ZiC→BC Clearly, the obtained grammar has the form specified by (1)–(3). Now we prove that L(G) = L(G0).The inclusion L(G)⊆L(G0) is obvious as for every derivation in Gwe use its counterpart derivation in G0replacing i-th context-free production from PCF by the sequence of productions r.i.1, r.i.2, r.i.3, r.i.4, r.i.5: wAw0r.i.1 =⇒wXiw0r.i.2 =⇒wXiYiw0r.i.3 =⇒ wZiYiw0r.i.4 =⇒wZiCw0r.i.5 =⇒wBCw0. In order to prove that L(G0)⊆L(G) we show that for every terminal derivation in G0we can construct a derivation in Gso that they both produce the same word. We use the counterpart productions from P\PCF to mimic analogous production. For the productions PCF 0we show that any deviation from the above defined sequence does not produce any new terminal derivation. First, we mention that the sequence of productions corresponding to i:A→BC starts by rewriting A on the new nonterminal, so, other productions not in r.i.∗cannot interfere the sequence. Yet, the production rule r.i.2 may generate extra Yi(for simplicity, we assume one extra Yigenerated). wAw0r.i.1 =⇒wXiw0(r.i.2)2 =⇒wXi(Yi)2w0r.i.3,r.i.4 =⇒ r.i.5 =⇒wBCYiw0. As we need to consider only terminal derivations we may assume that Yi will be necessary rewritten. The only rule to rewrite Yiis r.i.4. In order to perform it the letter Zimust precede by the letter Yi. It implies that the letter A must appear adjacently left from the letter Yi. Then the sequence of productions r.i.1, r.i.3, r.i.4, r.i.5 results to the same sentential form as if r.i.2 is applied once per every rewriting of A. wBCYiw0=⇒∗w1AYiw0 1 r.i.1 =⇒ wXiYiw0 1 r.i.3,r.i.4,r.i.5 =⇒w1BCw0 1. Therefore, it can be produced by rules from PCF .We also mention that in r.i.1−r.i.5 we start rewriting letters from N\N0by corresponding letters from N only after the letter Xiis rewritten. This imply that after r.i.4, we cannot insert additional Yiadjacently left from Ci. So, Zican be rewritten unambiguously. Finally, consider the case when the production r.i.4 is followed by some production from P\PCF which rewrites C. . . . r.i.4 =⇒wZiCiw0=⇒∗w2ZiCiw0 2 r.i.5 =⇒w2BCiw0 2. On the Power of Insertion P Systems of Small Size 39 As Zican be rewritten only if Ciappears to the right we may consider the equivalent derivation with the production r.i.5 applied directly after r.i.4.So the derivation is equivalent to the derivation with i:A→BC ∈PCF . Hence we proved L(G0)⊆L(G), and hence L(G0) = L(G). Every grammar in the normal form has the following property: every production can rewrite/add at most one (nonterminal) letter. Now we increase the maximal size of the context of insertion rules to two letters. It is known from [17] that the class INS2,2 2contains non-semilinear languages. Considering these systems with membrane regulation we get Theorem 9 hLSP3(ins2,2 2) = RE. Proof. We prove the theorem by simulating a 0-type grammar in the normal form from Lemma 8. Let G= (N, T, S, P ) be such a grammar. Suppose that rules in P are ordered and n=card(P). Now consider the following insertion P system, Π= (V, [1[2[3]3]2]1,{S$},∅,∅, R1, R2, R3, h, ϕ),where V=T∪N∪F∪F∪ {#,#,$},F={FA,|A∈N},F={FA,|A∈N}. We include into R1the following rules: (AB, #C, α, here),if AB →AC ∈P; (A, #C, Bα, here),if AB →CB ∈P; (A, C, α, here),if A→AC ∈P; (ε, C, Aα, here),if A→CA ∈P; (A, #δ, α, here),if A→δ∈P; ($, ε, ε, out), where α∈V\{#}.It may happen that the pair of letters AB subjected to be rewritten by a production AB →AC or AB →CB ∈Ris separated by letters that have been marked. We use two additional membranes to transfer a letter over marked ones. In order to transfer A∈Nwe add r.1.1 : (A, #FA, α, in2), α ∈V\{#} to the skin membrane. Then we add to the second membrane r.2.1 : (FA,#A, α0, out), r.2.2 : (FA,#A, α0, out), r.2.3 : (FAX, #FA,#, in3), r.2.4 : (FAFB,#FA,#, in3) r.2.5 : (FAX, #FA,#, in3), r.2.6 : (FAFB,#FA,#, in3) r.2.7 : (FA#, FA, α, in3), r.2.8 : (FA#, FA, α, in3), r.2.9 : (FA#, FA,#, in3), r.2.10 : (FA#, FA,#, in3), r.2.11 : (FA#, FA,#, in3), r.2.12 : (FA#, FA,#, in3), 40 A. Krassovitskiy for every X∈F∪N, FB∈F , α ∈V\{#,#}, α0∈ {ab |a∈N∪T, b ∈N∪T∪ {$}} ∪ {$}. Finally, we add to the third membrane the rules r.3.1 : (FA,#, α, out), α ∈V\{#}, r.3.2 : (FA,#, α, out), α ∈V\{#} The morphism his defined by h(a) = ½a, if a∈V\N, a# if a∈N. The weak coding ϕis defined by ϕ(a) = ½a, if a∈T, εif a∈V\T. We simulate the productions of Pin the skin membrane by marking nonterminals from Nand inserting corresponding letters of the productions. This is possible to do with insertion rules of weight (2,2,2) since the grammar has such a form so every production rewrite/add at most one letter. The simulation of the transfer is done in the second and the third membranes. The idea of the simulation is (1) to mark nonterminal we want to transfer, (2) jump over the marked letters with help of one special letter, at the end (3) mark the special letter and insert the original nonterminal. Since we use two letter contexts, in one step we can jump only over a single letter. Also we need to jump over the marking letter # as well as over marked nonterminals, and the letters inserted previously. We use letters FA∈Fand FA∈Fto keep information about the letter Awe want to transfer. In order to jump over # we introduce one additional marking symbol #.We mark letters from Fby #,and all other letters in V\{#,#} by #.E.g., in a words FA#,letter FAis unmarked. (1) The rule r.1.1 : (A, #FA, α, in2),specifies that every unmarked letter from Nmay be used for the transfer. (2) The rules r.2.3−r.2.12 in the second membrane specify that FAor FAis copied to the right in such a way that the inserted letter would not be marked. In order to do so, the appropriate rule chooses to insert either the overlined copy FAor the simple copy FA. The rules r.2.3−r.2.6 describe the possible jumps over one letter not in {#,#},and r.2.7−r.2.12 describe the possible jumps over the marking letters #, #.These rules send the sentential form to the third membrane. The rules in the third membrane mark one symbol FA∈For FA∈Fand send the sentential form back to the second membrane. (3) The rules r.2.1 and r.2.2 may terminate the transferring procedure and send the sentential form to the first membrane if letter $ or two letters from {ab |a∈N∪T, b ∈N∪T∪ {$}} appear in the the right context. On the Power of Insertion P Systems of Small Size 41 For example, consider the transfer of Ain the string AX#C$ (here, we underline inserting substrings) AX#C$r.1.1 =⇒A#FAX#C$r.2.3 =⇒$A#FAX#FA#C$r.3.1 =⇒ A#FA#X#FA#C$r.2.6 =⇒A#FA#X#FA#FAC$r.3.2 =⇒ A#FA#X#FA##FAC$r.2.1 =⇒A#FA#X#FA##FA#AC$ The sentential form preserves the following invariant: •The first membrane does not contain unmarked letters from F∪F . •There is exactly one unmarked letter F∪Fin the second membrane. •There are always two unmarked letters from F∪Fin the third membrane. We mention that the invariant is preserved by every derivation. Indeed, we start derivation from the axiom S$ that satisfies the invariant, then one unmarked symbol is inserted by r.1.1.Rules r.2.3−r.2.12 always add one more unmarked letter. And rules r.2.1, r.2.2, r.3.1, r.3.2 always mark one letter from F∪F . In order to verify that Πsimulates the same language as Gwe note that every reachable sentential form in Gwill be reachable also in Πby simulating the same production. Also we note that the derivation in Πmay terminate by the rule ($, ε, ε, out) only in the first membrane. Hence it guarantees that every transfer will be completed. It follows from the invariant that the simulation of the transfer is deterministic in the second membrane. There is a nondeterministic choice in the third membrane, where corresponding rules may mark one of the two unmarked letters. In the case the rule marks the rightmost letter, the derivation has to “jump” again over the inserted letter. The transfer satisfies the property that every terminating sequence replaces a nonterminal via arbitrary large string of marked letters, if it starts (by r.1.1) adjacently left from it. And in case the rule r.1.1 starts the transfer of a letter next to unmarked letter then it produces two marked symbols which do not affect on the result of the simulation. The output string wis in the language only if w0=ϕ(h−1(w)) is defined. This is the case only if the resulting output of Πdoes not contain unmarked nonterminals. On the other hand, every final derivation in Πhas its counterpart in G. By applying the inverse morphism h−1we filter out every sentential form with unmarked nonterminals from N. Hence, the corresponding derivation in Gis terminated. Finally, the weak coding ϕfilters away supplementary letters. Hence we have L(G) = L(Π). Remark 10 One may mention that there is a trade-off between the number of membranes and the maximal length of productions. By introducing additional nonterminals and fitting the grammar into the normal form we decrease the amount of used membranes. It is also the case in the other way: by growing the number of membranes we can simulate larger production rules. 42 A. Krassovitskiy 4 Conclusion This article investigates the expressive power of insertion P systems with encodings. The length of insertion rules and number of membranes are used as a measure of descriptional complexity of the system. In the article we use the fact that morphisms and weak codings are incorporated to insertion P systems. The obtained family hLS(t)P∗(ins1,1 ∗) serves to characterize matrix languages. When no membranes are used, the class hINS1,1 ∗equals the family of context-free languages. We proved the universality for the family hLSP∗(ins2,2 ∗). More precisely, for every recursively enumerable language we can construct an insertion P system of weight (2,2,2) and degree 3, so that applying an inverse morphism and a weak coding we generate the same language. Also, we want to mention that computational completeness of considered families indicates some trade-offs between descriptional complexity measures. Moreover, the descriptional complexity used in the paper may be extended by internal system parameters as, e.g., the size of alphabet, the number of rules per membrane, etc. Also, it seems quite promising to investigate decidable computational properties of the language family LS(t)P∗(ins∗,∗ ∗).We conjecture that it is incomparable with many known language families. We recall the open problem posed in [4], namely, whether hLSP1(ins2,2 ∗) equals RE One may see that in order to solve the problem by the technique used in the article, it is enough to find a concise way to transfer a letter over a marked context. In our case this can be reduced to the question whether it is possible to compute the membrane regulation in the skin membrane. Acknowledgments The author acknowledges the support PIF program of University Rovira i Virgili, and projectno. MTM2007-63422 from the Ministry of Science and Education of Spain. The author sincerely thanks the help of Yurii Rogozhin and Seghey Verlan without whose supervision the work would not been done. The author also warmly thanks Gheorghe P˘aun and Erzs`ebet Csuhaj-Varj´u who draw attention, during BWMC-09, to the possibility of using different normal forms of grammars to simulate P insertion-deletion systems. References 1. M. Daley, L. Kari, G. Gloor, R. Siromoney: Circular contextual insertions/deletions with applications to biomolecular computation. Proc. SPIRE’99 , 1999, 47–54. 2. M. Mutyam, K. Krithivasan, A. Siddhartha Reddy: On characterizing recursively enumerable languages by insertion grammars. Fundam. Inform., 64, 1-4 (2005), 317– 324 3. B.S. Galiukschov: Semicontextual grammars. Matematika Logica i Matematika Linguistika, Tallin University, 1981 38–50 (in Russian). On the Power of Insertion P Systems of Small Size 43 4. L. Kari, Petr Sos´ık: On the weight of universal insertion grammars. Theoretical Computer Sci., 396, 1-3 (2008), 264–270. 5. L. Kari: From micro-soft to bio-soft. Computing with DNA. Proc. on Biocomputing and Emergent Computations, 1997, 146–164. 6. L. Kari, Gh. P˘aun, G. Thierrin, S. Yu: At the crossroads of DNA computing and formal languages: characterizing RE using insertion-deletion systems. Proc. of 3rd DIMACS, 1997, 318–333. 7. L. Kari, G. Thierrin: Contextual insertion/deletion and computability. Information and Computation, 131, 1 (1996), 47–61. 8. A. Krassovitskiy, Yu. Rogozhin, S. Verlan: Further results on insertion-deletion systems with one-sided contexts. LNCS 5196, 2008, 333–344. 9. A. Krassovitskiy, Yu. Rogozhin, S. Verlan: One-sided insertion and deletion. Traditional and P systems case. Proc. CBM08, 2008, Austria, 53–64. 10. A. Krassovitskiy, Yu. Rogozhin, S. Verlan: Computational power of P systems with small size insertion and deletion rules. Proc. CSP08, Ireland, 2008, 137–148. 11. S. Marcus: Contextual grammars. Rev. Roum. Math. Pures Appl., 14 (1969), 1525– 1534. 12. M. Margenstern, Gh. P˘aun, Yu. Rogozhin, S. Verlan: Context-free insertion-deletion systems. Theoretical Computer Sci., 330 (2005), 339–348. 13. C. Martin-Vide, Gh. P˘aun, A. Salomaa: Characterizations of recursively enumerable languages by means of insertion grammars. Theoretical Computer Sci., 205, 1–2 (1998), 195–205. 14. A. Matveevici, Yu. Rogozhin, S. Verlan: Insertion-deletion systems with one-sided contexts. LNCS 4664, 2007, 205–217. 15. Gh. P˘aun: Marcus Contextual Grammars. Kluwer, Dordrecht, 1997. 16. Gh. P˘aun: Membrane Computing. An Introduction. Springer, Berlin, 2002. 17. Gh. P˘aun, G. Rozenberg, A. Salomaa: DNA Computing. New Computing Paradigms. Springer, Berlin, 1998. 18. G. Rozenberg, A. Salomaa, eds.: Handbook of Formal Languages. Springer, Berlin, 1997. 19. A. Takahara, T. Yokomori: On the computational power of insertion-deletion systems. LNCS, 2568, 2003, 269–280. 20. S. Verlan: On minimal context-free insertion-deletion systems. Journal of Automata, Languages and Combinatorics, 12, 1-2 (2007, 317-328. 50 M.A. Mart´ınez–del–Amor et al. The Tesla C1060 achieves 102 GB/sec of bandwidth to the off-chip memory (running at 800 MHz). This bandwidth is not enough for the big set of cores and the possibilities to saturate it are high. To obtain the maximum bandwidth available it is needed to coalesce accesses to the device memory. The coalesced accesses are obtained whenever the accesses are contiguous 16-word lines, otherwise a fraction of this bandwidth it is obtained. Coalesced accesses will be a critical point in the optimization process. In addition, the threads can use other memories like constant memory or texture memory. Reading from constant cache is as fast as reading from a registers, as long as all threads in the same warp read the same address. Texture Memory is optimized for 2D spatial locality (see Table 1). Table 1. Memory System on the Tesla C1060 Memory Location Size Latency Access Registers On-Chip 16384 32-bits Registers per SM ≃0 cycles R/W Shared Memory On-Chip 16 KB per SM ≃registers R/W Constant On-Chip 64 KB ≃registers R Texture On-Chip Up to Global >100 cycles R Local Off-Chip 4 GB 400-600 cycles R/W Global Off-Chip 4 GB 400-600 cycles R/W 3.2 Threading Model A SM is a hardware device specifically designed with multithreaded capabilities. Each SM manages and executes up to 1024 threads in hardware with zero scheduling overhead. Each thread has its own thread execution state and can execute an independent code path. The SMs execute threads in a Single-Instruction Multiple- Thread (SIMT) fashion [14]. Basically, in the SIMT model all the threads execute the same instruction on different piece of data. The SMs create, manage, schedule and execute threads in groups of 32 threads. This set of 32 threads is called Warp. Each SM can handle up to 32 Warps (1024 threads in total, see Table 2). Individual threads of the same Warp must be of the same type and start together at the same program address, but they are free to branch and execute independently. The execution flow begins with a set of Warps ready to be selected. The instruction unit selects one of them, which is ready for issue and execute instructions. The SM maps all the threads in an active Warp to the SP cores, and each thread executes independently with its own instructions and register state. Some threads of the active Warp can be inactive due to branching or predication, and this is also another critical point in the optimisation process. The maximum performance is achieved when all the threads in an active Warp takes the same path (the same execution flow). If the threads of a Warp diverge, the Warp serially executes each branch path taken, disabling threads that are not on that path, and when all the paths complete, the threads reconverge to the original execution path. Simulation of Recognizer P Systems by Using Manycore GPUs 51 Table 2. Major Hardware and Software Limitations programming on CUDA Configuration Parameters Limitation Threads/SM 1024 Thread Blocks/SM 8 32-bit Registers/SM 16384 Shared Memory/SM 16KB Threads/Block 512 Threads/Warp 32 Warps/SM 32 3.3 Parallel Computing with CUDA The GPU is, nowadays, a single-chip massively parallel system which is inexpensive and readily available. However, programming a highly-parallel system has historically been a domain of few experts [25]. The emergence of Compute Unified Device Architecture (CUDA) has helped develop highly-parallel applications easier than before. CUDA programming toolkit is an extension of ANSI C including several keywords and constructs. The GPU is seen as a coprocessor that executes data-parallel kernel functions. The user creates a program encompassing CPU code (Host code) and GPU code (Kernel code). These are separated and compiled by nvcc (Nvidia’s compiler for CUDA code). The host code is responsible for transfer data to and from the GPU memory (device memory) via API calls, to initiates the kernel code executed on the GPU. The threads executes the kernel code, and they are organized into a three-level hierarchy. At the highest level, each kernel creates a single grid that consists of many thread blocks. Besides, each thread block can contain up to 512 threads which can share data through Shared Memory and can perform barrier synchronization by invoking the --syncthreads primitive [25]. On the other hand, blocks can not perform synchronization. The synchronization across blocks can only be obtained by terminating the kernel. Finally, the threads within the block are organized into warps of 32 threads. Each block within the grid have their own identifier[18]. This identifier can be one, two or three dimensions depending on how the programmer has declared the grid. In the same way, each thread within the block have their own identifier which can be one, two or three dimensions as well. Combining thread and block identifiers, the threads can access to different data address and also select the work that they have to do. 52 M.A. Mart´ınez–del–Amor et al. 4 A Design of the Simulator for the Class of Recognizing P Systems 4.1 Algorithm Design Whenever we design algorithms in the CUDA programming model, our main effort is dividing the required work into processing pieces, which have to be processed by TB thread blocks of Tthreads each. Using a thread block size of T=256, we have empirically determined to obtain the overall best performance on the Tesla C1060. Each thread block access to one different set of input data, and assigns a single or small constant number of input elements to each thread. Each thread block can be considered independent to the other, and it is at this level at which internal communication (among threads) is cheap using explicit barriers to synchronize, and external communication (among blocks) becomes expensive, since global synchronization only can be achieved by the barrier implicit between successive kernel calls. The need of global synchronization in our designs requires successive kernel calls even to the same kernel. 4.2 P System Simulator with Active Membranes The simulator simulates a recognizer P system with active membranes, i.e Π= (O, H, µ, ω1, . . . , ωm, R) according to the notation described in section 2. The simulator is executed into two main stages: selection stage and execution stage. The selection stage consists of the search for the rules to be executed in each membrane. Once the rules have been selected, the execution stage consists of the execution of these rules. The selection stage takes the major part of the simulation time in the sequential code, since this part of the algorithm implies to check all the rules of the system in every membrane. So we have parallelized the selection stage on the GPU, and the execution stage is still executed on the CPU at this point of the implementation. The input data for the selection stage consists of the description of the membranes with their multisets (strings over O, labels associated with the membrane in H, etc...) and the set of rules Rto be selected. The output data of this stage will be the set of selected rules per membrane. Only the execution stage changes the information of the configuration. Besides, we have identified each membrane as a thread block where each thread represents an element of the alphabet O. Each thread block runs in parallel looking for the set of rules that has to execute, and each individual thread is responsible for identifying if there are some rules associated with the element that it represents, and if so, send it back to the execution stage. Finally, the CPU takes the control and executes the rules previously selected. As result of the execution stage, the membranes can vary including news elements, dissolving membranes, dividing membranes, etc. Therefore, we have to Simulation of Recognizer P Systems by Using Manycore GPUs 53 modify the input data for the selection stage with the newest structure of membranes, and then call the selection again. It is an iterative process until a system response is reached. Our simulator presents two restrictions: it can handle only two levels of membrane hierarchy for simplicity (the skin and the rest of elementary membranes), what is enough for solving lots of NP-complete problems; moreover, the number of objects in the alphabet must be divisible by a number smaller than 512 (the maximum thread block size), in order to distribute the objects among the threads equally. 5 A Case Study: Implementing a Solution to the N-Queens Problem In this section, we present a solution to the N-Queens problem, given by Miguel A. Guti´errez–Naranjo et al [8], using our simulator. The N-Queens problem is expressed as a formula in conjunctive normal form, in such way that one truth assignment of the formula is considered as N-Queens solution. A family of recognizer P system for the SAT problem [22] can state whether exists a solution to the formula or not sending yes or no to the environment. However, the yes ot no answer from the recognizer P system is not enough. Besides, the system needs to give us the way to encode the state of the N-Queens problem. The P system designed for solving the N-Queens problem is a modification of the P system for the SAT problem. It is an uniform family of deterministic recognizer P system which solves SAT as a decision problem (i.e., the P system sends yes or no to the environment in the last computation step), but it also stores the truth assignments that makes true the formula encoded in the elementary membranes of the halting configuration. 5.1 Implementation P-Lingua 1.0[4] is a programming language useful for defining P system models with active membranes. We use P-Lingua to encode a solution to the N-Queens problem, and also to generate a file that our simulator can use as input. Figure 2 shows the P-Lingua process to generate the input for our simulator. P-Linga 2.0[5] translates a model written in P-Lingua language into a binary file. A binary file is a file whose information is encoded in Bytes and bits (not understandable by humans like plain text), which is suitable for trying to compress the data. This binary file contains all the information of the P system (Alphabet, Labels, Rules, . . . ) which is executed by our simulator. In our tests, we use the recognizer P system for solving the 3-Queens problem. This problem creates 512 membranes and up to 1300 different objects. For the 4-Queens problem, the system would create 65536 membranes and up to 8000 54 M.A. Mart´ınez–del–Amor et al. Fig. 2. Generation of the simulator’s input different objects. Currently, our simulator can not handle this example due to memory space limitation (requires up to 8GB in device memory). This problem can be solved with overlays of data and subsequent calls to the GPU. We are working on this solution, and also in other solutions to make possible execute problems with bigger memory size constraints on our simulator. On the other hand, note that 2-Queens is a system with only 4 membranes, what is not enough for exploiting the parallelism in P systems. 6 Performance Analysis We now examine the experimental performance of our simulator. Our performance test are based on the solution to 3-Queen problem previously explained 5.1. Although this problem does not cover all the NP-complete problems that we want to simulate in our simulator, it states an example of how a NP-complete problem can be solved on the P system with active membranes simulator. We report the selection stage time which is executed on the GPU, and compare it with the selection stage for the sequential code. We do not include the cost of transferring input data from host CPU memory across the PCI-Express bus to the GPU’s on board memory. Selection is one building block of larger-scale computation. Our aim is to get a full implementation of the simulator on the GPU. In such case, the transfers across PCI-Express bus will be close to zero. The selection stage on the GPU takes about 195 msec. This is 12 times faster than the selection stage on the CPU which takes 2345 msec. We have used the NVIDIA GPU Tesla C1060 which has 240 execution cores and 4GB of device memory, plugged in a computer server with a Intel Core2 Quad CPU and 8GB of RAM, using the 32bits ubuntu server as Operating System. Our experimental results demonstrate the results we expect to see: a massivelyparallel problem such as selection of the rules in a P system with active membranes achieves faster running times on a massively-parallel architecture such as GPU. Simulation of Recognizer P Systems by Using Manycore GPUs 55 7 Conclusions and Future Work In this paper, we have presented a simulator for the class of recognizer P systems with active membranes using CUDA. The membrane computation has double parallel nature. The first level of parallelism is presented by the objects inside the membranes, and the second one is presented between membranes. Hence, we have simulated these P systems in a platform which provides those levels of parallelism. This platform is the GPU, with parallelism between thread blocks and threads. Besides, we have used a programming language called P-Lingua to generate a solution to the N-Queens problem, in order to test our simulator. Using the power and parallelism that provides the GPU to simulate P systems with active membranes is a new concept in the development applications for membrane computing. Even the GPU is not a cellular machine, its features help the researches to accelerate their simulations allowing the consolidation of the cellular machines as alternative to the traditional machines. The first version of the simulator is presented for P systems with active membranes, specifically, we have developed the selection stage of the simulator. In forthcoming versions, we will try to include the execution version in the GPU. This issue allows a completely parallel execution on the GPU, avoiding CPU-GPU transfers in every step, which degrades system performance. On the other hand, we shall adapt our simulator to use the resources available on the GPU at maximum. To develop general purpose programs on the GPU is easier than several years ago with tools such as CUDA. However, extracting the maximum performance on the GPU is still hard, so we need to make a deep analysis to obtain the maximum performance available for our simulator. It is also important to point out that this simulator is limited by the resources available on the GPU as well as the CPU (RAM, Device Memory, CPU, GPU). This limits the size of the instances of NP-complete problems whose solutions can be successfully simulated. In the following version of the simulator, we will try to reduce the memory requirements for the simulator in order to be able to simulate bigger instances of NP-complete problems. Moreover, it would be interesting to design heuristics to accelerate the computations of our simulator. Although, the massively parallel environment that provides the GPUs is good enough for the simulator, we need to go beyond. The newest cluster of GPUs provides a higher massively parallel environment, so we will attempt to scale to those systems to obtain better performance in our simulated codes. Acknowledgement The first three authors acknowledge the support of the project TIN2006–13425 of the Ministerio de Educaci´on y Ciencia of Spain, cofinanced by FEDER funds, and the support of the “Proyecto de Excelencia con Investigador de Reconocida Val´ıa” of the Junta de Andaluc´ıa under grant P08–TIC04200. The last three authors acknowledge the support of the project from the Spanish MEC and Euro- 56 M.A. Mart´ınez–del–Amor et al. pean Commission FEDER funds under grants “Consolider Ingenio-2010 CSD2006- 00046” and “TIN2006-15516-C04-03”, as well as by the EU FP7 NoE HiPEAC IST-217068. References 1. A. Alhazov, M.J. P´erez–Jim´enez: Uniform solution of QSAT using polarizationless active membranes. In J. Durand-Lose, M. Margenstern, eds., Machines, Computations, and Universality. LNCS 4664, Springer, 2007, 122–133. 2. I. Buck, T. Foley, D. Horn, J. Sugerman, K. Fatahalian, M. Houston, P. Hanrahan: Brook for GPUs: stream computing on graphics hardware. SIGGRAPH ’04, ACM Press (2004), 777–786. 3. G. Ciobanu, Gh. P˘aun, M.J. P´erez–Jim´enez, eds.: Applications of Membrane Computing. Springer, 2006. 4. D. D´ıaz–Pernil, I. P´erez–Hurtado, M.J. P´erez–Jim´enez, A. Riscos–N´u˜nez: P–Lingua: A programming language for membrane computing. In D. D´ıaz–Pernil, M.A. Guti´errez–Naranjo, C. Graciani–D´ıaz, Gh. P˘aun, I. P´erez–Hurtado, A. Riscos–N´u˜nez, eds., Proceedings of the 6th Brainstorming Week on Membrane Computing, Sevilla, F´enix Editora, 2008, 135–155. 5. M. Garc´ıa–Quismondo, R. Guti´errez–Escudero, I. P´erez–Hurtado, M.J. P´erez– Jim´enez: P–Lingua 2.0: added features and first applications. In this volume. 6. M. Garland, S.L. Grand, J. Nickolls, J. Anderson, J. Hardwick, S. Morton, E. Phillips, Y. Zhang, V. Volkov: Parallel computing experiences with CUDA. IEEE Micro, 28, 4 (2008), 13–27. 7. N.K. Govindaraju, D. Manocha: Cache–efficient numerical algorithms using graphics hardware. Parallel Comput., 33, 10–11 (2007), 663–684. 8. M.A. Guti´errez–Naranjo, M.A. Mart´ınez–del–Amor, I. P´erez–Hurtado, M.J. P´erez– Jim´enez: Solving the N−Queens puzzle with P systems. In this volume. 9. M.A. Guti´errez–Naranjo, M.J. P´erez–Jim´enez, A. Riscos–N´u˜nez: Available membrane computing software. In G. Ciobanu, Gh. P˘aun, M.J. P´er´ez–Jim´enez (eds.) Applications of Membrane Computing, Springer–Verlag, 2006. Chapter 15 (2006), 411–436. 10. M.A. Guti´errez–Naranjo, M.J. P´erez–Jim´enez, A. Riscos–N´u˜nez: Towards a programming language in cellular computing. Electronic Notes in Theoretical Computer Science, 123 (2005), 93–110. 11. M. Harris, S. Sengupta, J.D. Owens: Parallel prefix sum (Scan) with CUDA. GPU Gems, 3 (2007). 12. T.D. Hartley, U. Catalyurek, A. Ruiz, F. Igual, R. Mayo, M. Ujaldon: Biomedical image analysis on a cooperative cluster of GPUs and multicores. ICS ’08: Proceedings of the 22nd annual international conference on Supercomputing, ACM, 2008, 15–25. 13. M.D. Lam, E.E. Rothberg, M.E. Wolf: The cache performance and optimizations of blocked algorithms. ASPLOS-IV: Proceedings of the fourth international conference on Architectural support for programming languages and operating systems, ACM, 1991, 63–74. 14. E. Lindholm, J. Nickolls, S. Oberman, J. Montrym: NVIDIA Tesla. A unified graphics and computing architecture. IEEE Micro, 28, 2 (2008), pp. 39–55. Simulation of Recognizer P Systems by Using Manycore GPUs 57 15. W.R. Mark, R.S. Glanville, K. Akeley, M.J. Kilgard: Cg – a system for programming graphics hardware in a C–like language. SIGGRAPH ’03, ACM, 2003, 896–907. 16. J. Michalakes, M. Vachharajani: GPU acceleration of numerical weather prediction. IPDPS, 2008, 1–7. 17. J. Nickolls, I. Buck, M. Garland, K. Skadron: Scalable parallel programming with CUDA. Queue, 6, 2 (2008), pp. 40–53. 18. J D. Owens, M. Houston, D. Luebke, S. Green, J.E. Stone, J.C. Phillips: Gpu computing. Proceedings of the IEEE, 96, 5 (2008), 879–899. 19. J.D. Owens, D. Luebke, N. Govindaraju, M. Harris, J. Krger, A.E. Lefohn, T.J. Purcell: A survey of general–purpose computation on graphics hardware. Computer Graphics Forum, 26, 1 (2007), pp. 80–113. 20. Gh. P˘aun. Computing with membranes. Journal of Computer and System Sciences, 61, 1 (2000), 108–143, and Turku Center for Computer Science-TUCS Report No 208. 21. M. J. P´erez–Jim´enez, A. Romero–Jim´enez, F. Sancho–Caparrini: Complexity classes in models of cellular computing with membranes. Natural Computing, 2, 3 (2003), 265–285. 22. M. J. P´erez–Jim´enez, A. Romero–Jim´enez, F. Sancho–Caparrini: A polynomial complexity class in P systems using membrane division. In E. Csuhaj–Varj´u, C. Kintala, D. Wotschke, G. Vaszil, eds., Proceedings of the 5th Workshop on Descriptional Complexity of Formal Systems, DCFS 2003, Budapest, 2003, 284–294. 23. A. Ruiz, M. Ujaldon, J.A. Andrades, J. Becerra, K. Huang, T. Pan, J.H. Saltz: The GPU on biomedical image processing for color and phenotype analysis. BIBE, 2007, 1124–1128. 24. S. Ryoo, C. Rodrigues, S. Baghsorkhi, S. Stone, D. Kirk, W. mei Hwu: Optimization principles and application performance evaluation of a multithreaded GPU using CUDA. Proceedings of the 13th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming, 2008, 73–82. 25. S. Ryoo, C.I. Rodrigues, S.S. Stone, J.A. Stratton, Sain-Zee Ueng, S.S. Baghsorkhi, W.W. Hwu: Program optimization carving for GPU computing. J. Parallel Distrib. Comput., 68, 10 (2008), 1389–1401. 26. NVIDIA CUDA Programming Guide 2.0, 2008: http://developer.download. nvidia.com/compute/cuda/2 0/docs/NVIDIA CUDA Programming Guide 2.0.pdf 27. GPGPU organization. World Wide Web electronic publication: www.gpgpu.org 28. NVIDIA CUDA. World Wide Web electronic publication: www.nvidia.com/cuda Sleep-Awake Switch with Spiking Neural P Systems: A Basic Proposal and New Issues Jack Mario Mingo Computer Science Department University Carlos III of Madrid Avda. de la Universidad Carlos III, 22, Colmenarejo (Madrid), Spain [email protected] Summary. Spiking Neural P Systems are a kind of Membrane Systems developed with the aim of incorporating ideas from biological systems, known as spiking neurons, in the computational field. Initially, these systems were designed to take concepts of the neural science based on action potentials with the purpose of testing its possibilities from a computational point of view and not to be used as neurological models. In this work, a basic approach in the opposite sense is reviewed by means of the application of such systems on a well-known biological phenomenon. This phenomenon refers to the fluctuations among neural circuits which are responsible for swapping between awakeasleep states. This basic approach is analyzed and new issues are exposed. 1 Introduction Sleep is a highly organized and actively induced cerebral state with different stages [11]. It has been observed two kinds of sleep: REM-sleep and non-REM sleep, each one with a number of specific features. Specifically, non-REM sleep comprises four stages. Each stage is defined according to its activity in the electroencephalogram (EGG). Stage 1 contains alternative periods of alpha activity (8-12Hz), irregular speed activity and theta activity (3.5-7.5Hz). Stage 2 does not show alpha activity in the EGG although in this stage appears a phenomenon called sleep spindles (bursts of 12-14Hz sinusoidal waves) and, sometimes, high-voltage biphasic waves called K complexes. Stage 3 consists of delta activity (less than 3.5Hz) during part of its time (20-50%) and stage 4 consists of delta activity during the most time (more than 50%). Approximately 90 minutes after sleep starts human beings fall in non-REM sleep which is characterized by rapid eye movements, unsynchronize EGG, a nearly complete inhibition of skeletal muscle tone (atonia). The brain temperature and metabolic rate are high, equal to or greater than during the waking state and some sexual activity can be present as well. From a neural point of view, although sleep is controlled by the great part of the encephalon [4], there is a particularly important zone. It is the ventrolat- 66 J.M. Mingo •The ascending arousal system in the brainstem and the forebrain are necessaries to accumulate adenosine as consequence of a long activity of their neurons. This means that when the organism is awoken the adenosine is accumulated. The neural groups which are responsible to accumulate this necessity are called Accumulator System. •The dorsomedial nucleus of the hypothalamus (DMH) is the most active part in the suprachiasmatic nucleus (SCN) and it is responsible to control sleep-wake transitions following a 24-hour cycle. This neural group is called DMH System and its goal is to provide a motivation to awake. •The Acumulator System and the DMH system act as activators for the asleep and awake states but once they have activated their neural groups they stop to fire and neurons in the activation system (neurons of the groups LHA, TMN, LC and DR) and neurons in the regulator-VLPO combined system start to fire while the system remains sleeping or awakening. •Bearing in mind the previous suppositions the DMH system and the LHA are the neural groups that control the awake state and the Acumulator system and the Regulator System are the neural groups that control the asleep state. Figure 5 shows the possible connections that control sleep-wake switch from a neurological viewpoint. The connections were explained in Section 2. Fig. 5. Schematic drawing of the flip-flop circuit In Figure 5, LHA, DMH, TMN, LC, DR, VLPO (VLPO dense group), eVLPO (VLPO extended group) are neural groups defined previously in section two. Acumulator represents the region where the necessity to sleep is accumulated and Sleep-Awake Switch with Spiking Neural P Systems 67 Regulator would be a component of the Acumulator and its goal is to feed neurons in the VLPO. Both, accumulator and regulator groups are suppositions in the model because, from a neural point of view, regions that are responsible of accumulating adenosine are not known. Starting from the previous figure and the SN P system with inhibitory rules a basic model for sleep-wake system is shown in Figure 6 where the inhibitory rules are represented as discontinuous lines and the excitatory rules as continuous lines. Fig. 6. Sleep-Wake basic switch with a Spiking P System 4 Results of the Basic Model To analyze conveniently the system several suppositions must be taken into account: •The system supposes that an individual is awoken during 16 hours and is sleeping during 8 hours. Another configurations, for example, 18 hours awake/6 hours sleep, are possible but always a fixed period is maintained. •Each step of time in the system represents one hour of real time. The basic model described previously can be applied starting from an initial configuration. Table 7 shows results when the system starts in awake state. First row on each square shows initial spikes and the second one represents the rule 68 J.M. Mingo applied. An exclamation symbol close to spikes means that the spike gets out from the network to the environment and it is useful to identify two states: on or off. When neurons in the LHA are firing the system represents wakefulness and if neurons in the regulator system are actives the system falls asleep. A brief explanation shows these situations. In wakefulness the neurons in the activation system (LHA, TMN, LC and DR) are firing and they select and execute a rule (a+/a →a; 0 is selected in LHA and a2/a →a; 0 in TMN, LC and DR). The activated rule in LHA, TMN and LC locks during a step of time neurons in the VLPO, both the dense group and the extended group, because of inhibitory rules among them. A lock is shown in the table with the Bsymbol. DMH system acts like the system’s clock. In t= 0 DMH contains mspikes (m >> 16 is a necessary supposition if system does not finish) and it applies its only rule (a16 →a16; 23). This rule sends 16 spikes after 23 step of time because a daily cycle lasts 24 hours and the system tries to simulate that situation. This way, a complete cycle in the system lasts 24 steps of time (from 0 to 23). The only rule in DMH serves out to activate neurons in the LHA each 24 hours (steps of time). For its part, neurons in LHA are actives during 16 hours (steps of time). Along this period the necessity of sleep is accumulated in the accumulator system. When 15 steps of time have been consumed a state switch is started, the system gets to sleep and the accumulator system sends 8 spikes to the regulator system (it uses the rule a15 →a8; 0). To use this rule means that 7 spikes are losts and sleep state only lasts 8 hours (steps of time). Now, the regulator system controls sleep state executing its rule a+/a →a; 0. Once the system is sleeping, neurons of the dense and extended group are actives and execute the rule a2/a →a; 0. Besides, this rule locks neurons in TMN, LC and DR groups during a step of time. When t= 23, DMH is open again and it emits 16 spikes to LHA. Then, the system comes back to wakefulness and the process is repeated. Notation:n: neural groups, t: steps of time, ACU: accumulator system, REG: regulator system, DMH: dorsomedial nucleus of the hypothalamus, VPO: dense group of the ventrolateral preoptic nucleus, eVLPO: extended group of the ventrolateral preoptic nucleus, LHA: neurons in the lateral hypothalamus, TMN: neurons in the tuberomammillary nucleus, LC: neurons in the locus coeruleus and DR: neurons in the dorsal and median raphe nuclei 5 Comments and New Issues A computational basic model for sleep-wake switch is examined in this paper. The goal was to model a neurophysiological process by means of a computational device such as Spiking Neural P Systems. Traditionally, this type of computational mechanisms include biological concepts in order to test its possibilities from a computational point of view but this paper tries to apply them in a well-known neural process. In order to apply SN P Systems to this process a new definition of extended rules was necessary and this work proposes a definition for inhibitory rules with the aim of modeling the inhibitory connections among neural groups Sleep-Awake Switch with Spiking Neural P Systems 69 HHHH H n t0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 ACU -a a2a3a4a5a6a7a8a9a10 a11 a12 a13 a14 a15 - - - - - - - - - a... RA1 REG ----------------a8!a7!a6!a5!a4!a3!a2!a! - - . . . RG1RG1RG1RG1RG1RG1RG1RG1 DMH amam−16 ... RD1B B B B B B B B B B B B B B B B B B B B B B RD1B B VLPO a a2a2a2a2a2a2a2a... BBBBBBBBBBBBBBBBBRV1RV1RV1RV1RV1RV1RV1 EVLPO a a2a2a2a2a2a2a2a... BBBBBBBBBBBBBBBBBRE1RE1RE1RE1RE1RE1RE1 EVLPO a a2a2a2a2a2a2a2a... BBBBBBBBBBBBBBBBBRE1RE1RE1RE1RE1RE1RE1 LHA a16!a15!a14 !a13!a12 !a11!a10!a9!a8!a7!a6!a5!a4!a3!a2!a! - - - - - - - - a16!a15! . . . RH1RH1RH1RH1RH1RH1RH1RH1RH1RH1RH1RH1RH1RH1RH1RH1RH1RH1 TMN a2a2a2a2a2a2a2a2a2a2a2a2a2a2a2a2a2a a ... RT1RT1RT1RT1RT1RT1RT1RT1RT1RT1RT1RT1RT1RT1RT1RT1RT1B B B B B B B LC a2a2a2a2a2a2a2a2a2a2a2a2a2a2a2a2a2a a ... RL1RL1RL1RL1RL1RL1RL1RL1RL1RL1RL1RL1RL1RL1RL1RL1RL1B B B B B B B DR a2a2a2a2a2a2a2a2a2a2a2a2a2a2a2a2a2a a ... RR1RR1RR1RR1RR1RR1RR1RR1RR1RR1RR1RR1RR1RR1RR1RR1RR1B B B B B B B Fig. 7. Results of SN P System simulating the sleep-awake switch 70 J.M. Mingo which are involved in the sleep-awake transitions. Starting from an SN P System with inhibitory rules a basic model is built with two important constraints: 1. the system supposes that an individual is slept or is awoken during fixed periods previously established. 2. the system splits the sleep control in two subsystems: accumulator and regulator. This supposition produces a delay in transition between states. Specifically, a step of time is lost from sleep to wakefulness and vice versa because accumulator and regulator subsystems can not start to perform at the same time. This is the reason why the rule a15 →a8; 0 in the accumulator uses 15 as number instead of 16. In spite of this constraints the system can represent in a clear way a complex biological process and it defines formally such process. This formally definition would let possible implementations once an appropiate software was developed. Clarity of reading is also an advantage feature of SN P Systems in front of other alternatives like mathematical models based on equations, for example. However, a number of interesting questions arises about the natural process and the basic model that it would be precise to analyze and solve in order to come near the model and the biological reality. •The sleep and wakefulness states with fixed period are a very simplification because, obviously, neither human beings nor animals maintain along their life an established amount of time on each activity. Moreover, an individual does not sleep at all a day or more if really he or she cannot do it. This kind of flexibility is difficult to implement in the basic model even in SN P Systems as they are currently defined. •Because the system considers one hour as step of time, the transition between sleep and wakefulness is completed slower than it occurs really in the biological process. Shorter steps can be defined in the model but then the model would be larger and tedious to show results in a table format. A possibility would be implement steps of time with different magnitudes but this question is not considered by currents SN P Systems. •The inhibitory mechanism proposed by the model has several constraints and their effects are easily visible looking at the table. When a transition sleepawake is produced either the regulator system or the LHA neuron groups fires (an exclamation symbol appears close to spikes). That way, the system represents that a state is active but, during one or two steps, the other main neuron groups in the dominant state does not fire because they still are locked by neurons in the opposite control system. For example, when t= 16 and t= 17, which is a sleep state, neurons in VLPO (dense and extended group) are locked (a B symbol appears in the square) because neurons in TMN, LC and DR are executing their rules yet. A similar situation occurs when t= 24 and t= 25 for the wakefulness state. This is a problem originated in the definition of the SN P System because a step of time is consumed since a spike gets out of a Sleep-Awake Switch with Spiking Neural P Systems 71 neuron and comes in another one. Maybe, a best definition of inhibitory rules should be attended. •Sleep or wakefulness states in the biological processes are a behavior but in the model they are simulated as firing in neuron groups. The model associates LHA neurons and regulator system as the output neurons in a way that is different of usual SN P systems. A best approach would be if a special cell or item in the SN P system performed this function. •SN P systems generally send and receive signals from a state to another one but the system performs like a closed system in the sense that everything must be in the system. This is a drawback in this case because, for example, the suprachiasmatic nucleus does not receive only inputs of the considered neurons but other items. The basic model has replaced this fact with the delay property that SN P systems incorporate and the axiom (m >> 16) to simulate a non finite execution but better ideas on this matter would be a great contribution. •SN P systems usually work with individual neurons and the basic model represents neuron groups as a single cell. A complex and interesting question to solve is concerned with how the model would be modified to work with several neurons on each group, for example, several LHA, TMN, LC, DR and VLPO cells. •Besides the sleep-wake switch another biological process that has been studied from a neurological point of view are the transitions between non-REM and REM sleep. This transition involves practically the same neuron groups but it adds more complexity and connections among items. Contributions about this topic would also be interesting. References 1. A. Alhazov, R. Freund, M. Oswald, M. Slavkovik: Extended Spiking Neural P Systems Generating Strings and Vectors of Non-Negative Integers. In H.J. Hoogeboom, Gh. P˘aun, G. Rozenberg, eds., Workshop on Membrane Computing, WMC7, Leiden, the Netherlands 2006, LNCS 4361, Springer, 2007, 123–134. 2. I.I. Ardelean, D. Besozzi: On Modelling Ion Fluxes Across Biological Membranes with P Systems. Proceedings of the third brainstorming week on Membrane Computing, 2005. 3. J.H. Benington, S.K. Koldali, H.C. Heller: Monoaminergic and cholinergic modulation of REM-sleep timing in rats. Brain Research, 681 (1995), 141–146. 4. N. Carlson: Physiology of the Behavior. 8th Edition. Addison Wesly, 2006. 5. W. Gerstner, W. Kistler: Spiking Neuron Models. Single Neurons, Populations, Plasticity. Cambridge Univ. Press, 2002. 6. M. Ionescu, Gh. P˘aun, T. Yokomori: Spiking Neural P Systems. Fundamenta Informaticae, 71, 2-3 (2006), 279–308. 7. W. Mass, C. Bishop, eds.: Pulsed Neural Networks. MIT Press, Cambridge 1999. 8. J.M. Mingo: Una Aproximacin al Interruptor del Sueo Mediante Spiking Neural P Systems (only in Spanish). Not published. 72 J.M. Mingo 9. M. Nakao, A. Karashima, N. Katayama: Mathematical models of regulatory mechanisms of sleep-wake rhythms. Cellular and Molecular Life Science, 64 (2007), 1236– 1243. 10. Gh. P˘aun: Computing with membranes. Journal of Computer and System Sciences, 61 (2000), 108–143. 11. A. Rechtschaffen, J. Siegel: Sleep and Dreaming. Principles of Neuroscience. Fourth Edition, Edited by E.R. Kandel, J.H. Schwartz and T.M. Jessel, McGraw-Hill, New York, 2000, 936–947. 12. C.B. Saper, T.C. Chou, T.E. Scammell: The sleep switch: hypothalamic control of sleep and wakefulness. Trends in Neurosciences, 24, 12 (2001), 726–731. 13. C.B. Saper, T.E. Scammell, J. Lu: Hypothalamic regulation of sleep and circadian rhythms. Nature, 437 (2005), 1257–1263. The Computational Complexity of Uniformity and Semi-uniformity in Membrane Systems Niall Murphy1and Damien Woods2 1Department of Computer Science, National University of Ireland, Maynooth, Ireland http://www.cs.nuim.ie/∼nmurphy/ [email protected] 2Department of Computer Science and Artificial Intelligence, University of Seville, Spain http://www.cs.us.es/∼dwoods/ [email protected] Summary. We investigate computing models that are presented as families of finite computing devices with a uniformity condition on the entire family. Examples include circuits, membrane systems, DNA computers, cellular automata, tile assembly systems, and so on. However, in this list there are actually two distinct kinds of uniformity conditions. The first is the most common and well-understood, where each input length is mapped to a single computing device that computes on the finite set of inputs of that length. The second, called semi-uniformity, is where each input is mapped to a computing device for that input. The former notion is well-known and used in circuit complexity, while the latter notion is frequently found in literature on nature-inspired computing models, from the past 20 years or so. Are these two notions distinct or not? For many models it has been found that these notion are in fact the same, in the sense that the choice of uniformity or semi-uniformity leads to characterisations of the same complexity classes. Here, we buck this trend and show that these notions are actually distinct: we give classes of uniform membrane systems that are strictly weaker than their semi-uniform counterparts. This solves a known open problem in the theory of membrane systems. 1 Introduction In his famous 1984 paper on DNA computing [1], Adleman mapped a specific instance of the travelling salesman problem (TSP) to a set of DNA strands, and then used well-known biomolecular techniques to solve the problem. To assert generality for his algorithm, one would define a (simple) mapping from arbitrary TSP instances to sets of DNA strings. Then, in order to claim that this mapping is not doing the essential computation, it would have to be easily computable 74 N. Murphy, D. Woods (e.g. logspace computable). Circuit uniformity provides a well-established framework where we map each input length n∈Nto a circuit cn∈C, with a suitably simple mapping. However, Adleman did something different, he mapped a specific instance of the problem to a computing device. We call this notion semi-uniformity, and in fact a large number of computation models use semi-uniformity. This raises the immediate question of whether the notions of uniformity and semi-uniformity are equivalent. It has been shown in a number of models that whether one chooses to use uniformity or semi-uniformity does not affect the power of the model. However, in this paper we show that these notions are not equivalent. We prove that choosing one notion over another gives characterisations of completely different complexity classes, including known distinct classes. We prove this result for a computational model called membrane systems (also known as P-systems) [15]. Membrane computing is a branch of natural computing which defines computation models that are inspired by the structure and function of living cells. The membrane computing model is sufficiently formal that this question can be clearly stated, e.g. it is stated as Open Problem C in [16]. Why is this result surprising? Every class of problems solved by a uniform family of devices is contained in the analogous semi-uniform class, since one is a restriction of the other. However, in all membrane system models studied to date, the classes of problems solved by semi-uniform and uniform families turned out to be equal [4, 10, 19]. Specifically, if we want to solve some problem, by specifying a family of membrane systems (or some other model), it is often much easier to first use the more general notion of semi-uniformity, and then subsequently try to find a uniform solution. In almost all cases where a semi-uniform family was given for some problem [3, 11, 13, 19], at a later point a uniform version of the same result was published [2, 4, 13]. Here we prove that this improvement is not always possible. Since our main result proves something general about families of finite devices we would hope that, in the future, it can be applied to other computational models, besides membrane systems. Why? Firstly, our results are proved by converting the membrane system into a directed acyclic graph. Input acceptance is then rephrased as a graph reachability problem and this gives a very general tool that can be applied to other computational models (where we can find analogous graph representations). Secondly, the result concerns a general concept (uniformity/semiuniformity) that is independent of particular formalisms. Besides membrane systems and circuits, some other models that use notions of uniformity and semiuniformity include families of neural networks, molecular and DNA computers, tile assembly systems and cellular automata [6, 8, 12, 17, 18]. Our results could conceivably be applied to these models. We now briefly observe what happens when we relate the notion of semiuniformity to circuit complexity. We can easily define semi-uniformity for circuits. If the complexity class of the semi-uniformity function contains the prediction problem for circuits in the resulting family, then the semi-uniformity condition Uniformity and Semi-uniformity in Membrane Systems 75 Characterisation Power of (semi-)uniformity condition AC0 AC0 NC1 NC1 L L NL NL NC2 NC2 P P NP NP PSPACE PSPACE Fig. 1. Complexity classes that are characterised by the membrane systems studied in this paper. Characterisations by uniform systems are denoted by , and semi-uniform by . For example, Theorem 1 is illustrated by the fact that AC0-uniform systems characterise AC0, and that AC0-semi-uniform systems characterise NL. The previously known P[10] (indicated by ) and PSPACE [4, 20] (indicated by ) results, where semi-uniform and uniform classes have the same power are also shown. characterises the power of the model. If the semi-uniformity function is computable in the class that is characterised by the prediction problem for circuits in the resulting family, then we get the known characterisations for the analogous uniformity condition. However, in the uniform case it is not obvious what happens when we increase the uniformity beyond the power of the circuit, for example P-uniform AC0=AC0is an open problem [5]. Furthermore we should note that the uniformity condition in membrane systems preprocesses the input (as well as creating the device) and so is a seemingly different notion than circuit uniformity. If we add an analogous preprocessing step to circuits we see similar results as proven here for membrane systems: as soon as the preprocessing goes beyond the power of the circuit, we can ignore the circuit and let the preprocessing solve the problem. With preprocessing below the power of the circuit, the answer depends on the particular circuit model. In fact, if we restrict ourselves to polynomially sized circuits with only OR gates, we would see analogous results to those presented here, (i.e. our work shows that this circuit model is computationally equivalent to the AM0 −d membrane systems discussed in this work). 82 N. Murphy, D. Woods Π0that is equivalent in all respects except that there are exactly two input objects some extra rules. In particular, Π0accepts input xiff Πdoes. The previous argument permits us to consider only those systems that have two input objects {a,b}. Thus we restrict attention to the case that the input encoding function is of the form e:X→ {a,b}. We say that eis a characteristic function with range {a,b}. Let (e, f) –Pbe a (e, f)-uniform family of confluent recogniser membrane systems from the class AM0 −dthat solves problem X. We claim that there exists a family (e, f0)–Pthat also solves Xbut uses a uniformity function f0that produces membrane systems of a very restricted form. Consider the dependency graph of the membrane system f(n) = ΠX(n). In terms of reachability from {a,b}to {yes,no}in this graph (which corresponds to accepting/rejecting in the membrane system), the essential property is whether path(a,yes) is true or path(a,no) is true. However, this essential property is captured by the following (extremely simple) pair of dependency graphs: G1= (VG, EG1) and G2= (VG, EG2) where VG={a,b,yes,no},EG1={(a,yes),(b,no)}, and EG2={(a,no),(b,yes)}. Therefore if there is a family (e, f) –Pthat solves X, where frepresents valid, but arbitrary, dependency graphs, then there is another family of the form (e, f0) –P that also solves Xand is identical in every way except that f0represents dependency graphs of the restricted form just described. Now we prove the upper bound (e, f)–PMCAM0 −d⊆E, where e,fare respectively computable in E,F, with F⊆E, and the classes E,Fare as given in the statement. As we have just shown, for each problem Xin the class (e, f)–PMCAM0 −dthere is a family of the restricted form (e, f0) –Pthat solves X. To simulate (e, f0) –Pwith a given input x∈Xwe compute the pair e(x), f0(1|x|). The range of the function f0is two membrane systems, which correspond to the two (restricted) dependency graphs G1and G2from above, and whose reachability problems are (trivially) in the weakest Ethat we consider (E=AC0). Furthermore, since we only consider the case where F⊆Ethen we know that f0itself is computable in Eand we have (e, f)–PMCAM0 −d⊆E. The lower bound E⊆(e, f)–PMCAM0 −dis easy to show. We use the fact, shown above, that eis a characteristic function with access to the input word. Thus the following simple family computes any problem from E: function e(x) outputs aif xis a positive instance of Xand bif xis a negative instance of X, and fsimply maps ato yes and bto no.ut Acknowledgements Niall Murphy is funded by the Irish Research Council for Science, Engineering and Technology. Damien Woods is supported by a Project of Excellence from the Junta de Andaluca, grant number TIC-581. Uniformity and Semi-uniformity in Membrane Systems 83 References 1. Leonard Adleman. Molecular computation of solutions to combinatorial problems. Science, 266:1021–1024, 1994. 2. Artiom Alhazov, Carlos Mart´ın-Vide, and Linqiang Pan. Solving a PSPACE- complete problem by recognizing P Systems with restricted active membranes. Fundamenta Informaticae, 58(2):67–77, 2003. 3. Artiom Alhazov and Linqiang Pan. Polarizationless P Systems with active membranes. Grammars, 7:141–159, 2004. 4. Artiom Alhazov and Mario J. Prez-Jimnez. Uniform solution to QSAT using polarizationless active membranes. In Jrme Durand-Lose and Maurice Margenstern, editors, Machines, Computations and Universality (MCU), volume 4664 of LNCS, pages 122–133, Orlans, France, September 2007. Springer. 5. Eric Allender. Applications of time-bounded Kolmogorov complexity in complexity theory. In Osamu Watanabe, editor, Kolmogorov Complexity and Computational Complexity, chapter 1, pages 4–22. Springer, 1992. 6. Martyn Amos. Theoretical and Experimental DNA Computation. Natural Computing Series. Springer, 2005. 7. David A. Mix Barrington, Neil Immerman, and Howard Straubing. On uniformity within NC1.Journal of Computer and System Sciences, 41(3):274–306, 1990. 8. Matthew Cook, David Soloveichik, Erik Winfree, and Jehoshua Bruck. Programmability of chemical reaction networks. Technical report, Caltech Parallel and Distributed Systems Group [http://caltechparadise.library.caltech.edu/perl/oai2] (United States), 2008. In submission. 9. Miguel A. Gutirrez-Naranjo, Mario J. Prez-Jimnez, Agustn Riscos-Nez, and Francisco J. Romero-Campero. Computational efficiency of dissolution rules in membrane systems. International Journal of Computer Mathematics, 83(7):593–611, 2006. 10. Niall Murphy and Damien Woods. Active membrane systems without charges and using only symmetric elementary division characterise P. 8th International Workshop on Membrane Computing, LNCS, 4860:367–384, 2007. 11. Niall Murphy and Damien Woods. A characterisation of NL using membrane systems without charges and dissolution. Unconventional Computing, 7th International Conference, UC 2008, LNCS, 5204:164–176, 2008. 12. Ian Parberry. Circuit complexity and neural networks. MIT Press, 1994. 13. Mario J. Prez-Jimnez, Alvaro Romero-Jimnez, and Fernando Sancho-Caparrini. Complexity classes in models of cellular computing with membranes. Natural Computing, 2(3):265–285, 2003. 14. Gheorghe Pun. P Systems with active membranes: Attacking NP-Complete problems. Journal of Automata, Languages and Combinatorics, 6(1):75–90, 2001. 15. Gheorghe Pun. Membrane Computing. Springer-Verlag, Berlin, 2002. 16. Gheorghe Pun. Further twenty six open problems in membrane computing. In Proceedings of the Third Brainstorming Week on Membrane Computing, Sevilla (Spain), pages 249–262. Fnix Editoria, January 2005. 17. David Soloveichik and Erik Winfree. The computational power of Benenson automata. Theoretical Computer Science, 344:279–297, 2005. 18. David Soloveichik and Erik Winfree. Complexity of self-assembled shapes. SIAM J. Comput., 36(6):1544–1569, 2007. 19. Petr Sosk. The computational power of cell division in P systems: Beating down parallel computers? Natural Computing, 2(3):287–298, August 2003. 84 N. Murphy, D. Woods 20. Petr Sosk and Alfonso Rodrguez-Patn. Membrane computing and complexity theory: A characterization of PSPACE. Journal of Computer and System Sciences, 73(1):137–152, 2007. Appendix A A.1 Constructing dependency graphs in AC0 We are given a binary string xthat encodes a membrane system, Π. To make a dependency graph from a membrane system requires a constant number of parallel steps that are as follows. First, a row of circuits identifies all communication (type (b) and (c)) rules and uses the (static) membrane structure to determine the correct parent membranes, then writes out (a binary encoding of) edges representing these rules. Next, a row of circuits writes out all edges representing division (type (e)) rules, for example a rule [ a]h→[b][ c] becomes the edges (ah, bh),(ah, ch) in the dependency graph. In the final step we deal with evolution (type (a) rules) where it is possible to have polynomially many copies of polynomially many distinct objects on the right hand side of a rule (e.g. [ a]h→[bcbbcdee · · · z]h). To write out edges for these rules in constant time we take advantage of the fact that we require at most one edge for each object-membrane pair in O×H. We have a circuit for each element of {oh|o∈O, h ∈H}. The circuit for ohtakes as input (an encoding of) all rules in Rwhose left hand side is of the form [o]h. The circuit then, in a parallel manner, masks (an encoding of) the right hand side of the rule (for example [bbcdc]h) with the encoding of each object in O, (in the example, masking for (encoded) bwould produce (encoded) bb000). All encoded objects in the string are then ORed together so that if there was at least one copy of that object in the system we obtain a single instance of it. The circuit being unique for a specific left hand side [o]hnow writes out an encoding of the edge (oh, bh) and an encoding of all other edges for objects that existed on the right hand side of this rule in parallel. Structured Modeling with Hyperdag P Systems: Part A Radu Nicolescu, Michael J. Dinneen, Yun-Bum Kim Department of Computer Science, University of Auckland Private Bag 92019, Auckland, New Zealand [email protected] Summary. P systems provide a computational model based on the structure and interaction of living cells [9]. A P system consists of a hierarchical nesting of cell-like membranes, which can be visualized as a rooted tree. Although the P systems are computationally complete, many real world models, e.g., from socio-economic systems, databases, operating systems, distributed systems, seem to require more expressive power than provided by tree structures. Many such systems have a primary tree-like structure completed with shared or secondary communication channels. Modeling these as tree-based systems, while theoretically possible, is not very appealing, because it typically needs artificial extensions that introduce additional complexities, nonexistent in the originals. In this paper we propose and define a new model that combines structure and flexibility, called hyperdag P systems, in short, hP systems, which extend the definition of conventional P systems, by allowing dags, interpreted as hypergraphs, instead of trees, as models for the membrane structure. We investigate the relation between our hP systems and neural P systems. Despite using an apparently less powerful structure, i.e., a dag instead of a general graph, we argue that hP systems have essentially the same computational power as tissue and neural P systems. We argue that hP systems offer a structured approach to membrane-based modeling that is often closer to the behavior and underlying structure of the modeled objects. Additionally, we enable dynamical changes of the rewriting modes (e.g., to alternate between determinism and parallelism) and of the transfer modes (e.g., the switch between unicast or broadcast). In contrast, classical P systems, both tree and graph based P systems, seem to focus on a statical approach. We support our view with a simple but realistic example, inspired from computer networking, modeled as a hP system with a shared communication line (broadcast channel). In Part B of this paper we will explore this model further and support it with a more extensive set of examples. 86 R. Nicolescu, M.J. Dinneen, Y.-B. Kim 1 Introduction P systems provide a distributed computational model, based on the structure and interaction of living cells, first introduced by G. P˘aun in 1998 [8]. The model was initially based on transition rules, but was later expanded into a large family of related models. Essentially, all versions of P systems have a structure consisting of cell-like membranes and a set of rules that govern their evolution over time. Many of the “classical” versions use a structure where membranes correspond to nodes in a rooted tree. Such a structure is often visualized as Venn diagram where nesting denotes a parent/child relationship. For example, Figure 1 [10] shows the same P system structure with 9 membranes, labeled as 1, . . . , 9, both as a rooted tree and as a Venn diagram. 2 3 4 5 6 7 8 9 1 1 2 3 45 6 7 8 9 Fig. 1. A P system structure represented as a tree and as a Venn diagram. More, recently, tissue P systems [7] and neural P systems [9], here abbreviated as tP and nP systems, respectively, have been introduced, partially to overcome the limitations of the tree model. Essentially, these systems organize their cells in an arbitrary digraph. For example, ignoring for the moment the actual contents of cells (states, objects, rules), Figure 2 illustrates the membrane structure of a simple tP or nP system, consisting of 3 cells, σ1, σ2, σ3, where cell σ1is designated as the output cell. A large variety of rules have been used to describe the operational behavior of P systems, the main ones being: multiset rewriting rules, communication rules and membrane handling rules. Essentially, transition P systems and nP systems use multiset rewriting rules, P systems with symport/antiport operate by communicating immutable objects, P systems with active membranes combine all three type rules. For a comprehensive overview and more details, we refer the reader to [9, 10]. Structured Modeling with Hyperdag P Systems 87 σ1 σ2σ3 Fig. 2. A tP/nP system structure represented as a digraph. Besides theoretical computer science and biology, P systems have been applied to a variety of other domains, ranging from linguistics [5] to theoretically efficient solutions of NP-complete problems [14], or to model distributed algorithms [3, 6]. The underlying tree structure provides good support for reasoning and formal verification, good potential for efficient implementation on multi-core architectures, and an excellent visualization, very appealing to practitioners. Although the P systems are computationally complete, many real world models seem to require more expressive power, essentially trees augmented by shared or secondary communication channels. For example, the notion of a processing node having an unique parent is not true for (a) computer networks where a computer could simultaneously be attached to several subnets (e.g., to an Ethernet bus and to a wireless cell), (b) living organisms may be the result of multiple inheritance (e.g., the evolutionary “tree” is not really a tree, because of lateral gene transfer [4]) and (c) socio-economic scenarios where a player is often connected to and influenced by more than one factors [11, 12, 13]. Modeling these as tree-based systems, while theoretically possible, is not very appealing. Simulating shared or secondary channels requires artificial mechanisms that will ripple data up and down the tree, via a common ancestor. This could of course limit the merits of using a formal model. Tissue and neural P systems have been introduced to model such cases [7, 9]; details on neural P systems, in short, nP systems, are given in Section 3. However, these extensions are based on general graphs and, while allowing any direct communications, they also tend to obscure the structures already present in the modeled objects, limiting the advantages that a more structured approach could provide. Verification is more difficult without a clear modularization of concerns, practical parallel implementation could be less efficient, if the locality of reference is not enforced, and visualizations are not very meaningful, unless the primary structure is clearly emphasized. We do not think that we have to choose between structure and flexibility. We propose a solution that seems to combine both, i.e., flexibility without sacrificing the advantages of a structured approach. 88 R. Nicolescu, M.J. Dinneen, Y.-B. Kim Our main contribution in this paper is to propose a new model for P systems, called hyperdag P systems, in short, hP systems, that allows more flexible communications than tree-based models, while preserving a strong hierarchical structure. This model, defined in Section 4, (a) extends the tree structure of classical P systems to directed acyclic graphs (dags), (b) augments the operational rules of nP systems with broadcast facilities (via a go-sibling transfer tag), and (c) enables dynamical changes of the rewriting modes (e.g., to alternate between determinism and parallelism) and of the transfer modes (e.g., to switch between unicast or broadcast). In contrast, classical P systems, both tree and graph based P systems, seem to focus on a statical approach. We investigate the relation between our hP systems and neural P systems. Despite using an apparently less powerful structure, we show in Section 5 that our simple dag model has the same computational power as graph-based tissue and neural P systems. We argue that hP systems offer a structured approach to membrane-based modeling that is often closer to the behavior and underlying structure of the modeled objects. Because our extensions address the membrane topology, not the rules model, they can be applied to a variety of P system flavors, including transition systems and symport/antiport systems. We support our view with a realistic example (see Examples 8 and 9), inspired from computer networking, modeled as a hP system with a shared communication line (broadcast channel). Classical P systems allow a “nice” planar visualization, where the parent/child relationships between membranes are represented by Venn-like diagrams. We show in Section 6 that the extended membrane structure of hP systems can still be visualized by hierarchically nested planar regions. In this article we will restrict ourselves to P systems based on multiset rewriting rules, such as used by transition P systems and nP systems. However, because our extensions address the membrane topology, not the rules model, they can be applied to a variety of other P system flavors. 2 Preliminaries A (binary) relation Rover two sets Xand Yis a subset of their Cartesian product, R⊆X×Y. For A⊆Xand B⊆Y, we set R(A) = {y∈Y| ∃x∈A, (x, y)∈R}, R−1(B) = {x∈X| ∃y∈B, (x, y)∈R}. Adigraph (directed graph) Gis a pair (X, A), where Xis a finite set of elements called nodes (or vertices), and Ais a binary relation A⊆X×X, of elements called arcs. For an arc (x, y)∈A,xis a predecessor of yand yis a successor of x. A length n−1path is a sequence of ndistinct nodes x1, . . . , xn, such that {(x1, x2), . . . , (xn−1, xn)} ⊆ A. A cycle is a path x1, . . . , xn, where n≥1 and (xn, x1)∈A. Adag (directed acyclic graph) is a digraph (X, A) without cycles. For x∈ X,A−1(x) = A−1({x}) are x’s parents,A(x) = A({x}) are x’s children, and Structured Modeling with Hyperdag P Systems 89 A(A−1(x))\{x}=A(A−1({x}))\{x}are x’s siblings (siblings defines a symmetric relation). A node x∈Xis a source iff |A−1(x)|= 0, and x∈Xis a sink iff |A(x)|= 0. The height of a node xis the maximum length of all paths from xto a sink node. An arc (x, y) is transitive if there exists a path x1, . . . , xn, with x1=x, xn=yand n > 2. dags without transitive arcs are here called canonical. A (rooted unordered) tree is a dag with exactly one source, called root, and all other nodes have exactly one parent (predecessor). Sinks in a tree are also called leaves. A topological order of a dag is a linear reordering of vertices, in which each vertex xcomes before all its children vertices A(x). Dags and trees are typically represented with parent-child arcs on the top-down axis, i.e., sources/roots up and sinks/leaves down. Figure 3 shows a simple dag. Fig. 3. A simple dag. The parent-child axis is up-down. Here, plain lines indicate parentchild relations and dashed lines indicate siblings. We consider a variant hypergraph definition, based on multisets, as an extension of the classical definition [1], which is based on sets. A hypergraph is here a pair (X, E), where Xis a finite set of elements called nodes (or vertices), and Eis a finite multiset of subsets of X, i.e., e∈E⇔e⊆X. By using a multiset of edges, instead of a more conventional set of edges, we introduce an intensional element, where two extensionally equivalent hyperedges (i.e., hyperedges containing the same nodes) are not necessarily equal. A graph is a set based hypergraph where hyperedges are known as edges and contain exactly two nodes. Alternatively, a graph (X, E) can be interpreted as a digraph (X, A), where A={(x, y)| {x, y} ∈ E}. Hypergraphs (set or multiset based) can be represented by planar diagrams, where hyperedges are represented as regions delimited by images of Jordan curves (simple closed curves) [2]. With the above hypergraph definition, a height 1 dag (X, A) can be interpreted as a hypergraph (X, E), where Eis the multiset E={A(x)| |A−1(x)|= 0}. For example, Figure 4 represents, side by side, the dag D= ({a, b, c, d, e, f}, {(d, a),(d, b),(d, c),(e, b),(e, c),(f, b),(f, c)}) and its corresponding hypergraph H= ({a, b, c},{d, e, f}), where d={a, b, c}, e ={b, c}, f ={b, c}. Note that the apparently empty differences of regions are needed in the case of multiset based hypergraphs, to support the intensional (as opposed to the extensional) aspect: 90 R. Nicolescu, M.J. Dinneen, Y.-B. Kim here e6=f, despite containing the same nodes, band c, and neither enor fis included in d. d e a b c a b c d e f f Fig. 4. A simple height 1 dag and its corresponding hypergraph representation. Generalizing the above hypergraph definition, a height ngeneralized hypergraph is a system (X, E), recursively built via a sequence of nhypergraphs (X1, E1), . . . , (Xn, En) where X1=X,Xi∩Ei=∅, Xi+1 =Xi∪Ei,e∩Ei6= ∅for ∀e∈Ei+1 and E=Si∈{1,...,n}Ei. An arbitrary height ndag can be represented by a height ngeneralized hypergraph, where the hypergraph nodes correspond to dag sinks, and height ihyperedges correspond to height idag nodes, for i∈ {1, . . . , n}. We will later see that any generalized hypergraph that corresponds to a nontransitive dag can also be represented by hierarchically nested planar regions delimited by Jordan curves, where arcs are represented by direct nesting. For example, Figure 5 shows a height 2 dag and its corresponding height 2 hypergraph (X, E), where X=X1={a, b, c, d, e}, E1={f, g, h}, E2={i}, E ={f, g, h, i}. a b c d e a b c d e fgh i f g h i Fig. 5. A height 2 dag and its corresponding height 2 hypergraph. An alphabet Ois a finite non-empty sets of objects. We will assume that these alphabets are implicitly ordered. Multisets over an alphabet Oare represented as strings over O, such as on1 1. . . onk k, where oi∈O,ni≥0, and, in the canonical form, letters appear in sorted order, i.e., o1<· · · < ok, and ni≥1. The set of all multisets is denoted by O∗. For this representation, two strings are equivalent if they become equal after sorting, e.g., a2cbd0aand a3bc are equivalent representations of the same multiset {a, a, a, b, c}. Under this convention, the empty string λ Structured Modeling with Hyperdag P Systems 91 represents the empty multiset, and string concatenation represents multiset union, e.g., (a2c)·(ab) = a3bc. 3 Neural P Systems In this paper we present the definition of neural P systems as given in [9], that coincide with an early definition of tissue P systems as given in [7]. We define the following sets of tagged objects: Ogo ={(a, go)|a∈O},Oout ={(a, out)|a∈ O}, and we set Otot =O∪Ogo ∪Oout. For simplicity, we will use subscripts for these tagged objects, such as ago for (a, go) and aout for (a, out). We also define projection homomorphisms, here denoted in postfix notation: |O,|go,|out :O∗ tot → O∗, by o|O=o, ogo|go =o, oout|out =ofor o∈O, and otherwise λ. For example, a2a3 gob4bgo|go =a3b. Definition 1 (Neural P systems [9, 7]). Aneural P system (of degree m≥1) is a system: Π= (O, σ1, . . . , σm, syn, iout), where: 1. Ois an ordered finite non-empty alphabet of objects; 2. σ1, . . . , σmare cells, of the form σi= (Qi, si,0, wi,0, Pi),1≤i≤m, where: •Qiis a finite set (of states), •si,0∈Qiis the initial state, •wi,0∈O∗is the initial multiset of objects, •Piis a finite set of multiset rewriting rules of the form sx →s0x0ygozout, where s, s0∈Qi,x, x0∈O∗,ygo ∈O∗ go and zout ∈O∗ out, with the restriction that zout =λfor all i∈ {1, . . . , m}\{iout}. 3. syn is a set of digraph arcs on {1, . . . , m}, i.e., syn ⊆ {1, . . . , m}×{1, . . . , m}, representing unidirectional communication channels between cells, known as synapses; 4. iout ∈ {1, . . . , m}indicates the output cell, the only cell allowed to send objects to the “environment”. Example 1. To illustrate the operational behavior of nP systems, consider again the example of Figure 2, expanded now with states, rules and objects, see Figure 6. For simplicity, in this example only cell σ1provides rules. More formally, this nP system can be defined as the system Π1= (O, σ1, σ2, σ3, syn, iout), where: •O={a, b, c, d}; •σ1= ({s, t}, s, a2,{sa →sdbgocgo, sa →sd, s →t, td →tdout}); •σ2= ({s}, s, λ, ∅); •σ3= ({s}, s, λ, ∅); •syn ={(1,2),(1,3),(2,3),(3,1)}; •iout = 1. 98 R. Nicolescu, M.J. Dinneen, Y.-B. Kim •sx ⇒min s0yiff sw →s0w0∈Pi,w⊆x, and y= (x−w)∪w0; •sx ⇒par s0yiff sw →s0w0∈Pi,wk⊆x,wk+1 *x, for some k≥1, and y= (x−wk)∪w0k; •sx ⇒max s0yiff sw1→s0w0 1, . . . , swk→s0w0 k∈Pi,k≥1, such that w1. . . wk⊆x, y = (x−w1. . . wk)∪w0 1. . . w0 k, and there is no sw →s0w0∈Pi, such that w1. . . wkw⊆x(note that rules can be combined only if they start from the same state sand end in the same state s0). Definition 11 (Transition steps). Given two configurations C1= (s1w1, . . . , smwm)and C2= (s0 1w00 1, . . . , s0 mw00 m), we write C1⇒α,β C2, for αand β(as defined in Definition 9) if the conditions below are met. First, we apply rewriting steps (as defined in Definition 10) on each cell, i.e., siwi⇒α(i,si)s0 iw0 i,1≤i≤m. Secondly, we define z↑ j,k,z↓ j,k,z↔ j,k, the outgoing multisets from jto k, where j∈ {1, . . . , m}and, respectively, k∈δ−1(j),k∈δ(j),k∈δ(δ−1(j))\{j}: •If β(j, sj) = repl, then ◦z↑ j,k =w0 j|↑, for k∈δ−1(j); ◦z↓ j,k =w0 j|↓, for k∈δ(j); ◦z↔ j,k =w0 j|↔, for k∈δ(δ−1(j))\{j}. •If β(j, sj) = one, then ◦z↑ j,kj=w0 j|↑, for an arbitrary kj∈δ−1(j), and z↑ j,k =λfor k∈δ−1(j)\{kj}; ◦z↓ j,kj=w0 j|↓, for an arbitrary kj∈δ(j), and z↓ j,k =λfor k∈δ(j)\{kj}; ◦z↔ j,kj=w0 j|↔, for an arbitrary kj∈δ(δ−1(j))\{j}, and z↔ j,k =λfor k∈ δ(δ−1(j))\{j, kj}. •If β(j, sj) = spread, then ◦ {z↑ j,k}k∈δ−1(j)is an arbitrary multiset partition of w0 j|↑; ◦ {z↓ j,k}k∈δ(j)is an arbitrary multiset partition of w0 j|↓; ◦ {z↔ j,k}k∈δ(δ−1(j))\{j}is an arbitrary multiset partition of w0 j|↔. Finally, we set w00 i=w0 i|O∪[ j∈δ−1(i) z↑ j,i ∪[ j∈δ(i) z↓ j,i ∪[ j∈δ(δ−1(i))\{i} z↔ j,i, for i∈ {1, . . . , m}. Definition 12 (Halting and results). If no more transitions are possible, the hP system halts. For halted hP system, the computational result is the multiset that was cumulatively sent out (to the “environment”) from the output cells Iout. The numerical result is the set of vectors consisting of the object multiplicities in the multiset result. Example 7. As examples, consider two hP systems, Π3and Π4, both functional equivalent a functional equivalent of the earlier Π2nP system. Π3= (O, σ1, σ2, σ3, δ, Iout), see Figure 9, where: Structured Modeling with Hyperdag P Systems 99 •O={a}; •σ1= ({s}, s, a, {sa →sa↓, sa →saout}); •σ2= ({s}, s, λ, {sa →sa↑}); •σ3= ({s}, s, λ, {sa →sa↑}); •δ={(1,2),(1,3)}; •Iout ={1}. s, a sa →sa↓ sa →saout s, λ sa →sa↑ s, λ sa →sa↑ σ1 σ2σ3 Fig. 9. Π3, a simple hP system (equivalent to the Π2nP system of Figure 7). Π4= (O, σ1, σ2, σ3, σ4, σ5, δ, Iout), see Figure 10, where: •O={a}; •σ1= ({s}, s, a, {sa →sa↔, sa →saout}); •σ2= ({s}, s, λ, {sa →sa↔}); •σ3= ({s}, s, λ, {sa →sa↔}); •σ4= ({s}, s, λ, ∅); •σ5= ({s}, s, λ, ∅); •δ={(4,1),(4,2),(5,1),(5,3)}; •Iout ={1}. 5 Relations Between P Systems, Neural P Systems and Hyperdag P Systems Theorem 1 (Hyperdag P systems include non-dissolving P systems). Any non-dissolving transition P system can be simulated by a hP system, with the same number of steps and object transfers. 100 R. Nicolescu, M.J. Dinneen, Y.-B. Kim s, a sa →sa↔ sa →saout s, λ sa →sa↔ s, λ sa →sa↔ σ1 σ2σ3 σ4σ5 s, λ s, λ Fig. 10. Π4, another simple hP system (equivalent to the Π2nP system of Figure 7). Proof. Given a non-dissolving, transition P system Π[10], we build a functionally equivalent hP system Hby the following transformation f. Essentially, we use the same elements, with minor adjustments. As the underlying structure, we can reuse the rooted tree structure of the P systems, because any rooted tree is a dag. Π= (O, C, µ, w1, . . . , wm, R1, . . . , Rm, io), f(Π) = (O0, σ0 1, . . . , σ0 m, δ, Iout). •O0=O; •σ0 i= (Q0 i, s0 i,0, w0 i,0, P0 i), 1 ≤i≤m, where: ◦Q0 i={s}, where sis any symbol /∈O; ◦s0 i,0=s; ◦w0 i,0=wi; ◦P0 i={su →sv0|u→v∈Ri}, where v0is a translation of vby the following homomorphism: (O∪O×Tar)∗→O∗, such that a→a, (a, here)→a, (a, out)→a↑, (a, in)→a↓; •δ=µ; •Iout ={io}; •The object rewrite mode is the max constant function, i.e., α(i, s) = max, for i∈ {1, . . . , m}, s ∈Qi; •The object transfer mode is the spread constant function, i.e., β(i, s) = spread, for i∈ {1, . . . , m}, s ∈Qi. Tags go-child(↓), go-parent(↑) correspond to P system target indications in, out, respectively. An empty tag corresponds to P system target indication here. Object rewrite and transfer modes of hP systems are a superset of object rewrite and transfer mode of P systems. We omit here the rest of the proof which is now straightforward but lengthy. Remark 1. We leave open the case of dissolving P systems, which can be simulated, but not properly subsumed by hP systems. Proving that hP systems also cover nP systems appears more daunting. However, here we will use a natural interpretation of hP systems, where the bulk of the computing will be done by the sink nodes, and the upper nodes (parents) will function mostly as communication channels. Structured Modeling with Hyperdag P Systems 101 Remark 2. The combination of go-sibling (↔) with repl object transfer mode enable the efficient modeling of a communication bus, using only one hyperedge or, in the corresponding dag, narcs. In contrast, any formal systems that use graph edges (or digraph arcs) to model 1:1 communication channels will need n(n−1) separate edges (or 2n(n−1) arcs) to model the associated complete subgraph (clique). It is expected that this modeling improvement will also translate into a complexity advantage, if we use the number of messages measure. In hP systems, a local broadcast needs only one message to siblings, while needing n−1 messages in graph/digraph based systems. Example 8. Figure 11 shows the structure of an hP system that models a computer network. Four computers are connected to “Ethernet Bus 1”, the other four computers are connected to “Ethernet Bus 2”, while two of the first group and two of the second group are at the same time connected to a wireless cell. In this figure we also suggest that “Ethernet Bus 1” and “Ethernet Bus 2” are themselves connected to a higher level communication hub, in a generalized hypergraph. Example 9. Figure 12 shows the computer network of Figure 11, modeled as a graph (if we omit arrows) or as a digraph (if we consider the arrows). Note that the graph/digraph models, such as nP, do not support the grouping concept, i.e., there is no direct way to mark the nodes a, b, c, d as being part of the “Ethernet Bus 1”, etc. We can now sketch the proof of the theorem comparing hP systems and nP systems. Theorem 2 (Hyperdag P systems can simulate bidirectional nP systems). Any bidirectional nP system can be simulated by a hP system, with the same number of steps and object transfers. Proof. Given a bidirectional nP system Π, we build a functionally equivalent hP system Hby the following transformation f. As the underlying structure, we use a dag of height 1, where the cells are sink nodes, and the syn arcs are reified as height 1 nodes. Without loss of generality, we assume that in the nP systems synapses are distinct from cells. Π= (O, σ1, . . . , σm, syn, iout), f(Π) = (O0, σ0 1, . . . , σ0 m+|syn|, δ, Iout). •O0=O; •σ0 i=f(σi), for i∈ {1, . . . , m}, where: ◦Q0 i=Qi; ◦s0 i,0=si,0; ◦w0 i,0=wi,0; ◦P0 i={u→v0|u→v∈Pi}, where v0is a translation of vby the following homomorphism: O∗ tot →O∗, such that a→a,ago →a↔,aout →aout ; 102 R. Nicolescu, M.J. Dinneen, Y.-B. Kim Ethernet Bus 1 Ethernet Bus 2 Wireless Bus Ethernet Bus 1 Ethernet Bus 2 a b c d e f g h a b c d e f g h Wireless Bus Fig. 11. A computer network and its corresponding hP/hypergraph representation. a b c d e f g h (Ethernet Bus 1 ) (Ethernet Bus 2 )(Wireless Bus) Fig. 12. The graph/digraph representations of the computer network of Figure 11. •σ0 m+1, . . . , σ0 m+|syn|is an arbitrary ordering of elements in syn; •δ={(e, u),(e, v)|e= (u, v)∈syn}; •Iout ={iout}; •The object rewrite mode is a constant function, i.e., α(i, s) = α0, for i∈ {1, . . . , m +|syn|}, s ∈Qi, where α0∈ {min, par, max}; Structured Modeling with Hyperdag P Systems 103 •The object transfer mode is a constant function, i.e., β(i, s) = β0, for i∈ {1, . . . , m +|syn|}, s ∈Qi, where β0∈ {repl, one, spread}. Here the nodes corresponding to synapses are inactive as they just link neighboring cells. Essentially, the cells keep their original nP rules, with minor adjustments. A go tag in nP rules corresponds to a sibling(↔) tag in hP rules. Object rewrite and transfer modes of hP systems are a superset of object rewrite and transfer modes of nP systems. We omit here the rest of the proof which is now straightforward but lengthy. Remark 3. We leave open the case of non-bidirectional nP systems, which can be simulated, but not properly subsumed by hP systems. 6 Planar Representation of hP Systems Classical tree-based P systems allow a “nice” planar representation, where the parent/child relationships between membranes are represented by Venn-like diagrams. Can we extend this representation to cover our dag-based hP systems? In this section we will show that any hP system structurally based on a canonical dag can still be intensionally represented by hierarchically nested planar regions, delimited by Jordan curves (simple closed curves). Conversely, we also show that any set of hierarchically nested planar regions delimited by Jordan curves can be interpreted as a canonical dag, which can form the structural basis of a number of hP systems. We will first show how to represent a canonical dag as a set of hierarchically nested planar regions. Algorithm 3 (Algorithm for visually representing a canonical dag) Without loss of generality, we consider a canonical dag (V, δ) of order n, where vertices are topologically ordered according to the order implied by the arcs, by considering parents before the children, i.e., V={vi|i∈ {1, . . . , n}}, where (vi, vi+1)∈δ. Figure 13 shows side by side a simple height 1 canonical dag and its corresponding hypergraph representation. Note the intensional representation (as opposed to the extensional one): v2is not totally included in v1, although all vertices included in v2, i.e., v4and v5, are also included in v1. A possible topological order is v1, v2, v3, v4, v5. For each vertex vi, we associate a distance ψi=1 2(n−i+1) , for i∈ {1, . . . , n}. For Figure 13, ψi=1 32 ,1 16 ,1 8,1 4,1 2, for i∈ {1, . . . , n}. We process the vertices in reverse topological order vn, . . . , v1, at each step i representing the current vertex viby a planar region Ri. First, we set parallel horizontal axis Xoand Xp, vertically separated by distance 3(n−1). Secondly, we set points o1, . . . , onon Xo, such that oiand oi+1 are separated by distance 3, for 1 ≤i≤n−1. We define oias the origin point of vi, 104 R. Nicolescu, M.J. Dinneen, Y.-B. Kim v3v4v5 v1 v2 1 2 345 Fig. 13. A simple canonical dag and its corresponding hypergraph representation. and write oi=origin(vi). Finally, we set points p1, . . . , pnon Xp, such that piand pi+1 are separated by distance 3, for 1 ≤i≤n−1. We define pias the corridor point of vi. Figure 14 shows the construction of Xo, Xp, oiand pi, for the dag of Figure 13, where n= 5. o3o4o5 o1o2 p1p2p3p4p5 Xo Xp Fig. 14. Construction of Xo, Xc, oiand pi, for the dag of Figure 13, where n= 5. If the current vertex viis a sink, then Riis a circle with with radius 1 2centered at oi. If the current vertex viis a non-sink, then Riis constructed as follows: Assume that the children of viare w1, . . . , wni, and their (already created) regions are S1, . . . , Sni. Consider line segments l0, l1, . . . , lni, where l0is bounded by oiand pi, and ljis bounded by piand origin(wj), for j∈ {1, . . . , ni}. Let L0,L1, . . . , Lni, T1, . . . , Tnibe the regions enclosed by Jordan curves around l0,l1, . . . , lni, S1, . . . , Sni, at a distance ψi, and let R0 i=L0∪Sj=1,...,niLj∪Sj=1,...,niTj. We define Rias the external contour of R0 i. This definition will discard all internal holes, if any, without introducing any additional containment relations between our regions. The details of our construction guarantee that no internal hole will ever contain an origin point. ut Figure 15 shows an intermediate step (left) and the final step (right) of applying Algorithm 3 on Figure 13. The representation of Figure 15 is topologically Structured Modeling with Hyperdag P Systems 105 equivalent to the hypergraph representation of Figure 13 (right). Figure 16 shows the side by side, another dag and its corresponding planar region representation; internal holes are represented by dotted lines. Our objective here was not to create “nice” visualizations, but to prove that it is possible to represent an arbitrary canonical dag, i.e., an arbitrary hP system structurally based on a canonical dag, by hierarchically nested planar regions. Xo Xp R1R2R3R4R5 Xo Xp R2R3R4R5 Fig. 15. An intermediate step and the final step of applying Algorithm 3 on Figure 13. 1 2 3 4 R1R2R3R4 Xo Xp Fig. 16. A height two dag and its corresponding representation, built by Algorithm 3. We will next show that, for any finite set of hierarchically nested planar regions, we can build a corresponding canonical dag (i.e., the underlying structure of a hP system). 106 R. Nicolescu, M.J. Dinneen, Y.-B. Kim Algorithm 4 (Algorithm for building a canonical dag from finite set of hierarchically nested planar regions) Assume that we have nhierarchically nested planar regions, 1. Label each planar region by Ri,i∈ {1, . . . , n}, 2. If Ridirectly nests Rjthen draw an arc from a vertex vito a vertex vj, i, j ∈ {1, . . . , n},i6=j. ut We now show that a canonical graph produced from Algorithm 4 does not contain any cycles. Our proof is by contradiction. Let us assume a directed graph Gproduced from Algorithm 4 contains a cycle vi, . . . , vk, . . . , vi. Then every vertex in a cycle has an incoming arc. If vertex vkis a maximal element in a cycle, with respect to direct nesting, then its corresponding planar region Rkhave the largest region area among planar regions in a cycle. Since no other planar region in a cycle can contain Rk, there are no arc incident to vertex vk. Hence, there is no cycle in G. Remark 4. We leave open the problem of representing dags (that is hP systems) that contain transitive arcs. 7 Summary We have proposed a new model, as an extension of P systems, that provides a better communication structure and we believe is often more convenient for modeling realworld applications based on tree structures augmented with secondary or shared communication channels. We have shown that hP systems functionally extends the basic functionality of transition P systems and neural P systems, even though the underlying structure of hP systems is different. In the dag/hypergraph model of hP systems we can have a natural separation of computing cells (sinks) from communication cells (hyperedges). This model also allows us to easily represent multiple inheritance and/or to distribute computational results (as specified by a dag) amongst several different parts of a membrane structure. We note that the operational behavior of hP systems is separate from the topological structure of a membrane system. In this paper, we illustrated hP systems using the computational rules of nP systems, where multisets of objects are repeatedly changed within cells, by using a fixed set of multiset rewriting rules, or transferred between cells, using several possible transfer modes. Finally, we provided a intuitive visualization of hP systems, by showing that any set of hierarchically nested planar regions (which represents any set of cells ordered by containment) is equivalent to, or modeled by, a dag without transitive arcs. We provided simple algorithms to translate between these two interpretations. Part B of this paper will explore this model further and support it with a more extensive set of examples. Structured Modeling with Hyperdag P Systems 107 References 1. C. Berge: Hypergraphs. Combinatorics of Finite Sets. Elsevier Science Publishers, 1989. 2. C. Carath´eodory: Theory of Functions of a Complex Variable. Vol. 1, Chelsea Publishing Company, 1954. 3. G. Ciobanu: Distributed algorithms over communicating membrane systems. BioSystems, 70, 2 (2003), 123–133. 4. W.F. Doolittle: Uprooting the tree of life. Scientific American, 282 (2000), 90–95. 5. T.-O. Ishdorj, M. Ionescu: Replicative-distribution rules in P systems with active membranes. Proceedings of the First International Colloquium on Theoretical Aspects of Computing (ICTAC 2004), 68–83. 6. C. Li: Master Thesis. The University of Auckland, Supervisor: R. Nicolescu, 2008. 7. C. Mart´ın-Vide, G. P˘aun, J. Pazos, A. Rodr´ıguez-Pat´on: Tissue P systems. Theoretical Computer Science, 296, 2 (2003), 295–326. 8. G. P˘aun: Computing with membranes. Journal of Computer and System Sciences, 61, 1 (2000), 108–143 (and Turku Center for Computer Science, TUCS Report 208, November 1998). 9. G. P˘aun: Membrane Computing–An Introduction. Springer-Verlag, 2002. 10. G. P˘aun: Introduction to membrane computing. Proceeding of the First Brainstorming Workshop on Uncertainty in Membrane Computing, 2004, 17–65. 11. G. P˘aun, R.A. P˘aun: Membrane computing as a framework for modeling economic processes. Seventh International Symposium on Symbolic and Numeric Algorithms for Scientific Computing (SYNASC 2005), 11–18. 12. M. Slikker, A. Van Den Nouweland: Social and Economic Networks in Cooperative Game Theory. Kluwer Academic Publishers, 2001. 13. V.I. Voloshin: Coloring Mixed Hypergraphs: Theory, Algorithms and Applications. American Mathematical Society, 2002. 14. C. Zandron, C. Ferretti, G. Mauri: Solving NP-complete problems using P systems with active membranes. Proceedings of the Second International Conference on Unconventional Models of Computation, 2000, 289–301. 114 A. ObtuÃlowicz 26. Newman, M. E. J., The structure and function of complex networks, SIAM Review 45 (2003), pp. 167–256. 27. ObtuÃlowicz, A., On mathematical modeling of anatomical assembly, spatial features, and functional organization of cortex by application of hereditarily finite sets, in: Proceedings of 9th Workshop on Membrane Computing, Edinburgh, July 28 – July 31, 2008, ed. P. Frisco et al., Technical Report HW-MASC-TR-0061, School of Mathematical and Computer Sciences, Heriot–Watt University, Edinburgh, UK, 2008, pp. 371–382. 28. P˘aun, Gh., Perez-Jimenez, M. J., Spiking neural P systems. Recent results, research topics, in: 6th Brainstorming Week on Membrane Computing, Sevilla 2008, web page. 29. Ravasz, E., Barab´asi, A.-L., Hierarchical organization in complex networks, Physical Review E 67 (2003), 026112. 30. Sporns, O., Small–world connectivity, motif composition, and complexity of fractal neuronal connections, BioSystems 85 (2006), pp. 55–64. 31. Sporns, O., Chialvo, D. R., Kaiser, M., Hilgetag, C. C., Organization, development and function of complex brain networks, Trends in Cognitive Sciences 8 (2004), pp. 418–425. 32. Sporns, O., Tononi, G., Edelman, G. M., Theoretical neuroanatomy:relating anatomical and functional connectivity in graphs and cortial connection matrices, Cerebral Cortex 10 (2000), pp. 127–141. 33. Wuchty, S., Ravasz, E., Barab´asi, A.-L., The Architecture of Biological Networks. The Discovery of Initial Fluxes of Metabolic P Systems Roberto Pagliarini, Vincenzo Manca Verona University Computer Science Department Strada Le Grazie 15, 37134 Verona, Italy. {roberto.pagliarini,vincenzo.manca}@univr.it Summary. A central issue in systems biology is the study of efficient methods to infer fluxes of biological reactions starting from experimental data. Among the different techniques proposed in the last years, in the theory of Metabolic P systems Log-Gain principles have been introduced, which prove to be helpful for deducing biological fluxes from temporal series of observed dynamics. However, crucial tasks remain to be performed for a complete suitable application of these principles. In particular the algebraic systems introduced by the Log-Gain principles require the knowledge of the initial fluxes associated with a set of biochemical reactions. In this paper we propose an algorithm for estimating initial fluxes, which is tested in two case studies. 1 Introduction In the last years, the problem of reverse-engineering of biological phenomena from experimental data has spurred increasing interest in the scientific communities. For these reasons, many computational models inspired from biology have been given. Among these models, the Metabolic P systems [9, 10], shortly MP systems, proved to be relevant in the analysis of dynamics of biochemical processes [4, 12, 14, 13]. MP systems intend to model metabolic systems, that is, structures where matter of different types is subject to reactions, or transformations of various types. The importance of these computational models is their potential applicability to the reverse-engineering problem of biological phenomena. In fact, the MP systems introduce a theory, called Log-Gain [8], intrinsically related to the structure of these computational models. This theory provides a method for constructing MP models of real phenomena from time-series of observed dynamics. In fact, given a real system which can be observed in its evolution, then almost all the elements occurring in the definition of MP system can be, in principle, deduced by macroscopic observations of the system [9]. 116 R. Pagliarini, V. Manca The only component which cannot be directly measured is a set of regulation functions which state the reaction fluxes, that is, the amount of reactants transformed by the reactions at any state of the system. These functions depend on internal microscopic processes on which molecules are involved. However, Log- Gain theory provides a way for deducing them from time-series of the states of a given system along a sufficient number of observation instants. A key point for achieving this task consists in the discovery of the fluxes associated to the passage of a metabolic system from two initial observation instants. In this paper an algorithm is proposed for estimating these initial fluxes from few steps of observation. The present paper is organized as follows. Section 2 is devoted to the definition of Metabolic P Systems. Section 3 briefly recalls the Log-Gain theory. In Section 4 we describe the algorithm that solves our problem. Section 5 reports some experimental results obtained by the new framework. Some further remarks and directions for future researches are discussed in the last section. 2 Metabolic P Systems MP systems are a special class of dynamical systems (the reader can find some details concerning dynamical aspects of MP systems in [11]), based on P systems [3, 16, 17, 18], which are related to metabolic processes. MP systems are essentially constituted by multiset grammars where rules are regulated by specific functions depending on the state of the system. From a Membrane Computing point of view, MP systems can be seen as deterministic mono-membrane P systems where the transitions between states are calculated by a suitable recurrent operation. In an MP system the overall variation, in a macroscopic time interval, of the whole system under investigation is considered. In this manner, the evolution law of the system consists in the knowledge of the contribution of each reaction in the passage between any two instants separated by such an interval. Therefore, dynamics is given at discrete steps, and at each step, it is ruled by a partition of matter among the reactions transforming it. The principle underlying the partitioning is called mass partition principle. This principle replaces the mass action law1of ODE systems. The mass partition principle defines the transformation rate of object populations rather than single objects, according to a suitable generalization of chemical laws [9]. 1The foundation of this law is the theory of molecular collisions. The first formulation of this law, formulated by Waage and Guldberg [21], is the following: “the rate of any given chemical reaction is proportional to the product of the concentrations of the reactants”. The Discovery of Initial Fluxes of Metabolic P Systems 117 2.1 MP systems: a formal definition The following definition introduces the MP systems in a formal way (N,Z, and R respectively denote the sets of natural, integer, and real numbers). Definition 1 (MP system) An MP system Mis specified by the following construct: M= (X, R, V, H, Φ, ν, µ, τ, δ) where X,Rand Vare finite disjoint sets, and moreover the following conditions hold, with n, m, k ∈N: •X={x1, x2, . . . , xn}is a finite set of substances. This set represents the types of molecules; •R={r1, r2, . . . , rm}is a finite set of reactions. A reaction ris a pair of type αr→βr, where αridentifies the multiset of the reactants (substrates) of r and βridentifies the multiset of the products of r(λrepresents the empty multiset). The stoichiometric matrix Aof a set Rof reactions over a set X of substances is A= (Ax,r |x∈X, r ∈R)with Ax,r =|βr|x− |αr|x, where |αr|xand |βr|xrespectively denote the number of occurrences of xin αrand βr. Of course, a reaction rcan be seen as the vector r= (Ax,r |x∈X)of Zn. We also set Rα(x) = {r∈R|x∈αr},Rβ(x) = {r∈R|x∈βr}, and R(x) = Rα(x)∪Rβ(x); •V={v1, v2, . . . , vk}is a finite set of parameters (such as pressure, temperature,. . . ); •H={hv|v∈V}is a set of parameters evolution functions. The function hv:N→Rstates the value of parameter v, and H[i] = (hv(i)|v∈V); •Φ={ϕr|r∈R}is the set of flux regulation maps, where, for each r∈R, ϕr:Rn+k→R. Let q∈Rnbe the vector of substances values and s∈Rkbe the vector of parameters values. Then (q, s)∈Rn+kis the state of the system. We set by U(q, s) = (ϕr(q, s)|r∈R)the flux vector in the state (q, s); •νis a natural number which specifies the number of molecules of a (conventional) mole of M; •µis a function which assigns, to each x∈X, the mass µ(x)of a mole of x (with respect with to some measure units); •τis the temporal interval between two consecutive observation steps; •X[i] = (x1[i], x2[i], . . . , xn[i]), for each i∈N, is the vector of substances values at the step i, and X[0] are the initial values of substances. The dynamics δ:N→Rnof the system is completely identified by the following recurrent equation, called Equational Metabolic Algorithm shortly EMA: X[i+ 1] = A×U(X[i], H[i]) + X[i] (1) where Ais the stoichiometric matrix of reactions having dimension n×m, while ×,+, are the usual matrix product and vector sum. We denote by EMA[i] the system (1) at the step i. By using the formulation introduced above it is simple to note that we can obtain the vector X[i+ 1] from vectors X[i]and U(X[i], X[i]). 118 R. Pagliarini, V. Manca 3 Log-Gain Theory: A Brief Recall The starting point of the Log-Gain theory for MP systems [20] is the Allometry Law [2, 6] which assumes a proportion between the relative variations of the fluxes of a reaction rand the sum of relative variations of tuners of r, that is, magnitude influencing r. The relative variation of a variable xis given, in differential notation and with respect to the time variable t, by d(lg x)/dt. This explains the term “Log-Gain”. Given a dynamics of an MP system, we will use the following simplified notations, for i∈N, and r∈R: ur[i] = ϕr(X[i], H[i]) and U[i] = (ur[i]|r∈R) (2) Assuming to know the vectors X[i] and X[i+1], the equation (1) can be rewritten in the following form, which we called ADA[i] (Avogadro and Dalton Aggregation [10]): X[i+ 1] = A×U[i] + X[i] (3) The formula (3) expresses a system of nequations and mvariables (nis the number of substances and mthe number of reactions) which is assumed to have maximal rank. This supposition is not restrictive. In fact, if it does not hold the rows which are linearly depend on other rows are removed. Formally ADA[i] is the same to system EMA[i] introduced in Section 2. However, these two systems have dual interpretations. In fact, in EMA[i], the vectors U[i] and X[i] are known, and the vector X[i+ 1] is computed by means of them, while in ADA[i], the vector X[i+ 1] −X[i] is known and U[i] is computed by solving the system, as we will see by formula (6). Usually, in a biochemical phenomenon, the number of reactions is greater than the number of substances, and this means that the system (3) has more than one solution. Therefore, fluxes cannot be univocally deduced by means of ADA. The following principle [8] allows us to add more equations to the above system in order to get a univocally solvable system which could provide the flux vector. Definition 2 (Discrete Log-Gain) Let (z[i]|i∈N)a real valued sequence. Then, the discrete log-gain of zis given by the following equation: Lg(z[i]) = z[i+ 1] −z[i] z[i](4) Principle 1 (Log-Gain regulation) Let U[i], for i≥0, be the vector of fluxes at step i. Then the Log-Gain regulation can be expressed in terms of matrix and vector operations: (U[i+1]−U[i])/U[i] = B×L[i] + C⊗P[i+ 1] (5) where: The Discovery of Initial Fluxes of Metabolic P Systems 119 •B= (pr,z |r∈R, z ∈X∪V)where pr,z ∈ {0,1}with pr,z = 1 if zis a tuner of rand pr,z = 0 otherwise; •L[i] = (Lg(z[i]) |z∈X∪V)is the column vector of substances and parameters log-gains ; •P[i+ 1] is a column vector of values associated with the reactions and called (Log-Gain) offsets at step i+ 1; •C= (cr|r∈R), where cr= 1 if r∈R0, while cr= 0 otherwise, and R0is a subset of reactions having the Covering Offset Log Gain Property, that is, it is a set of nlinear independent vectors of Zn; • × denotes the usual matrix product; •+,−,/,⊗denote the component-wise sum, subtraction, division and product2. If we assume to know the flux unit vector at step iand put together the equations (5) and (3) at steps iand i+ 1 respectively, we get the following linear system called Offset Log-Gain Adjustment module at step i, shortly OLGA[i], in which the number of variables (here reported in bold font) is equal to the number of equations: A×U[i+1] = X[i+ 2] −X[i+ 1] (6) (U[i+1]−U[i])/U[i] = B×L[i] + C⊗P[i+1] Now, if the vectors X[i] and V[i], for 0 ≤i≤l, where l∈N, are obtained by experimental measures, then it is possible to solve OLGA[i] for i= 0, . . . , l −1, obtaining the vector U[i] for i∈[1, l −1]. 4 An Algorithm for the Estimation of Initial Metabolic Fluxes The method described in the previous section assumes the knowledge of the initial values of fluxes. Problem 1 (Initial Fluxes Problem) Given X[0] and H[0], find a flux vector U[0] such that it satisfies the initial dynamics, that is: X[1] u A ×U[0] + X[0] where umeans approximate equality. The algorithm given below circumvents the Initial Fluxes Problem by using the knowledge about the dynamics in the first evolution steps in order to approximate the amount of substances which is not transformed, we call inertia of the system (at a given step). 2Given two n×mmatrices Aand B, the operation A⊗Binvolves the action of multiplying component-wise each element of Aby the corresponding element of B. 120 R. Pagliarini, V. Manca 4.1 The proposed algorithm Our approach is based on the assumption that if the inertia of each substance is known, then only a part of substances has to be partitioned among the reactions which require to consume them. The main steps of the algorithm are described in the following of this section. Step 1. The goal of the first step is to evaluate grossly the initial fluxes at the step 0 by assuming that they are proportional to the reactants, that is, for all r∈R: ˆur[i] = kryr[i] (7) where kr∈R, and yr[i] is the product of all substance quantities, at the step i, which are reactants for r. We suppose that if αr=λthen yr= 1, and we set ˆ U[i] = (ˆur[i]|r∈R) (8) For example, in a metabolic system with three kinds of substances, a,b,c, and with a set of reactions given in the first column of the Table 1, the relationships among the fluxes of these reactions and their reactants are reported in the second column of the Table 1. Let us consider the following system, called Local-Stoichiometric Module at the Reactions Maps r1:a→bc kr1a r2:b→a kr2b r3:c→ab kr3c r4:c→cc kr4c Table 1. Reactions and their flux regulation maps of the Local-Stoichiometric Module. step i: x[i+ 1] −x[i] = X r∈R(x) Ax,r ˆur[i]∀x∈X(9) If we assume that the constants kr, with r∈R, do not sensibly change in few steps, then by applying the system (9) for a sufficient number of steps we can obtain a square linear system of dimension mhaving maximum rank. In the example reported in Table 1, we have a Local-Stoichiometric Module of 3 equations which initially has 4 unknowns. It has rank 3. At the second iteration of this module we get other 3 equations and the rank of Local-Stoichiometric Module is maximum. Thus, we can obtain a system of equation having unique solution. In general, if we start with the Local-Stoichiometric Module at the step 0 then we can compute the vector ˆ U[0] = (ˆur[0] |r∈R) by applying the Local-Stoichiometric module a suitable number of steps. The Discovery of Initial Fluxes of Metabolic P Systems 121 Step 2. The aim of this step is to approximate the inertia of the system. We split this step in two sub-steps. In the first one we take nlinear independent reactions, obtaining a set R0,according to the Covering Offset Log Gain Property. Then, we use the set R0to obtain an OLGA[1] module, with U[0] = ˆ U[0], where ˆ U[0] is the vector of fluxes computed in the previous step. We will indicate with U∗[1] = (u∗ r[1] |r∈R) the solution of this system. However, if some elements of this vector is a negative real value, then we choose others nlinear independent reactions and reapply the procedure above describe (it easy to prove that a positive vector must exist). In the second sub-step we compute, for each x∈X, the inertia, indicated by ¯x, by applying the following equation: ¯x[1] = x[1] −X r∈Rα(x) u∗ r[1],∀x∈X(10) Step 3 In the last step we obtain the vector of fluxes at the evolution step 1 by solving an optimization problem. In fact, the vector U◦= (u◦ r|r∈R) we search has to be a strictly positive vector of Rm(positive in each component) which satisfies the following nequations: x[1] −¯x[1] = X r∈Rα(x) u◦ r[1],∀x∈X(11) and it is bounded, for each component, by the following constraint: u◦ r[1] ≤   min xj[1]−¯xj[1] |αr|xj |xj∈αrif αr6=λ krif αr=λ (12) and such that U◦= min ξ∈RmkA×ξ−(X[2] −X[1])k(13) 5 Experiments In this section, in order to evaluate the performance of our algorithm, we apply it to two case studies: i) a synthetic oscillatory metabolic system, and ii) the Belousov-Zhabotinsky reaction [1, 7, 19, 22]. 122 R. Pagliarini, V. Manca Reactions Flux regulation maps r1:a→aa ϕ1=k1a/(k1+k2c+k4b+k6) r2:a→b ϕ2=k2ac/(k1+k2c+k4b+k6) r3:b→λ ϕ3=k3b/(k3+k6) r4:a→c ϕ4=k4ab/(k1+k2c+k4b+k6) r5:c→λ ϕ5=k5c/(k5+k6) X[0] = (100 100 1) k1=k3=k5= 4, k2=k4= 0.02, k6= 100 Table 2. Sirius’ reactions and maps. 5.1 A synthetic metabolic system Let us consider the synthetic non-cooperative metabolic system called Sirius and given in Table 2 [9]. Firstly, we generate the dynamics of this model for 1000 steps by using the flux regulation maps of Sirius. Then, we use our algorithm to approximate the vector of fluxes U◦[1] at the evolution step 1. Starting from U◦[1], by applying OLGA[i] for i= 1,2, . . . , 900, we deduce the vectors U[i], for i= 2,3, . . . , 899, according to the Log-Gain theory. Figure 1 shows the fluxes relative to the dynamics of Sirius initially generated, while the Figure 2 shows the inferred fluxes. These results show an almost complete accordance. Fig. 1. The values of Sirius’ fluxes calculated by using the flux regulation maps given in Table 2. The Discovery of Initial Fluxes of Metabolic P Systems 123 Fig. 2. The values of Sirius’ fluxes calculated by applying the Log-gain theory and the initial vector of fluxes inferred by the proposed algorithm. 5.2 A biochemical case study In this subsection the application of the algorithm to approximate the initial fluxes of the Belousov-Zhabotinsky reaction, also known as BZ reaction, is discussed. This reaction represents a famous example of a biochemical oscillatory phenomenon. Its importance is that it is the first evidence of a biochemical clock. Although the stoichiometry of the BZ reaction is quite complicated, several simplified mathematical models of this phenomenon have been proposed. In particular, Prigogine and Nicolis [15] proposed a simplified formulation of the dynamics of the BZ reaction, called Brusselator, whose oscillating behavior is represented by only two substances, xand yrespectively, and it is governed by the following system of differential equations: dx dt =k1−k2x+k3x2y−k4x(14) dy dt =k2x−k3x2y where k1= 100, k2= 3, k3= 10−4and k4= 1 represent constant rates. The numerical solution of the system system (14) from initial conditions x= 1 and 130 L. Pan, Gh. P˘aun The rules of type (1) are firing (we also say spiking)rules, and they are applied as follows. If the neuron σicontains kspikes, and ak∈L(E), k ≥c, then the rule E/ac→a;dcan be applied. The application of this rule means removing cspikes (thus only k−cremain in σi), the neuron is fired, and it produces a spike after dtime units (a global clock is assumed, marking the time for the whole system, hence the functioning of the system is synchronized). If d= 0, then the spike is emitted immediately, if d= 1, then the spike is emitted in the next step, etc. If the rule is used in step tand d≥1, then in steps t, t +1, t +2, . . . , t +d−1 the neuron is closed (this corresponds to the refractory period from neurobiology), so that it cannot receive new spikes (if a neuron has a synapse to a closed neuron and tries to send a spike along it, then that particular spike is lost). In the step t+d, the neuron spikes and becomes again open, so that it can receive spikes (which can be used starting with the step t+d+ 1). The rules of type (2) are forgetting rules and they are applied as follows: if the neuron σicontains exactly sspikes, then the rule as→λfrom Rican be used, meaning that all sspikes are removed from σi. If a rule E/ac→a;dof type (1) has E=ac, then we will write it in the following simplified form: ac→a;d. In each time unit, if a neuron σican use one of its rules, then a rule from Rimust be used. Since two firing rules, E1/ac1→a;d1and E2/ac2→a;d2, can have L(E1)∩L(E2)6=∅, it is possible that two or more rules can be applied in a neuron, and in that case, only one of them is chosen non-deterministically. Note however that, by definition, if a firing rule is applicable, then no forgetting rule is applicable, and vice versa. Thus, the rules are used in the sequential manner in each neuron, but neurons function in parallel with each other. The initial configuration of the system is described by the numbers n1, n2, . . . , nm, of spikes present in each neuron. During a computation, the “state” of the system is described by both by the number of spikes present in each neuron, and by the open/closed condition of each neuron: if a neuron is closed, then we have to specify when it will become open again. Using the rules as described above, one can define transitions among configurations. Any sequence of transitions starting in the initial configuration is called a computation. A computation halts if it reaches a configuration where all neurons are open and no rule can be used. With any computation (halting or not) we associate a spike train, the sequence of zeros and ones describing the behavior of the output neuron: if the output neuron spikes, then we write 1, otherwise we write 0. An SN P system can be used in various ways. In the generative mode, we start from the initial configuration and we define the result of a computation (i) either as the number of steps between the first two spikes sent out by the output neuron, or (ii) as the number of spikes present in neuron σout when the computation halts (note that in the first case we do not request that the computation halts after sending out two spikes). We denote by N2(Π) the set of numbers computed in the first way by an SN P system Πand by Ngen(Π) the set of numbers generated New Normal Forms for Spiking Neural P Systems 131 by Πin the second case. We can also use Πin the accepting mode: a number nis introduced in the system in the form of a number f(n) of spikes placed in neuron σin, for a well-specified mapping f, and the number nis accepted if and only if the computation halts. (Alternatively, we can introduce the number to be recognized as the distance in time between two spikes entering neuron σin from the environment.) We denote by Nacc(Π) the set of numbers accepted by Π. In the generative case, the neuron (with label) in is ignored, in the accepting mode the neuron out is ignored (in most cases below, we identify the neuron σi with its label i, so we say “neuron i” understanding that we speak about “neuron σi”). We can also use an SN P system in the computing mode, introducing a number in neuron in and obtaining a result in neuron out, but we do not consider this case here. A neuron σi(in the initial configuration of an SN P system) is characterized by ni, the number of spikes present in it, and by Ri, its associated set of rules. An SN P system is said to be in the kR-normal form, for some k≥1, if there are at most kdifferent sets R1, . . . , Rkof rules used in the mneurons of the system. An SN P system is said to be in the knR-normal form, for some k≥1, if there are at most kdifferent pairs (n1, R1),...,(nk, Rk) describing the mneurons of the system. We denote by NαSNP(kβ, forg, dley) the families of all sets Nα(Π) computed by SN P systems in the kβ-normal form, for α∈ {2, gen, acc},β∈ {R, nR}, and k≥1; if no forgetting rules are used, then we remove the indication forg from the notation; if all rules have delay d= 0, then we remove the indication dley from the notation. 4 A 3R-Normal Form Result We are going now to prove the main result mentioned in the Introduction: SN P systems with only three different sets of rules are universal when generating numbers encoded in the first two spikes of the spike train. Theorem 1. NRE =N2SNP (3R, dley). Proof. We show that NRE ⊆N2SNP (3R, dley); the converse inclusion is straightforward (or we can invoke for it the Turing-Church thesis). Let us consider a register machine M= (m, H, l0, lh, I) with the properties specified in Section 2. We construct an SN P system Πwhich simulates Min the way somewhat standard already when proving that a class of SN P systems is universal. Specifically, we construct modules ADD and SUB to simulate the instructions of M, as well as an output module FIN which provides the result (in the form of a suitable spike train). Each register rof Mwill have a neuron rin Π, and if the register contains the number n, then the associated neuron will contain 2nspikes. The modules will be given in a graphical form, indicating their initial configuration, the synapses, and, for each neuron, the associated set of rules; all neurons 132 L. Pan, Gh. P˘aun are initially empty, with the exception of the neuron associated with the initial label, l0, of M, which contains one spike, and with exception of a few other neurons, as shown in the following figures. ' & $ % ' & $ % ' & $ % ' & $ % ' & $ % ' & $ %' & $ % ' & $ % # " !# " ! # " !# " ! # " !         ' & $ %   - ZZZZZ Z~ HHHH Hj  * ?         QQQQQQQQQ Qs ? @@ @R      @@@@ @R      ? ? ? ? -  QQQ Qs Q Q Qk  -    = > li a(aa)∗/a →a; 0 a→a; 0 li1 a(aa)∗/a →a; 0 a→a; 0 r li2 a(aa)∗/a →a; 0 a→a; 0 li3a(aa)∗/a →a; 0 a→a; 0 li6a(aa)∗/a →a; 0 a→a; 0 li4a(aaa)∗/a →a; 0 a→a; 1 li5 a(aa)∗/a →a; 0 a→a; 0 li7 a(aa)∗/a →a; 0 a→a; 0 li8 a(aa)∗/a →a; 0 a→a; 0 li9 a2 a(aaa)∗/a →a; 0 a→a; 1 li10 a2 a(aaa)∗/a →a; 0 a→a; 1 li11 a(aa)∗/a →a; 0 a→a; 0 li12 a(aa)∗/a →a; 0 a→a; 0 ljlk li13 a(aa)∗/a →a; 0 a→a; 0 li14 a(aa)∗/a →a; 0 a→a; 0 Fig. 1. Module ADD, simulating li: (ADD(r), lj, lk) We consider the following three sets of rules: R1={a(aa)∗/a →a; 0, a →a; 0}, R2={a(aa)∗/a3→a; 0, a →a; 1}, R3={a(aaa)∗/a →a; 0, a →a; 1}. New Normal Forms for Spiking Neural P Systems 133 The ADD module used to simulate an addition instruction li: (ADD(r), lj, lk) is indicated in Figure 1. No rule in R1can be applied in the presence of an even number of spikes. If a spike enters the neuron (with the label) li, then this neuron starts using its rules; initially, this is the case with neuron l0. Neuron lispikes and one spike is sent to both neuron li1and neuron li2, which also spike in the next step. In this way, two spikes are sent to neuron r, and this represents the increment of register rby one. Neuron li2also sends a spike to neurons li3,li4, and li5. Neurons li3and li5spike immediately, while neuron li4can non-deterministically choose either rule to use as both of them are enabled by the existence of a single spike – this ensures the non-deterministic passage to one of the instructions ljor lk.. Assume that σli4uses the rule a(aa)∗/a →a; 0. This means that in the next step σli8receives two spikes, hence no rule here can be used. Simultaneously, neurons li6and li7receive one spike each, and both of them spike. In this way, σli9receives one spike and σli7continues having one spike. Neuron li9contains now a number of spikes of the form 3n+ 3, for some n≥0 (initially we had two spikes here, hence n= 0) and no rule is enabled. In the next step, this neuron receives one further spike, and the first rule is fired (the number of spikes is now 3(n+ 1) + 1). All neurons ljand li11, li12 receive one spike. The last two neurons send back to σli9one spike each, hence the number of spikes in this neuron will be again congruent with 2 modulo 3, as at the beginning. Thus, the neuron associated with the label ljhas been activated. If neuron li4uses the rule a→a; 1, then σli7receives two spikes at the same (after one time unit) time and this branch remains idle, while neurons li8, li10, li13, li14 behave like neurons li7, li9, li11, li12, and eventually σlkis activated and the number of spikes from σli10 returns to the form 3s+ 2, for some s≥0. The simulation of the ADD instruction is correctly completed. The SUB module used to simulate a subtraction instruction li: (SUB(r), lj, lk) is shown in Figure 2. Because the reader has the experience of examining the work of the ADD module, this time we do not write explicitly the rules, but the sets R1, R2, R3as defined above. Like in the case of the ADD module, the SUB module starts to work when a spike enters the neuron with the label li. The functioning of each neuron is similar to the previous case (the rules to be used are chosen in the same way and eventually the neurons remain with a number of spikes like that in the starting configuration). The neuron lisends a spike to the neurons li1, and r. If register ris not empty, then the rule a(aa)∗/a3→a; 0 of R2will be applied. Assume that this is the case. This means that σrspikes immediately, hence σli2 receives two spikes (one from σli1) and is doing nothing, while neuron li3receives one spike and it fires. A spike is sent to each of the three neurons li7, li8, and li9. This last neuron will send a spike to σli11 , which will spike, thus activating the neuron associated with label lj. The two spikes sent by σli7, σli8to σli10 wait here, as no rule is enabled for a number of spikes of the form 3n+ 2. In the next step, a spike comes from neuron li11, hence σli10 ends with a number of spikes which is a multiple of 3, hence no rule is activated. 134 L. Pan, Gh. P˘aun          # " !                                 @@@R ?        + SSSS Sw Z Z~ BB BN     / JJJ J^ @@@@@@@@ @R HHHHHHHHHHH Hj   +    SSSS Sw    /                    ? ? liR1 li1 R1r R2 li2R1 li4 R1li5 R1 li6 R1 li3 R1 li7 R1 li8 R1 li9 R1 li10 a3 R3 li11 a3 R3 lklj -  Fig. 2. Module SUB, simulating li: (SUB(r), lj, lk) If the register ris empty, then in σrwe have to use the rule a→a; 1. The neuron li2receives a spike from σli1and in the next step it fires, at the same time with the move of spikes produced at the previous step in σrand kept there because of the delay. In this moment, all neurons li2, li3, li4, li5, and li6contains one spike. Neurons li4, li5, li6send their spikes to σli10 and σli11 , but they immediately receive one spike from σli2. This also happens with neurons li7, li8, li9, which receive one spike each from σli3. Neuron li10 spikes and activated lk, sending at the same time one spike to σli11 , thus completing here the number of spikes to a multiples of three. Similarly, in the next step, σli10 (resp., σli11 ) receives three spikes each, from neurons li4, li7, li8(respectively, li5, li6, li9). The simulation of the SUB instruction New Normal Forms for Spiking Neural P Systems 135 is correctly completed, with the neurons containing numbers of spikes of the same parity as in the beginning. The modules have different neurons, precisely identified by the label of the respective instruction of M. Modules ADD do not interfere. However, a problem appears with modules SUB: when simulating an instruction li: (SUB(r), lj, lk), neuron σrsend one spike to all neurons ls2, ls3from modules associated with instructions ls: (SUB(r), lu, lv) (that is, subtracting from the same register r). However, no undesired effect appears: the spikes arrive simultaneously in neurons ls2, ls3, hence they send one spike to each of the three neurons “below” them, which, in turn, send their spikes to neurons ls10, ls11; each of these neurons gets three spikes, hence no rule can be used here, the spikes are just accumulated (in a number which continues to be a multiple of 3).                         / JJJJ^ - ?      - !! !HH Hj    H H HY lh R1 1 R2 lh1 R1 lh2 R1 R1 out Fig. 3. The FIN module The addition and subtraction modules simulate the computation of M. In order to produce the number generated by Mas the distance between the first two spikes sent out by the system Πwe use the module FIN from Figure 3. It is triggered when Mreaches the lh:HALT instruction. At this point a single spike is sent to neuron 1, and at the same time to σlh1. Neuron σ1sends a spike to each neuron lh1, lh2, and out. The output neuron spikes (for the first time). Neurons σ1and σlh1 continuously exchange spikes, hence at each step from now on neuron σ1contains an odd number of spikes and fires. Neuron out gets two spikes in each step, one from σ1and one from σlh2, hence nothing happens. When the content of σ1is exhausted, the rule a→a; 1 must be used here. The neuron is closed, the spike of σlh1is lost, σout receives only one spike, from σlh2, and spikes for the second time. The work of the system continues forever, because of the interchange of spikes between σlh1and σ1, but we are interested only in the distance between the first 136 L. Pan, Gh. P˘aun two spikes emitted by σout, and this distance is equal to the number stored in register 1 in the end of the computation of M. Consequently, N(M) = N2(Π) and this concludes the proof. 2 In the previous figures one can see that the set R1appears in neurons having zero or one spike (the case of σl0) in the initial configuration, R2only with zero spikes, and the set R3appears in neurons with zero, two, or three spikes. This means that, if we also consider the number of spikes present in a neuron in the initial configuration when defining the type of a neuron, then the previous 3R- normal form becomes a 6nR-normal form. Corollary 1. NRE =N2SNP (6nR, dley). When considering the generated number encoded in the number of spikes present in the output neuron, then several simplifications of the previous constructions are possible. First, the module FIN is no longer necessary; moreover, when a spike is sent to neuron lh, the computation will halt. Because no instruction is performed in the register machine after reaching the instruction lh:HALT, we provide no outgoing synapse for neuron lh, so it does matter which rules are present in this neuron, no change is implied on the result of the computation. Then, we only write to register 1, hence to neuron 1 we do not have to apply SUB operations; this means that we can only add spikes to this neuron, namely, one at a time, while any rule can be used inside because no outgoing synapse is present for this neuron. However, always we have the same number of types of neurons, because three types are necessary in modules ADD and SUB. (Similarly, we can use for neuron 1 a FIN module which halves the number of spikes and sends them to another neuron, which leads back to a construction as above.) The results are again as above (the details are left to the reader): Corollary 2. NRE =NgenSNP (3R, dley) = NgenSNP(6nR, dley). 5 Final Remarks The accepting case brings further simplifications: the ADD instructions are deterministic (hence only one type of neurons is necessary – see Figure 4, where R1 is as above), and the FIN module is no longer necessary (we consider the input given as the number of spikes initially present in neuron σin, without taking into account this number when defining the types of neurons). However, if we also take into consideration the number of spikes present in the neurons, then we get the following types: (0, R1),(1, R1),(0, R2),(3, R3), that is, we have the next result: Corollary 3. NRE =NaccSNP(3R, dley) = NaccSNP(4nR, dley). New Normal Forms for Spiking Neural P Systems 137                     ? @@ @R ? @@@ @R     li R1 li1 R1 li2 R1 r lj Fig. 4. Module ADD, simulating li: (ADD(r), lj) We do not have a construction for module SUB using only two types of neurons. As mentioned in the Introduction, there are several open problems and research topics suggested by the previous results. We conclude by mentioning a few basic ones. Is the result in Theorem 1 optimal, or a 2R-normal form or even a 1R-normal form result is valid? Extend this study to other classes of SN P systems, with other types of rules or with other modes of using the rules. Important remark: The proof of Theorem 1 implicitly shows that the universality of SN P systems can be obtained without using forgetting rules. The result was first stated in [2], but the construction of the SUB module as given in that paper has a bug: the interaction between neurons from different SUB modules acting on the same register is not examined and the undesired interactions avoided. This is done in the construction from Figure 2, as explicitly mentioned above. Note that our construction uses the delay feature; in [1] it is proved that both the forgetting rules and the delay feature can be avoided without losing the universality, but also there it seems that the interaction of neurons in different SUB modules is not carefully checked. Whether or not the idea in our module SUB can be used to obtain such a stronger normal form remains to be seen. Acknowledgements The work of L. Pan was supported by National Natural Science Foundation of China (Grant Nos. 60674106, 30870826, 60703047, and 60803113), Program for New Century Excellent Talents in University (NCET-05-0612), Ph.D. Programs Foundation of Ministry of Education of China (20060487014), Chenguang Program of Wuhan (200750731262), HUST-SRF (2007Z015A), and Natural Science 138 L. Pan, Gh. P˘aun Foundation of Hubei Province (2008CDB113 and 2008CDB180). The work of Gh. P˘aun was supported by Proyecto de Excelencia con Investigador de Reconocida Val´ıa, de la Junta de Andaluc´ıa, grant P08 – TIC 04200. References 1. M. Garc´ıa-Arnau, D. P´erez, A. Rodr´ıguez-Pat´on, P. Sos´ık: Spiking neural P systems: Stronger normal forms. Proc. Fifth Brainstorming Week on Membrane Computing (M.A. Guti´errez-Naranjo et al., eds.), Fenix Editora, Sevilla, 2007, 157–178. 2. O.H. Ibarra, A. P˘aun, Gh. P˘aun, A. Rodriguez-Paton, P. Sosik, S. Woodworth: Normal forms for spiking neural P systems. Theoretical Computer Science, 372, 2-3 (2007), 196–217. 3. M. Ionescu, Gh. P˘aun, T. Yokomori: Spiking neural P systems. Fundamenta Informaticae, 71, 2-3 (2006), 279–308 4. Gh. P˘aun, G. Rozenberg, A. Salomaa, eds.: Handbook of Membrane Computing. Oxford University Press, 2010 (in press). 5. The P Systems Website: http://ppage.psystems.eu. Spiking Neural P Systems with Anti-Spikes Linqiang Pan1,2, Gheorghe P˘aun2,3 1Department of Control Science and Engineering Huazhong University of Science and Technology Wuhan 430074, Hubei, China [email protected], [email protected] 2Department of Computer Science and Artificial Intelligence University of Sevilla Avda. Reina Mercedes s/n, 41012 Sevilla, Spain 3Institute of Mathematics of the Romanian Academy PO Box 1-764, 014700 Bucure¸sti, Romania [email protected], [email protected] Summary. Besides usual spikes employed in spiking neural P systems, we consider “anti-spikes”, which participate in spiking and forgetting rules, but also annihilate spikes when meeting in the same neuron. This simple extension of spiking neural P systems is shown to considerably simplify the universality proofs in this area: all rules become of the form bc→b0or bc→λ, where b, b0are spikes or anti-spikes. Therefore, the regular expressions which control the spiking are the simplest possible, identifying only a singleton. A possible variation is not to produce anti-spikes in neurons, but to consider some “inhibitory synapses”, which transform the spikes which pass along them into antispikes. Also in this case, universality is rather easy to obtain, with rules of the above simple forms. 1 Introduction The spiking neural P systems (in short, SN P systems) were introduced in [4], and then investigated in a large number of papers. We refer to the respective chapter of [7] for general information in this area, and to the membrane computing website from [9] for details. In this note, we consider a variation of SN P systems which was suggested several times, i.e., involving inhibitory impulses/spikes or inhibitory synapses and investigated in a few papers under various interpretations/formalizations – see, e.g., [1], [2], [5], [8]. The definition we take here for such spikes – we call them anti-spikes (somewhat thinking to anti-matter) – considers having, besides usual “positive” spikes denoted by a, objects denoted by ¯a, which participate in spiking or forgetting rules as usual spikes, but also in implicit rules of the form a¯a→λ: if an anti-spike meets a spike in a given neuron, then they annihilate each other, 242 J. Wang, T.-O. Ishdorj, L. Pan In the next section, in particular, we aim to show the fact that the assumed precomputed work space used in [6] to solve SAT can be pre-computed practically in advance in polynomial time by SN P systems with budding rules. Then, a solution to SAT problem is given by the systems with already pre-computed work space. 4 Uniform solution to sat by (dendritic) SN P systems Our SN P system with budding rules is composed of two subsequent subsystems: construction of a SN P system structure which meant to solve sat problem uniformly and the SN P systems family, [6], which solves the sat problem efficiently – for the sake of simplicity, we avoid the neuron budding and the spike firing rules are used at the same time in each subsystem. Π= (O, Σ, H, syn, R, soma, out) where: 1. O={a}is the singleton alphabet; 2. His a finite set of labels for neurons, H⊇H0={soma, out, e0, e1, e2, e3, b1, b2, b3, c, s, +,−} is the labels for neurons initially given; 3. Σ={σi|i∈H0}is the set of initial neurons; 4. syn ⊆H×His a finite set of synapses, with (i, i)/∈syn for i∈H), syn ⊇syn0={(e, ei)|0≤i≤3, e ∈ {+,−}} ∪ {(e0, bi)|1≤i≤3} ∪ {(b3, c),(s, +),(+,−),(−,+), λ}is the set of synapses initially in use; 5. Ris a set of neuron budding and extended spiking rules specified as follows. Fig. 4. The initial topological structure (new born dendrite) of the system Π:soma and out neurons, generator. Efficiency of Spiking Neural P Systems 243 Constructing the system structure The system initially contains an input neuron σsoma, an output neuron σout, and a sub-structure so-called generator block Gcomposed of the set of neurons Σ and the set of synapses syn0,|Σ|=|syn0|= 13, the corresponding topological structure is illustrated in Figure 4. The generation mechanism is governed by only neuron budding rules and controlled by the labels of budding neurons and the created synapses. The labels of each neuron in a subsystem in layer nencodes an associated truth assignment. The system construction algorithm consists of two main parts: A. To generate the dendritic-tree sub-structure (the layers 0 – nin Figure 3, exponentially large in n) and the truth assignments for nBoolean variables. The process starts from the initial neuron σsoma (the root node). B. To complete the network structure. The subsystems in nth layer of the system establish connections to the generator block according to the truth assignments represented in those subsystems, and they are expanded by further three layers, finally converged to the output neuron σout. A. The dendritic-tree generation process, controlled by the labels of neurons as well as the synapses, starts from the initial neuron σsoma (cell body). It is noteworthy that since the truth assignments associated with the subsystems in nth layer are encoded in the labels of those neurons compose each subsystem, the truth assignments are being generated while the dendritic-tree has been constructed. The label of a neuron σcis a sequence of the form c= (k, j, x(p) k) = (k, j, xk(1) = p) = (k, j, p, xk2,...,xkk), p∈ {t, f}, where the first pair (k, j) indicates the location of the neuron on the dendritic-tree: kis the layer number, jis the place where the neuron is in its subsystem, the subsequence x(p) krepresents a string of length kformed by Boolean values tand fbeing generated. Whereas pin x(p) kindicates that the first entry of the subsequence is exactly p– which is later importantly used in the budding rules to distinguish the being generated truth assignments from a same neuron. Moreover, we stress again that the labels of neurons of a chain of length kin a layer krepresents a truth-assignment vof length k, precisely, vis a sequence formed by xk(j), 1 ≤j≤k, of the neuron labels c= (k, j, xk) of a chain. However, hence each chain or subsystem of a layer structure is a separate unit and associates with a truth assignment, all the truth assignments represented in a layer are distinguishable from each other. In other words, the truth assignments are encoded in both the labels of neurons of the chains and its layer structure of composing units too. We do not care which assignment is associated with which subsystem of the layer. In this phase of computation three types of budding rules are performed for the role: budding rules of type a0) applied to the neuron σsoma which initiates the generation of the structure; the dendritic-tree structure is constructed from the layer 0 towards the layer n, for each layer two types of rules such as a1)n−1 times and a2) once, are alternated, total n×nsteps needed to complete. 244 J. Wang, T.-O. Ishdorj, L. Pan a0) [ ]csoma →(csoma, c(1,1,t))[ ]c(1,1,t), [ ]csoma →(csoma, c(1,1,f))[ ]c(1,1,f ), where (csoma, c(1,1,t)),(csoma, c(1,1,f))∈syn. The initial neuron σcsoma buds two new neurons as the rules a0) apply to it simultaneously. The newly produced neurons are: σc(1,1,t)with a synapse (csoma, c(1,1,t)) coming from the father neuron and σc(1,1,f )connected with the father neuron by a synapse (csoma, c(1,1,f)), respectively. Where the symbols tand fin the neuron labels indicate truth values t(rue) and f(alse), respectively, hence the two truth assignments (t) and (f) of length 1 for a single Boolean variable y1are formed. Note that the left hand side of each rule a0) (where λ∈syn0is omitted) requires its interaction environment is empty i.e no synapse exists connected to the neuron σcsoma . Once the rules have applied, the interaction environment of the neuron σcsoma has been evolved having two new synapses going out are created, which makes those rules are not applicable to this neuron anymore. Thus, the base of the first layer of the dendritic-tree has been established, at the first step of the computation. An almost complete system structure for sat(2, m) is depicted in Figure 5, which is worth to follow during the construction. To complete the established layer 1 (in general, i, 1 ≤i≤n), the rules of type a1) generate the 2 (in general 2inumber of) subsystems or the chains of n neurons. a1) (c(k,j−1,x(p) k), c(k,j,x(p) k))[ ]c(k,j,x(p) k) →(c(k,j,x(p) k), c(k,j+1,x(p) k))[ ]c(k,j+1,x(p) k) , p∈ {t, f}, 1 ≤j≤n−1, 1 ≤k≤n,c(k,0,x(p) k)=c(k−1,n,xk−1),x(p) k= (p, xk−1)∈ {t, f}k. The chains composed of nneurons in a layer kare generated by iterative applications of the rules of type a1) in n−1 steps. This rule can be applied in a neuron of type σc(k,j,x(p) k), 1 ≤j≤n−1, when its interaction environment is provided in which exists a single synapse (c(k,j−1,x(p) k), c(k,j,x(p) k)) coming to the neuron. Then each rule buds a single neuron σc(k,j+1,x(p) k)with a synaptic connection (c(k,j,x(p) k), c(k,j+1,x(p) k)), where the second entry (j+ 1) of the neuron label differs from the father neuron as its corresponding label entry as (j), otherwise the rest of the labeling sequence is inherited from the father neuron’s label; xkis a truth assignment of length kover {t, f}. The newly created synapse changes the interaction environment of the father neuron, which prevents another application of the rule. As soon as the last neurons, whose second entry of the label is n, of the layer are produced, the next two types of rules are enabled to apply to those neurons as follows. Efficiency of Spiking Neural P Systems 245 a2) (c(k,n−1,x(p) k), c(k,n,x(p) k))[ ]c(k,n,x(p) k) →(c(k,n,x(p) k), c(k+1,1,t,x(p) k))[ ]c(k+1,1,t,x(p) k) , (c(k,n−1,x(p) k), c(k,n,x(p) k))[ ]c(k,n,x(p) k) →(c(k,n,x(p) k), c(k+1,1,f,x(p) k))[ ]c(k+1,1,f,x(p) k) , p∈t, f, and 1 ≤k≤n−1. Those two rules of type a2) apply simultaneously to each last neuron of type σc(k,n,x(p) k)of each chain in the current layer k, the interaction environments must satisfy the rule condition. As a result, each neuron buds two new neurons with respective synapses. The next layer of the system is thus established. Hence the interaction environment of each father neuron extended by two new synapses, none of these rules is possible to apply again to those neurons. We shall look at the labels of newly produced pairs of type σc(k+1,1,t,x(p) k)and σc(k+1,1,f,x(p) k), the labels are formed as follows: first of all the pair (k+ 1,1) corresponds to the neuron location where k+1 indicates the new layer number while 1 says the neuron is the very first one in its corresponding chain of length nin the new layer; the rest of the labeling sequence as (t, x(p) k) or (f, x(p) k) represents a new truth assignment for Boolean variables x1, x2, . . . , xk+1, where the newly inserted symbol tor fassociates with a truth value t(rue) or f(alse), respectively, while x(p) kis an heritage from the father neuron. Thus, all the possible 2k+1 different truth assignments are generated in layer k+ 1. The truth assignment generation steps for two Boolean variables y1, y2can be observed as described in Figure 5. The rules of type a1) are enabled in turn to complete the newly established layer by continued generation of the chains of length n. By the alternated applications of the rules of types a1) (n−1 times) and a2) (once), in n2steps the layers from 0 to nare, the dendritic-tree, constructed by means an exponential work space and all the truth assignments of length nare generated. Now, we come to the part B of the algorithm. B. The pre-computation to construct the SN P system structure continues until it converges to the output neuron in a further few steps. The main function of this part of the algorithm is to design the substructure which is devoted to the test of the satisfiability of truth assignments against the clauses and to the exploration of the possibility whether any solution to the clauses of the propositional formula exists. The very first task in part Bis to connect the layer nto the generator block appropriately according to the truth assignments formed in this layer. We recall here that, in layer n, there are 2nnumber of subsystems each one is composed of a sequence of nneurons (chains). However, each subsystem injectively corresponds to a different truth assignment of length n. More precisely, taking the label of a neuron σc(n,j,xn)in layer n, where the subsequence xn= (xn1, xn2,...,xnn)∈ {t, f}nrepresents a truth assignment. We associate jth entry of xnwith the jth neuron of considering subsystem, thus, 246 J. Wang, T.-O. Ishdorj, L. Pan Fig. 5. An almost complete structure of Πsystem for sat(2, m) (maturated dendrite tree). The neuron budding rules used in each computation step are indicated by their labels in the corresponding neurons, while the spiking rules are presented too. each neuron is indicated by an abstract triple (n, j, xn(j)), where 1 ≤j≤n, and associated with a truth value xnj. This way, a truth assignment of length nis represented by the nneurons (labels) of a subsystem. For instance, in a case n= 2 as described in Figure 5, 22= 4 different truth assignments of length 2 have been generated for two Boolean variables y1and y2 and each one is associated with a subsystem of layer n= 2. Technically, the first subsystem is composed of two neurons with labels c(2,1,f,t)and c(2,2,f,t), respectively. Whereas the former one associates with Boolean t(rue) value as x2= (f, t) and x2(1) = f, while the later one with f(alse) value as x2(2) = t, and then altogether forms an assignment (f, t); the case with other subsystems are the same where (t, t), (f, f), (t, f), respectively, are generated; one can see that the four truth assignments are well distinguished from each other by the layer structure of four subsystems (chains). Efficiency of Spiking Neural P Systems 247 The next synapse creation (budding) rules establish three synapses coming to the neurons which associate with a Boolean t(rue) value while four synapses to the neurons associated with f(alse) value, from the generator block. b1) (c(n,j−1,xn), c(n,j,xn))[ ]c(n,j,xn(j)=p)→(c(n,j,xn(j)=p), cei)[ ]ei, 1≤j≤n,p∈ {t, f}and s≤i≤3, where s= 1 if p=t,s= 0 if p=f, c(n,0,xn)=c(n−1,n,xn). Those neurons σc(n,j,xn(j)=t)whose interaction environment satisfies the condition of the rules b1) create three synapses coming from the neurons σcei, 1≤i≤3, while the neurons σc(n,j,xn(j)=f)establish synapses coming to it from the four neurons σcei, 0 ≤i≤3, of the generator block. The synapse creation rules of type b1) and the neuron budding rules of type a1) are applied to the same neurons in layer nat the same time in a consequent n−1 steps as their interaction environments coincide. Again looking at Figure 5, neuron σc(2,1,f,t)associates with f(alse) value gets 4 synapses from the generator as neuron σc(2,1,f,t)gets 3 cause its identity of t(rue) value. b2) (c(n,n−1,xn), c(n,n,xn))[ ]c(n,n,xn)→(c(n,n,xn),1)[ ]1. The rule of type b2) applies parallel to the last neurons of the layer nand produce the neurons σ1forming a new layer n+ 1. Meantime the rules of type b1) create synapses from the same neurons of layer nto the generator block at last. b3) (c(n,n,xn),1)[ ]1→(1,2)[ ]2, b4) (c(n,n,xn),1)[ ]1→(bi,1)[ ]bi, 1 ≤i≤3. As rules of type b3) apply to the neurons σ1and bud neurons σ2, rules of type b4) apply too and create three synapses coming from the neurons σbi, 1≤i≤3, to each neuron σ1. Thus, layer n+ 2 is formed. b5) (1,2)[ ]2→(2,3)[ ]3, b6) (1,2)[ ]2→(c, 2)[ ]c. The rules of types b5) and b6) apply simultaneously to a neuron σ2with a synapse (1,2). As a result, the former one buds a new neuron σ3, while the later one makes a new connection coming from the neuron σcas (c, 2). All other neurons σ2get the same effect by the rules as the maximal parallel applications of the rules. b7) (2,3)[ ]3→(3, out)[ ]out. The pre-computation of the SN P system structure construction is completed by forming the converged connections from the neurons σ3to the output neuron σout, by means the rules of type b7) are applied parallel to all neurons of layer n+ 3. Thus, the SN P system device structure totally empty of spikes which is to solve all the instances of sat(n, m), has been (pre-)computed in a polynomial time. The next computation stage (post-computation) to solve sat(n, m) is plugged-in as follows. 248 J. Wang, T.-O. Ishdorj, L. Pan Solving sat Any given instance γnof sat(n, m) is encoded in a sequence of spikes. Each clause Ciof γncan be seen as a disjunction of at most nliterals: for each j∈ {1,2,...,m}, either yjoccurs in Ci, or ¬yjoccurs, or none of them occurs. In order to distinguish these three situations we define the spike variables αij , for 1 ≤i≤mand 1 ≤ j≤n, as variables whose values are amounts of spikes, and we assign to them the following values: αij =     aif yjoccurs in Ci a2if ¬yjoccurs in Ci λotherwise. In this way, clause Ciwill be represented by the sequence αi1αi2···αin of spike variables; in order to represent the entire formula γnwe just concatenate the representations of the single clauses, thus obtaining the sequence α11α12 ···α1nα21α22 ···α2n···αm1αm2···αmn. As an example, the representation of γ3= (y1∨ ¬y2)∧(y1∨y3) is aa2λaλa. The spiking rules residing in the neurons of the system which perform for solving the introduced problem are listed below with a brief description for each. But we do not go detailed explanation of each rule functions here, we prefer it refer to Section 3 and the paper [6], also the neuron budding rules are out of usage in this stage. A given instance γn∈sat(n, m) encoded in a spike sequence is introduced into the system structure and will be processed by the spiking rules according to their roles in each step of the computation. c1) [a→a]csoma ; [ a2→a2]csoma ; [a→a;n2−n−1]s. At each computation step of introducing the input, we insert 0, 1 or 2 spikes into the system through the input neuron σsoma, according to the value of the spike variable αij we are considering in the representation of γn. Meantime we insert a single spike ainto neuron σsonce, which excites the generator block. c2) [a→a]c(k,j,xk); [ a2→a2]c(k,j,xk) 1≤k≤n−1,1≤j≤n, xk∈ {t, f}k. Each spike inserted into the input neuron is duplicated here and transmit along the first layer of the system towards next layers. When a spike passes a touching point – neuron with label of type c(k,n,xk), it is duplicated and enter into next layer, etc., finally 2ncopies of them will take place at the layer n. Once the copies of a clause are taken place on the neurons of the chains of length nin layer n, the combined functioning of the generator block and the layer ntests the assignments against each copy of the clause in consideration. For this purpose, the rules c3)−c5) are used. c3) [a→a]ei; 0 ≤i≤3, [a→a;n−1]+; [ a→a;n−1]−. Efficiency of Spiking Neural P Systems 249 The generator block and its spiking rules. The generator block provides 3 and 4 spikes, respectively, in each nsteps to the neurons associated with truth values tand f, of layer n, in order to test the satisfiability of the truth assignments against a clause which has been taken place through the layer. c4) [a→a]tt; [ a3→λ]tt; [ a2→a2]t1; [a4→a]tt; [a5→λ]tt; [ a2→a]t0; tt=c(n,j,xn(j)=t),1≤j≤n, t1=c(n,j,xn(j)=t),1≤j≤n−1, t0=c(n,n,xn(n)=t), xn∈ {t, f}n. The spiking rules residing in the neurons which associate Boolean t(rue) value in layer n. The rules a2→a2,a2→a, and a→aused to transmit the spike variables a,a2along the chains. Once a clause placed, each neuron associated with t(rue) value contains either of a spike aor a2or empty. As a spike variable arepresents a truth variable y, to which a spike true value a3sent by the generator is assigned and it results an yes answer as a4, then it passes to the neuron σ1along the chain with a saying that a truth variable of the clause is satisfied by true value of a truth assignment or simply the clause is satisfied by a truth assignment of the corresponding chain. Meanwhile, a true value a3 is assigned to the spike variables a2associates to truth variable ¬yand empty wherever, which give a result no by means the rules a3→λand a5→λare performed. c5) [a→a]ff; [ a4→λ]ff; [ a2→a2]f1; [a5→λ]ff; [a6→a]ff; [ a2→a]f0; ff=c(n,j,xn(j)=f),1≤j≤n, f1=c(n,j,xn(j)=f),1≤j≤n−1, f0=c(n,n,xn(n)=f), xn∈ {t, f}n. The spiking rules residing in the neurons which associate with Boolean f(alse) value in layer n. The functioning of the rules is similar as rules c5). c6) [a→a;n−1]bi; 1 ≤i≤3, [a2/a →a]1; [a3→λ]1; [a4→a]1; [a5→a]1. Whether an assignment satisfies the considered clause or not is checked by a combined functioning of the neurons with label 1 in layer n+1 and the neurons with label bi,1≤i≤3. c7) [a→λ]2; [ a2→a]2; [a→a]c. With a combined function of neuron σc, neuron σ2emits a spike into neuron σ3 if the corresponding assignment satisfies the under consideration clause here, otherwise no spike is emitted. c8) [am→a]3; [a+/a →a]out. Neurons with label 3 count how many clauses of the instance γnare satisfied by the associated truth assignments. If any of those neurons get mspikes, which fire, hence the number of spikes that reach neuron out is the number of 250 J. Wang, T.-O. Ishdorj, L. Pan assignments that satisfy all the clauses of γn. Thus, the output neuron fires if it has got at least one spike by means the problem has a positive solution, otherwise there is no assignment satisfies the instance γn. This stage of the computation ends at the n2+nm + 4th step. The entire system halts in total at most in 2(n2+nm + 4) number of computation steps. Thus, we got a full (deterministic, polynomial time and uniform) solution to sat(n,m) in the framework of SN P systems. 5 Discussion The present paper concerns the efficiency of SN P systems, we proposed a way to solve NP-complete problems, particularly sat, in polynomial time. Specifically, the neuron budding rule is introduced in the framework of SN P systems, which a new feature enhances the efficiency of the systems to generate necessary work space. Neuron budding rules drive the mechanism of neuron production and synapse creation according to the interaction of a neuron with its environment (described by its synapses connected to other neurons). A very restricted type of rule of neuron budding, at most one synapse is involved in an environment, is used, but it is powerful enough to solve the considered problem, sat. The solution to sat by SN P systems with neuron budding contains two computation stages: first, constructing the device structure which has no spikes inside, second, introducing the considered problem to be solved encoded in spikes into the device. The system works in deterministic and maximally parallel manner. The whole mechanism we considered here for solutions to computational intractable problems is elegant from computational complexity theory point of view as the designed algorithm can be computed by a deterministic Turing machine in polynomial time; the operation of neuron budding is well motivated by neural biology. We believe that SN P systems use the restricted budding rules can be an efficient computing tool to solve other NP hard problems. The SN P systems with neuron budding rules can be extended by introducing more general rules, which in some sense capture the dynamic interaction of neurons with their environment. One possible form of such general rules is as follows: Ai[ ]iBi→Cj[ ]jDj, where Ai, Biand Cj, Djare the set of synapses coming to and going out from, respectively, the specified neurons σiand σj. Clearly, in such general rules, more than one synapses are involved in the environment of the considered neuron. Acknowledgments The work of Tseren-Onolt Ishdorj was supported by the project “Biotarget”, it is a joint project between U. of Turku and ˚ Abo Akademi University funded by the Academy of Finland. The work of L. Pan was supported by National Natural Science Foundation of China (Grant Nos. 60674106, 30870826, Efficiency of Spiking Neural P Systems 251 60703047, and 60533010), Program for New Century Excellent Talents in University (NCET-05-0612), Ph.D. Programs Foundation of Ministry of Education of China (20060487014), Chenguang Program of Wuhan (200750731262), HUST-SRF (2007Z015A), and Natural Science Foundation of Hubei Province (2008CDB113 and 2008CDB180). References 1. H. Chen, M. Ionescu, T.-O. Ishdorj: On the efficiency of spiking neural P systems. Proc. 8th Inter. Conf. on Electronics, Information, and Communication, Ulanbator, Mongolia, June 2006, 49–52. 2. H. Chen, M. Ionescu, T.-O. Ishdorj, A. P˘aun, Gh. P˘aun, M.J. P´erez-Jim´enez: Spiking neural P systems with extended rules. Fourth Brainstorming Week on Membrane Computing (M.A. Guti´errez-Naranjo, Gh. P˘aun, A. Riscos-N´u˜nez, F.J. Romero- Campero, eds.), vol. I, RGNC Report 02/2006, Research Group on Natural Computing, Sevilla University, F´enix Editora, 2006, 241–266. 3. H. Chen, R. Freund, M. Ionescu, Gh. P˘aun, M.J. P´erez-Jim´enez: On string languages generated by spiking neural P systems. Fundamenta Informaticae 75 (2007), 141–162. 4. V. Danos, J. Feret, W. Fontana, R.Harmer and J. Krivine, Investigation of a Biological Repair Scheme, In Proceedings of 9th International Workshop on Membrane Computing, Edinburgh, UK, July 28-31, 2008, LNCS 5391, 1–12, 2009. 5. M. Ionescu, Gh. P˘aun, T. Yokomori: Spiking neural P systems. Fundamenta Informaticae, 71, 2–3 (2006), 279–308. 6. T.-O. Ishdorj, A. Leporati: Uniform solutions to SAT and 3-SAT by spiking neural P systems with pre-computed resources. Natural Computing, 7, 4 (2008), 519–534. 7. T.-O. Ishdorj, A. Leporati, L. Pan, X. Zhang, X. Zeng: Uniform solutions to QSAT and Q3SAT by spiking neural P systems with pre-computed resources. In the present volume. 8. M.A. Guti´errez-Naranjo, A. Leporati: Solving numerical NP-complete problems by spiking neural P systems with pre-computed resources. Proceedings of Sixth Brainstorming Weeks on Membrane Computing, Sevilla University, Sevilla, February 2-8, 2008, 193–210. 9. A. Leporati, C. Zandron, C. Ferretti, G. Mauri: Solving numerical NP-complete problem with spiking neural P systems. Membrane Computing, International Workshop, WMC8, Thessaloniki, Greece, 2007, Selected and Invited Papers (G. Eleftherakis, P. Kefalas, Gh. P˘aun, G. Rozenberg, A. Salomaa, eds.), LNCS 4860, Springer-Verlag, Berlin, 2007, 336–352. 10. A. Leporati, G. Mauri, C. Zandron, Gh. P˘aun, M. J. P´erez-Jim´enez: Uniform solutions to SAT and SUBset SUM by spiking neural P systems. Submitted. 11. A. Leporati, C. Zandron, C. Ferretti, G. Mauri: On the computational power of spiking neural P systems. International Journal of Unconventional Computing, 2007, in press. 12. Gh. P˘aun: Membrane Computing – An Introduction. Springer-Verlag, Berlin, 2002. 13. Gh. P˘aun: Twenty six research topics about spiking neural P systems. Fifth brainstorming week on membrane computing, Fenix Editora, Sevilla, 2007, 263-280. 14. M. Sipser: Introduction to the Theory of Computation. PWS Publishing Company, Boston, 1997. 15. Think and Grow Toys: http://www.tagtoys.com/dendrites.htm