Simulation of Spatial P system models

Simulation of Spatial P system models

Theoretical Computer Science 529 (2014) 11–45 Contents lists available at ScienceDirect Theoretical Computer Science www.elsevier.com/locate/tcs Si...

2MB Sizes 0 Downloads 11 Views

Theoretical Computer Science 529 (2014) 11–45

Contents lists available at ScienceDirect

Theoretical Computer Science www.elsevier.com/locate/tcs

Simulation of Spatial P system models Roberto Barbuti, Andrea Maggiolo-Schettini, Paolo Milazzo ∗ , Giovanni Pardini Dipartimento di Informatica, Università di Pisa, Italy

a r t i c l e

i n f o

Article history: Received 16 August 2012 Received in revised form 12 July 2013 Accepted 6 August 2013 Communicated by T. Yokomori Keywords: Membrane computing Simulation algorithms Spatiality

a b s t r a c t Spatial P systems are an extension of the P systems formalism in which objects and membranes are embedded into a two-dimensional discrete space. Spatial P systems are characterised by the distinction between ordinary objects and mutually exclusive objects, with the constraint that any position can accommodate any number of ordinary objects, and at most one mutually exclusive object. The presence of mutually exclusive objects makes the simulation of Spatial P system models more complex than that of standard P systems. In this paper, we present a polynomial-time algorithm for the simulation of a restricted version of Spatial P systems where the restriction consists in considering only mutually exclusive objects and rules having exactly one reactant and one product. This version of Spatial P systems, although very restricted, is expressive enough to model interesting biological systems. In particular, we show how it can be used to simulate two models describing different dynamics of fish populations, namely the dynamics of territorial fish and the formation and movement of herring schools. In addition, the simulation methodology we propose can be adapted to simulate richer versions of Spatial P systems. © 2013 Elsevier B.V. All rights reserved.

1. Introduction In the field of Membrane Computing, the P systems formalism [1] has been proposed as a distributed and parallel computing model inspired by the structure and the functioning of the living cell. A P system is composed of a hierarchy of membranes, each of them containing a multiset of objects, which are processed by evolution rules. Evolution rules allow the description of the behaviour of a model, for example by representing chemical reactions. The basic P systems model associates with each membrane a set of rules to be applied to the (only) objects contained in the membrane. A rule specifies reactants and products: when a rule is applied, the reactants are removed from the membrane and the products are either (i) added inside the membrane, or (ii) sent into an inner membrane, or (iii) sent outside the membrane. Many variants and extensions of P systems exist that include features which increase their expressiveness and which are based on different evolution strategies. In particular, we cite P systems with promoters and inhibitors [2], and P systems with priority among evolution rules [1]. In P systems with promoters and inhibitors, the applicability of each rule can be constrained by the presence or absence of other objects in the membrane. In particular, the multiset of promoters represent the objects which are needed for the rule to be applicable. On the contrary, the multiset of inhibitors represent those objects whose presence effectively blocks the application of the rule. With regard to P systems with priority among evolution rules, a rule is applicable in a membrane only if no higher priority rule is applicable, namely only if each higher priority rule

*

Corresponding author. E-mail addresses: [email protected] (R. Barbuti), [email protected] (A. Maggiolo-Schettini), [email protected] (P. Milazzo), [email protected] (G. Pardini). 0304-3975/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.tcs.2013.08.002

12

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

requires a multiset of reactants which are not present in the membrane. See [3] for the definition of variants of P systems and [4] for a complete bibliography. In [5] we proposed Spatial P systems, an extension of P systems in which objects and membranes are embedded in a two-dimensional discrete space. Objects are associated with precise positions inside membranes, and evolve by means of the application of evolution rules. As for standard P systems, rules specify the objects which are consumed and the ones which are produced. In particular, rules specify the relative positions of the various objects involved, with respect to the position in which the rule is applied. A feature of Spatial P systems is the distinction between ordinary objects and mutually exclusive objects. Every position inside a membrane can accommodate an arbitrary number of ordinary objects, but at most one mutually exclusive object. In this paper we consider a revised version of Spatial P systems that extends the one given in [5] by including more general evolution rules, rule promoters, a notion of priority, and arbitrary shapes of membranes. These extensions allow a more natural modelling of systems and phenomena of interest. In particular, Spatial P systems can be profitably used to model and analyse the dynamics of populations (of biological or artificial entities) at the individual level and by taking into account spatial and environmental aspects of such systems. The aim of the paper is to investigate simulation algorithms for Spatial P systems. Simulation is probably the most effective analysis approach for individual-based population models since it can give precise model-based descriptions of the population dynamics. The definition of an efficient simulation algorithm for Spatial P systems is not an easy task. In particular, the presence of mutually exclusive objects makes it rather complex to simulate an evolution step in accordance with the semantics of the formalism. We present a polynomial-time algorithm for the simulation of a restricted version of Spatial P systems where the restriction consists in considering only mutually exclusive objects and rules having exactly one reactant and one product (and an arbitrary number of promoters). We call these rules singular. Focusing on singular rules will allow us to work on a simulation strategy tailored to the handling of mutually exclusive objects in space. This is the critical aspect of Spatial P systems simulation and the algorithm we develop for the restricted case of singular rules will be the basis for the definition of algorithms dealing with richer forms of evolution rule. In this paper we discuss how the algorithm can be extended to more general cases and, in particular, we hypothesise that a simulation algorithm for the general formalism would be NP-complete. However, a deep investigation of these aspects is left as future work. Though restricted, Spatial P systems including singular rules only are expressive enough to model interesting systems. As case studies we consider the dynamics of territorial fish and the schooling behaviour of herring. We construct Spatial P system models of such phenomena by only using singular rules. A simulation tool for both the model has been developed by implementing the proposed polynomial-time simulation algorithm. Such a tool will allow us to test the efficiency of the simulation algorithm in practice on a model of relevant size. In addition, we exploit the possibility of Spatial P system models of representing membranes of any shape in order to simulate the behaviour of herring in the presence of obstacles and of spatial constraints. Related work. Spatial modelling allows the description of the position of the elements of a system, with respect to some more or less abstract notion of space. Various formalisms proposed in the literature allow different levels of representation of the space. The simplest form of spatial modelling is given by the ability of representing compartments, which allow the representation of different locations for the elements of the system. In the biological setting, compartments are often entailed by membranes. Since membranes play an important role in many biological processes, most computer science formalisms for biology allow some form of compartmental modelling (e.g. [6,7]). In order to faithfully describe some biological systems, less abstract representations of space than the one provided by compartments, are needed. Beyond computer science, one of the earliest approaches in this direction involved the use of differential equations models, which allow the description of many different spatially-aware biological processes, such as reaction diffusion systems. Models of this kind and their extensions have been used, for example, as the basis for models of insect dispersal and chemotaxis [8]. Other kinds of spatial models involve spatial pattern formation [8,9]. Since these models are based on continuous variables, they may not be adequate to describe discrete biological entities. In computer science, a formalism including spatial features are Cellular Automata [10], initially devised as a computational formalism inspired by biological behaviours. A cellular automaton is composed of a finite grid of cells, where each cell has an associated state taken from a finite set of different states. Time is discrete, and at each step of the evolution all the cells change state in accordance with a rule, which is characteristic of the particular cellular automaton model. The rule is deterministic and “local”, in the sense that the new state of a cell is determined only on the basis of the previous states of the cell itself and of nearby cells. Many variants of Cellular Automata have been defined. Cellular Automata have been developed as a computational tool inspired by biological behaviours. However, even if their focus has not been in the modelling of biological systems, they have been later used for this purpose. In this view, Cellular Automata are particularly suitable for describing the evolution of populations of many similar entities whose behaviour is based on local interactions. For example, they have been used to describe tumour growth [11], showing how a simple rule driving the evolution can lead to a complex behaviour of the whole population. Other extensions of cellular automata, for example including stochasticity, have also been proposed [11]. Spatial P systems use a representation of the space analogous to that of Cellular Automata. However, in Spatial P systems, cells contain objects, and interactions occur between objects. For this reason, each position in a Spatial P system can accommodate any number of objects, thus representing a possibly infinite number of different “states”.

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

13

The Spatial Calculus of Looping Sequences (Spatial CLS) [12] is an extension of the Calculus of Looping Sequences (CLS), a formalism geared towards the modelling of cellular systems. Spatial CLS is based on term rewriting and allows a direct representation of membranes in the syntax of the calculus. The evolution of a system is modelled by rewrite rules, which can be used to describe both biochemical reactions and more complex behaviours, such as rearrangements in the membrane structure of the system. Spatial CLS models are based on two/three-dimensional continuous space, and allow an accurate description of the motion of the elements and of their size. In particular, Spatial CLS allows the description of the space occupied by elements and membranes, which can change their sizes dynamically as the system evolves. Space conflicts which may occur can be resolved by performing a rearrangement of the elements and the membranes. The Calculus of Wrapped Compartments (CWC) [13] is another formalism for the description of biological systems based on a form of term rewriting. CWC has been used in [14] to model spatial aspects of the plant-fungus interactions that happen in the context of the Arbuscular Mycorrhizal symbiosis. Moreover, a spatial extension of CWC has been proposed in [15] by embedding CWC in a surface language. Such a language can be used to model both spatial transformations (in a two-dimensional discrete space) and biochemical transformations. Language constructs can be translated into standard CWC models to enable stochastic simulation. SpacePI [16] is an extension of the π -calculus with (continuous) space and time. In SpacePI, positions in a continuous space, such as R2 , are associated with processes, and processes can move autonomously according to a movement function. In SpacePI, processes move uniformly, namely on a constant direction with constant speed, during fixed time intervals, at the end of which they can change direction and speed according to their movement function. The calculus is deterministic and processes may, and are required to, communicate only when they are close enough, with respect to their Euclidean distance. In particular, an interaction radius is associated with each action and co-action. Two processes communicate as soon as the sum of the interaction radii of their complementary actions are equal to the distance between processes. SpacePI does not allow the direct modelling of compartments and membranes. Another extension of the π -calculus with spatial features is the 3π process algebra [17], where processes are embedded in a 3D space. The calculus is based on affine geometry, and processes can interact by exchanging channel names or geometric data. The movement of a process is realised by applying a frame shift operation (described by an affine map) to the process itself. A peculiarity of the calculus is that a process has no visibility of its location in the global space, but can nevertheless communicate its position to other processes which can be used, for example, to compute the distance between two processes. Membranes and compartments cannot be described explicitly. In the field of ecological modelling, the Process Algebra with Locations for Population Systems (PALPS) [18] has been proposed as a domain-specific process calculus for the modelling of population systems. In PALPS, each individual is modelled by a process, each possessing a species and a location. Locations are related by means of a neighbourhood relation, and individuals can either move non-deterministically through nearby locations, or be involved in other actions, such as preying. The behaviour of individuals can be determined by local conditions, for example by taking into account the number of individuals present in a location. Also P systems have been applied in the field of Ecology [19–21]. The extension of P systems proposed in [19], based on a probabilistic extension of P systems introduced in [22], is used to describe the behaviour of metapopulations. Metapopulations are local populations living in spatially separated areas (called patches), where populations can interact, and individuals can disperse from a patch to nearby patches. Models for metapopulations aim at discovering how the fragmented habitat influences local and global population behaviour. In the proposed model, objects are used to model different species (predators and preys), and patches are represented as elementary membranes in a flat membrane structure. Patches form the nodes of an undirected weighted graph, with edges describing a neighbourhood relation between patches, modelling spatial proximity. Costs associated with edges model the effort of individuals for migrating from a patch to another. Another extension of P systems, called multienvironment P systems [23], has been applied to the modelling of ecosystems [24,25]. A multienvironment P system model is composed of a set of environments, each containing a P system, from a particular class of probabilistic P systems. All the P systems in environments have the same skeleton, namely they use the same alphabet of symbols, membrane structure, and sets of evolution rules. Objects in an environment, i.e. outside the skin membrane of a P system, can move from an environment to another by means of the application of special communication rules. The movement of objects to other environments, as well as the internal behaviour of P systems, is driven by time-varying probabilistic values associated with the rules. In Tissue-like P systems [26,27] the membrane hierarchy of ordinary P systems is replaced by a graph whose nodes are elementary membranes, called cells, and the edges describe the communication channels between cells. Such a device is useful to describe systems composed of many entities living together, such as in tissues, organs and organisms. Tissue-like P systems have been extended to Population P systems [28], which provide mechanisms to dynamically modify the structure of the communication graph among cells. Lattice Population P systems [29,30] are a variant of Population P systems in which space is represented as a finite regular lattice where each cell of the lattice contains a distinct stochastic P system describing the behaviour of an individual. Translocation rules describe the possible movements of objects from a cell to another. All the above mentioned variants of P systems provide a space representation which is more abstract than that of Spatial P systems, since they are ultimately based on a graph of cells among which objects are able to move. More importantly, they do not provide any ability to model the space occupied by objects, as instead it is possible in Spatial P systems by exploiting mutually-exclusive objects. At the end of the paper we show, by means of the two case studies presented, how Spatial P systems can be used effectively for the modelling of ecosystems, by exploiting the features provided.

14

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

An important field of application of P systems is that of the formal modelling of multi-agent systems, namely highlydynamic and parallel systems formed by “agents” which collaborate and communicate to achieve specific global goals. The applicability of variants of P systems as the formal basis behind the modelling of multi-agent systems has been investigated in [31]. In particular, it has been observed that Membrane Computing systems (which include P systems and their variants) are well-suited to model reactions to external stimuli, while they are somewhat less suited to model the internal logic of agents. Being an extension of P systems, sharing all the basic features with them, such observations also extends to Spatial P systems. Structure of the paper. The paper is structured as follows. In Section 2 we recall the definition of P systems. In Section 3 we present the revised version of Spatial P systems we consider for the definition of the simulation algorithms. In Section 4 we propose a simulation algorithm for Spatial P systems with singular rules, we prove its correctness with respect to the semantics of the formalism and we show that the execution of a step of the algorithm, for a given Spatial P system model, is polynomial with respect to the size of the skin membrane, and to the number of objects, rules and promoters. In Section 5 we discuss extensions of the simulation algorithms to less restricted variants of Spatial P systems. In Section 6.2 we model and simulate the case studies. Finally, in Section 7 we draw our conclusions and discuss future work. 2. P systems In this section we recall the definition of P systems [3]. In the following, we often denote multisets over a finite alphabet as strings. More precisely, let V ∗ be the set of all strings over an alphabet V , including the empty one, denoted by λ. In a string representing a multiset we use exponents to represent multiple instances of the same elements in a concise way. For instance, with a2 b3 we denote a multiset consisting of two instances of a and three instances of b. Moreover, given a multiset M, the multiplicity of an element x in M is denoted as m(M, x). For instance, m(a2 b3 , b) = 3. We use the following notations for operations on multisets, where the multiplicity is significant. The union operator between multisets is denoted by , while the difference operator is denoted by \\. Inclusion between two multisets A, B is denoted by AFB , meaning that A is a submultiset of B . Proper inclusion, namely if AFB and A = B , is denoted by D. A P system is composed of a hierarchy of membranes. Each membrane is labelled univocally by a natural number and contains a multiset of objects and a set of evolution rules. Evolution rules describe how the objects of the system evolve, for example they can be used to describe chemical reactions, i.e. rules in which some objects interact and, as a result, they are transformed into some other objects. Given a membrane m, its evolution rules in the set R m can be applied only to the objects contained in the same membrane, and not in any other membrane. Formal semantics of different versions of P systems are presented in [32–36]. An evolution rule is of the form u → v, where u and v are multisets whose elements are called reactants and products, respectively. When a rule is applied, the reactants are removed from the membrane and the products are added to the target membrane, which could be a different membrane than the one in which the rule is applied. Given a membrane m, the products of a rule associated with m are described by a multiset of tuples of the following forms: (a, here), usually written as ahere , meaning that the object a is added to the same membrane m; (a, out), usually written as aout , meaning that the object a is to be sent out of the membrane; (a, inx ), usually written as ainx , meaning that the object a is to be sent into the child membrane labelled by x. Formally, let TAR be the set of object targets {here, out} ∪ {ini | i ∈ N}, and ∗ . The target indication here is often omitted, let V tar = V × TAR. An evolution rule u → v is such that u ∈ V ∗ and v ∈ V tar ∗ therefore a rule of the form u → v where v ∈ V is to be interpreted as if all the objects v have target indication here. The number of reactants of a rule, that is, the length of u, is called the radius of such a rule. An evolution rule is said to be cooperative if its radius is greater than 1, otherwise the rule is called non-cooperative. This naming is also extended to P system models, that is, a non-cooperative P system is such that all its rules are non-cooperative, otherwise it is a cooperative P system. A formal definition of P systems follows. Definition 2.1. A P system is a tuple Π = ( V , μ, w 1 , . . . , w n , R 1 , . . . , R n ) where:

• V is a finite alphabet whose elements are called objects; • μ ⊂ N × N describes the tree-structure of membranes, where (i , j ) ∈ μ denotes that the membrane labelled by j is contained in the membrane labelled by i;

• w i , with 1  i  n, are strings from V ∗ representing multisets over V associated with membranes 1, 2, . . . , n of μ; • R i , with 1  i  n, are finite sets of evolution rules associated with membranes 1, 2, . . . , n of μ. A characteristic of P systems is the way in which rules are applied in each step, namely with maximal parallelism. In each step, evolution rules are applied in a maximal non-deter ministic way in all membranes, that is, in each membrane, a multiset of rules is selected non-deterministically to consume the membrane objects, in such a way that no other rule can be applied to the objects which are not involved in any rule application. A configuration of a P system is given by an association of its membranes with multisets of objects. The multisets of objects w 1 , . . . , w n associated with membranes in the definition of a P system Π represent the initial configuration

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

15

Fig. 1. An example of P system model.

of Π . A computation is a sequence of transitions between configurations of a given P system Π starting from the initial configuration. Each transition of a computation describes a maximally parallel step. A computation is successful if and only if it reaches a configuration in which no rule is applicable. The result of a successful computation can be defined in different ways. Unsuccessful computations are those computations which never halt, thus yielding no result. Example 2.1. Fig. 1 depicts a P system with two membranes, labelled 1 and 2. The rules r1 and r2 are associated with 2 causes a copy of membrane 1, while membrane 2 has no rules associated with it. An application of rule r1 = a → ab in2 c in 2 object b and two copies of object c to be sent into the inner membrane 2. The object a is still present after the application, 3 , instead, can be applied to a pair of objects a, and since it appears in the right-hand part of the rule. Rule r2 = a2 → c in 2 results in sending three copies of the object c into membrane 2. The initial state, as depicted, contains two copies of object a in membrane 1, and no objects in membrane 2. At the first step, either rule r1 or r2 is applied. In fact, both the rules are enabled, since their reactants are present in the membrane. Actually, if r1 is applied to an object a, then the maximality requires it to be applied also to the other copy of a. This application sends the objects b2 c 4 into membrane 2. The objects contained in membrane 1 remain aa after the application, therefore the double application of rule r1 can be repeated in the subsequent step. Whenever rule r2 is applied, it causes the two copies of a in membrane 1 to disappear, thus terminating the computation. In such a case, the objects ccc are sent into membrane 2. Therefore, any computation of this P system is composed of a sequences of steps in which only r1 is applied (twice per step), followed by a last step in which rule r2 is applied once. Therefore, whenever the P system terminates, membrane 2 contains a multiset of objects bk c 2k+3 , for some k  0. P systems are often equipped with a priority relation among evolution rules [3]. The priority among rules is described by the partial order relations ρi ⊆ R i × R i , for each i ∈ {1, . . . , n}. Given the partial order ρ , for some membrane, a pair (r1 , r2 ) ∈ ρ , with r1 = r2 , means that the rule r2 has higher priority than r1 , and is usually written as r1 < r2 . The applicability of rules is constrained as follows. A rule r1 can be applied only if no other higher priority rule r2 > r1 is enabled in the same step, that is, if the reactants for r2 are not available at the beginning of the step. Note that, even if the reactants of rule r2 are present, it can happen that such a rule is not applied, for example because there is another non-comparable rule r 3 (i.e. neither r2 < r3 nor r3 < r2 ) competing with r2 for the same reactants. Even in such a case, any lower priority rule r1 < r2 still cannot be applied. Example 2.2. Consider the rules r1 = ab → c, r2 = a → d, r3 = b → e, with relative priorities r1 > r3 , namely ρ = {(r1 , r1 ), (r2 , r2 ), (r3 , r3 ), (r3 , r1 )}. Given a multiset of objects bb, the only possible state that can be obtained after one step is ee, resulting from the application of rule r3 twice. Note that r3 could be used since rule r1 could not, as the reactants ab are not present in the given state. Instead, given the objects aab, the following states can be obtained: cd, by applying both r1 and r2 once; and ddb, by applying r2 twice. Note that rule r3 cannot be applied, even if the object b is not used by any rule, since the higher priority rule r1 is enabled in the initial state aab. Other features of P systems that are often considered are promoters and inhibitors [2,3]. In this case evolution rules are extended to specify, for each rule, two multisets of objects denoted as promoters and inhibitors. They are used to constrain the applicability of a rule on the basis of the objects which are present in the membrane. Precisely, a rule can be applied only if all the promoters, with their multiplicities, are present in the membrane, and not all of the inhibitors are present in the membrane [2]. A rule with promoters and inhibitors is usually denoted as u → v |x1 ,...,xk ,¬ y 1 ,...,¬ yh , where u ∈ V ∗ and ∗ denote the reactants and the products, respectively, while each x ∈ V denotes a promoter, and each y denotes an v ∈ V tar i j inhibitor. Formally, given a multiset of objects w contained in a membrane, the rule u → v |x1 ,...,xk ,¬ y 1 ,...,¬ yh can be applied / w \\u. Note that objects corresponding to promoters and inhibitors only if ∀i = 1, . . . , k. xi ∈ w \\u and ∃ j ∈ {1, . . . , h}. y j ∈ can evolve independently, as they can be consumed by the application of other rules. For example, a rule a → b|c means that an object a can be transformed to b only if there is at least one object c present in the membrane. Conversely, the rule aa → b|¬b means that two copies of a can be transformed to b only if no b is present.

16

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

Fig. 2. The spatial structure of membranes and objects of a Spatial P system, with the reference coordinate system shown.

3. Spatial P systems Spatial P systems, originally proposed in [5], are a spatial extension of P systems in which membranes and objects are embedded in a two-dimensional discrete space with natural coordinates Z2 . Each object in a Spatial P system model is associated with a position inside a membrane. Evolution rules are associated with membranes, and are also extended to allow specifying the positions of reactants and products with respect to the actual position in which the rule is to be applied. A distinctive feature of Spatial P systems are the mutually exclusive objects, which are subject to the constraint that there can be at most one of them in each position. There are no constraints on the other (ordinary) objects. The distinction between ordinary and mutually exclusive objects can be thought of as an abstraction of the real-world space occupied by the objects. As with standard P systems, at each step rules are applied with maximal parallelism. From the theoretical point of view, the feature of mutually-exclusive objects provided by Spatial P systems is a powerful feature. In fact, as shown in [5], Spatial P systems using only non-cooperative evolution rules (namely rules in which there is only one object in their left-hand part) achieve Turing universality, by exploiting the semantics of mutually exclusive objects. In this paper we consider a revised and extended version of Spatial P systems. The most important differences with respect to the original definition given in [5] are: (i) membranes are no more restricted to have rectangular shape, but instead can have more liberal shapes; (ii) the format of evolution rules is a generalisation of the one presented previously, with also the ability to denote promoters; (iii) priority among rules can be described by using a partial order relation. Example 3.1. We present a simple model of bacteria which eat food dispersed in an environment. The aim of this example is just to help the reader gain some insights on the functioning of Spatial P systems, before presenting the formalisation of the structure and semantics in the following sections. Fig. 2a depicts both (i) the membrane structure, and (ii) the position of the objects contained in membranes, which form a particular configuration. Similarly to standard P systems, the outer skin membrane is labelled with 1, and contains all other membranes and objects. In the example, membrane 1 contains membrane 2, besides a number of objects of kind b, each denoting a bacterium, and f , representing single food molecules. Moreover, we consider objects b as mutually-exclusive objects, while f is an ordinary object. The idea is that bacteria constitute much bigger entities than food, therefore no more than one bacterium can be present in a position at any time. In this way, interpenetration of bacteria is automatically prevented by the semantics of Spatial P systems. On the other hand, an object f is the formal counterpart of some small molecule ingested by bacteria as food, thus they are abstracted as ordinary non-mutually-exclusive objects, for which many of them can occur in the same position. The figure also shows the reference coordinate systems used, in which grid cells correspond to the possible (discrete) positions for the objects; the bottom-left cell corresponds to position (0, 0). The position of objects in a configuration is always expressed with respect to the reference coordinate system. Therefore, the depicted configuration is composed of three copies of object b, at positions (1, 3), (2, 2), (6, 2), and six copies of object f , one for each position (2, 0), (3, 0), (4, 0), (6, 2), and two copies in position (6, 3). Finally, note that membrane bounds have no thickness, thus they conceptually occur in between adjacent cells. The version of Spatial P systems we consider here includes both promoters and priorities. Promoters are as in standard P systems and are useful for describing rules in which the presence of some objects enables the rule to be applied. Priorities among evolution rules, instead, have a different meaning with respect to the usual semantics of P systems. Recall that, in

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

17

standard P systems, a rule can be applied in a membrane only if there is no higher priority rule which is enabled, that is, the reactants of each higher priority rule are not present in the membrane region at the beginning of the step. Instead, in Spatial P systems we permit lower priority rules to be applied even if there are higher priority rules whose reactants are present. This difference is motivated by the fact that, in Spatial P systems, because of the presence of mutually-exclusive objects, it can happen that a rule for which the reactants are present cannot be actually applied, since that would create a conflict among mutually-exclusive objects. Precisely, the semantics of Spatial P systems requires that there does not exist a set of positions in which it is possible to replace the multisets of rules selected in each of them with “greater” multisets, and still obtain a selection of rules which is applicable with respect to the reactants, and which does not create conflicts among mutually-exclusive objects. In this context, a greater multiset means increasing the multiplicity of application of any rule r, given the possibility of cancelling the application of any lower priority rule r < r. Example 3.2. (Continues from Example 3.1.) In order to keep the model simple, we constrain bacteria to move only one cell at a time, in either of the directions up/down/left/right. When there are no f in adjacent positions, then the bacterium moves randomly in either one of the four directions. Otherwise, it moves over the adjacent cell which contains an f (or one of them non-deterministically in case of multiple adjacent f ), and then eats it. Fig. 2b shows the set of evolution rules associated with membrane 1 from Fig. 2a. These rules are just a minimal set of rules which model the intended behaviour. However, a more complete model would require further rules to deal with corner cases, which we do not consider in this example. (For instance, when two bacteria want to move to a same cell with food, since only one of them succeeds, what is the expected behaviour of the other? should it stay in the same position or wander around?) The form of rules is u Jπ K → v, where u denote reactants, π are the promoters, and v are the products. Each object appearing in a rule is specified as a p , where p is the relative coordinate of the object with respect to the actual position in which the rule is going to be applied. An important aspect is that a rule can be applied any number of times to any position inside a membrane. Albeit not used in this example, keep in mind that a rule could allow arbitrary long jumps for an object from one position to another. In this example, for the sake of simplicity, the movement of bacteria is limited to one cell at a time, thus we do not exploit such an ability. Rules from Fig. 2b can be considered as partitioned into three sets, ordered with respect to their priority. The lowest priority rules {r6 , r7 , r8 , r9 } model the random movement in case adjacent cells are free from food. In particular, rules r6 , r8 capture the left/right movement, while rules r7 , r9 capture the upward/downward movement, respectively. Note that the positions appearing in the rules is relative to the actual position in which the rule is applied. For example, rule r9 can be applied to object b in position (1, 3), by causing it to move into position (1, 3) + (0, −1) = (1, 2). Such a rule could also be applied to the other b in (2, 2), moving it to (2, 1), but not to object b in position (6, 2), since the resulting position (6, 1) would be outside the membrane. As for standard P systems, in Spatial P systems membrane crossing is possible only by using specific operators in/out (not present in this example). The set of higher priority rules {r2 , r3 , r4 , r5 }, instead, captures the behaviour when there is food in one (or more) of the adjacent cells. These rules use promoters to check for the presence of the movement of food in an adjacent cell, and allow moving the bacterium into such a cell. Since f appears only as a promoter and not as a reactant, such an object is not deleted upon rule application. Note that, since bacteria are modelled as mutually-exclusive objects, if there are more bacteria willing to move into a same cell containing food, only one of them succeeds. Finally, the highest priority rule r1 models the eating of a single f by a bacterium occurring in the same position. For example, as regards the configuration depicted in Fig. 2a, rule r1 is the highest priority rule applicable to position (6, 2), causing its content to be transformed from bf to b. 3.1. Formal definition Membranes. In order to formally define the shape of a membrane, we consider the bottom-left vertex of each cell of the discrete two-dimensional space. Two vertexes are adjacent iff they correspond to bottom-left vertices of adjacent cells. In this way, a membrane can be described by the circular list of vertices through which it passes, where each vertex has to be adjacent to both the previous and the next position in the list. Membranes are defined as simple closed curves composed only of horizontal and vertical edges, and need to be properly nested. Sibling membranes must not overlap, and membranes cannot exceed the bounds of their parent membrane. The formal description of membranes in a Spatial P system model is composed of the following elements, where we assume that membranes are labelled by {1, . . . , n}:

• a tree structure describing the containment hierarchy among membranes, represented by the set μ ⊂ {1, . . . , n} × {1, . . . , n}, where (i , j ) ∈ μ denotes that the membrane labelled by j is contained in the membrane labelled by i; • a function σ : {1, . . . , n} → Z2 , which describes membrane bounds by associating to each membrane i ∈ {1, . . . , n} a sequence of distinct vertices σ (i ) = p 1 , . . . , p s  through which the membrane edges pass. In order to ensure that the membrane is composed only of horizontal and vertical edges, each vertex p i must be adjacent to the next one p i +1 . In addition, the last vertex p s must be adjacent to the first one p 1 . Adjacency of vertices and the

18

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

condition that all of them must be distinct ensures that the membrane defines a simple closed curve. We also require each membrane to be properly contained in its parent membrane (except for the outer membrane), and that vertices forming membrane bounds of different membranes are distinct: ∀i 1 , i 2 . i 1 = i 2 ⇒ σ (i 1 ) ∩ σ (i 2 ) = ∅. This last constraint ensures that different membranes do not overlap. The structure of membranes can be seen as a partition of the space bounded by the skin membrane, where a position belongs to the region of a membrane if and only if it is contained within its bounds and not contained in any other child membrane. The set of all positions belonging to a membrane is called a region. Let Extents(i ) denote the set of positions inside the closed curve denoted by σ (i ). The set of positions forming the region of a membrane is defined by the function  Region : {1, . . . , n} → P (Z2 ) as Region(i ) = Extents(i ) \ (i , j )∈μ Extents( j ). The constraint that each membrane is properly contained in its parent membrane formally corresponds to requiring that ∀(i , j ) ∈ μ. Extents( j ) ⊂ Extents(i ). For the sake of simplicity, we also define Region for the special value 0, representing all the positions outside the skin membrane, as Region(0) = Z2 \ Extents(1). Example 3.3. As regards membranes depicted in Fig. 2a, their containment hierarchy is formally defined by the tree structure μ = {(1, 2)}, while their bounds are described by the following sequences of vertices: σ (1) = (2, 0), (2, 1), (1, 1), (1, 2), (0, 2), (0, 4), (1, 4), (1, 5), (2, 5), (2, 6), (5, 6), (5, 5), (6, 5), (6, 4), (7, 4), (7, 2), (6, 2), (6, 1), (5, 1), (5, 0); σ (2) = (4, 2), (4, 3), (4, 4), (5, 4), (5, 3), (5, 2). Moreover, the extents of membrane 1 correspond to the positions inside its bounds, including those of membrane 2. In particular, Extents(2) = {(4, 3), (4, 2)} corresponds to the positions coloured in grey in Fig. 2a, and it is such that Extents(2) ⊂ Extents(1). On the other hand, their regions are distinct, since Region(2) = Extents(2) and Region(1) = Extents(1) \ Region(2). Objects and rules. Each object in a Spatial P system model is associated with a position in the region of a membrane. The two kinds of objects, ordinary objects and mutually exclusive (ME) objects, are represented by two disjoint sets V and E, respectively. Evolution rules are composed of reactants, promoters and products. A set of evolution rules is associated with each membrane, and each rule from such a set can be applied any number of times to any position inside membrane region. Each object appearing in a rule carries a coordinate which is interpreted with respect to the actual position in membrane region where the rule is going to be applied. Therefore, in order to apply a rule in a position, the relative positions appearing in the rule have to be instantiated to actual positions inside membrane region. In particular, in order for the rule to be applied to a position in membrane region, both its reactants and promoters must be present in their relative positions. Similarly, the resulting position of products refers to the actual position where rule is applied. Formally, evolution rules have the following form:

q y (u 1 ) p 1 . . . (uk ) pk (π1 )q1 . . . (πγ )qγ → ( v 1 )t1 . . . ( v h )th where {u i }i , {πi }i , { v i }i denote reactants, promoters, and products, respectively. In particular, (u 1 ) p 1 , . . . , (uk ) pk denote the reactants as strings of objects u i , each of them specifying a relative position p i ∈ Z2 . Reactants are consumed by the application of the rule, and each object in a configuration can be consumed only by one rule application. Positions p 1 , . . . , pk are assumed to be distinct. On the other hand, (π1 )q1 , . . . , (πγ )qγ denote the promoters of the rule, whose presence is necessary for the rule to be applicable, but which are not consumed by its application. Similarly to the previous case, each πi is a string of objects with relative position q i ∈ Z2 , and q1 , . . . , qγ are assumed to be distinct. Finally, ( v 1 )t1 , . . . , ( v h )th represent the products of the rule, where each v i is a string of objects, and each t i can take one of the following forms:

• t i ∈ Z2 , denoting a relative movement within membrane region; • outd , with d ∈ Z2 , denoting outwards membrane crossing towards a position in the region of the parent membrane, if any, while objects sent out of the skin membrane disappear from the system;

• in j ,d , with j ∈ {1, . . . , n} and d ∈ Z2 , denoting the movement towards a position in the region of the child membrane j. Also in this case, t 1 , . . . , th are assumed to be distinct. A generic list of reactants (u 1 ) p 1 . . . (uk ) pk , promoters (π1 )q1 . . . (πγ )qγ , and products ( v 1 )t1 . . . ( v h )th may be denoted by

k

a symbol u, π , or v, respectively. The radius of a rule is defined as the number of reactants, namely i =1 |u i |. An evolution rule without promoters is said to be cooperative if its radius is greater than 1, otherwise the rule is called non-cooperative. Moreover, a rewrite rule r is called singular iff it contains exactly one reactant and one product. Note that the definition of evolution rules we have given does not include inhibitors. General definition. In order to present the formal definition of Spatial P systems, we need to formally introduce the notion of configuration of a Spatial P system, which extends the notion of configuration of P systems with positions. The presence of the mutually exclusive objects motivates also the definition of a notion of conflict-free configuration, namely a configuration satisfying the mutual exclusion constraints these objects are subjected to. In the following of the paper, positions which refer to the global reference coordinate system are termed absolute positions, to distinguish them from relative positions used only in rules, which instead refer to the actual position of application.

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

19

Definition 3.1 (Configuration). A configuration W : Extents(1) → ( V ∪ E )∗ of a Spatial P system is a mapping associating each position in the extents of the skin membrane of the considered Spatial P system with a multiset of objects. Given a position (x, y ) ∈ Extents(1), the multiset of objects associated with it by W is denoted w (x, y) . A configuration describes the contents of each position within the extents of the skin membrane, with a mapping from absolute positions in Extents(1), i.e. relative to the reference coordinate system, to a multiset of objects. Definition 3.2 (Conflict-free). A configuration W of a Spatial P system is conflict-free, denoted conflict-free( W ), iff each posi tion contains at most one mutually exclusive object, namely ∀ p ∈ Extents(1). x∈ E m( w p , x)  1. Definition 3.3 (Spatial P system). A Spatial P system is a tuple ( V , E , μ, σ , W 1 , ( R 1 , 1 ), . . . , ( R n , n )) where:

• • • • •

V and E are disjoint alphabets, denoting ordinary objects and mutually exclusive objects respectively;

μ ⊂ {1, . . . , n} × {1, . . . , n} describes the tree-structure of membranes, as explained before; σ : {1, . . . , n} → Z2  describes the shape and position of each membrane i ∈ {1, . . . , n}; W 1 is the initial configuration of the system, which is assumed to be conflict-free;

( R i , i ), for i ∈ {1, . . . , n}, describe the finite sets of evolution rules R i associated with the membranes, together with a partial order relation i on R i describing the priorities among the rules; in case the priorities are not used, namely if, for every membrane i, the partial order relation is the identity relation I R i , then its indication can be omitted.

As for P systems, computations of Spatial P systems consist of steps in which evolution rules are applied with maximal parallelism. The (non-deterministic) choice of the rules to be applied at each step takes into account the priorities among rules. On the other hand, since each rule can be applied to different positions, the coordinates syntactically appearing in it are to be interpreted with the actual absolute position in which the rule is applied. In other words, positions occurring in rules are relative to the actual absolute position of application. It is worth noting that priorities in Spatial P systems are not as in standard P systems. In Spatial P systems priorities act as constraints for the application of rule instances, rather than of rules. This means that an instance of a low-priority rule can be applied if it does not compete for available reactants with instances of the higher priority rules. As an example, let us consider two evolution rules r1 = (ab)(0,0) → (c )(0,0) and r2 = (a)(0,0) → (d)(0,0) , and a priority relation such that r1 > r2 . Moreover, let us assume that a position of the Spatial P system contains the multiset of objects a4 b2 . The interpretation of priorities in Spatial P systems is such that the only possible maximally parallel step that can be performed consists in applying rule r1 twice (thus consuming reactants a2 b2 ) and rule r2 twice (thus consuming the remaining reactants a2 ). The lower priority rule r2 can be applied even if r1 is applied since the latter cannot consume all of the a’s. On the other hand rule r2 cannot be applied more than twice since this would require the higher priority rule r1 to be applied one time less. 3.2. Semantics In order to define the semantics, some auxiliary definitions are needed. We define the notion of p-enabled rule, which allows us to determine whether a given rule can be considered for application in a given absolute position p with respect to the reference coordinate system. In particular, it takes into account the relative positions of the objects, and ensures that their corresponding target positions are contained in the correct membrane regions. Definition 3.4. Given a membrane m ∈ {1, . . . , n}, and a position p = (x, y ) ∈ Region(m), an evolution rule r = (u 1 ) p 1 . . . (uk ) pk J(π1 )q1 . . . (πγ )qγ K → ( v 1 )t1 . . . ( v h )th ∈ R m is p-enabled iff the following conditions are satisfied: 1. ∀i ∈ {1, . . . , k}. p + p i ∈ Region(m) 2. ∀ j ∈ {1, . . . , γ }.⎧p + q j ∈ Region(m) ⎨ p + t i ∈ Region(m) if t i ∈ Z2 3. ∀i ∈ {1, . . . , h}. p + d ∈ Region(parent(m)) if t i = outd , with d ∈ Z2 ⎩ p + d ∈ Region(m ) if t i = inm ,d , with (m, m ) ∈ μ and d ∈ Z2 where parent( j ) = i, if (i , j ) ∈ Enabled( R , p ).

μ; and parent(1) = 0. The set of all p-enabled rules of a set of rules R is denoted by

Condition 1 ensures that each relative position p i ∈ Z2 of a reactant refers to a position p + p i which is contained in the membrane region, i.e. p + p i is inside membrane bounds and does not belong to the region of any inner membrane. Condition 2 is the analogous of the previous one as regards promoters. Condition 3 deals with products. The first case, analogously to conditions 1 and 2, ensures that for any target position t i ∈ Z2 appearing in the rule, the resulting position p + t i with respect to the current position p is contained in the

20

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

membrane region. If the rule specifies an outd target, then the second case ensures that the resulting position p + d is in the region of the parent membrane. Note that, if the current membrane is the skin membrane, then the resulting position must be outside it. Finally, if the rule specifies a target inm ,d , then m must be a child membrane of m and the resulting position must belong to the region of m . In each step of the evolution of a Spatial P system, some evolution rules are chosen and applied to the system state by removing all the reactants and adding all the products. In particular, in each step, for every membrane m and for every position p in membrane region, a multiset of p-enabled evolution rules is chosen non-deterministically, with some constraints needed to ensure that mutual exclusivity among objects, priority, and maximality are handled correctly. For each position of a Spatial P system model, the multisets of rules selected for application are described by a selection, defined as follows. Definition 3.5 (Selection). A selection of rules is described by a function S from positions to multisets of rules, such that for each position p ∈ Extents(1), S ( p ) = R is a multiset of p-enabled rules from R m , where membrane m is such that p ∈ Region(m). We assume that ∀ p ∈ Extents(1). ∃R. S ( p ) = R. Finally, an empty selection S is a selection such that ∀ p ∈ Extents(1). S ( p ) = ∅. We introduce some auxiliary functions which deal with the reactants, promoters and products of the rules described by a selection. Given a rule r and a relative position p, let η(r , p ), Π(r , p ), and T (r , p ) represent the multisets of reactants, promoters, and products, respectively, appearing in the rule in the relative position p. Formally:



  η (u 1 ) p1 . . . (uk ) pk Jπ K → v , p = u i if ∃i . p = p i ∅





otherwise

πi if ∃i . p = p i Π u J(π1 ) p 1 . . . (πγ ) p γ K → v , p = ∅

otherwise if ∃i . t i = p ∨ t i = out p ∨ t i = in j , p T u Jπ K → ( v 1 )t1 . . . ( v h )th , p = ∅ otherwise.



These functions are extended to multisets of rules as follows: η(R, p ) = r ∈R η(r , p ), Π(R, p ) = r ∈R Π(r , p ), and

T (R, p ) = r ∈R T (r , p ). Finally, given a selection S and an absolute position q ∈ Extents(1), we denote:







vi

1. by η( S , q) the multiset of objects which are required as reactants in q; 2. by Π( S , q) the objects required as promoters in q; 3. by T ( S , q) the objects produced in q by the application of the rules. Note that in the definition of functions η , Π and T to be applied to selections we make a little abuse of notation. In fact the functions applied to selections are not an extension of the functions applied to multisets of rules. The difference is that in the case of multisets of rules, η , Π and T require as second argument a relative position, whereas in the case of selections η, Π and T for selections is the following: the second

argument is an absolute position.

The formal definition of functions

η( S , q) = ( p,R)∈ S η(R, q − p ), Π( S , q) = ( p,R)∈ S Π(R, q − p ), and T ( S , q) = ( p,R)∈ S T (R, q − p ). A selection represents a potential step of a Spatial P system, since it specifies for each position a multiset of rules that could be applied. A selection is valid if it can correspond to an actual step of the system, namely if all of its rules can be applied in the positions they are associated with, and the application of such rules agrees with maximal parallelism, priorities and the constraints on the mutual exclusive objects. We now define notions of applicability and mutual exclusivity of a selection that will be used later to formally define the notion of valid selection. Definition 3.6 (Applicability). Let W be a configuration; a selection of rules S satisfies the property of applicability iff the reactants and promoters required by each rule are present in the system, i.e.:

∀ p ∈ Extents(1). η( S , p )F w p ∧ Π( S , p )F w p . Definition 3.7 (Mutual exclusivity). Let W be a configuration; a selection of rules S satisfies the property of mutual exclusivity of objects iff each position p in the configuration resulting from the application of the S to W contains at most one mutually exclusive object, i.e.:

∀ p ∈ Extents(1).





m( w p , x) − m







η( S , p ), x + m T ( S , p ), x  1.

x∈ E

(ζ )

Let us denote by S W the set of all selections satisfying a set of properties ζ ⊆ {appl, excl}, for a given configuration W , where appl and excl refer to the properties of applicability and mutual exclusivity, respectively. In order to formalise the notion of valid selection we need to introduce the following relations on multisets of rules, and on selections.

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

21

Definition 3.8. The priority relation  on multisets of rules, and its extension to selections, are defined as follows:

R  R SS

       ∀r ∈ R . m R , r < m(R, r ) ⇒ ∃r > r . m R , r > m R, r ;

⇐⇒



⇐⇒



∀ p . S ( p )  S ( p ).

(1) (2)

Given two multisets of rules R, R , they are in the priority relation R  R iff whenever R contains a rule r whose multiplicity is smaller than the multiplicity it has in R, there exists some other rule r in both R and R such that r > r and its multiplicity in R is greater than the multiplicity it has in R. In other words, R  R means that either R and R are the same multiset, or the highest priority rules (there can be more than one of these since priority is a partial order) the multiplicity of which is not the same in R and R must have a higher multiplicity in R than in R. The priority relation is extended to selections by considering the multisets of rules selected for each position. Now we are ready to formally define the notion of valid selection, namely a selection such that there does not exist a higher priority selection which satisfies the properties of both applicability and mutual exclusivity. (appl,excl)

Definition 3.9 (Valid selection). Let W be a configuration; a selection S ∈ S W of priority and maximality, which is formally defined as follows:

is valid for W iff it satisfies the property

∀S  S. S ∈ / S (appl,excl) . (ζ )

valid valid The set of all valid selections is denoted by S W . In the notations S W and S W , the indication of the configuration W can be omitted if there is no risk of ambiguity. Now we give a formal definition of the semantics of Spatial P systems as a Labelled Transition System. In the definition of the semantics (and in the rest of the paper) we use an auxiliary function apply that from a configuration W and a selection S gives a new configuration W that is the result of application to W of the rules in S. Function apply is formally defined as follows:

apply( W , S ) = W

where ∀ p ∈ Extents(1). w p = w p \\η( S , p )  T ( S , p ).

Definition 3.10 (Semantics). The semantics of a Spatial P system P = ( V , E , μ, σ , W 1 , ( R 1 , ρ1 ), . . . , ( R n , ρn )) is defined as the Labelled Transition System (W , S , →) where W denotes the set of all possible configurations, S denotes the set of S all possible selections, and the transition relation →: W × S × W is such that W −→ W iff (i) W , W are configurations valid having domain Extents(1), (ii) S ∈ S W , and (iii) W = apply( W , S ). A computation of a Spatial P system from the initial configuration W 1 is a (possibly empty) sequence of transitions S1

S2

S k−1

valid W 1 −−→ W 2 −−→ · · · −−−→ W k , for some valid selections S i ∈ S W , i = 1, . . . , k − 1. A configuration W k is said to be i reachable if there exists a computation having it as its last configuration. A successful computation is a computation such that no transitions are possible in the last configuration. Since we assume that the initial configuration W 1 does not contain conflicts among mutually exclusive objects, also any reachable W does not contain any conflict. This is formalised by the following lemma and corollary.

Lemma 3.1. Given a configuration W of a Spatial P systems P , it holds

   S valid ∀ S ∈ SW . ∀W . W − → W ⇒ conflict-free W where the labelled transition relation is obtained in accordance with the semantics of P . Proof. Selection S is valid. By Definition 3.9 a valid selection satisfies the mutual exclusivity property (Definition 3.7). Such a property implies that the result of the application of the rules of the selection on any configuration is conflict-free. 2 Corollary. Given a Spatial P system P , all its reachable configurations are conflict-free. Proof. Follows from Lemma 3.1 and from the fact that the initial configuration of a Spatial P system is always assumed to be conflict-free (Definition 3.3). 2 4. Simulation of Spatial P systems We develop a polynomial time algorithm for the simulation of Spatial P systems containing only mutually-exclusive objects and singular rules (see Section 3.1), and in which the priority relation among the rules with the same reactant, w.r.t. each membrane, defines a total ordering. By considering such a limited kind of models, we are able to highlight the

22

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

functioning of our algorithm in the basic case, and prove its correctness, before discussing the algorithms for some more complex cases as extensions of the basic one. Let us consider singular rules (see Section 3.1), that is, rules of the form (ax )(0,0) Jπ K → (a y ) p for some promoters π , (a)

position p, and ax , a y ∈ E. Given a membrane i, we denote by R i ⊆ R i the subset of rules in which ax appears as the reactant of the rule. We present a simulation algorithm which efficiently simulates Spatial P system models constrained as follows:

• there are n mutually-exclusive objects: E = {a1 , . . . , an }, and no ordinary objects (V = ∅); • only singular rules are allowed; • the priority relation on rules in R (i ax ) defines a total ordering, for all membranes i. 4.1. The idea behind the algorithm Recall that a selection describes which rules are to be applied to all the positions contained in the extents of the skin membrane. In other words, it describes, all at once, the rule applications for all the positions inside all membranes. The algorithm is based on the idea that, given a configuration of a Spatial P system, a valid selection can be constructed incrementally, by starting from an empty selection and by either adding new rules one-by-one or replacing existing rules with higher-priority ones. Then, at each iteration of the incremental construction, it is sufficient to perform the following steps sequentially: 1. considering the configuration, randomly pick an object contained in some position (in any membrane); 2. either try to include in the selection a new rule for it (if not previously done) or replace the already included rule with another with higher priority (this is made easy by the fact that priority is assumed to be a total order); this step is performed as long as the resulting selection is conflict-free. The two steps above are iterated in a cycle until a fixpoint selection has been obtained, namely a selection in which no rule, in any position, can be further added or replaced without creating conflict. A selection like this will be called filling selection. This process of incremental construction results in a selection in which rules are applied in a maximal way since no further rule can be added to it. The filling selection computed so far, however, might be not valid. In particular, there could exist rules with higher priority that could be successfully applied to the considered configuration. Let us take, for example, the situation in which the considered configuration includes two different mutually-exclusive objects in two adjacent positions. Moreover, let r1 and r2 be rules of the Spatial P system that, when applied together, swap such objects. If these are the only available objects and rules, the execution of the algorithm terminates immediately by giving an empty selection as result. This is because the selection computed at the first step of the algorithm (which includes either r 1 or r2 ) is not conflict-free. The algorithm described so far terminates with a selection such that no higher-priority selection can be obtained by replacing only one rule with a higher-priority one. However, a higher-priority selection could be obtained by replacing all at once a number of rules with the same number of higher-priority rules. Computing a valid selection out of a filling one requires a second algorithm to be given that identifies suitable groups of rules in the filling selection and replaces them with higher-priority ones. This can be done by constructing a sort of dependency graph and by looking for cycles in it. 4.2. Auxiliary definitions Before presenting the algorithm, we need to introduce some concepts. Let us denote by S [R/ p ] the selection S such that ∀q = p . S (q) = S (q) and S ( p ) = R. The first definition we give is that of improved multisets of rules. (appl)

Definition 4.1. Given a configuration W and a selection S ∈ S W for the position p is defined as follows:

R ∈ improved( W , S , p )

⇐⇒

, the set of improved multisets of rules improved( W , S , p )

  (appl) R = r ∧ ∀r ∈ S ( p ). r < r ∧ S R / p ∈ S W .

Let R = S ( p ) be the currently selected multiset of rules for a position p, that is either R = {} or R = {r } for some r. Informally, the set of improved multisets of rules improved( W , S , p ) corresponds to all the possible singletons R = {r } such that: (i) the priority of r is higher than the one of rule r selected in R, if any; and (ii) the selection resulting from the replacement of R with R in position p still satisfies the constraint of applicability (Definition 3.6). The definition of improved multisets is the basis for the definition of a filling selection, which will be used in the definition of the simulation algorithm. Let target be a function that gives the resulting (relative) position associated with a rule, that is, target(a(0,0) Jπ K → a p ) = p . Such a function is also extended to singleton sets as target({r }) = target(r ). Moreover, in case no rule is applied to an object, since the object continues to occupy its current position, we can assume

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

23

that the resulting target position coincides with the original position. In other words, this corresponds to considering the application of a rule of the form a(0,0) → a(0,0) . For this reason we assume target(∅) = (0, 0). (appl,excl)

Definition 4.2. Given a configuration W , a selection S ∈ S W



   ∀ p . ∀R ∈ improved( W , S , p ). target S ( p ) = target R filling

We denote by S W

is said to be filling for W iff the following property holds:

⇒





/ SW S R / p ∈

(excl)

.

the set of all filling selections for a given configuration W .

A filling selection is such that the resulting position of each rule r having higher priority than the currently selected rule r, if any, is already occupied. If no rule is selected for the position p, then the resulting position of each p-enabled rule is occupied. This implies that, given a rule r applied to an object a in a position p ∈ Region(i ) for some i, it is not possible (a) to replace the application of r with any other p-enabled rule r ∈ R i such that r > r and target(r ) = target(r ), since the resulting position is already occupied. Note that the definition of filling selection considers, for possible replacements, each position singularly. Moreover, the condition target( S ( p )) = target(R ) is necessary to prevent the replacement of a rule r with another rule r having the same target position. As it will be discussed later, such a constraint is needed to ensure the correctness of the algorithm (in particular, to ensure that the algorithm could possibly generate any valid selection). 4.3. The simulation algorithm The simulation algorithm for the restricted case performs a sequence of steps, where each step corresponds to a transition of the Spatial P system and is simulated by the algorithm SimulateStep. Given a conflict-free configuration W , SimulateStep( W ) non-deterministically generates a configuration W simulating a transition of the Spatial P system, according to the semantics given in Definition 3.10. A step is composed of the following two phases:

• (phase 1) determine a filling selection S F ∈ S filling which can possibly be not valid. • (phase 2) find a valid selection S  ∈ S valid such that S F ≺ S  . This is done by computing a sequence of selections S F = S 1 ≺ S 2 ≺ · · · ≺ S m = S  in which each S i ∈ S (appl,excl) . Note that, if selection S F is valid, then S  = S F . Finally, the computed valid selection S  is applied to the current configuration W giving W = apply( W , S  ) as result. Formal definition. The algorithm SimulateStep is formally defined in Algorithm 1, which uses the algorithms GenerateSelection and ValidateSelection implementing phase 1 and phase 2 of the algorithm, respectively. Algorithm 1 SimulateStep( W ) 1: 2: 3: 4:

let S ⊥ be such that ∀ p ∈ Extents(1). S ⊥ ( p ) = ∅ S F = GenerateSelection( W , S ⊥ ) S  = ValidateSelection( W , S F ) return apply( W , S  )

Algorithm 2 GenerateSelection( W , S ) 1: while ∃ p . ∃R ∈ improved( W , S , p ).

target( S ( p )) = target(R ) ∧ S [R / p ] ∈ S W 2: S = S [R / p ] 3: end while 4: return S

(appl,excl)

do

Algorithm GenerateSelection, defined in Algorithm 2, starts from the empty selection S ⊥ ∈ S (appl,excl) and then iteratively constructs greater selections by performing, in each iteration, a replacement of a rule in a position. Precisely, in each iteration, the current multiset S ( p ) selected in a position p is replaced with an improved multiset R ∈ improved( W , S , p ) which causes a change in the set of positions occupied by the mutually exclusive objects, and such that we have the updated selection S = S [R / p ] ∈ S (appl,excl) . That is, selection S has to satisfy the constraints on applicability (Definition 3.6) and mutual exclusivity (Definition 3.7). Finally, the algorithm stops generating a filling selection S F ∈ S filling . Formally, the requirement that the updated selection S be such that the set of occupied positions changes is expressed by the condition target( S ( p )) = target(R ). Such a condition is actually needed to prove the correctness of the algorithm (see also Definition 4.2). Algorithm ValidateSelection, defined in Algorithm 3, implements phase 2 of the algorithm. It computes a sequence of selections S F = S 1 ≺ S 2 ≺ · · · ≺ S m = S  where each selection S i is determined by finding a set of occupied positions Q in

24

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

Algorithm 3 ValidateSelection( W , S ) 1: repeat W ,S 2: let G dep = (V , E ) be the dependency graph defined by Eqs. (3) and (4)

3: if there exists a cycle C = p 1 , p 2 , . . . , pk  in the graph, with p 1 = pk then 4: let S = S 5: for i = 1 to k − 1 do 6: let (R , q) ∈ replacements( p i , p i +1 ) 7: S = S [R /q] 8: end for 9: S=S 10: end if 11: until there are no more cycles 12: return S

the preceding selection S i −1 where it is possible to replace all together the applied rules with higher-priority ones without violating mutual exclusivity, namely S i ∈ S (appl,excl) for all i. W ,S The implementation of ValidateSelection uses, at each iteration, a directed graph G dep = (V , E ), where E ⊆ V × V , describW ,S

ing the dependencies between applications of rules. By construction, graph G dep is such that each directed cycle describes a possible replacement of rules with respect to a set of occupied positions Q . The ValidateSelection algorithm iteratively W ,S constructs the graph G dep , searches for a cycle C , and performs the replacement of rules described by C by updating the current selection S. This procedure is iterated until there are no more cycles, terminating with a valid selection. The conW ,S struction of the dependency graph G dep , with respect to a configuration W and a selection S, is formally detailed in the following. Let S be a selection of a given a configuration W . The set of nodes V ⊆ Extents(1) of the graph, corresponding to non-empty positions in the configuration W = apply( W , S ), is defined as:

      V = q  W = apply( W , S ) ∧ m w q , x = 1

(3)

x∈E

where function apply is defined as in Definition 3.10. As regards the edges, consider each non-empty position q in W = apply( W , S ). Then, the set of edges E of the dependency graph is the minimal set such that:

R ∈ improved( W , S , q)

⇒









q + target S (q) , q + target R



∈ E.

(4)

In this way, there is an edge in the graph iff the current multiset S (q), selected in some position q, can be replaced by a multiset R ∈ improved( W , S , q), in which either (a) an application of a rule r is replaced by the application of another rule r with higher priority; or (b) no rule is currently selected in S (q) and one becomes applied. In both cases, the replacement results in the position p = q + target( S (q)) being freed and p = q + target(R ) being occupied, thus the edge ( p , p ) is present in the graph. Note that p and p may actually be the same position. Note that different replacements can be represented by the same edge ( p , p ) ∈ E and, consequently, there can be different replacements which correspond to the same cycle C . We denote by replacements( p , p ) the set of replacements corresponding to the edge ( p , p ) ∈ E , represented as pairs of the form (R , q). The case of different replacements which correspond to the same edge occurs when there are different rule applications which would put an object in the same position, among the p-enabled rules of priority higher than that of the current rule. Since the algorithm does not take into account the actual object present in a position, all those rules can be considered equivalent for the purpose of the replacement. In fact, selecting any of those rules would produce a replacement which can be applied. Moreover, there is no need of selecting the rule with the highest priority, since the highest priority rule could be found and applied in some subsequent iteration. filling Once phase 2 is completed, the algorithm has obtained a maximal filling selection S  ∈ max(S W ), which is also a valid selection according to Definition 3.9. That is, it is not possible to apply, in any position, any other higher priority mutually exclusive rule without violating either the constraint on applicability or the constraint on mutual exclusivity from Definitions 3.6 and 3.7, respectively. Example 4.1. This example illustrates how the algorithm works in the case in which a valid selection is obtained at the end of phase 1, by the algorithm GenerateSelection. Consider a Spatial P system composed of a single rectangular skin membrane having size 3 × 2, defined as Π = ({ }, {a, b}, { }, σ , W 1 , ( R 1 , ρ1 )) where σ = {1 → (0, 0), (0, 1), (0, 2), (1, 2), (2, 2), (3, 2), (3, 1), (3, 0), (2, 0), (1, 0)}, rules R 1 = {r1 : a(0,0) → a(2,0) , r2 : a(0,0) → a(1,0) , r3 : b(0,0) → b(0,1) } with priorities ρ1 such that r1 > r2 . As regards the initial state, we consider the state W 1 depicted in Fig. 3. From such a state, there are three possible transitions that the Spatial P system can perform, resulting in one of the possible configurations W 2 , W 3 , W 4 , as shown in Fig. 3. Let us consider the initial empty selection S 1 = S ⊥ ∈ S (appl) . According to the definition of improved, the set of improved multisets of rules improved( W 1 , S 1 , q) for each non-empty position q ∈ {(0, 1), (1, 0), (2, 0)} are as follows: (0, 1) → {{r1 }, {r2 }}, (1, 0) → {{r3 }} and (2, 0) → {{r3 }}.

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

25

Fig. 3. The initial configuration W 1 and the configurations reachable after one step, for Example 4.1.

Fig. 4. Tentative configurations obtained during algorithm execution (Example 4.1).

Fig. 5. The dependency graph of Example 4.1.

Suppose that the algorithm selects p = (0, 1) with R = {r2 }. Note that, since no rule is currently applied in p, target( S 1 ( p )) = target({ }) = (0, 0). Moreover, since target(R ) = (1, 0) = (0, 0), and the resulting selection satisfies the constraint on mutual exclusivity, the current selection S is updated to S 2 = S 1 [{r2 }/ p ] ∈ S (appl,excl) . The current tentative configuration U 2 = apply( W 1 , S 2 ) is shown in Fig. 4. At the subsequent iteration improved( W 1 , S 2 , (0, 1)) = {{r1 }}, while the improved sets of positions (1, 0) and (2, 0) are unchanged. Assume that the algorithm selects again p = (0, 1), with R = {r1 }. Since the constraints are satisfied also in this case, the new selection obtained is S 3 = S 2 [{r1 }/ p ], corresponding to the configuration U 3 = apply( W 1 , S 3 ) (Fig. 4). The improved set for position (0, 1) now becomes empty, since there is no other rule with higher priority than r2 that can be applied in such a position. Assuming that position p = (1, 0) is selected, then rule r3 is applied, obtaining S 4 = S 3 [{r3 }/ p ], corresponding to the configuration U 4 . The last position considered is p = (2, 0), with R = {r3 }. However, in this case such a rule cannot be applied, since S 4 [{r3 }/ p ] ∈ / S (excl) . In fact, it would cause a conflict between the mutually exclusive objects a, b in position (2, 1), as depicted by configuration  U 5 in Fig. 4. Since there are no more improved sets to be considered, GenerateSelection terminates with selection S F = S 4 , corresponding to configuration U 4 . W ,S As regards phase 2, algorithm ValidateSelection considers the dependency graph G dep1 F = (V , E ), where V = {(1, 1), (2, 0), (2, 1)} are the non-empty positions of configuration U 4 . As regards the edges E , the improved sets are empty for each position except q = (2, 0), for which improved( W , S F , q) = {{r3 }}. Therefore, the only edge ( p 1 , p 2 ) ∈ E is from position p 1 = q + target( S F (q)) = (2, 0) to position p 2 = q + target({r3 }) = (2, 1). Note that such an edge correctly describes that, if rule r3 is applied the object b in position (2, 0), then position p 1 would be freed and p 2 occupied. The dependency graph is shown in Fig. 5. Since there are no cycles, algorithm ValidateSelection terminates with the valid selection S  = S F . Finally, it is easy to see that algorithm GenerateSelection, by selecting different positions and rules at each iteration, can generate each of the possible resulting configuration W 2 , W 3 , and W 4 . Example 4.2. This example shows a case in which algorithm GenerateSelection produces a selection which is not valid, as it violates the property of priority and maximality from Definition 3.9. Consider a Spatial P system with objects V = ∅, E = {a, b, c }, having only one membrane, and with initial configuration W 1 as shown in Fig. 6. The evolution rules and their relative priorities are shown in the following: (2 )

ra

(1 )

: a(0,0) → a(1,1) > ra : a(0,0) → a(2,1)

(2 )

(1 )

(2 )

: c (0,0) → c (−2,1) > rc : c (0,0) → c (−1,1) .

rb : b(0,0) → b(0,1) > rb : b(0,0) → b(−1,1) rc

(1 )

Let p 1 = (0, 1), p 2 = (1, 1), p 3 = (2, 1) denote the positions above objects a, b, and c in W 1 , respectively. We have that a (1) (2) (1) (2) can be moved either to p 3 by ra or in p 2 by ra . Moreover, b can be moved either to p 1 by rb or to p 2 by rb . Finally, c (1)

(2)

can be moved either to p 2 by rc or to p 1 by rc . (2) (1) (2) (1) At the beginning, the improved sets for each non-empty position in W 1 are: (0, 0) → {ra , ra }, (1, 0) → {rb , rb }, and (2)

(1)

(1)

(2, 0) → {rc , rc }. Suppose that position (0, 0) is considered at the first iteration, and that ra (1)

obtaining the tentative configuration V 2 depicted in Fig. 6. Let us suppose that rb (1)

configuration V 3 . Finally, suppose that rc

is selected and applied, is applied, giving the conflict-free

is applied in position (2, 0) giving the configuration V 4 .

26

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

Fig. 6. Configurations of Example 4.2.

Fig. 7. The dependency graph of Example 4.2.

(2)

(2)

(2)

With respect to configuration V 4 , the improved sets are (0, 0) → {ra }, (1, 0) → {rb }, and (2, 0) → {rc }, meaning that it could be possible to replace, in each position from Q = {(0, 0), (1, 0), (2, 0)}, the currently selected rule with another higher priority rule. Nevertheless, by considering each position in Q singularly, the replacement of a rule with a higher priority one would cause a conflict between mutually exclusive objects in each case. For example, the application of rule (2) ra would cause a conflict in position p 2 , which is already occupied by a c object, and similarly for the other cases. filling

Therefore, phase 1 terminates by producing the selection S F ∈ S W (1)

(1)

(1)

describing the application of the rules ra , rb , and

rc , and which corresponds to configuration V 4 . (1) (2) (1) (2) The selection of rules S F is not valid, since it is possible to replace the application rb with rb and rc with rc , obtaining a greater selection (according to the priority relation defined by Eq. (6)) that still satisfies the properties of both W ,S applicability and mutual exclusivity. The dependency graph G dep1 F = (V , E ), used by algorithm ValidateSelection, is such (1)

(1)

that V = { p 1 , p 2 , p 3 }. As regards the edges E , note that the improved sets are such that (0, 0) → {ra }, (1, 0) → {rb }, and (1)

(1)

(2, 0) → {rc }. Therefore, considering position (0, 0), there is an edge ( p 3 , p 2 ) ∈ E corresponding to the replacement of ra (2) (1) (2) with ra . If we consider position (1, 0), there is an edge ( p 1 , p 2 ) ∈ E corresponding to the replacement of rb with rb . (1) (2) Finally, the replacement of rc with rc in position (2, 0) gives the edge ( p 2 , p 1 ) ∈ E . The dependency graph is shown in Fig. 7. Note that the only cycle C = p 1 , p 2 , p 1  in the dependency graph corresponds to the simultaneous replacements of rules in positions (1, 0) and (2, 0). Algorithm ValidateSelection then performs the replacement, giving the valid selection of (1)

(2)

(2)

rules ra , rb , rc , which corresponds to the final configuration V 5 shown in Fig. 6. 4.4. Correctness of the simulation algorithm The following theorem proves that simulation algorithm SimulateStep behaves correctly with respect to the semantics of Spatial P systems. Theorem 4.1. S ∀S. W − → W

⇐⇒

W can be generated by SimulateStep( W ).

Proof. Case ⇐) Phase 1 of the algorithm may produce a selection S F which is not valid because it violates property on priority and maximality (Definition 3.9). By definition of the priority relation on selections  (Eq. (6)), this means that there exists a valid selection S   S F , with S  ∈ S valid . In phase 2, such a valid selection S  is determined. Let S i be a selection satisfying both the properties of applicability (appl) and mutual-exclusivity (excl), but not valid, i.e. S i ∈ S (appl,excl) \ S valid . In order to find a selection S i +1  S i satisfying both appl and excl, it is necessary to find a set of positions in which some rule is applied in each position, and such that either (i) the multiplicity of a rule is increased (m(R, r ) < m(R , r )), or (ii) there is a rule with higher priority whose multiplicity is increased (∃r > r . m(R , r ) > m(R, r )). Since in each position there is at most one object, case (i) cannot occur. For the same reason, in case (ii), increasing the multiplicity of a rule r > r, with respect to the rule r currently applied in a position, actually means replacing the application of r with an application of r . Recall that phase 1 produces a filling selection, namely such that for any position p ∈ Region(i ), for some i, in which a rule r is applied to an object a, all the positions which would be occupied by the application of any other p-enabled rule (a) r ∈ R i , with r > r, are already occupied. This implies that, for each replacement of a rule that results in a position x being freed and a position y being occupied (which can be the same as x), there are both (i) a replacement which occupies x and (ii) a replacement which frees position y. Formally, for any position p in which a rule r ( p ) = a(0,0) Jπ K → a p 1 is replaced by a rule r ( p ) = a(0,0) Jπ K → a p 2 , there exists a position p in which a rule r ( p ) = b(0,0) Jπ K → bq 1 is replaced by some rule r ( p ) such that p + p 2 = p + q1 . This corresponds to finding a cycle in the dependency graph G dep constructed in algorithm ValidateSelection. By applying such a replacement, we obtain a selection S i +1  S i . Note that there exists a selection S i +1  S i iff there is such a replacement of rules. Since there are no infinite chains of selections S 1 ≺ S 2 ≺ · · · ≺ S i ≺ · · · in which

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

27

all the selections satisfy properties appl and excl but are not valid, by iterating this procedure we eventually reach a valid selection S  . Case ⇒) Let us denote by −−−→1 the relation describing an iteration of the algorithm GenerateSelection, implementing phase 1. That is, S 1 −−−→1 S 2 iff, starting from selection S 1 , S 2 is reachable after one iteration of algorithm GenerateSelection. Let ⇒1 be the transitive closure of −−−→1 . As regards phase 2, we denote by −−−→2 the relation describing an iteration of the ValidateSelection algorithm, and by ⇒2 the transitive closure of −−−→2 . Moreover, given a configuration W , and (appl,excl)

, let us denote by occ( S ) the set of occupied positions in the configuration resulting from the a selection S ∈ S W application of S to W , that is, given W = apply( W , S ), then occ( S ) = { p | w p = ∅}. Let min(S filling ) be the set of minimal elements of S filling . That is, min(S filling ) ⊆ S filling is such that: (i) ∀ S ∈ min(S filling ).  S ∈ S filling . S ≺ S; and (ii) ∀ S ∈ S filling . ∃ S ∈ min(S filling ). S  S . In order to prove the thesis, it is sufficient to prove that: 1. ∀ S ∈ min(S filling ). ⊥ ⇒1 S, that is, any minimal filling selection can be obtained from GenerateSelection; 2. ∀ S ∈ S valid . ∃ S ∈ min(S filling ). S ⇒2 S, that is, for any valid selection S there is a minimal filling selection S such that algorithm ValidateSelection can obtain S starting from S . As regards part 1, we prove that, given a minimal filling selection S, for all selections S ≺ S, GenerateSelection algorithm can perform an iteration from S to some other S which is still less than or equal to selection S. That is, given a minimal filling selection S ∈ min(S filling ):

∀ S ∈ S (appl,excl) . S ≺ S

⇒

∃ S ∈ S (appl,excl) . S −−→1 S ∧ S ≺ S  S .

This implies that GenerateSelection can generate any minimal filling selection S. Consider a selection S ∈ min(S filling ), and a selection S ∈ S (appl,excl) such that S ≺ S. Let us consider the sets of occupied position for selections S and S , namely occ( S ) and occ( S ), respectively. There are two possible cases: either occ( S ) = occ( S ), or not. On one hand, assume that occ( S ) = occ( S ). Then there must exist a position p in which the currently applied rule r, i.e. S ( p ) = {r }, can be replaced by the algorithm with a rule r  r such that target(r ) ∈ / occ( S ), and which is smaller than the rule r selected for p in S, i.e. S ( p ) = {r } ⇒ r  r . In fact, suppose that such a position does not exist. Since S ∈ S filling , also S ∈ S filling and, therefore, S would not be minimal, contradicting the hypothesis. On the other hand, assume that occ( S ) = occ( S ). Note that the number of occupied positions in S and S is the same, namely |occ( S )| = |occ( S )|. Therefore, there exists a position q ∈ occ( S ) \ occ( S ). Consider the selection S, and let p be the position in which there is a rule r applied in p that puts an object in q, that is, S ( p ) = {r } ∧ target(r ) = q. It is sufficient that the algorithm selects position p to search for a greater rule r to replace the currently applied rule r, with S ( p ) = {r }. In fact, assume that the algorithm always selects the lowest applicable rule r such that r  r  r and target(r ) = target(r ). In the worst case, if no rules rˆ such that r ≺ rˆ ≺ r can replace the currently applied rule r, then the algorithm selects rule r . In such a case, the replacement could be forbidden only if target(r ) = target(r ). However, this is absurd, since it would imply that S is not minimal, contradicting the hypothesis. As regards part 2, consider a valid selection S ∈ S valid . In order to obtain a minimal filling selection, it is sufficient to construct a sequence of filling selections S = S 1  S 2  · · ·  S k = S , where each S i +1 is obtained from S i by replacing the rules selected for some positions with lower priority rules, ensuring that the set of positions occupied remains the same, i.e. occ( S i +1 ) = occ( S ). This construction corresponds to the converse of algorithm ValidateSelection, and terminates as soon as a minimal filling selection S = S k ∈ min(S filling ) is obtained. Therefore S ⇒2 S. 2 4.5. Algorithm complexity In order to make simulation of Spatial P system models feasible it is important to ensure that the considered simulation algorithm has a reasonable complexity. In particular, the simulation of a single evolution step of the model should be rather efficient since in many cases it will be repeated thousands or millions of times in every single simulation run. The complexity of the structure and of the semantics of Spatial P systems requires several computations to be performed at each simulation step. However, in the restricted case of Spatial P systems with singular rules and priority as a total order we show that, for a given Spatial P system model, a simulation step has polynomial complexity with respect to (i) the number of cells in the extents of the skin membrane, (ii) the number of objects, (iii) the number of rules, and (iv) the number of promoters in rules. We separately discuss the cost of algorithms GenerateSelection and ValidateSelection. GenerateSelection. We briefly discuss a simple concrete implementation of algorithm GenerateSelection (Algorithm 2), which starts from the empty selection S ⊥ . The idea of the algorithm is to maintain, for each non-empty position, the list of rules which can be applied in such a position, ordered with respect to the their priorities. Each step consists in searching for a position in which the currently applied rule (if any) can be replaced with a higher priority one, which can be done by examining the list of rules for such a position.

28

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

The algorithm works as follows: 1. associate, with each non-empty position p in some membrane i, the list of rules l( p ) = [r1 , . . . , rk ] ordered by priority, ( W ( p )) where {r1 , . . . , rk } = R i and ∀1  i < k. r i < r i +1 ; 2. set Q to be the set of initially occupied positions, i.e. Q = { p | W ( p ) = ∅}; 3. find a position p ∈ Q to “improve”, i.e. in which the currently applied rule in l( p ) can be replaced with a higher priority one, and update Q accordingly; 4. iterate the last step until no more replacements can be performed. Let N = |Extents(1)| be the total number of positions of the Spatial P system. Let K = max{| R i | | i ∈ {1, . . . , n}} be the (a) maximum number of rules of a membrane. Let  K = max{| R i | s.t. i ∈ {1, . . . , n}, a ∈ E } be the maximum number of rules applicable to any object, in any membrane. Note that  K is an upper bound of the length of each list l( p ) for any position p. Moreover, let m = |{ p | W ( p ) = ∅}| denote the number of positions occupied by some object. As regards the initialisation performed in step 1, let H denote the maximum number, among all rules, of different positions in which the promoters of a rule are contained, i.e. H = max{γ | ∃i . ∃(ax )(0,0) J(π1 )q1 . . . (πγ )qγ K → (a y ) p ∈ R i }. Initialising the lists costs O ( N K H ), since it demands checking the applicability of each membrane rule in each position, by checking whether all the promoters are present in their relative positions. In step 2, set Q can be initialised in O ( N ). As regards step 3, note that the number of rules not yet considered (i.e. the sum of the lengths of the lists of each position) is initially not greater than m K , since each object has at most  K possible moves. Initially, the worst case corresponds to visiting all the lists for each occupied position, and finding a replacement only for the last rule considered. Since, once a rule is applied, it is never considered again, the algorithm performs at most m K steps, where each step means visiting all the rules not yet considered. Since the number of rules not yet considered decreases by (at least) one in every step, and updating Q costs O (1), the total cost of step 3 is O ( N K H + m2  K 2 ).

m K

i =0 ( O (i )

+ O (1)) = O (m2  K 2 ). The total cost of the algorithm is thus

ValidateSelection. The dependency graph, constructed at the beginning of the algorithm, has at most κ = m K edges, since the are m objects and at most  K possible rule applications for each object. Thus, constructing the graph has complexity O (m K ). The complexity of searching a cycle in the graph is equal to the number of edges in the graph. After each iteration, some edges are removed, and thus the complexity of finding other cycles decreases. In the worst case, a cycle composed of only two nodes is found in each iteration, therefore the algorithm performs at most κ /2 iterations. The dependency graph, instead of being recreated at each iteration, can be incrementally updated at each iteration to reflect the application of a set of rules. Let p i denote the position occupied by the i-th rule in the list. Let ( p , p ) be an edge of the cycle, representing the application of a rule. The graph can be updated as follows: (i) remove all the edges corresponding to the rule itself and the other lower priority rules; (ii) replace each edge of the form ( p , z) (for some z) corresponding to some higher priority rule with ( p , z). The replacing of edge ( p , z) with ( p , z) means that higher priority rules can still be used for the object in p , which would cause it to move to another position z. The cost of the i-th iteration is O (i ), since the prevalent cost comes from updating the graph. In fact, the removal of an edge costs O (1), and each edge is removed at most once. Moreover, the cost of replacing edges is inversely proportional to the step, as the maximum number of edges exiting from a node decreases after each iteration. Therefore, the total cost of κ /2 ValidateSelection algorithm is O (i ) = O (κ 2 ) = O (m2  K 2 ). i =0 Given the cost of the algorithms composing the SimulateStep algorithm, as detailed in the previous paragraphs, and that the application of the selection of rules costs O ( N ), the total time complexity of the SimulateStep algorithm is O ( N K H + m2  K 2 ). 5. Extensions The simulation algorithm presented in the previous section handles only the quite restricted form of Spatial P system models specified in Section 4.3. This has allowed us to focus on the core aspects of Spatial P systems that make its simulation difficult, which stems from the presence of mutually-exclusive objects along with the maximally-parallel semantics for rule applications. Nevertheless, some of the restrictions introduced can easily be lifted, and dealt with by means of simple modifications of the proposed algorithm, as we discuss in this section. These modifications can all be used at the same time. We actually present some simple cases which are useful for the modelling of real-world systems, without claiming to be complete in regards to the possible extensions. Most of these extensions are then exploited in the modelling examples presented in Section 6.2. It is important however to keep in mind the core restriction that we keep: each rule must always contain exactly one mutually-exclusive object among the reactants, and one among the products. We speculate that, if we relaxed this constraint, then the problem of determining a step in the evolution of Spatial P systems would become NP complete. We informally discuss this aspect at the end of this section.

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

29

5.1. The role of ordinary objects Let us consider singular rules, namely rules composed of exactly one reactant and one product, both of which are mutually-exclusive objects. It is easy to see that the presence of ordinary objects among the promoters and/or products of rules does not substantially affect the complexity of the simulation algorithm. In fact, as regards promoters, their presence can be easily checked at the beginning of a step, while having some ordinary objects as products only affects the applicability of the rule as regards the resulting target positions, which must occur inside the membrane. This last aspect can be easily checked just once for each rule application considered in a step. On the other hand, the situation in case of ordinary objects among reactants is diverse. A simple case is that of ordinary objects among reactants referring to the position of application (i.e. to the relative position (0, 0)). In this case, we expect that the algorithm complexity is not influenced substantially. In fact, since due to the presence of a mutually-exclusive object at most one rule can be applied to a position, then whenever a rule is to be replaced with a higher priority rule it would be sufficient to also check whether the reactants are available at that precise time. Other cases, such as when there are also rules containing only ordinary objects among both reactants and products, need further investigations, that we leave as future work. 5.2. Priorities as a partial order In the following we discuss how the algorithm can be extended to deal with a priority relation on rules which is a partial order, thus removing the constraint of being a total order. First of all, instead of a single priority relation i for each membrane i, we discuss the possibility of having a different priority relation  p for each position p ∈ Region(i ) inside a membrane. In order to make algorithms GenerateSelection and ValidateSelection work in this case, it is necessary to introduce pointwise priority relations, and to redefine accordingly the function improved from Definition 4.1. Definition 5.1. Let p ∈ Region(i ) be a position, for some membrane i, and ρ p be a total order among the rules R i . The pointwise priority relation on multisets of rules and its extension to selections are defined as follows: R ext p R

S

ext

S



⇐⇒ ⇐⇒

      ∀r ∈ R . m R , r < m(R, r ) ⇒ ∃r >ext p r . m R , r > m R, r ; ∀ p. S ( p)

ext p



S ( p)

(5) (6)

 r. where r >ext p r iff (r , r ) ∈ ρ p and r = (appl)

Definition 5.2. Let W be a configuration and S ∈ S W be a selection. Consider a position p ∈ Region(i ), for some membrane i, and let ρ p be a total order among the rules R i . Then the set of improved multisets of rules improved( W , S , p ) for the position p is such that:

R ∈ improved( W , S , p )

⇐⇒

  (appl) R = r ∧ ∀r ∈ S ( p ). r
where r
The idea of this extension is to provide for an initial phase (called phase 0) which precedes phase 1, with respect to the high level description of the algorithm presented in Section 4.3, whose purpose is to compute, for each position p, a linear ordering  p of the p-enabled rules. A linear ordering, which is formally defined in the following, is a total order constructed from the partial order i defining priorities among rules R i , with p ∈ Region(i ). This enables the subsequent phase 1 and phase 2 to continue working as before, except for the fact that they need to use, for each position, a possibly different total order as the priority relation of the rules. Definition 5.3. Given a set of rules R with partial order , a linear ordering of the rules in R is a total order ext corresponding to a sequence r1 , r2 , . . . , rk  such that: 1. k = | R | ∧ ∀r ∈ R . ∃i . r = r i ; 2. ∀i , j . r i > r j ⇒ i < j. The first constraint ensures that the sequence r1 , . . . , rk  is a permutation of the elements from R, while the second constraint ensures that it respects the relative priority ordering between the rules. Note that a sequence r1 , . . . , rk  represents the inverse of a linear extension of the partially ordered set of rules R. Each time the algorithm SimulateStep is executed, it has to randomly generate, for each non-empty position p ∈ Region(i ), for some i, a linear ordering l = r1 , . . . , rk  of the p-enabled rules, namely of Enabled( R i , p ) = {r1 , . . . , rk }. Then it uses such a linear ordering as the pointwise priority relation ρ p for the position p.

30

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

Theorem 5.1. The extension of the SimulateStep algorithm, using a linear ordering of the rules enabled in a position as the ordering of the rules, behaves correctly with respect to the semantics of Spatial P systems. Proof. Any selection generated by the linear algorithm is valid with respect to the original partial order of the rules. In fact, given a selection of rules S generated by the algorithm and a position p, the property (b) of Definition 5.3 implies that any higher priority rule than the rule applied in p has been considered by the algorithm. In particular, the selection produced in phase 1 of the algorithm is also filling for the original partial order of the rules. As regards phase 2, if there are no replacements which can be performed, then this implies that there are no possible replacements also if one considers the original partial order of the rules. However, note that a replacement performed either in phase 1 or phase 2 may not correspond to obtaining a greater selection of rules with respect to the original partial order of rules. For example, consider the rules r1 , r2 , r3 , such that r1 > r3 , r2 > r3 , and the linear ordering r1 , r2 , r3 , for some position. Suppose the currently selected rule is r2 , and that the algorithm replaces r2 with r3 . Then, in this case, r3 is greater than r2 w.r.t. the linear ordering, but not w.r.t. the original partial order. On the other hand, assuming that every possible linear ordering of the rules can be generated in any position ensures that every possible valid selection can be generated by the algorithm. Indeed, it is sufficient that algorithm GenerateSelection considers a linear ordering of rules r1 , . . . , r x−1 , r x , r x+1 , . . . , rk , with r x ∈ S ( p ), such that the rules r x+1 , . . . , rk which follow r x in the linear ordering are the all and only rules having lower priority than r x , i.e. ∀i > x. r i < r x . 2 5.3. Extensions with richer evolution rules We have shown that it is possible to efficiently simulate Spatial P system models in which the form of evolution rules is quite restricted. As we are going to discuss in the following, we argue that this is not always the case, and in particular that the problem in the unrestricted case is NP complete. Therefore, since the actual complexity of simulation of a Spatial P system model lies in the form of the used evolution rules, this calls for a thorough investigations of the trade-off between simulation complexity and features. For example, as we have shown in the previous section, a simple extension of the proposed algorithm allows using a partial order for rule priority. We believe that the algorithm for the restricted case can be extended to other interesting cases, using different forms of evolution rules and other features, while still retaining a polynomial time complexity. Let us consider a Spatial P system model containing only mutually-exclusive objects, no restrictions on the form of rules, and priority as a total order. Note that the difference with respect to the restricted model considered for the development of the simulation algorithm in Section 4 lies only in the form of the rules, which are not required to be singular. Models of this kind can be simulated with an algorithm which iteratively performs two phases analogous to the two phases of algorithm SimulateStep for the restricted case (Algorithm 1), until a fixed point is reached. The two phases are the following: 1. improve the selection of rules in each position considered singularly, until no more replacements are possible; 2. search for groups of replacements in different positions which can be applied all at once. As regards phase 1, it can be implemented by using an algorithm analogous to GenerateSelection (Algorithm 2), which has polynomial time complexity. As regards phase 2, the underlying problem can be reduced to searching a hypercycle in a directed hypergraph, which is proved to be NP complete in [37]. The fact that we have to deal with hypergraphs instead of simple graphs is because a replacement of a rule with a higher priority one may free and occupy any number of positions, instead of having always one freed and one occupied as in the restricted case. Provided that the algorithm is correct, and that, given any hypergraph H it is possible to construct a Spatial P system such that the hypergraph constructed in phase 2 is exactly H , this suffices to prove that simulation problem in this case is NP complete. Thus, also the generic case, which includes that case, is also NP complete. 6. Applications In this section we discuss different applications of Spatial P systems, and present the results of their simulations. The first model describes the behaviour of territorial fish, such as Garibaldi (Hypsypops rubicundus), which protects a small area around the nest from any other fish coming near, by attacking any fish which is too close to the nest. Such a Spatial P system model is composed of territorial fish, each guarding a different area, and wandering fish, which move around randomly. The description is quite simple, since it uses only four kinds of rules for describing the behaviour of both fish. The second model captures the schooling behaviour of fish, such as herring, in which many single individuals aggregate to form schools, and move in a coordinated manner. This model is quite more complicated than the first one, since the movement of single herring is described at a lower level than the case of territorial fish. In particular, once the direction of movement of a herring is determined from the relative positions of nearby individuals, all the possible moves are described by a rule for each position in between its current position and the determined target position. Once a formal model is available, the modifications to the model can be implemented more easily than by modifying a concrete simulator; in our opinion, dealing with evolution rules (and their formal syntax) is easier than having to deal

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

31

Fig. 8. Description of relative positions included in sets Q i = {positions labelled by i }, for i ∈ {1, 2, 3, 4}, with respect to the origin (0, 0) in the central cell.

with concrete implementations in some programming language (such as, in our case, C#). This is the main motivation for developing a formal Spatial P systems model for the Boids model An important aspect of both models lies in the modelling of fish as mutually-exclusive objects. This allows us to analyse the behaviour of our models in case of high local density of fish, since the semantics of such objects prevent them from interpenetrating. This demonstrates an advantage of the Spatial P systems with respect to other modelling formalisms, in which such a boundary condition must be explicitly taken into account in the behavioural description. Moreover, having a formal model allows modifications to the models to be performed easily, since Spatial P systems provide a formal syntax which is easy to read and modify. In other words, dealing with high-level evolution rules is easier than having to deal with concrete implementations in some programming language, thus improving the ordinary cycle towards the development of a model, i.e., (i) making a hypothesis on the functioning of a systems, (ii) formalising the behaviour, (iii) performing simulations to test the validity of the hypothesis. 6.1. Territorial fish The model. We consider a number of territorial fish tf(1) , . . . , tf(n) , each protecting an area of radius 4 units around a their respective nests. The position of each nest is represented by an ordinary object s(1) , . . . , s(n) . Finally, wandering fish are denoted by the mutually-exclusive object w ∈ E. Therefore, the ordinary objects V and mutually-exclusive objects V of this model are:



E = tf(1) , . . . , tf(n) , w ;



V = s(1) , . . . , s(n) . The territory around a nest controlled by a territorial fish can be represented by a set of positions relative to the nest itself. Let us consider Fig. 8, in which central position denoted (0, 0) represents the centre of the coordinate system. We denote by Q 0 , Q 1 , Q 2 , Q 3 , Q 4 the sets of (relative) positions having distance 0, 1, 2, 3, 4, from the nest, respectively. For example, Q 1 = {(−1, 0), (0, 1), (1, 0), (0, −1)}, Q 2 = {(−2, 0), (−1, 1), (0, 2), (1, 1), (2, 0), (1, −1), (0, −2), (−1, −1)}, which correspond to the positions labelled 1 and 2 in Fig. 8, respectively. Note that Q 0 = {(0, 0)}. Formally, sets Q i can be defined as Q l = { p ∈ Z2 | l − 1 <  p 2  l}, for all l ∈ {0, 1, 2, 3, 4}. Moreover, we denote the entire set of positions of the territory covered by a territorial fish as Q ∗ = Q 0 ∪ Q 1 ∪ Q 2 ∪ Q 3 ∪ Q 4 (note that this also includes the origin). As regards movement of fish (both territorial and wandering ones), they are able to move only one cell at a time, in either one of the directions in D = {(−1, 0), (0, 1), (1, 0), (0, −1)}. That is, a fish in position p can move, in one step of the evolution of the Spatial P system, only to one of the positions in Z = { p + (−1, 0), p + (0, 1), p + (1, 0), p + (0, −1)}. Given a set X ⊆ Z2 and a position p ∈ Z2 , we define the translated set of positions in X with respect to p by X ( p ) = { p + q | q ∈ X }. Moreover, we assume function nearest( p , X ), giving the positions in X closest to p, using the Euclidean distance. Formally nearest( p , X ) = {q ∈ X | ∀q ∈ X . q − p 2  q − p 2 }. Using such notations, we define moves( p , q), where p is the position of a territorial fish moving towards q, which gives the positions, among those in Z and within territory Q ∗ , closest to target position q, i.e. moves( p , q) = nearest(q, D ( p ) ∩ Q ∗ ). Informally, the behaviour of fish is as follows.

• Territorial fish. When its territory is free from other fish, a territorial fish stays in the nest. We assume that the fish is able to see all its territory at once, therefore as soon as a wandering fish enters its territory, it moves towards it in order to chase it away. In case of multiple wandering fish, the territorial fish chases the one closest to the nest. When no other wandering fish is present, the territorial fish moves back to the nest. • Wandering fish. When there are no territorial fish nearby, a wandering fish either moves non-deterministically in one of the directions in D = {(−1, 0), (0, 1), (1, 0), (0, −1)}, or stays in its current position. Whenever a territorial fish attacks it, precisely if a territorial fish comes close to it within distance 2, the wandering fish runs away from the attacking fish by moving in the opposite direction from the territorial fish.

32

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

(i )

Evolution rules for territorial fish. The behaviour of a generic territorial fish i ∈ {1, . . . , n} is formalised by sets of rules R l , (i )

(i ) (i ) (i ) R l , for l ∈ {0, 1, 2, 3, 4}, and R ⊥ . Rules in both R l and R l , for l ∈ {0, 1, 2, 3, 4}, capture the behaviour in which the (i ) territorial fish chases a wandering fish in the territory, while R ⊥ captures the case in which no wandering fish is present.

Formally, for i ∈ {1, . . . , n}, they are defined as follows:



 p ∈ Q ∗ , q ∈ Q l , p ∈ moves( p , q) l ∈ { 0, 1 , 2 , 3 , 4 }

(i ) q (i ) y (i ) (i )  R l = tf p s(0,0) , wq → tf p p ∈ Q ∗ , q ∈ Q l l ∈ { 0, 1 , 2 , 3 , 4 }    q y (i ) (i ) (i ) (i ) R ⊥ = tf p s(0,0) → tf p  p ∈ Q ∗ , p = (0, 0), p ∈ moves p , (0, 0) (i )



y (i ) q (i ) (i ) s(0,0) , wq → tf p

R l = tf p

(i )

(i )

(i )

(i )

(i )

(i )

(i )

(i )

(i )

with priority R 1 > R 1 > R 2 > R 2 > R 3 > R 3 > R 4 > R 4 > R ⊥ . (i ) (i ) (i ) (i ) Set R l contains rules of the form tf p Js(0,0) , wq K → tf p . Such a rule is to be applied only to the position of the nest (i )

s(i ) , since s(0,0) occurs as a promoter. The rule reads as follows: if there is wandering fish w in position q ∈ Q l , then the territorial fish tf(i ) moves from position p to p ∈ moves( p , q). In other words, if there is a wandering fish at distance l from the nest, the territorial fish moves towards it. (i ) (i ) (i ) (i ) (i ) On the other hand, the lower priority rule set R l < R l contains rules of the form tf p Js(0,0) , wq K → tf p . Therefore, (i )

(i )

(i )

(i )

(i )

(i )

(i )

(i )

for each rule tf p Js(0,0) , wq K → tf p ∈ R l there exists a corresponding rule tf p Js(0,0) , wq K → tf p ∈ R l , having the same reactants and promoters as the former, and in which the territorial fish is not moved from its position. Since the former (i ) has lower priority than the latter, such a rule can be applied only if the first one cannot. Rule set R l is used to avoid (i )

unwanted behaviours in the case a rule from R l cannot be applied because of the constraint on mutually-exclusive objects, as discussed later. (i ) (i ) (i ) (i ) Rules from R ⊥ are a simplified form of the previous case: tf p Js(0,0) K → tf p , with p ∈ moves( p , (0, 0)). This rule can also be applied only to the position of the nest s(i ) , and describes that the territorial fish tf(i ) moves back towards the nest. (i ) (i ) (i ) (i ) (i ) (i ) (i ) (i ) (i ) Priority R 1 > R 1 > R 2 > R 2 > R 3 > R 3 > R 4 > R 4 > R ⊥ among the rules captures the intuitive behaviour of

territorial fish introduced previously. First of all, note that, since object tf(i ) occurs as a reactant in all the rules from (i ) (i ) (i ) R l , R l , R ⊥ , then at most one rule can be applied. Thus, according to the priority, if there is any wandering fish at (i )

distance 1 (i.e. a rule in R 1

is applicable), then the territorial fish chases it, irregardless of any other wandering fish (i )

(i )

(i )

(i )

at greater distance (namely if any rule from R j , with j ∈ {2, 3, 4}, is applicable). The case for rules in R 2 , R 3 , R 4 is analogous. (i ) (i ) (i ) (i ) However, when a wandering fish is present in the territory, it may happen that a rule tf p Js(0,0) , wq K → tf p in R j , for

some j, is not applicable because the target position p of the territorial fish is occupied, preventing it to move. According to the semantics of Spatial P systems, in such a case, any lower priority rule can be applied. In this case, the presence of the (i ) (i ) (i ) (i ) (i ) (i ) lower priority rule tf p Js(0,0) , wq K → tf p in R j prevents any other rule from any set R j with j > j, or R ⊥ to be applied, thus avoiding that the territorial fish chases another wandering fish or goes back to the nest.

w Evolution rules for wandering fish. The behaviour of wandering fish is formalised by two sets of rules R w att. , and R rel. , where the former describes the behaviour when attacked by a territorial fish, while the latter describes the relaxed situation, when there are no territorial fish in the vicinity.





 Rw att. = w(0,0) Jtf p K → w p p ∈ Q 1 ∪ Q 2 , p ∈ nearest(− p , D )







.

 Rw rel. = w(0,0) → wδ δ ∈ D ∪ (0, 0)

w The priorities among rules are R w att. > R rel. . We consider a wandering fish attacked when there is a territorial fish within distance 2 from it. In fact, rule w(0,0) Jtf p K → w p from R w att. , whose application is centred in the position of the wandering fish w, is applicable whenever there is a territorial fish tf in relative position p ∈ Q 1 ∪ Q 2 (i.e. at a distance no more than 2). In this case, the territorial fish moves in the opposite direction from the territorial fish, namely to position p among D ( p ) closest to − p (i.e. furthest from p); formally p ∈ nearest(− p , D ). On the other hand, when there are no territorial fish nearby, the wan, dering fish either moves randomly in any direction, or stays in the same position, by applying any rule w(0,0) → wδ in R w rel. where δ ∈ D ∪ {(0, 0)}. Thus, assuming the current position of w in p, it arbitrarily moves to any position in D ( p ), where D = D ∪ {(0, 0)}.

Example 6.1. We illustrate the actual application of the rules with a small example. Let us consider the situation depicted in Fig. 9a, in which there is a nest s(i ) , a territorial fish tf(i ) in the same position as the nest, and two wandering fish nearby. Note that Fig. 9a shows only a fragment of a bigger membrane, and we assume both that the other (not shown) positions are empty, and that membrane bounds are far away from the depicted positions. This last condition allows us to avoid the corner case in which rule application is prevented due to positions falling outside of the membrane.

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

33

Fig. 9. (a) Pictorial view of the initial configuration of a fragment of the Spatial P system described in Example 6.1, with two wandering fish w within the territory of tf(i ) . (b,c) The evolution rules applicable to object s(i ) , and to the leftmost wandering fish, with the relative priorities specified. (Rules are labelled univocally for convenience.)

Fig. 10. (a) One of the possible configurations reached after one step from the configuration depicted in Fig. 9(a). (b) The rules applicable to s(i ) . Rules from the same set are incomparable w.r.t. the priority relation; this intuitively means that they have the same priority. (See Example 6.1.)

Fig. 9b shows the sets of rules applicable to the position of the nest s(i ) . The rules are labelled univocally for convenience of the presentation. The leftmost wandering fish is at distance 3 from the nest, while the other fish is at distance 4, thus (i ) (i ) (i ) (i ) there is a rule applicable from each set R 3 , R 3 , R 4 , R 4 . Therefore, the territorial fish will try to chase the nearest wandering fish, by applying r1 , otherwise it will remain in its position by applying r2 . Note that rules r3 and r4 , albeit being enabled, can never actually be applied. As regards the leftmost wandering fish, the rules applicable to it are shown in Fig. 9c. In particular, rule r5 captures the behaviour in which the fish is attacked by the territorial fish tf(i ) in the relative position (0, −2), by moving the fish model the relaxed situation in which the fish wanders one cell above. The lower priority rules r6 , r7 , r8 , r9 , r10 from R w rel. around. Therefore, in case the escape rule r5 of the wandering fish cannot be applied, it just wanders around neglecting the presence of the territorial fish. This is done to keep the model simple; in fact, similarly to territorial fish, we could have provided a fallback behaviour in case the escape rule r5 by allowing the wandering fish to either try escape towards other nearby positions or, at worst, stay in the same position. Finally, as regards the other wandering fish, the rules applicable to it are just R w , since there are no territorial fish in the vicinity of it. rel. Fig. 10a shows the configuration reached after one step, in which rules r1 has been applied to the territorial fish, rule r5 has been applied to the leftmost wandering fish, and r10 to the other wandering fish. The rules applicable to the territorial fish are shown in Fig. 10b. Differently from the previous case, the distance of the two wandering fish from the nest is (i ) the same, namely 3. Thus there are two rules from R 3 applicable to tf(i ) , namely r11 in which it chases the leftmost wandering fish, and r12 in which it chases the other fish. The semantics prescribes that one of the two rules is selected (i ) non-deterministically. Note that there is also rule r15 from set R ⊥ enabled, used to ensure that the territorial fish goes back to the nest when its territory is free from wandering fish. However, since there are two wandering fish in the territory, such a rule can never be applied as there is always a higher priority rule applicable. As regards the wandering fish, the set of rules applicable to each of them are the same as in the previous step.

34

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

Fig. 11. Two possible configurations reached after the one shown in Fig. 10, by assuming that the rightmost wandering fish has not moved in the last step (i.e. rule r10 from Fig. 9 has been applied to it). Note that leftmost wandering fish has moved outside of the considered fragment, since it had a territorial fish chasing it (i.e. rule r5 from Fig. 9 has been applied to it).

Fig. 12. Simulation of a single territorial fish (depicted as a filled square) chasing a wandering fish (as circle) which has entered its territory. The territory guarded by the fish is green coloured (grey in black-and-white print).

Finally, Figs. 11a and 11b show the two possible configurations reached after one further step, assuming that rule r 6 is applied to the rightmost wandering fish, i.e. it did not move. On the other hand, the other wandering fish just moved outside the portion of the considered membrane, since the escape rule r5 has been applied to it. This depicts the two possible positions of the territorial fish resulting from the non-deterministic application of either r 11 or r12 to it. Discussion. We have implemented an ad hoc simulator1 for the territorial fish model, implementing the simulation algorithm presented in Section 4, together with the modifications discussed in Section 5 for dealing with ordinary objects and priority as a generic partial order. The model is composed of a single square membrane with side length 20. We report the results of the simulation of two variants of the same model, both containing one territorial fish guarding its territory. In the first simulation, there are three wandering fish, while the second simulation describes a more crowded environment, with 25 wandering fish. Initially the territorial fish is positioned in the nest, while the wandering fish are positioned randomly in the entire space. Fig. 12 shows a sequence of configurations reached during the first simulation. In this case, the nest is at position (9, 9). The territorial fish is depicted as a filled square, and its territory is shown as colour in green (grey in B/W print). Wandering fish are drawn as circles. Initially, at step (a), the territorial fish is in the nest, and the three wandering fish are outside the territory. Next (b), the wandering fish enters the territory, and the territorial fish starts moving towards it (c). Since the wandering fish keeps moving, at step (d) the territorial fish changes direction. Note that, at step (d), the wandering fish still does not considers itself attacked, since the territorial is at distance 3 from it (according to Fig. 8), thus it wanders around normally. At step (e), the two fish are adjacent. The wandering fish therefore moves upwards, trying to escape from the territorial fish, which moves upwards as well. Note that, due to the maximally-parallel semantics, both of them move simultaneously, thus in the next step they are still adjacent. At step (g) the wandering fish finally exits the territory, and the territorial fish goes back to the nest (steps (h–l)), since in the meantime no other wandering fish enters its territory. Fig. 13 shows a single territorial fish in a more crowded situation, with the aim of illustrating the behaviour in case of multiple wandering fish within the territory. In this simulation, the position of the nest is (10, 11) (the simulator places

1

Source code of the simulator is available upon request and from [38].

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

35

Fig. 13. Simulation illustrating the behaviour of a single territorial fish facing the presence of multiple wandering fish inside its territory. The model consists of 25 wandering fish.

the nest randomly in the space). We report a sequence of configurations reached during the simulation, where initially the territorial fish is chasing a wandering fish at distance 4 from the nest (relative position (0, 4)). Once the wandering fish exits the territory at step (b), the territorial fish starts to move back to the nest, but another wandering fish immediately enters the territory in position (2, −3). While it moves downwards, the situation changes again, since in step (g) another wandering fish enters the territory in position (4, 0). Some further steps are depicted in the figure; in the last step shown, the territory is free from wandering fish. Implementation. We have implemented the restricted simulation algorithm from Section 4, which only allows singular rules and assumes a total order priority relation. Since ordinary objects (denoting nests) occurs only as rule promoters, and are never modified, they can be dealt with easily as discussed in Section 5. Moreover, we have also exploited Theorem 5.1, which allowed us to handle the partial order priority relation with the restricted algorithm. The simulator has been implemented in the C# programming language, and simulations performed with a Mono virtual machine running on the Linux OS and using a 2.40 GHz CPU. Using this configuration, 1000 steps of the evolution are simulated in just a few seconds, for both of the models presented. This is made possible by the fact that the simulator is able to efficiently determine which rules are applicable to each position. In particular, the ad hoc simulator examines each fish individually, and dynamically generates all the sets of rules applicable to the fish (either territorial or wandering) ordered by their relative priority, as in the examples shown in Figs. 9, 10. For each positions, it then generates a random linear ordering of the rules, according to Section 5.2. All the lists of rules obtained, for all the position considered, are finally fed into the restricted resolution algorithm which finds a valid selection of rules, which is subsequently applied to the configuration to conclude the simulation of the step. 6.2. The schooling behaviour of herring In this section we present a Spatial P system model of the schooling behaviour of fish. Some kinds of fish, such as herring, swim together forming schools. This behaviour has different advantages, for example it is used against predators. One interesting problem with this kind of systems is that the single fish are able to organise themselves into schools by following local attraction rules. That is, this behaviour is not driven by some external entity which controls the behaviour of the single fish, but the common behaviour emerges as the result of local interactions between fish. Besides fish, similar behaviours also occur for other species, such as the flocking behaviour of birds, or the motion of herds of animals. The purpose of this model is to show that Spatial P systems can be used effectively to describe the behaviour of biological processes at a quite low level by specifying, in the setting of Spatial P systems, the Boids simulation algorithm for aggregate motion [39]. With respect to the territorial fish model presented in the previous section, the movement of fish in the schooling model is represented with higher accuracy; in fact, instead of moving one cell at a time, fish can move at a bigger distance (the actual maximum distance covered in a step is specified with a parameter). Furthermore, by modelling fish as mutually-exclusive objects, collisions between fish are prevented automatically, since such an aspect is already taken into account by the semantics, without requiring any other modelling effort from the modeller. In this manner, movement is represented more faithfully, since if the resulting position for a herring cannot be occupied because of another fish, then it is also allowed to move into intermediate positions, with decreasing priority associated with positions further away from the expected resulting position. As we will show in the section presenting the simulation results, this aspect is particularly important when simulating crowded situations which occur, for example, when there is an obstacle along the direction of movement of a school.

36

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

Fig. 14. The possible directions of a herring.

Fig. 15. The different forces, in the Boids model, determining the resulting direction of an individual.

6.2.1. Spatial P systems model The Spatial P systems model is based upon Boids model, a well-known algorithm for simulating this kind of aggregate motion proposed by Reynolds in 1987 [39]. In this model, the movement of a fish is driven by the sum of different forces, calculated with respect to distance and direction of movement of nearby fish. Forces are described by real-valued vectors in the continuous space R2 . Individuals can “see” only a small space around themselves, in order to determine the relative distance and behaviour of nearby individuals. In our model, a direction is associated with each herring, describing the direction where the herring is headed. There are eight possible directions, corresponding to the directions shown in Fig. 14. In the Spatial P system model, we use a different mutually exclusive symbol to denote each possible direction of a herring. In particular, the set of mutually exclusive objects is E = {1, 2, 3, 4, 5, 6, 7, 8}, where each symbol corresponds to a direction as shown in Fig. 14. There are no ordinary objects: V = ∅. Fig. 15 shows the forces affecting the behaviour of an individual (depicted as an upward arrow in the middle of the figure), with respect to the position and direction of other herring within its Field of View (grey-coloured cells). The bold black vector originating from the considered individual represents the resultant of the sum of different forces (depicted as dashed lines), which are different in the three cases. Firstly, there is a repulsive separation force that causes the fish to move away from nearby fish, whose magnitude is inversely proportional to their relative distances. Secondly, there is an attractive cohesion force which instead tends to keep the school together, by driving the fish towards the direction in which there are most of them. Thirdly, there is an alignment force, which tends to align the direction of the herring to the most common direction among nearby individuals. The sum of these three forces is used to determine the resulting direction of each fish. By using these rules, randomly positioned fish in a space are able to organise themselves into schools, in which individuals move in a coordinated manner. The Field of View (FOV) of a herring depends upon the actual direction of the herring, and the parameters VisibilityDistance and VisibilityAngle. Their role is exemplified in Fig. 16, where the cells forming the FOV are coloured in grey. Let us consider a VisibilityAngle α ∈ [0, 2π ], and a VisibilityDistance l ∈ R+ . Moreover, let us assume a reference coordinate system whose origin coincides with the position of the herring. The FOV of a herring having direction d ∈ {1, . . . , 8} is formally defined as the following set of positions w.r.t. the herring:









FOVlα (d) = p ∈ Z2   p 2  l, angle vect(d), p  α /2

where angle( p , q) ∈ [0, π ] denotes the angle between the two vectors p , q, and vect(d) denotes the unit vector corresponding to direction d. Formally, angle( p , q) = arccos(( p · q)/( p q)), where p · q denotes the dot product between p and q, and vect (d) = (cos θ, sin θ), where θ = (d − 1)2π /8. The model depends on a number of parameters, as detailed in the following:

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

37

Fig. 16. The Field of View (FOV, coloured in grey) of a herring having direction 3 (depicted in lower part of the figure as a black arrow pointing upwards). The FOV corresponds to cells forming an angle with respect to the direction of the herring of no more than α /2, and being within distance l from its position (α denotes the visibility angle, l denotes the visibility distance).

• Speed: a real value describing the speed of a herring (distance covered by a herring in one time unit) assumed to be the same for all of them;

• VisibilityDistance: the distance at which the herring can see, as depicted in Fig. 16; • VisibilityAngle: describes the Field of View (FOV) of the herring, namely the angular extent that the herring can see ahead of itself (Fig. 16);

• SeparationCoefficient: a real value greater than 1, used to scale down the repulsive effect of nearby individuals with respect to their relative distance; a greater value means that the repulsive effect decreases faster according to their distance; • CohesionCoefficient: how much the cohesion force influences the resulting direction; • AlignmentCoefficient: how much the alignment force influences the resulting direction. Fig. 17 summarises the parameters and the formulae used in our implementation of the Boids model to determine the movement of herring. Given a herring in a position p ∈ Z2 , the separation force v sep is computed as the sum of the vectors xi − p from each nearby herring to the current herring, multiplied by a value xi − p −s , where s is the SeparationCoefficient. Therefore the effect of the separation force decreases with the increase of the distance between the two herring. The cohesion force v coh corresponds to the vector from p to the “mean position” among the positions of nearby herring, multiplied by the CohesionCoefficient. The mean position is actually obtained by summing up all the vectors corresponding to the position of nearby herring, and dividing by the number of nearby herring. Finally, the alignment force v align (actually, a unit vector representing a direction) is obtained by summing up each unit vector corresponding to the direction of each nearby herring, multiplied by the AlignmentCoefficient. The vector force that determines the movement of a herring is computed as the resultant of its separation, cohesion and alignment forces, with the latter two multiplied by the CohesionCoefficient and the AlignmentCoefficient, respectively. Finally, the new position p reached by the herring after one time unit is obtained by adding to p the resulting displacement q. Formally, it corresponds to the resultant force v properly normalised, scaled by the Speed factor, and then discretised. Moreover, q also determines the resulting new direction d ∈ {1, . . . , 8} of the herring, which is defined as the direction among {1, . . . , 8} for which the corresponding vector vect(d ) has minimal angular distance to q. Evolution rules. In terms of Spatial P systems the model is described by using singular rules describing the movement of each individual herring and having the following form:

  (d)(0,0) J(d1 )x1 . . . (dk )xk K → d q where d ∈ {1, . . . , 8} is the current direction of the herring, d1 , . . . , dk ∈ {1, . . . , 8} and x1 , . . . , x8 ∈ Z2 are the directions and positions of nearby herring, d ∈ {1, . . . , 8} is the resulting direction, and q ∈ Z2 is the expected displacement of the herring. Promoters are used to see nearby herring in the FOV of the herring. A set of promoters (d1 )x1 . . . (dk )xk denotes a set of k nearby herring, each at position xi and moving towards direction di . For each possible configuration of nearby herring in the FOV there is a rule in the model describing how the herring is expected to move. Formally, given visibility angle α and visibility distance l, a configuration of herring in the FOV is described by a partial mapping Γ : FOVlα (d)  {1, . . . , 8}, which associates with each non-empty position in the FOV the direction of the herring contained in it. Such a mapping may be partial, in order to deal with empty cells; in fact, a position x for which Γ (x) is undefined means that such a position x is empty. We also define the set of positional objects specified by a mapping Γ as Γ = {d p | ( p , d) ∈ Γ }. By a slight abuse of notation, a rule can thus be expressed as (d)(0,0) JΓ K → (d )q . Priorities are used also to prevent conflicts that may arise as a consequence of the movements of herring. Given a configuration Γ of nearby herring there are different rules, with decreasing priority and the same left-hand side, that describe the possible movements of the herring. For example, Fig. 18 shows a herring, its resulting direction, and the

38

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

Fig. 17. Computation of the position reached by a herring in one step, by taking into account the position and direction of the herring in its Field of View. (All positions are expressed with respect to the reference coordinate system.)

possible resulting positions q1 , q2 , . . . , q6 , in decreasing order of priority. In particular, the highest priority rule tries to put an object in the most distant position along the chosen direction, i.e. q1 in Fig. 18. That is, the preferred position q1 is the one corresponding to a movement along the chosen direction, with the specified speed, which corresponds to displacement ( v / v ) · Speed , as computed in Fig. 17. If, as the result of the movement of other herring, position q1 happens to be occupied, then there is a rule with lower priority which tries to put an object in the immediately preceding position q2 , and so on until the least priority rule which tries to put the herring position q6 , adjacent to the current position. Therefore, given some promoters Γ and an object d, there are rules r1 > r2 > · · · > r6 , which try to put the object d in one of the resulting positions q1 , q2 , . . . , q6 , respectively. Formally, given a resulting displacement q = ( v / v ) · Speed (where v is defined in Fig. 17) the list of possible positions reachable(q) reachable by a herring, ordered by decreasing priority, is defined as q1 , . . . , qh  where: (i) q1 = q; (ii) for all j ∈ {1, . . . , h − 1}, it holds that q j , q j +1 are adjacent and that q j  > q j +1 ; (iii) qh is adjacent to (0, 0). By using the introduced notations, we are now ready to present the complete set of evolution rules of the model. Let Γ = ∅ be a non-empty configuration of nearby herring for a herring having direction d. Moreover, with reference to Fig. 17, let q be the resulting displacement, and d be the new resulting direction, and assume reachable (q) = q1 , . . . , qh . The following sets of rules are defined:



q

R dΓ = d(0,0) Γ

y

q y q y q y → dq 1 > d(0,0) Γ → dq 2 > · · · > d(0,0) Γ → dq h > d(0,0) Γ → d (0,0)

R d∅ = {d(0,0) → d , d(0,0) → d } where d = (d + 2) mod 8 and d = (d + 6) mod 8 model a rotation of the direction by either 90◦ or −90◦ , respectively. Γ

Γ

Moreover, the priority relation is such that ∀Γ1 , Γ2 : FOVlα (d)  {1, . . . , 8}, whenever Γ1 ⊇ Γ2 then R d 1  R d 2 .

Sets of rules R dΓ deal with the case in which the configuration of individuals in the FOV is non-empty, while R d∅ deal with the empty FOV. Rules in the first set move the herring along towards the resulting position. In order to correctly model the behaviour in which the fish first chooses a direction, and then tries to move along that direction, discarding any other direction, it is necessary to include in R d∅ a special “no-move” rule d(0,0) JΓ K → d (0,0) with lowest priority. Such a special rule describes the lack of movement for the object, and just captures the change of direction of the fish. The priority relation ensures that a rule which considers “more” promoters than another, i.e. for which more herring are specified in FOV, have higher priority. More precisely, if a rule r1 differs from another rule r2 only from the fact that the configuration of nearby individuals Γ2 of r2 is a strict subset of Γ1 from r1 , then it must hold r1 < r2 . This ensures that all of the herring that are present in the FOV are taken into account in the computation of the resulting position of the herring.

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

Fig. 18. An example of possible moves of a herring.

39

Fig. 19. Reference values of the parameters for the schooling model.

In case there are no nearby fish, rules in R d∅ model a random change of direction for the fish, by rotating of 45◦ either clockwise or counter-clockwise. Since the current definition of Spatial P systems does not provide probabilities among rules, it is sufficient to provide different rules with the alternative outcome. Note that the special “no-move” rule in R dΓ is

necessary to avoid falling back to rules R d∅ when there are higher priority rules enabled which are not applied because of conflicts among mutually exclusive objects. The rules are associated with the skin membrane of the system, having a square shape of 1000 × 1000 cells. In the simulations (presented in the next section), we use inner membranes to model obstacles. Note that the behaviour of herring both near the bounds of the skin membrane and near obstacles does not need to be explicitly modelled, since it is handled automatically by the semantics. In fact, recall that a rule is applicable to a position only if all its relative position fall inside the region of the membrane. For this reason, and due to the form of rules, any rule which tries to move a herring to a position outside membrane region is disabled. Moreover, the applicability of a rule also takes into account the promoters, which must occur inside membrane region. However, this also does not pose a problem, since there is a promoter only for positions occupied by herring, and not for empty positions. In this manner, on the one hand, any rule which refers to herring contained in the FOV but falling outside membrane region are not considered. On the other hand, lower priority rules which refers only to herring inside the membrane region can still be applied.

6.2.2. Simulation of the model The simulation of the schooling model has been performed by adapting the ad hoc simulator of the territorial fish model from Section 6.1. In particular, the simulator engine is the same, which implements the restricted simulation algorithm from Section 4; the only difference lies in the way the set of rules are determined.2 Nevertheless, similarly to the territorial fish, the simulator is able to efficiently generate all the rules applicable to each herring, ordered by their relative priority. In this manner, the algorithm efficiently determines which rules are applicable in each position. As for the territorial model, we also exploited the extensions presented in Section 5, in particular as regards the ability to use a priority relation which is a partial order. The model consists of a square grid of size 1000 × 1000 cells, which is initially filled with 1000 fish randomly spread in the whole grid. We performed simulations by varying all of the most relevant model parameters. Videos of the simulations performed are available online at [38]. The simulated time is 3000 time units and the reference parameter values used are shown in Fig. 19. A complete simulation of 3000 steps, running on the Linux OS and using a 2.40 GHz CPU, is carried out in just a few minutes, with some variability related to the number of herring simulated. Fig. 20 depicts the state at different times of a simulation in which the reference parameter values are used. The grid is represented as a bitmap image where each pixel represents a cell. Each non-empty position in the grid is represented by a black pixel. The directions of fish are not shown. Fig. 21 shows a graph representing the number of clusters of individuals (herring schools) and the covering of individuals by clusters (the fraction of individuals that belong to a cluster) over time. Clustering of herring positions has been performed by using the DBSCAN (Density-Based Spatial Clustering of Applications with Noise) algorithm [40], and in particular the implementation provided by the Python library Scikit-learn [41]. The DBSCAN algorithm is able to discover clusters of arbitrary shapes by taking into account the local density of points in the space. This ability is fundamental to our application, since herring schools can have different shapes while having more or less constant density throughout each school. Another important feature of this algorithm is the ability to filter out noise, namely lone points and small clusters, which we do not want to be counted as schools. The DBSCAN algorithm has two input parameters that govern its behaviour: Eps ∈ R+ and MinPts ∈ N. Informally, the meaning of these parameters is that each point in a cluster has to contain, within a radius of Eps (its neighbourhood), at least MinPts points. The actual values of the parameters used in our analyses are Eps = 0.5 and MinPts = 10, which have been empirically determined.

2

Source code of the simulator is available upon request and from [38].

40

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

Fig. 20. States of a simulation in which reference parameters are used. Images represent a portion of the simulation state corresponding to 200 × 200 cells at different times. Black points in the images are individuals of the modelled population of herring.

Fig. 21. Time course of the number of clusters and of the fraction of the herring population covered by clusters.

Figs. 20, 21 show that herring tend to form rather big schools. In particular, Fig. 21 shows that after a transient phase of about 1000 time units in which initially a high number of small clusters are formed, the population reaches a steady state in which, on average, 97.46% (on average) of the individuals are involved on average in 5.28 clusters. Average values are computed by excluding the initial transient phase. Moreover, Fig. 20c allows us to appreciate the shape of schools obtained by using the reference parameters; the school depicted is moving in the top-left direction. Fig. 22 shows simulation results obtained by varying parameter values. In the figure we have one panel for each parameter. Each panel reports, for a number of different parameter values, the average number of clusters and the average covering of a simulation. As before, average values are computed by excluding the initial transient period of 1000 time units. We remark that the average size of a cluster can be computed as follows:

Average cluster size =

Covering ratio × Population size Number of clusters

.

Panels 22a and 22b show that visibility is strongly related with the ability of herring of forming schools. The covering of the population by clusters increases with both VisibilityDistance and VisibilityAngle. On the other hand, the dynamics of the number of clusters and of their size is a bit different in the two cases. In the case of VisibilityDistance the number of clusters decreases for higher values of the parameter. In particular, for small values (i.e. 5) herring form many small clusters, and the coverage is quite low at less than 60%. For average values of the parameter (i.e. 10) they form more clusters, but still rather small; finally, for high values of the parameter (i.e. 16–30) they form a few big clusters. In the case of VisibilityAngle, for small values (i.e. 45–90 degrees) herring form a lot of small clusters. For average values of the parameter (i.e. 170–200 degrees) clusters become bigger, since the coverage in these cases is similar while the number of clusters decreases with the increase in the visibility angle; finally, for the highest value of the parameter (i.e. 360) they tend to split again into smaller clusters (10, on average) than for mid-range values.

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

41

Fig. 22. Simulation results (number of clusters and covering ratio) obtained by varying parameter values of the model.

The value of the VisibilityAngle parameter influences also the shape of clusters in a significant way. Fig. 23 shows examples of clusters for different values of the parameter. When the angle of visibility is narrow (e.g. 45 degrees) herring can see nearby herring only in a very small range in front of them. As a consequence, they are not able to form schools, as shown in Fig. 23a where only a few groups composed of a small number of individuals are formed. When VisibilityAngle reaches values greater or equal to 180 degrees herring can see nearby herring also on their sides. As a consequence, schools are formed successfully, and their shape is depicted in Fig. 23b. Finally, with full visibility (360 degrees) herring can see nearby herring also behind them. This makes school formation very fast but also undermines the stability of big clusters, therefore their size is smaller, as in the cluster depicted in Fig. 23c. Panel 22c shows results of simulations in which CohesionCoefficient is varied. This parameter has limited influence on the number of clusters and on the covering ratio. On the other hand, as shown in Fig. 24, the CohesionCoefficient has a very strong influence on the shape of clusters. In fact, even if the average size of clusters (in terms of number of involved

42

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

Fig. 23. Shapes of clusters by varying VisibilityAngle. The cluster in (b) is moving downwards and the one in (c) is moving towards the top left corner.

Fig. 24. Shapes of clusters by varying CohesionCoefficient. Both clusters are moving to the right.

Fig. 25. Shapes of clusters by varying AlignmentCoefficient. The cluster in (a) is moving from right to left and the cluster in (b) is moving towards the top-left corner.

individuals) does not change much among all of the parameter values, the density of the individuals in a cluster (and consequently the area covered by a cluster) changes significantly. Panel 22d shows results of simulations in which SeparationCoefficient is varied. In this case, what is influenced by the parameter is only the number of clusters, which decreases with higher values of the parameter. (Recall that a higher value of SeparationCoefficent means a more reduced repulsive effect.) Finally, Panel 22e shows results of simulations in which AlignmentCoefficient is varied. This parameter influences how fast a herring joins a school that is passing nearby. For small values of the parameter (i.e. 2) we have that herring are not very fast in joining schools, and, as a consequence, the number of schools remains rather small (herring may align too late to join a school). Another consequence of this slow reactivity of herring is that clusters are small, and have tapered shape (see Fig. 25a). Increasing the AlignmentCoefficient to 5 allows the formation of bigger clusters, and their number remains more or less the same for values 10–18. Fig. 25b shows the shape of a cluster with alignment coefficient 8.

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

43

Fig. 26. Snapshots at different times of a simulation in which reference parameters and a small obstacle are considered.

Fig. 27. Snapshots after 75 time units of two simulations with reference parameters and different obstacles.

In addition to the simulations in which we varied the model parameters, we performed also some simulations to study the dynamics of a school that encounters an obstacle during its movement. We considered a school consisting of 100 herring moving from left to right. In addition we considered three different obstacles: a small one (10 × 10), a medium one (10 × 20) and a big one composed by two parts (10 × 20 + 40 × 40). Obstacles are described as suitably shaped membranes. Simulations lasted 200 time units each. Snapshots from a simulation in which reference movement parameter values and a small obstacle are used are depicted in Fig. 26. In this situation the school passes the obstacle without difficulties. In Fig. 27 we show snapshots of simulations in which the same parameters are used, but where a medium obstacle (panel 27a) and a big obstacle (panel 27b) are considered. In these cases, the school passes the obstacles but then it continues as two separate schools. This behaviour is mainly due to the size of the obstacles compared with the value of the VisibilityDistance parameter. In fact, the value of the parameter is 14 while the distance between the two groups obtained after separation by the obstacle is 20 for the medium obstacle and 40 for the big obstacle (these values correspond to the size of the obstacles). Consequently, herring of one of the two schools cannot see any herring of the other school, and this prevents the reunion of the two groups. Panel 28a shows a snapshot of simulation in which VisibilityDistance has been set to 32. In this case, even with a medium size obstacle the school does not split. On the contrary, we can force a school to split by significantly increasing the AlignmentCoefficient. In Panel 28b we show a snapshot of a simulation in which this parameter is set to 40. Such a high value forces herring to continue to go straight without getting together even after passing a small obstacle. The variation of other parameter values either does not influence significantly the response of the school to the obstacle, or makes the dynamics of the school less regular. 7. Conclusions Spatial P systems, defined in [5], extend ordinary P systems by embedding objects and membranes in a two-dimensional discrete space. A peculiar feature of Spatial P system is the distinction among ordinary objects and mutually-exclusive objects. The former can be considered the formal counterpart of small molecules, while the latter can be used to capture, in an abstract way, the space occupied by objects. In fact, mutually-exclusive objects are subject to the constraint that, at any time, at most one of such a kind can be contained in a cell.

44

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

Fig. 28. Snapshots after 75 time units of two simulation in which parameters are varied.

In this paper, we have presented a new version of Spatial P systems which extends the original definition by allowing (i) a more general format for evolution rules, also including the specification of promoters; (ii) the specification of a priority relation among rules; and (iii) arbitrary membrane shapes. We have considered the problem of simulating Spatial P system models, and we have given a polynomial-time algorithm for a restricted version of models, composed of only mutually exclusive objects and in which rules have exactly one reactant and one product (called singular rules). Moreover, the priority relation is assumed to be a total order. Extensions of the simulation algorithm to deal with other variants have also been presented. Finally, we have discussed the general case, for which we hypothesise that a simulation algorithm would be NP-complete. We have left the investigation of these aspects as a future work. The applicability of the formalism to model real-world ecological systems have been demonstrated by two case studies: the first models the behaviour of territorial fish, which guard a small territory around their nests from other fish; the second captures the movement of fish school, such as herring, which are able to organise themselves into schools moving in a coordinated manner. These models show that Spatial P systems can be used effectively to capture spatial behaviours at different levels of abstraction. In particular, the territorial fish model is composed of a small number of rules, where priority plays a crucial role in obtaining a small model composed of easily readable rules. The simulations of such a model show the high expressivity of the formalism. On the other hand, the schooling model shows how Spatial P systems can be used at a lower level of abstraction, by encoding the fundamental Boids simulation algorithm for coordinated movement as a Spatial P system model. Finally, we remark the natural way of modelling fish as mutually-exclusive objects, which captures the inherent constraint of non-interpenetration between physical objects. In conclusion, the restricted kind of Spatial P systems that the proposed simulation algorithm deals with is still expressive enough to model interesting systems, which can then be simulated efficiently. As a future work, we plan to investigate the applicability of Spatial P systems to other kinds of systems. On the one hand, Spatial P systems can be used to develop epidemiology models in spatially dispersed populations. Another promising application context consists in the modelling of the heterogeneous distribution of biomolecules and organelles within the cell. The importance of taking into account spatial heterogeneity for the description of intracellular processes has been discussed in [42], while aspects of simulations have been investigated in [43]. On the other hand, we plan to implement a simulator for Spatial P systems models, by also taking advantage of the ad hoc simulator already developed for the case studies. This would allow an easier construction of models, without requiring any low-level programming. For this purpose, we also plan to develop an input language with the ability to express rule schemata, which can be very useful to specify bigger models composed of many similar rules, such as in the case of the schooling behaviour model we have presented. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]

G. P˘aun, Computing with membranes, J. Comput. Syst. Sci. 61 (1) (2000) 108–143. P. Bottoni, C. Martín-Vide, G. P˘aun, G. Rozenberg, Membrane systems with promoters/inhibitors, Acta Inform. 38 (2002) 695–720. G. P˘aun, Membrane Computing. An Introduction, Springer-Verlag, Berlin, 2002. P systems web page, http://ppage.psystems.eu. R. Barbuti, A. Maggiolo-Schettini, P. Milazzo, G. Pardini, L. Tesei, Spatial P systems, Nat. Comput. 10 (1) (2011) 3–16. A. Regev, E.M. Panina, W. Silverman, L. Cardelli, E.Y. Shapiro, BioAmbients: An abstraction for biological compartments, Theor. Comput. Sci. 325 (1) (2004) 141–167. L. Cardelli, Brane calculi – Interactions of biological membranes, in: V. Danos, V. Schachter (Eds.), Computational Methods in Systems Biology, in: Lecture Notes in Computer Science, vol. 3082, Springer, Berlin/Heidelberg, 2005, pp. 257–278. J.D. Murray, Mathematical Biology: I. An Introduction, 3rd edition, Springer-Verlag, 2002. J.D. Murray, Mathematical Biology: II. Spatial Models and Biomedical Applications, 3rd edition, Springer-Verlag, 2003. J.V. Neumann, Theory of Self-Reproducing Automata, University of Illinois Press, 1966. M. Patel, S. Nagl, Mathematical models of cancer, in: S. Nagl (Ed.), Cancer Bioinformatics: From Therapy Design to Treatment, John Wiley & Sons, 2006, pp. 59–93, Ch. 4. R. Barbuti, A. Maggiolo-Schettini, P. Milazzo, G. Pardini, Spatial calculus of looping sequences, Theor. Comput. Sci. 412 (43) (2011) 5976–6001.

R. Barbuti et al. / Theoretical Computer Science 529 (2014) 11–45

45

[13] M. Coppo, F. Damiani, M. Drocco, E. Grassi, A. Troina, Stochastic calculus of wrapped compartments, in: Proc. Quantitative Aspects of Programming Languages 2010, QAPL 2010, in: EPTCS, vol. 28, 2010, pp. 82–98. [14] C. Calcagno, M. Coppo, F. Damiani, M. Drocco, E. Sciacca, S. Spinella, A. Troina, Modelling spatial interactions in the arbuscular mycorrhizal symbiosis using the calculus of wrapped compartments, in: Proc. Computational Models for Cell Processes 2011, CompMod’11, in: EPTCS, vol. 67, 2011, pp. 3–18. [15] L. Bioglio, C. Calcagno, M. Coppo, F. Damiani, E. Sciacca, S. Spinella, A. Troina, A spatial calculus of wrapped compartments, in: Proc. Membrane Computing and Biologically Inspired Process Calculi 2011, MeCBIC’11, 2011, pp. 25–39. [16] M. John, R. Ewald, A.M. Uhrmacher, A spatial extension to the π -calculus, in: Proceedings of the First Workshop From Biology To Concurrency and Back, FBTC 2007, Electron. Notes Theor. Comput. Sci. 194 (3) (2008) 133–148. [17] L. Cardelli, P. Gardner, Processes in space, in: F. Ferreira, B. Löwe, E. Mayordomo, L. Mendes Gomes (Eds.), Programs, Proofs, Processes, in: Lecture Notes in Computer Science, vol. 6158, Springer, Berlin/Heidelberg, 2010, pp. 78–87. [18] M. Antonaki, A. Philippou, A process calculus for spatially-explicit ecological models, in: G. Ciobanu (Ed.), Proceedings 6th Workshop on Membrane Computing and Biologically Inspired Process Calculi, in: EPTCS, vol. 100, 2012, pp. 14–28. [19] D. Besozzi, P. Cazzaniga, D. Pescini, G. Mauri, Modelling metapopulations with stochastic membrane systems, Biosystems 91 (3) (2008) 499–514. [20] M. Cardona, M.À. Colomer, M.J. Pérez-Jiménez, D. Sanuy, A. Margalida, Modeling ecosystems using p systems: The bearded vulture, a case study, in: D.W. Corne, P. Frisco, G. P˘aun, G. Rozenberg, A. Salomaa (Eds.), Membrane Computing: 9th International Workshop, in: Lecture Notes in Computer Science, vol. 5391, 2009, pp. 137–156. [21] M. Cardona, M.À. Colomer, A. Margalida, I. Pérez-Hurtado, M.J. Pérez-Jiménez, D. Sanuy, A p system based model of an ecosystem of some scavenger birds, in: G. P˘aun, M.J. Pérez-Jiménez, A. Riscos-Núñez, G. Rozenberg, A. Salomaa (Eds.), Membrane Computing, in: Lecture Notes in Computer Science, vol. 5957, Springer, Berlin/Heidelberg, 2010, pp. 182–195. [22] D. Pescini, D. Besozzi, G. Mauri, C. Zandron, Dynamical probabilistic P systems, Int. J. Found. Comput. Sci. 17 (1) (2006) 183–204. [23] F.J. Romero-Campero, M.J. Pérez-Jiménez, A model of the quorum sensing system in vibrio fischeri using P systems, Artif. Life 14 (1) (2008) 95–109. [24] M. Cardona, M.À. Colomer, A. Margalida, A. Palau, I. Pérez-Hurtado, M.J. Pérez-Jiménez, D. Sanuy, A computational modeling for real ecosystems based on P systems, Nat. Comput. 10 (1) (2011) 39–53. [25] M.À. Colomer, A. Margalida, D. Sanuy, M.J. Pérez-Jiménez, A bio-inspired computing model as a new tool for modeling ecosystems: The avian scavengers as a case study, Ecol. Model. 222 (1) (2011) 33–47. [26] C. Martín-Vide, G. P˘aun, J. Pazos, A. Rodríguez-Patón, Tissue P systems, Theor. Comput. Sci. 296 (2) (2003) 295–326. [27] R. Freund, G. P˘aun, M.J. Pérez-Jiménez, Tissue P systems with channel states, Theor. Comput. Sci. 330 (1) (2005) 101–116. [28] F. Bernardini, M. Gheorghe, Population P systems, J. Univers. Comput. Sci. 10 (5) (2004) 509–539. [29] F.J. Romero-Campero, J. Twycross, M. Cámara, M. Bennett, M. Gheorghe, N. Krasnogor, Modular assembly of cell systems biology models using P systems, Int. J. Found. Comput. Sci. 20 (3) (2009) 427–442. [30] F.J. Romero-Campero, N. Krasnogor, An approach to the engineering of cellular models based on P systems, in: Mathematical Theory and Computational Practice, in: Lecture Notes in Computer Science, vol. 5635, Springer, Berlin/Heidelberg, 2009, pp. 430–436. [31] P. Kefalas, I. Stamatopoulou, Modelling of multi-agent systems: Experiences with membrane computing and future challenges, in: P. Milazzo, M. de J. Pérez-Jiménez (Eds.), Proceedings First Workshop on Applications of Membrane Computing, Concurrency and Agent-based Modelling in POPulation Biology, in: EPTCS, vol. 33, 2010, pp. 71–82. [32] N. Busi, Using well-structured transition systems to decide divergence for catalytic P systems, Theor. Comput. Sci. 372 (2–3) (2007) 125–135. [33] O. Andrei, G. Ciobanu, D. Lucanu, A rewriting logic framework for operational semantics of membrane systems, Theor. Comput. Sci. 373 (3) (2007) 163–181. [34] R. Barbuti, A. Maggiolo-Schettini, P. Milazzo, S. Tini, Compositional semantics and behavioral equivalences for P systems, Theor. Comput. Sci. 395 (1) (2008) 77–100. [35] R. Barbuti, A. Maggiolo-Schettini, P. Milazzo, S. Tini, Compositional semantics of spiking neural P systems, J. Log. Algebr. Program. 79 (6) (2010) 304–316. [36] R. Barbuti, A. Maggiolo-Schettini, P. Milazzo, S. Tini, An overview on operational semantics in membrane computing, Int. J. Found. Comput. Sci. 22 (1) (2011) 119–131. [37] C. Özturan, On finding hypercycles in chemical reaction networks, Appl. Math. Lett. 21 (9) (2008) 881–884. [38] Research group on modelling, simulation and verification of biological systems – supplemental material, http://www.di.unipi.it/msvbio/wiki/ SpatialPSystems. [39] C.W. Reynolds, Flocks, herds and schools: A distributed behavioral model, SIGGRAPH Comput. Graph. 21 (1987) 25–34. [40] M. Ester, H.-P. Kriegel, J. Sander, X. Xu, A density-based algorithm for discovering clusters in large spatial databases with noise, in: E. Simoudis, J. Han, U. Fayyad (Eds.), Proceedings of the Second International Conference on Knowledge Discovery and Data Mining, 1996, pp. 226–231. [41] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, D. E., Scikit-learn: Machine learning in python, J. Mach. Learn. Res. 12 (2011) 2825–2830. [42] J. Elf, M. Ehrenberg, Spontaneous separation of bi-stable biochemical systems into spatial domains of opposite phases, Syst. Biol. IEE Proc. 1 (2) (2004) 230–236. [43] J. Hattne, D. Fange, J. Elf, Stochastic reaction–diffusion simulation with MesoRD, Bioinformatics 21 (12) (2005) 2923–2924.