Modeling regenerative processes with Membrane Computing

Information Sciences 381 (2017) 229–249 Contents lists available at ScienceDirect Information Sciences journal homepage: www.elsevier.com/locate/ins Modeling regenerative processes with membrane computing Manuel García-Quismondo a , ∗, Michael Levin b , Daniel Lobo c a National Zoological Park, Smithsonian Institution, 3001 Connecticut Ave. NW, Washington, DC 20 0 08, USA b Department of Biology, Tufts University, 200 Boston Ave, Medford, MA 02155, USA c Department of Biological Sciences, University of Maryland, Baltimore County, 10 0 0 Hilltop Circle, Baltimore, MD 21250, USA a r t i c l e i n f o Article history: Received 5 January 2016 Revised 11 October 2016 Accepted 25 November 2016 Available online 27 November 2016 Keywords: Planarian Modeling Regeneration P systems Membrane computing Bioinformatics a b s t r a c t Understanding the remarkable ability of some organisms to restore their anatomical shape following the amputation of large parts of their bodies is currently a major unsolved question in regenerative biology and biomedicine. Despite rapid advances in the molecular processes required for regeneration, a systems level, algorithmic understanding of this process has remained elusive. For this reason, the field needs new computational paradigms to help model the flow of information during regeneration. Membrane computing is a branch of natural computing that studies the properties and applications of theoretical computing devices known as P systems. These systems are an abstraction of the structure and functioning of a living cell, as well as its organization in tissues. Here, we propose a model of regenerative processes in planarian worms based on P systems, which recapitulates several aspects of regenerative pattern regulation. Our results demonstrate that it is possible to apply a novel computational framework to help understand pattern regulation in regenerative biology. ©2016 Elsevier Inc. All rights reserved. 1. Introduction Membrane computing is a branch of natural computing introduced by P ˘ aun in 1998 [63] . This discipline studies the properties and applications of theoretical computing devices known as P systems, which are an abstraction of the structure and functioning of a living cell, as well as its organization in tissues and other higher-order structures [27,45,46,55–57,64,79,81– 83] . A P system defines a structure consisting of a graph of interconnected compartments , which are symbolic entities that could represent biological cells or group of cells in an organism. Each compartment contains a multiset of objects, which may be molecules with potentially different electrical charges. There are rules that dictate how objects are created, removed, or migrate across compartments and during a computation, objects are processed by means of rewriting rules . P system rules are abstractions of the biochemical and electrical reactions that occur inside living cells. Objects can move between two compartments if there exists an edge or link that connects the two compartments. This represents the biological signaling that occurs between two cellular entities. Spurred by the success of membrane computing as a modeling framework, we propose herein a new class of P systems called regenerative P systems . This computational framework provides a mathematical formalism to model regenerative processes in biological organisms. One of the major unsolved questions facing basic biology and biomedicine today is understanding the remarkable ability of some organisms to restore their anatomical shape following amputation of large portions of their bodies [1,73] . For ∗Corresponding author. E-mail addresses: mgarcia.univ[email protected] (M. García-Quismondo), Michael.Le[email protected] (M. Levin), [email protected] (D. Lobo). http://dx.doi.org/10.1016/j.ins.2016.11.017 0020-0255/© 2016 Elsevier Inc. All rights reserved. 230 M. García-Quismondo et al. / Information Sciences 381 (2017) 229–249 example, salamanders can regenerate their limbs, eyes, jaws, hearts, and portions of the brain [7] . Understanding shape homeostasis and restoration is critical, not only for advances in regenerative medicine, but also for fundamental issues in developmental biology and evolution [26,29,32,62,69] . While much progress has been recently made in identifying the molecular signaling required for regeneration [6] , we are still very far from understanding the control of shape. The molecular details are becoming clearer, but the information flow, computation, and control policies for organizing cell behaviors towards large-scale anatomical outcomes are largely unknown [33] . Planaria are a powerful model system in this field, because they can regenerate every part of their bodies [34,70] . We need to develop computational models that can explain how regeneration creates exactly what is missing, in the right location, and stops when the target morphology has been achieved. Several types of models have been proposed [11,49,80] , yet constructive models are the important to explore. These models show the sufficient steps to restore pattern, not only pathways of events necessary for regeneration [2,35,38] . This is largely an unexplored field, and we are still groping for the correct formalism and appropriate data representation for such models and the algorithms they implement [3,5,50,68,74] . To enrich the field of possibilities, and introduce biologists to another way to think about the controls of pattern regulation, we present here a model of planarian regeneration based on P systems. A number of modeling formalisms have been attempted in planaria, including [15,35,72,74] and others described in [14] ; however, P systems have not been investigated in this context. We decided to select membrane computing as our modeling framework for regenerative processes because of its biologically–relevant properties. The explicit spatial structure of compartments in P systems facilitates the abstraction of morphological regions in regenerative biological organisms. Furthermore, the intercellular communication and signaling through proteins, small molecules, or biophysical phenomena essential to regeneration can be mathematically modeled with signals and objects between the compartments. Our model of regeneration can be easily extended by quantizing these signals in packages. This permits us to model the magnitude of these cellular signals as the difference in the number of objects of a given class. In addition, the modular property of P systems permits us to bundle any specific aspect of regenerative processes into a subset of rules. These rules can then be modified without affecting the overall dynamics of the model [13,66] . For example, the specific signaling mechanisms for regenerating the brain and eyes can be modeled with an isolated set of rules without interferences with other aspects of the model. In our proposed model, the flatworm morphology is represented as a rectangular compartment grid. Each compartment represents a morphological region, and each region belongs to one of the following parts in the morphology: head, trunk or tail. Signals can be sent between two compartments as long as there exists a link between them. We use objects to express the membership of a compartment to a worm region and to represent organs inside these regions. We also use links to denote communication between compartments. Our model aims to reproduce molecular and electrical signals sent between cells at a local level, which eventually configure emergent morphologies [33,75] . This modeling approach bears some resemblance to peer-to-peer networks. In both frameworks, a consensus system behavior emerges from the pairwise communication between nodes [16] . In fact, fault diagnosis of electrical and computer networks has been successfully modeled using P systems [58] . We have tested our modeling framework with a set of planarian in silico experiments that recapitulate the resultant phenotypes observed in vivo . Section 2 introduces regenerative P systems as a new class of membrane computing devices. In Section 3 , the model is described and its behavior is explained. Section 4 simulates the model for several scenarios, encompassing interesting case studies. Section 5 compares the model predictions in these scenarios with experimental results. Finally, our conclusions are reported in Section 6 . 2. Regenerative P systems In this section, we propose a novel framework named regenerative P systems to model biological regeneration processes. This framework is mainly inspired by kernel P systems, a previous approach in membrane computing [20–22,28] . Prior to defining the framework, we describe some preliminary concepts as defined in [18] : An alphabet is a non-empty set whose elements are called symbols . Likewise, a multiset w over an alphabet is a pair w = (, f) where f : → N is a mapping. For each x ∈ we say that f ( x ) is the multiplicity of the symbol x in w . If w = (, f) is a multiset, then its support is defined as supp(w ) = { x ∈ | f(x ) > 0 } . A multiset is finite if its support is a finite set. A set is a multiset such that the multiplicity of each element of its support is greater or equal to 1, that is, the multiset can contain more than one object of the same class. If w = (, f) is a finite multiset over , and supp(w ) = { a 1 , . . . , a k } then it will be denoted as w = a f(a 1 ) 1 . . . a f(a k ) k (here the order is irrelevant), and we say that f(a 1 ) + ···+ f(a k ) is the cardinal of w , denoted by | w |. The empty multiset is denoted by ∅ . We also denote by M ( ) the set of all finite multisets over . Consider w 1 = (, f 1 ) and w 2 = (, f 2 ) multisets over . We define the following concepts: •The union of w 1 and w 2 , denoted by w 1 + w 2 is the multiset ( , g ), where g = f 1 + f 2 , that is, g(x ) = f 1 (x ) + f 2 (x ) for each x ∈ . Likewise, when w 1 is updated as w 1 ← w 1 + w 2 , we say that w 2 is generated into w 1 . •The relative complement of w 2 in w 1 , denoted by w 1 ࢨw 2 is the multiset ( , g ), where g(x ) = f 1 (x ) −f 2 (x ) if f 1 ( x ) ≥f 2 ( x ) and g(x ) = 0 otherwise. Likewise, when w 1 is updated as w 1 ← w 1 −w 2 , then we say that w 2 is consumed from w 1 . We also say that w 1 is a submultiset of w 2 , denoted by w 1 ⊆w 2 , if f 1 ( x ) ≤f 2 ( x ) for each x ∈ . M. García-Quismondo et al. / Information Sciences 381 (2017) 229–249 231 Let w = (, f) a multiset over and A a set. We define the intersection w ∩ A as the multiset ( , g ), where g(x ) = f(x ) for each x ∈ ∩ A , and g(x ) = 0 otherwise. We also must define some other concepts inspired by the idea of guard found in [22] . Please note that some definitions are slightly different from those provided in that reference: Let be an alphabet and let = { a | a ∈ } , provided ∩ = ∅ . We say that Co() = ∪ is the alphabet of object constraints over . We also recursively define a guard over an alphabet as a predicate over M ( ). To do so, let w ∈ M ( ) be a multiset over at a given instant. When a guard does not evaluate to true, then it evaluates to false: •An empty string is a guard over , and we say that it always evaluates to true . •An object constraint a ∈ Co ( ) is a guard over . We say that a evaluates to true for w if (a ∈ ∧ a ∈ supp(w )) ∨ (a ∈ ∧ a / ∈ supp(w )) . •Let g 1 , g 2 be two guards over . Then, g 1 ∧ g 2 and g 1 ∨ g 2 are guards. g 1 ∧ g 2 evaluates to true if both g 1 and g 2 evaluate to true. Likewise, g 1 ∨ g 2 evaluates to true if either g 1 or g 2 (or both) evaluate to true. In addition, the concept of script (see below) has been inspired by Probabilistic Guarded Scripted P systems [19] . When a script is generated, it triggers a procedure on the current configuration of the system. The execution of this procedure may have different effects on the ongoing computation. As an example, the multisets associated with each one of the compartments can be updated or the computation can be halted. A regenerative P system of degree ( m, n ), m ≥1, n ≥1 is defined as: = (H, G, , , { w h i | 1 ≤i ≤q } , R , { (s k , τk ) | 1 ≤k ≤n } ) , where: •H = { h 1 , . . . , h m } is an alphabet of labels . •G = (V, S) is a non-directed graph without loops where V = { h 1 , . . . , h q } , q ≤m is the initial set of nodes labeling compartments in the system and S is the initial set of edges (links) over V . •and are finite alphabets such that ∩ = ∅ and σ∈ . Elements in are called objects , and elements in are called scripts . •For each i, 1 ≤i ≤q, w h i ∈ M() is the initial multiset associated with compartment h i . •R is a finite set of rules of the following types: −[ u ] h i → [ v , s ] h j with u, v ∈ M ( ), 1 ≤i, j ≤m and s ∈ . If s = σ, then we omit it in the notation. −[ u ] h i ⇒ [ v ] h j or [ u ] h i ⇒ [ v ] h j with u, v ∈ M ( ), 1 ≤i, j ≤m , and h i  = h j . −[ u ] h i ⇒ [ v ] h i [ w ] h j with u, v, w ∈ M ( ), 1 ≤i, j ≤m and h i  = h j . •Each rule from R has an associated priority, that is, there exists a function p R from R into Z . We denote p R (r) in r as (p R (r)) prior to r , e.g., r ≡(p R (r))[ u ] h i → [ v , s ] h j . If p R (r) is equal to 0, then we omit it in the notation. •Each rule r ∈ R has an associated guard g r . We denote g r in r as { g r } prior to r and to p R (r) , e.g., r ≡ { g r } (p R (r))[ u ] h i → [ v , s ] h j . If g r is the empty string, then we omit it in the notation. •{( s k , τk ) | 1 ≤k ≤n } is a set of pairs which contains the pair ( σ, F ) such as s k ∈ and τk is a computable predicate over the set of configurations of (see below). F is a predicate which always evaluates to false. For each two pairs ( s i , τi ), ( s j , τj ), if i  = j then s i  = s j . We define a configuration C t = (V t , S t , { w h i ,t | h i ∈ V t } ) of a regenerative P system of degree ( m, n ) at any given time t ≥0 as a tuple whose components are the following: a) V t is the (possibly equal) subset of compartment labels from H ;b) S t the set of links between labels; c) w h i ,t | h i ∈ V t are the multisets of objects associated with each of the compartments of the system. The initial configuration of this regenerative P system is C 0 = (V, S, { w h i | h i ∈ V } ) Next, let w h i ,t be the multiset associated with compartment h i at any given time t ≥0. Likewise, let S t be the set of links in configuration C t . Then, we define the applicability and application effects of rules in the regenerative P system as follows: •A rule of the form r ≡{ g r } (p R (r))[ u ] h i → [ v , s ] h j is applicable on compartment h i if g r evaluates to true, u ⊆w h i ,t and (h i = h j ∨ (h i , h j ) ∈ S) . Upon its application, u is consumed from w h i ,t and v is generated into w h j ,t . The role of script s is explained later in this section. •A rule of the form r ≡{ g r } (p R (r))[ u ] h i ⇒ [ v ] h j is applicable on compartment h i if g r evaluates to true, u ⊆w h i ,t , h j ∈ V and ( h i , h j ) ∈ S . Upon its application, u is consumed from w h i ,t and v is generated into w h j ,t . Moreover, S ← S ∪ {( h i , h j )}. •A rule of the form r ≡{ g r } (p R (r))[ u ] h i ⇒ [ v ] h j is applicable on compartment h i if g r evaluates to true, u ⊆w h i ,t , h j ∈ V and ( h i , h j ) ∈ S . Upon its application, u is consumed from w h i ,t and v is generated into w h j ,t . Moreover, S ← S ࢨ{( h i , h j )}. •A rule of the form r ≡{ g r } (p R (r))[ u ] h i ⇒ [ v ] h i [ w ] h j is applicable on compartment h i if g r evaluates to true, u ⊆w h i ,t and h j ∈ V . Upon its application, u and v are respectively consumed from and generated into w h i ,t . Also, a new compartment h j with associated multiset w is created, updating V t as V t ← V t ∪ { h j }. Finally, S ← S ∪ {( h i , h j )}. Given a set of rules R , we say that a multiset of applicable rules at time t is mutually applicable if they altogether do not consume more objects than available. Moreover, we say that a mutually applicable multiset of rules is maximal in the case that, if we add any rule, then the resulting multiset is not mutually applicable. 232 M. García-Quismondo et al. / Information Sciences 381 (2017) 229–249 Next, the semantics of the model is defined. At every time t and until the computation halts, a transition step consists of the following stages: 1) all applicable rules are selected. If there exists none of such rules, then the current configuration is a halting configuration and the computation halts; 2) from these rules, those with the maximum priority are kept; 3) from the remainder, a maximal multiset of rules is non-deterministically chosen. The following restriction is applied: if any rule from this set does not consume any object, then it must appear once in the multiset; 4) all rules from such a multiset are applied; 5) for every rule of the form r ≡{ g r } (p R (r))[ u ] h i → [ v , s ] h j which has been applied in the current transition step, the predicate associated with s is evaluated for the current configuration. Objects consumed and generated in the current transition step are considered for this application. If any of these predicates evaluates to true, then the current configuration is a halting configuration and the computation halts. Finally, we define a computation C = { C 0 , C 1 , . . . , C n } as a sequence of configurations such that: a) the first term of the sequence is the initial configuration of ; b) each remaining configuration is obtained by applying the rules defined in  on the previous configuration in a transition step; c) the last term of the sequence is a halting configuration. In summary, a regenerative P system can be viewed as a set of q interconnected compartments , labeled by h 1 , ... , h q such that: a) w h 1 , ... , w h q are finite multisets over representing the objects initially placed in the compartments of the system; b) G is a non-directed graph whose edges (links) specify explicit connections among compartments; c) R is the set of rules that define the transition steps between configuration of the system and each rule r is associated with an integer number p R E (r) in Z describing a priority amongst rules, in such a way that only those applicable rules with the maximum priority can be applied; d) is a set of symbols triggering halting predicates, in such a way that, if one of the triggered predicates is complied at a given time, then the system halts. To illustrate the computations and the application of rules in a regenerative P system, let us consider a system such as: = (H = { 1 , 2 , 3 } , G = {{ 1 , 2 } , { (1 , 2) }} , = { a, b, c, d, e } , = { s } , { w 1 = { a 4 , c, d} , w 2 = { a 2 , b 2 }} R = { r 1 ≡{ ¯ d } (−1)[ a 2 , e ] 1 → [ b, c 2 , s ] 2 , r 2 ≡{ ¯ e } (0)[ a 2 , d] 1 ⇒ [ b, e ] 1 , r 3 ≡{ d} (0)[ c] 1 ⇒ [ c] 2 , r 4 ≡{ d} (0)[ b 2 ] 2 ⇒ [ b] 3 , r 5 ≡(1)[ a 2 ] 2 ⇒ [ b, c, d] 2 [ c] 3 } , { (s, τ) } ) where τ( C ), being C a configuration of , causes the system to halt. Let us consider the initial configuration of this regenerative P system. In this configuration, only rules r 2 , r 3 and r 5 are applicable. Although there exists a link between compartments 1 and 2, r 1 cannot be applied because compartment 1 contains one d and zero e  s, and r 1 requires that compartment 1 contains zero d  s and one e . Moreover, rule r 2 is associated with the same compartment and has a higher priority (1). In addition, since there exists no link between compartments 2 and 3 and r 5 has a higher priority than r 4 . By applying r 2 and r 3 , two a  s, one c and one d are consumed and one b and one e are generated in compartment 1. Consequently, the multiset associated with compartment 1 becomes { a 2 , b, e } in the next configuration. Plus, one object c is sent to compartment 2. By applying rule r 5 , two a  s are consumed and one b , one c and one d are generated in compartment 2. In addition, compartment 3 is created and linked with compartment 2. This compartment contains now one c . Therefore, the multiset associated with compartment 2 becomes { b 3 , c 2 , d } and the multiset associated with compartment 3 becomes { c }. In the next configuration, only r 1 and r 4 can be applied. Both r 2 and r 3 require one d in compartment 1, but this compartment contains zero d  s in the present configuration. Likewise, r 5 cannot be applied because there are zero a  s in compartment 2 and compartment 3 has already been created. Upon applying r 1 , one e and two a  s are consumed from compartment 1 and one b and two c  s are generated in compartment 2. Moreover, the script s is generated. Likewise, the application of r 4 consumes two b  s from compartment 2 and generates one b in compartment 3. After all rules in this transition step have been applied, the link between compartments 2 and 3 is removed. Thus, the multiset associated with compartment 1 becomes { b }, the multiset associated with compartment 2 becomes { b 2 , c 4 , d }, and the multiset associated with compartment 3 becomes { b, c }. Once these rules have been applied, the predicate associated with the script s is triggered and the system halts. Table 1 describes this computation. 3. Regenerative P system model of planarian regeneration During planarian regeneration, precise intercellular signaling mechanisms are responsible for the growth and repatterning of the body regions (head, trunk, and tail) and organs with correct proportions and in the correct locations [34] . M. García-Quismondo et al. / Information Sciences 381 (2017) 229–249 233 Table 1 An example of a regenerative P system computation. Initial configuration [ a 4 , c, d] 1 −[ a 2 , b 2 ] 2 Transition step 1 r 2 ≡{ ¯ e } (0)[ a 2 , d] 1 ⇒ [ b, e ] 1 r 3 ≡{ d }(0)[ c ] 1 ⇒ [ c ] 2 r 5 ≡(1)[ a 2 ] 2 ⇒ [ b, c, d ] 2 [ c ] 3 Configuration 2 [ a 2 , b, e ] 1 −[ b 3 , c 2 , d] 2 −[ c] 3 Transition step 2 r 1 ≡{ ¯ d } (−1)[ a 2 , e ] 1 → [ b, c 2 , s ] 2 r 4 ≡{ d} (0)[ b 2 ] 2 ⇒ [ b] 3 Configuration 3 [ b] 1 −[ b 2 , c 4 , d] 2 [ b, c ] 3 The system halts. Table 2 Model parameters used in the experiments. Parameter Description Value length Maximum distance from head to tail in any 4 morphology generated by the model width Maximum distance from left to right in any 3 morphology generated by the model lp Maximum number of time steps that any part 7 allocation stage can take bclkmax Maximum number of time steps that all part 31 allocation stages can take, plus the required overhead to discard unused signals oclkmax Maximum number of time steps that all organ 8 allocations can take, plus the required overhead to discard unused signals Although this process is not fully understood, it is thought to involve proteins and small molecules that are secreted by cells in different regions. These molecules travel by diffusion or active transport to other regions of the worm. Once there, they act as signals to instruct the fate of the new regenerated tissue [4,31,76] . For example, the family of proteins in the Wnt pathway are involved in the regeneration of a tail [60] , whereas notum is involved in the signaling of the head region [61] . Furthermore, bioelectrical signals related to head or tail regeneration can also travel through gap junctions in the cells and act as regeneration signals [54] . Understanding the exchange and processing of these signals is an essential step for inferring the mechanisms responsible for the regeneration in the worm [36] . To this end, we present here a model within the framework of regenerative P systems that captures, in an abstract fashion, the cell signaling processes that enable a planarian worm to regenerate. Our model repeats iteratively the same cycle (sequence of steps), which consists of two sub-cycles. On the first sub-cycle, the model orients itself: it locates the head, the tail, the left-hand side and the right-hand side. If necessary, it also grows at a rhythm of one compartment per cycle. This procedure is illustrated in Figs. 2 –5 . The second half-cycle involves the allocation of organs inside the morphology. This process is performed case by case, depending on the organ that is being allocated. Some organs must be allocated in the center of the morphology, such as the pharynx. In contrast, other organs must be on the outside regions of the morphology, such as the eyes and brains. In total, the model takes 39 time steps to complete a cycle. Fig. 1 represents a regenerative P systems model of a planarian worm of length 4 and width 3 (see Table 2 ). The compartments labeled in make up the trunk, those labeled in and compose the left-hand and right-hand sides, respectively; and those labeled in and represent the tail and the head, respectively. Compartments with two colors belong to two regions simultaneously. The circles in each compartment represent the organs in each morphology part. Moreover, a black line between two compartments represents a link between them, through with objects (some proteins, small molecules, etc.) can pass. The model of planarian regeneration consists of a regenerative P system defined as follows: = (H, G, , , { w h i | 1 ≤i ≤m } , R , { (s k , τk ) | 1 ≤k ≤n } ) , where: •H = { (1 , 1) , (1 , 2) , . . . , (1 , np + 1) , (2 , 1) , (2 , 2) , . . . , (length, width ) } , where length and width are the maximum distance from head to tail and left to right, respectively, in any possible morphology generated by the model. •G = (V, S) is a non-directed graph representing the initial morphology regions and their connections. Each node in G represents a region in the initial morphology, and each edge represents an initial link between two regions. The specific nodes in V depend on the simulated experiment. 234 M. García-Quismondo et al. / Information Sciences 381 (2017) 229–249 Fig. 1. A visual representation of a regenerative P system model of a planarian worm. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article). We say that two compartments are adjacent if the first or the second index (but not both) of their labels differ by one unit. For instance, [ ] 1,1 , [ ] 1,2 are adjacent, but not [ ] 2,1 , [ ] 1,2 . Then, there exists an edge in S for each pair of adjacent compartments whose labels belong to V . That is to say: S = { (i 1 , j 1 ) , (i 2 , j 2 ) } , 1 ≤i 1 , i 2 ≤length, 1 ≤j 1 , j 2 ≤width , | i 1 −i 2 | + | j 1 −j 2 | = 1 •Alphabet is composed of objects representing signals, markers and clocks. The meaning of bclkmax and oclkmax is explained in Table 2 . = { fet m,n , rep m,n , rec m,n , 1 ≤m ≤length, 1 ≤n ≤width } ∪ { l ss i 1 ,i 2 , l ssb i 1 ,i 2 , rss i 1 ,i 2 , rssb i 1 ,i 2 , hes i 1 ,i 2 , hesb i 1 ,i 2 , tas i 1 ,i 2 , tasb i 1 ,i 2 , 1 ≤i 1 ≤max (lengt h, widt h ) , 0 ≤i 2 ≤max (lengt h, widt h ) } ∪ { les, res, less, ress, lesb, resb} ∪ { he, ta, le, re, lb, rb, phar, pharsl, phars v } ∪ { bclk bclki , 1 ≤bclki ≤bclkmax } ∪ { sclk sclki , 1 ≤sclki ≤4 } ∪ { oclk oclki , 1 ≤oclki ≤oclkmax } −Objects fet m, n , rep m, n and rec m, n , 1 ≤m ≤length , 1 ≤n ≤width , represent fetch, reply and received signals which assist each compartment in knowing to which morphology part it belongs. −Objects l ss i 1 ,i 2 , l ssb i 1 ,i 2 , rss i 1 ,i 2 , rssb i 1 ,i 2 , hes i 1 ,i 2 , hesb i 1 ,i 2 , t as i 1 ,i 2 and t asb i 1 ,i 2 , 1 ≤i 1 ≤max (lengt h, widt h ) , 0 ≤i 2 ≤ max (lengt h, widt h ) determine the location of the pharynx in the body. Specifically, the model tries to allocate a single pharynx in the center of the trunk. This task is accomplished by sending a signal towards the left-hand side ( l ss i 1 ,i 2 and l ssb i 1 ,i 2 ), right-hand side r ss i 1 ,i 2 and r ssb i 1 ,i 2 , head ( hes i 1 ,i 2 and hesb i 1 ,i 2 ) and tail ( t as i 1 ,i 2 and t asb i 1 ,i 2 ), so as to measure the distance from each compartment to each one of these regions. When these signals bounce back and the differences in longitudinal (head to tail) and ventral (left to right) distances are off by at most one unit towards the left-hand side and the tail, then the model allocates the pharynx in that compartment. −Objects les, less, lesb, res, ress and resb are signals used to allocate the left ( les, less and lesb ), right ( res, ress and resb ) eyes and brains and left ( lvncs, lvncss and lvncsb ) and right ( rvncs, rvncss and rvncsb ) ventral nerve cords. To do so, M. García-Quismondo et al. / Information Sciences 381 (2017) 229–249 235 Fig. 2. Initial (left) and final (right) morphologies of a model of tail regeneration. The top diagrams are cartoonlike illustrations of the morphologies and the bottom diagrams are visual representations of their P system configurations. each compartment in the head sends a signal to the left ( les and less ) and right ( res and ress ). If there is no response ( lesb and resb , respectively) in the next step, then the model allocates a left eye and brain. The same applies for the right eye and brain with objects res, ress and resb . −Objects he and ta are markers used to determine which compartments belong to the head and tail. Objects le and lb represent the left eye and brain and re and rbrvnc represent their right counterparts. Finally, object phar represents the pharynx and objects pharsl and pharsv are auxiliary signals to allocate it. −Objects bclk bclki , 1 ≤bclki ≤bclkmax , sclk sclki , 1 ≤sclki ≤5 and oclk oclki , 1 ≤oclki ≤oclkmax are clocks which assist on determining to which morphology part each compartment belongs ( bclk bclki and sclk sclki ) and in which compartment each organ is ( oclk oclki ). •Alphabet is composed of scripts representing computable predicates over configurations from : = { sc 1 , σ} . The predicate τ1 associated to sc 1 is computed as follows. At the beginning of the computation, a variable conf storing the latest configuration is set to a null value. Then, when sc 1 is generated, if conf is set to null then τ1 evaluates to false. Otherwise, τ1 evaluates to true at a certain instant if the configuration of the system at that instant is equal to the value in conf . We say that two configurations C 1 and C 2 are equal if both of the following holds: −For each compartment i in C 1 there exists a single compartment i  in C 2 with the same label and associated multiset than i . Likewise, for each compartment i in C 2 there exists a single compartment i  in C 1 with the same label and associated multiset than i . −For each two compartments i and j which are linked in C 1 , there exist two linked compartments i  , j  in C 2 such as i  has the same label as i and j  has the same label as j . Likewise, for each two compartments i and j which are linked in C 2 , there exist two linked compartments i  , j  in C 1 such as i  has the same label as i and j  has the same label as j . •The initial multisets are: −M i,j = { bclk 1 } , (i, j) ∈ V, that is, each compartment contains a region clock set at time 1. •The rule set R is composed of the following rules: −First, the halting condition is checked. To do so, the model generates a script which triggers the halting predicate. Note that this is only necessary for a single compartment. (2)[ bclk 1 → fet 1 , 1 , bclk 2 , sc 1 ] 2 , 1 236 M. García-Quismondo et al. / Information Sciences 381 (2017) 229–249 If (2, 1) ∈ V , then such a rule is replaced by: (2)[ bclk 1 → fet i −1 ,j , bclk 2 , sc 1 ] i,j provided that ( i, j ) ∈ V . −Next, the model associates each compartment with a morphology part. This procedure consists of four stages, each one for a region marker. Each compartment sends a signal towards the direction in which the region must be. For instance, in the left side stage each compartment sends a signal fet m, n towards the left-hand side. If it does not receive a signal rep m, n in the next step, it means that there exists no compartment at its left. Therefore, the compartment generates a marker ls to allocate itself in the left-hand side. If the compartment does not have the minimum possible index (1), it must divide towards the left and pass on its left side marker. If there exists any compartment to the right of the newly created compartment, then these compartments are linked. The following rules model this procedure for the leftside (1), head (2), right side (3), and tail (4). lp is experimentally tuned so that each stage consists of enough steps to allocate a region. • (1) [ bclk 1 → fet m,n −1 , bclk 2 ] m,n [ fet m,n −1 ] m,n → [ rep m,n ] m,n −1 (−1) [ r ep m,n −1 ] m,n → [ r ec m,n ] m,n −1 { g 1 , 1 } (−1) [ → ls ] m,n { g 1 , 2 } [ bcl k lp−1 , l s ] m,n ⇒ [ bcl k lp ] m,n [ scl k 1 , l s ] m,n −1 { g 3 } [ ] m,n ⇒ [ ] m,n −1 [ sclk 1 → bclk lp+1 ] m,n −1 1 ≤m ≤length 2 ≤n ≤width where g 1 , 1 ≡bclk 1 ∧ sclk 1 ∧ rec m,n −1 ∧ ls ∧ bclk lp−2 , g 1 , 2 ≡bclk 1 ∧ sclk 1 ∧ rec m,n −1 and g 3 ≡ 4 h =1 bclk h . • (1) [ bclk lp+2 → fet m −1 ,n , bclk lp+3 ] m,n [ fet m −1 ,n ] m,n → [ rep m,n ] m −1 ,n (−1) [ r ep m −1 ,n ] m,n → [ r ec m,n ] m −1 ,n { g 2 , 1 } (−1) [ → ta ] m,n { g 2 , 2 } [ bclk 2 ·lp−1 , ta ] m,n ⇒ [ bcl k 2 ·lp ] m,n [ scl k 2 , ta ] m −1 ,n { g 3 } [ ] m,n ⇒ [ ] m −1 ,n [ sclk 2 → bclk 2 ·lp+1 ] m −1 ,n 2 ≤m ≤length 1 ≤n ≤width where g 2 , 1 ≡bclk 1 ∧ sclk 1 ∧ rec m,n −1 ∧ ta ∧ bclk 2 ·lp−2 and g 2 , 2 ≡bclk 1 ∧ sclk 1 ∧ rec m −1 ,n . • (1) [ bclk 2 ·lp+2 → fet m,n +1 , bclk 2 ·lp+3 ] m,n [ fet m,n +1 ] m,n → [ rep m,n ] m,n +1 (−1) [ r ep m,n +1 ] m,n → [ r ec m,n ] m,n +1 { g 3 , 1 } (−1) [ → rs ] m,n { g 3 , 2 } [ bclk 3 ·lp−1 , rs ] m,n ⇒ [ bcl k 3 ·lp ] m,n [ scl k 3 , rs ] m,n +1 { g 3 } [ ] m,n ⇒ [ ] m,n +1 [ sclk 3 → bclk 3 ·lp+1 ] m,n +1 1 ≤m ≤length 1 ≤n ≤width −1 where g 3 , 1 ≡bclk 1 ∧ sclk 1 ∧ rec m,n −1 ∧ rs ∧ bclk 3 ·lp−2 and g 3 , 2 ≡bclk 1 ∧ sclk 1 ∧ rec m,n +1 . • (1) [ bclk 3 ·lp+2 → fet m +1 ,n , bclk 3 ·lp+3 ] m,n [ fet m +1 ,n ] m,n → [ rep m,n ] m +1 ,n (−1) [ r ep m +1 ,n ] m,n → [ r ec m,n ] m +1 ,n { g 4 , 1 } (−1) [ → he ] m,n { g 3 , 2 } [ bclk 4 ·lp−1 , rs ] m,n ⇒ [ bcl k 4 ·lp ] m,n [ scl k 4 , he ] m +1 ,n { g 3 } [ ] m,n ⇒ [ ] m +1 ,n [ sclk 4 → bclk 4 ·lp+1 ] m +1 ,n 1 ≤m ≤length −1 1 ≤n ≤width where g 4 , 1 ≡bclk 1 ∧ sclk 1 ∧ rec m,n −1 ∧ he ∧ bclk 4 ·lp−2 and g 4 , 2 ≡bclk 1 ∧ sclk 1 ∧ rec m +1 ,n . −In case that none of the above rules have been applied at a certain step, there exists another auxiliary set of rules which increase the morphology part clock in each compartment. M. García-Quismondo et al. / Information Sciences 381 (2017) 229–249 237 • (−1) [ bclk i → bclk i +1 ] m,n ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 1 ≤m ≤length 1 ≤n ≤width 1 ≤i ≤bclkmax −1 i  = 2 i  = lp + 3 i  = 2 ·lp + 3 i  = 3 ·lp + 3 • [ bclk 2+ j → bclk 3+ j ] m,n [ bclk i + lp+3 → bclk i + lp+4 ] m,n [ bclk i +2 ·lp+3 → bclk i +2 ·lp+4 ] m,n [ bclk i +3 ·lp+3 → bclk i +3 ·lp+4 ] m,n [ bclk bclkmax → oclk 1 ] m,n  1 ≤m ≤length 1 ≤n ≤width 0 ≤j ≤1 −After the morphology part allocation stage has been completed, the objects rec m, n that remain unused are cleared. (1) [ rec m,n −1 , clk cl kit·l p−2 → c lk cl kit·l p−1 ] m,n (1) [ rec m −1 ,n , clk cl kit·l p−2 → c lk cl kit·l p−1 ] m,n (1) [ rec m,n +1 , clk cl kit·l p−2 → c lk cl kit·l p−1 ] m,n (1) [ rec m +1 ,n , clk cl kit·l p−2 → c lk cl kit·l p−1 ] m,n ⎧ ⎨ ⎩ 1 ≤m ≤length 1 ≤n ≤width 1 ≤clkit ≤4 −Our model does not permit the same compartment to simultaneously belong to the head and the tail regions, nor to the left-hand and the right-hand side. Therefore, if any of these pairs of object markers coincides in a compartment at the same time, they are cleared out. Cells often implement bi-stable negative feedback signaling to ensure that only one outcome dominates. The purpose of these rules is to avoid incorrect behaviors in the model resulting from having the same compartment associated with opposing parts. (−1) [ he, ta → ] m,n (−1) [ ls, rs → ] m,n 1 ≤m ≤length 1 ≤n ≤width −These rules initiate and locate the pharynx in the center of the trunk. The model considers that all compartments that do not belong either to the head and the tail are automatically part of the trunk. Therefore, head and tail compartments are excluded from the pharynx allocation. Likewise, compartments which belong either to the left-hand or the right-hand side are excluded from being assigned as ’center’. { he ∧ ta ∧ ls ∧ rs ∧ oclk 1 } [ → lss n, 0 , rss n, 0 , hes m, 0 , tas m, 0 ] m,n { he ∧ ta ∧ ls ∧ rs } [ pharsl, phars v → phar] m,n 2 ≤m ≤length −1 2 ≤n ≤width −1 −After all compartments have been associated with a morphology part, but before the organs are allocated, the model flushes organ markers. { oclk 1 } [ phar → ] m,n { oclk 1 } [ e → ] m,n { oclk 1 } [ b → ] m,n 1 ≤m ≤length 1 ≤n ≤width −Compartments that have been allocated to the head send a signal to the left and right to allocate the eyes and brains. If they do not receive any response from the compartment on the left in the next step, then the compartment is at the leftmost position in the head. Consequently, it generates an eye and a brain. The same applies for the right-hand side. { he ∧ oclk 1 } [ → les, res ] m,n { he } [ l es ] m,n → [ l ess ] m,n −1 { he } [ r es ] m,n → [ r ess ] m,n +1 { he } [ l ess ] m,n → [ l esb] m,n +1 { he } [ r ess ] m,n → [ r esb] m,n −1 { he ∧ oclk oclkmax } [ lesb → ] m,n { he ∧ oclk oclkmax } [ resb → ] m,n { he ∧ oclk oclkmax ∧ lesb } [ → e, b] m,n { he ∧ oclk oclkmax ∧ resb } [ → e, b] m,n 2 ≤m ≤length 1 ≤n ≤width −The organ clock advances one unit per computation step. 244 M. García-Quismondo et al. / Information Sciences 381 (2017) 229–249 Fig. 8. Initial (left) and final (right) morphologies of an experiment of regeneration of a posterior head in a worm exposed to octanol. In response to the drug, specific rules defined in the model are altered, resulting in the regeneration of a complete posterior head and the restoration of the original double–head morphology in which it was amputated. tablished by comparing the sequence indexes between the wounded region and the rest of the remaining body. During the decrementation stage, the regions corresponding to the missing indexes are sequentially regenerated. Although this model keeps track of the head–tail polarity, territories are not associated with the head, trunk or tail. Instead, they are placed in a sequence between 0 (tail) and a predefined number n (head). In contrast to ours, Slack’s model is unidimensional, modeling the flatworm morphology as sequential segments along the anterior–posterior axis. Meinhardt [48] proposed a regeneration framework based on dynamical systems of diffusing substances with local self– activation and lateral inhibition. Distinction between regions is regulated by a set of interacting products. In this way, the concentration levels of different products dictate the location of different regions, delimiting boundaries among them. Region allocation is more expressive in this model than in ours because the association of regions with body parts is implemented in a continuous manner. A real value is assigned to the membership of a body part to a region. Even if our current model does not implement this type of numerical tagging, it could be implemented in our framework, since objects inside compartments can have a cardinality greater than one. Thus, the membership of a compartment in a body part could be encoded as the cardinality of the body part marker in that compartment. Recently, a comprehensive dynamic genetic regulatory model of planarian regeneration was reverse–engineered directly from experimental data, including amputations, genetic knock–downs, and drugs [35,36] . Unlike Slack’s, this model can explicitly delimit the head, trunk and tail regions in the morphology. The authors used a machine learning algorithm to automatically generate the model from a dataset of experiments and validated its predictions by performing new experiments in the laboratory [39] . However, like Slack’s, this model does not implement organ allocation, a feature included in our model presented here. Finally, it is worth mentioning a cellular aggregation model based on P systems [44] ; in this work, the authors proposed a P system model to capture cellular aggregation in Dictyostelium discoideum . This organism alternates between two morphologies: a mass of cells spread in the environment that aggregates into a moving slug in response to specific environ- M. García-Quismondo et al. / Information Sciences 381 (2017) 229–249 245 mental conditions. The authors defined rules for biomechanical motion that were based on the proximity across cells and attraction/repulsion dynamics. However, since this work captures the behavior of a different biological process, we have not included a detailed comparison of this work with our model. After this comparative analysis, we conclude that one of the new features provided by our model is the allocation of organs inside body parts. This feature is absent in all the previously proposed frameworks that we have analyzed. Organs are easily abstracted in our model due to the hierarchical structure of P systems, in which compartments may contain objects and/or, as it is the case in other membrane computing frameworks [13] , other compartments. This gives place to a hierarchy of potentially infinite levels and facilitates the modeling of hierarchical structures such as organs and other body parts. 5.2. Validation To test our model, we reproduced five different experiments from the literature. For each of these experiments, our model reproduced the resulting morphologies observed in the in vivo experiments. Our model was able to regenerate an entire morphology even in the absence of a tail–head dipole. Fig. 5 displays a representative example of this feature. At the beginning of this simulation, there exists no compartment signaling the head nor the tail. However, the model generates a head and a tail at opposite sides of the morphology. Region planar proportions are also regenerated correctly by our model. For instance, in the experiment from Fig. 4 two new columns of compartments are created, yielding a morphology with a ventral size of three columns. These proportions are equivalent to those of a wild–type flatworm. Our model also succeeds in generating new organs and placing them in the correct region of the body. For example, the experiment from Fig. 4 correctly generates a new brain and a new eye in the right side of the head, and a new pharynx in the trunk. This mechanism works even when the final morphology differs from a wild–type flatworm. As an example, the experiment from Fig. 6 generates a double–head worm. In this morphology, each head slice contains two eyes and two brains, giving a final morphology with four brains and four eyes. To test our model further, we included two additional experiments. In these experiments, the original morphology is not the wild–type, but a double–head worm ( Figs. 7 and 8 ). For the first experiment, a small notch from the posterior head is removed. The idea of this experiment is to test the model predictions when the morphology before amputation is not wild–type. Interestingly, the same model that can regenerate the wild–type morphology now accurately predicts that the ectopic posterior head is regenerated, restoring the double–head morphology. This experiment proves that the model can yield correct predictions even for unconventional morphologies. In the second experiment, we considered a case in which the whole posterior head of a double–head worm is amputated. A set of experiments from the Planform database in which octanol was applied during regeneration [54] suggests that the model should regenerate a new ectopic head. To account for the presence of octanol in the media in this experiment, we altered a rule module in the model. This is necessary because the compartment structure retains no information about amputated parts. Thus, it is not possible to encode this information in the compartments themselves. Without this change, the model would have regenerated a posterior tail instead of a head, yielding a wild–type morphology. But with these new rules, the model regenerates the posterior head, explaining persistence of a double–head phenotype after amputation [54] . To the best of our knowledge, this feature is not explained by any other current computational model. 5.3. Future work Despite the success of our model in reproducing these experiments, there are still some important aspects that have not been adequately captured. For instance, it has been shown that chemical compounds [47,51] , radiation [9] and electromagnetic fields [8,10,52] influence regeneration by targeting specific endogenous components. These features, such as proliferating stem cells and bioelectrical signaling among cells, are absent in our model. We plan to incorporate these aspects in future work. The fundamental notion of P systems, encoding information on membranes and passing signals among subunits, are a good match for the known biology underlying regeneration. We believe that it is possible to extend this formalism to include both biochemical and biophysical signaling mechanisms. The current version of our model has several other limitations. For example, region sizes and adjacencies are respectively encoded on the number of membranes associated with each region and the links between membranes that belong to different regions. However, other relevant information, such as angles between regions and relative organ orientations [38] are not represented. The integration of these features into our model will be performed in the future. Another feature to integrate in this model is the ability to represent permanent changes to the target morphology. Recent experimental data have shown that bioelectrical connectivity among cells can be altered, changing a flatworm to a doublehead worm permanently. Amputated fragments from this worm can regenerate double–head worms [54] . Such editing of the target morphology represents an important feature for models of regeneration, since they force explicit statements of how complex target shape is encoded, or emerges, from the properties of cells or cell networks [40] . Yet another line of future work is the automatic design of regenerative models based on membrane computing. The current version of our model was designed manually. This is obviously impractical for large sets of experiments. A future line of work will include the development of machine learning algorithms to yield regeneration models of planarian flatworms 246 M. García-Quismondo et al. / Information Sciences 381 (2017) 229–249 from large experimental repositories [35,37,38] . Not only would these models be able to predict the resulting morphology from a new experiment, but they would also help elucidate the internal workings of regenerative processes in flatworms at a cellular level. This is possible because of the tremendous expressiveness of membrane computing. In particular, its semantics is based on rewriting rules that can be bundled into modules with associated functions [13,43,66] . In our model, these modules represent specific aspects of regenerative processes in planarian flatworms. In a framework for the automatic design of models based on membrane computing, the modules would behave as white boxes. The application of rewriting rules could be interpreted as bioelectrical and biochemical processes. Such modularity is a very attractive feature, as it is known to facilitate the development and evolvability of complex forms [30,41,42,44,71,77,78] . The expressive power of this modularity is also corroborated in this work. In the experiment shown in Fig. 8 , the desired behavior of the model could be adjusted by simply switching a module. The easiness with which the rules can be modified in order to achieve specific patterning results is one of the major assets of membrane computing as a modeling framework for intracellular [12,25,67] and inter–cellular [66] phenomena. This class of models facilitates the linkage between high–level patterning outcomes and mechanism–level rules. This is a challenging inverse problem [40] that these systems handle well, especially in comparison with more traditional frameworks such as differential equations [12] . 6. Conclusions In this work, we presented a new application of the rich computational formalism of membrane computing for describing models of regenerative biology. P systems are an excellent framework for modeling the signaling mechanisms in regenerative organisms due to their inherent capacity to model morphological regions as membranes and their algorithmic (generative) nature to model dynamic shape restoration. Unlike molecular pathway networks, P system models can be easily designed to recapitulate large-scale patterning behaviors. This facilitates the development of testable models without needing machine learning methods to search and optimize parameters. Based on this novel formalism, we presented as a proof-of-concept a novel model of planarian regeneration able to recapitulate several aspects of the extraordinary ability of these organisms to restore their regional configuration and organ locations after large-scale amputations. Future work will extend the complexity of this model, derive novel specific predictions from it, and validate them at the bench. Novel applications of computer science formalisms to the biological problem of pattern formation and control, such as the work presented here, is essential for leveraging the continuous progress and experimental results from molecular biology experiments and translating them into a mechanistic understanding leading to biomedical advances. Acknowledgments Manuel García-Quismondo acknowledges the support of the National FPU Grant Programme from the Spanish Ministry of Education. Michael Levin gratefully acknowledges an Allen Discovery Center award from The Paul G. Allen Frontiers Group, and support of the G. Harold and Leila Y. Mathers Charitable Foundation. References [1] J.A. Baddour, K. Sousounis, P.A. Tsonis, Organ repair and regeneration: an overview, Birth Defects Res. Part C: Embryo Today: Rev. 96 (1) (2012) 1–29, doi: 10.1002/bdrc.21006 . [2] N. Barkai, D. Ben-Zvi, Big frog, small frog maintaining proportions in embryonic development, FEBS J. 276 (5) (2009) 1196–1207, doi: 10.1111/j. 1742-4658.2008.06854.x . [3] D. Basanta, M. Miodownik, B. Baum, The evolution of robust development and homeostasis in artificial organisms, PLoS Comput. Biol. 4 (3) (2008) e10 0 0 030, doi: 10.1371/journal.pcbi.10 0 0 030 . [4] W.S. Beane, J. Morokuma, J.M. Lemire, M. Levin, Bioelectric signaling regulates head and organ size during planarian regeneration, Development 140 (2) (2013) 313–322, doi: 10.1242/dev.086900 . [5] N. Bessonov, M. Levin, N. Morozova, N. Reinberg, A. Tosenberger, V. Volpert, On a model of pattern regeneration based on cell memory, PLoS One 10 (2) (2015) e0118091, doi: 10.1371/journal.pone.0118091 . [6] K.D. Birnbaum, A.S. Alvarado, Slicing across kingdoms: regeneration in plants and animals, Cell 132 (2008) 697–710, doi: 10.1016/j.cell.2008.01.040 . [7] J.P. Brockes, A. Kumar, Appendage regeneration in adult vertebrates and implications for regenerative medicine, Science 310 (5756) (2005) 1919–1923, doi: 10.1126/science.1115200 . [8] F.A. Brown , Effects and after-effects on planarians of reversals of the horizontal magnetic vector, Nature 209 (5022) (1966) 533–535 . [9] F.A. Brown, Y.H. Park, Seasonal variations in sign and strength of Gamma-Taxis in Planarians, Nature 202 (4931) (1964) 469–471, doi: 10.1038/202469a0 . [10] H.M. Brown , T.E. Ogden , The electrical response of the planarian ocellus, J. Gen. Physiol. 51 (2) (1968) 237–253 . [11] O. Chara, E.M. Tanaka, L. Brusch, Mathematical modeling of regenerative processes, Curr. Top. Dev. Biol. 108 (2014) 283–317, doi: 10.1016/ b978- 0- 12- 391498- 9.0 0 011-5 . [12] S. Cheruku , A. P ˘ aun , F.J. Romero-Campero , M.J. Pérez-Jiménez , O.H. Ibarra , Simulating FAS-induced apoptosis by using P systems, Prog. Nat. Sci. 17 (2007) 424–431 . [13] M.A . Colomer, A . Margalida, M.J. Pérez-Jiménez, Population dynamics P system (PDP) models: a standardized protocol for describing and applying novel bio-inspired computing tools, PLoS One 8 (4) (2013) 1–13, doi: 10.1371/journal.pone.0060698 . [14] F. Durant, D. Lobo, J. Hammelman, M. Levin, Physiological controls of large-scale patterning in planarian regeneration: a molecular and computational perspective on growth and form, Regeneration 3 (2) (2016) 78–102, doi: 10.1002/reg2.54 . [15] M. Emmons-Bell, F. Durant, J. Hammelman, N. Bessonov, V. Volpert, J. Morokuma, K. Pinet, D. Adams, A. Pietak, D. Lobo, M. Levin, Gap junctional blockade stochastically induces different species-Specific head anatomies in genetically wild-Type girardia dorotocephala flatworms, Int. J. Mol. Sci. 16 (12) (2015) 27865–27896, doi: 10.3390/ijms161126065 . [16] G. Fox, S. Pallickara, X. Rao, Towards enabling peer-to-peer grids, Concurr. Comput.: Pract. Exp. 17 (7–8) (2005) 1109–1131, doi: 10.1002/cpe.863 . [17] M. García-Quismondo , R. Gutiérrez-Escudero , M.A. Martínez-del Amor , E.F. Orejuela-Pinedo , I. Pérez-Hurtado , P-Lingua 2.0: a software framework for cell-like p systems, Int. J. Comput. Commun. Control 4 (3) (2010) 234–243 . M. García-Quismondo et al. / Information Sciences 381 (2017) 229–249 247 [18] M. García-Quismondo, M. Martínez-del Amor, M. Pérez-Jiménez, Probabilistic guarded P systems, a new formal modelling framework, in: M. Gheorghe, G. Rozenberg, A. Salomaa, P. SosÃk, C. Zandron (Eds.), Membrane Computing, Lecture Notes in Computer Science, 8961, Springer International Publishing, 2014, pp. 194–214, doi: 10.1007/978- 3- 319- 14370- 5 _ 12 . [19] M. García-Quismondo, J.M. Reed, F.S. Chew, M.A. Martínez-del Amor, M.J. Pérez-Jiménez, Evolutionary response of a native butterfly to concurrent plant 2 invasions: simulation of population dynamics, PLoS Comput. Biol. Unpublished results . [20] M. Gheorghe, F. Ipate, C. Dragomir, A kernel P system, in: Proceedings of the Tenth Brainstorming Week on Membrane Computing, I, Fenix editora, 2012, pp. 153–170, doi: 10.1007/978- 3- 642- 11467- 0 _ 20 . [21] M. Gheorghe, F. Ipate, C. Dragomir, L. Mierla, L. Valencia-Cabrera, M. García-Quismondo, M.J. Pérez-Jiménez, Kernel P systems –version I, Proceedings of the Eleventh Brainstorming Week on Membrane Computing(2013a) 97–124. [22] M. Gheorghe, F. Ipate, R. Lefticaru, M.J. Pérez-Jiménez, A. Turcanu, L. Valencia Cabrera, M. García-Quismondo, L. Mierla, 3-col problem modelling using simple kernel P systems, Int. J. Comput. Math. 90 (4) (2013b) 816–830, doi: 10.1080/00207160.2012.743712 . [23] L.V. Goodrich, D. Strutt, Principles of planar polarity in animal development, Development 138 (10) (2011) 1877–1892, doi: 10.1242/dev.054080 . [24] R.S. Gray, I. Roszko, L. Solnica-Krezel, Planar cell polarity: coordinating morphogenetic cell behaviors with embryonic polarity, Dev. Cell 21 (1) (2011) 120–133, doi: 10.1016/j.devcel.2011.06.011 . [25] T. Hinze , S. Hayat , T. Lenser , N. Matsumaru , P. Dittrich , Hill kinetics meets P systems: a case study on gene regulatory networks as computing agents in silico and in vivo, in: Proceedings of the 8th International Conference on Membrane computing, WMC’07, Springer-Verlag, Berlin, Heidelberg, 2007, pp. 320–335 . [26] D.E. Ingber, M. Levin, What lies at the interface of regenerative medicine and developmental biology? Development 134 (14) (2007) 2541–2547, doi: 10.1242/dev.003707 . [27] M. Ionescu , G. P ˘ aun , T. Yokomori , Spiking neural P systems, Fundam. Inform. 71 (2,3) (2006) 279–308 . [28] F. Ipate, R. Lefticaru, L. Mierla, L. Valencia-Cabrera, H. Han, G. Zhang, C. Dragomir, M.J. Pérez-Jiménez, M. Gheorghe, Kernel P systems: Applications and implementations, in: Z. Yin, L. Pan, X. Fang (Eds.), Proceedings of the Eighth International Conference on Bio-Inspired Computing: Theories and Applications (BIC-TA), 2013, Advances in Intelligent Systems and Computing, 212, Springer, Berlin, Heidelberg, 2013, pp. 1081–1089, doi: 10.1007/978- 3- 642- 37502- 6 _ 126 . [29] J. Jaeger, J. Sharpe, On the Concept of Mechanism in Development, Oxford University Press, pp. 56–78, 10.1093/acprof:oso/9780199671427.001.0001 [30] C.P. Klingenberg, Studying morphological integration and modularity at multiple levels: concepts and analysis, Philos. Trans. R. Soc. Lond. B: Biol. Sci. 369 (1649) (2014), doi: 10.1098/rstb.2013.0249 . [31] R. Lander, C.P. Petersen, Wnt, Ptk7, and FGFRL expression gradients control trunk positional identity in planarian regeneration., eLife 5 (2016), doi: 10. 7554/eLife.12850 . [32] M. Levin, The wisdom of the body: future techniques and approaches to morphogenetic fields in regenerative medicine, developmental biology and cancer, Regen. Med. 6 (6) (2011) 667–673, doi: 10.2217/rme.11.69 . [33] M. Levin, Morphogenetic fields in embryogenesis, regeneration, and cancer: non-local control of complex patterning, BioSystems 109 (3) (2012) 243– 261, doi: 10.1016/j.biosystems.2012.04.005 . [34] D. Lobo, W.S. Beane, M. Levin, Modeling planarian regeneration: a primer for reverse-engineering the worm, PLoS Comput. Biol. 8 (4) (2012) e1002481, doi: 10.1371/journal.pcbi.1002481 . [35] D. Lobo, M. Levin, Inferring regulatory networks from experimental morphological phenotypes: a computational method reverse-engineers planarian regeneration, PLoS Comput. Biol. 11 (6) (2015) e1004295, doi: 10.1371/journal.pcbi.1004295 . [36] D. Lobo, M. Levin, Computing a Worm: Reverse-Engineering Planarian Regeneration, vol. 2, Springer International Publishing, p. 637–654, 10.1007/978- 3-319-33921-4_2 [37] D. Lobo, T.J. Malone, M. Levin, Planform: an application and database of graph-encoded planarian regenerative experiments, Bioinformatics 29 (8) (2013) 1098–1100, doi: 10.1093/bioinformatics/btt088 . [38] D. Lobo, T.J. Malone, M. Levin, Towards a bioinformatics of patterning: a computational approach to understanding regulative morphogenesis, Biol. Open 2 (2) (2013) 156–169, doi: 10.1242/bio.20123400 . [39] D. Lobo, J. Morokuma, M. Levin, Computational discovery and in vivo validation of hnf4 as a regulatory gene in planarian regeneration., Bioinformatics (2016), doi: 10.1093/bioinformatics/btw299 . [40] D. Lobo, M. Solano, G.A. Bubenik, M. Levin, A linear-encoding model explains the variability of the target morphology in regeneration, J. R. Soc. Interface 11 (92) (2014), doi: 10.1098/rsif.2013.0918 . [41] D. Lobo, F.J. Vico, Evolution of form and function in a model of differentiated multicellular or ganisms with gene regulatory networks, BioSystems 102 (2–3) (2010a) 112–123. http://dx.doi.org/10.1016/j.biosystems.2010.08.003 . [42] D. Lobo, F.J. Vico, Evolutionary development of tensegrity structures, BioSystems 101 (3) (2010b) 167–176. http://dx.doi.org/10.1016/j.biosystems.2010. 06.005 . [43] D. Lobo, F.J. Vico, J. Dassow, Graph grammars with string-regulated rewriting, Theor. Comput. Sci. 412 (43) (2011) 6101–6111. http://dx.doi.org/10.1016/ j.tcs.2011.07.004 . [44] V. Manca, G. Pardini, Morphogenesis through moving membranes, Nat. Comput. 13 (3) (2014) 403–419, doi: 10.1007/s11047- 013- 9407- 4 . [45] C. Martín-Vide, G. Paun, J. Pazos, A. Rodríguez-Patón, Tissue P systems, Theor. Comput. Sci. 296 (2) (2003) 295–326, doi: 10.1016/S0304-3975(02) 00659-X . [46] C. Martín-Vide, J. Pazos, G. P ˘ aun, R.-P. A., A new class of symbolic abstract neural nets: tissue P systems, in: O. Ibarra, L. Zhang (Eds.), Computing and Combinatorics, Lecture Notes in Computer Science, 2387, Springer, Berlin, Heidelberg, 2002, pp. 290–299, doi: 10.1007/3- 540- 45655- 4 _ 32 . [47] P. Mason, Chemo-klino-kinesis in planarian food location, Anim, Behav, 23 (1975) 460–469. http://dx.doi.org/10.1016/0 0 03-3472(75)90 095-0 . [48] H. Meinhardt, Models for the generation and interpretation of gradients, Cold Spring Harb. Perspect. Biol. 1 (4) (20 09) a0 01362, doi: 10.1101/ cshperspect.a0 01362 . 20 066097[pmid] [49] H. Meinhardt, A. Gierer, Pattern formation by local self-activation and lateral inhibition, Bioessays 22 (8) (20 0 0) 753–760, doi: 10.1002/ 1521-1878(20 0 0 08)22:8  753::AID-BIES9  3.0.CO;2-Z . [50] J. Miller, Evolving developmental programs for adaptation, morphogenesis, and self-repair, in: W. Banzhaf, J. Ziegler, T. Christaller, P. Dittrich, J. Kim (Eds.), Advances in Artificial Life, Lecture Notes in Computer Science, 2801, Springer, Berlin, Heidelberg, 2003, pp. 256–265, doi: 10.1007/ 978- 3- 540- 39432- 7 _ 28 . [51] S. Miyamoto , A. Shimozawa , Chemotaxis in the freshwater planarian, dugesia japonica japonica (behavior biology and ecology), Zool. Sci. 2 (3) (1985) 389–395 . [52] N.J. Murugan, L.M. Karbowski, R.M. Lafrenie, M.A. Persinger, Temporally-patterned magnetic fields induce complete fragmentation in planaria, PLoS One 8 (4) (2013) e61714, doi: 10.1371/journal.pone.0061714 . [53] T. Nogi, M. Levin, Characterization of innexin gene expression and functional roles of gap-junctional communication in planarian regeneration, Dev. Biol. 287 (2) (2005) 314–335, doi: 10.1016/j.ydbio.2005.09.002 . [54] N.J. Oviedo, J. Morokuma, P. Walentek, I.P. Kema, M.B. Gu, J.-M. Ahn, J.S. Hwang, T. Gojobori, M. Levin, Long-range neural and gap junction proteinmediated cues control polarity during planarian regeneration, Dev. Biol. 339 (1) (2010) 188–199. http://dx.doi.org/10.1016/j.ydbio.2009.12.012 . [55] L. Pan, M.J. Pérez-Jiménez, Computational complexity of tissue-like p systems, J. Complex. 26 (3) (2010) 296–315. http://dx.doi.org/10.1016/j.jco.2010. 03.001 . [56] L. Pan, G. P ˘ aun, M.J. Pérez-Jiménez, Spiking neural P systems with neuron division and budding, Sci. China Inf. Sci. 54 (8) (2011) 1596–1607, doi: 10. 1007/s11432-011-4303-y . 248 M. García-Quismondo et al. / Information Sciences 381 (2017) 229–249 [57] H. Peng, J. Wang, M.J. Pérez-Jiménez, A. Riscos-Núñez, An unsupervised learning algorithm for membrane computing, Inf. Sci. 304 (2015) 80–91, doi: 10.1016/j.ins.2015.01.019 . [58] H. Peng, J. Wang, M.J. Pérez-Jiménez, H. Wang, J. Shao, T. Wang, Fuzzy reasoning spiking neural P system for fault diagnosis, Inf. Sci. 235 (2013) 106–116, doi: 10.1016/j.ins.2012.07.015 . [59] I. Pérez-Hurtado, M. García-Quismondo, M.A.M. del Amor, R.M. Gutiérrez = Escudero, E.F. Orejuela-Pinedo, L. Valencia-Cabrera, L.F. Macías-Ramos, P- Lingua webpage, 2009, ( https://www.p-lingua.org ). [60] C.P. Petersen, P.W. Reddien, A wound-induced Wnt expression program controls planarian regeneration polarity, Proc. Natl. Acad. Sci. 106 (40) (2009) 17061–17066, doi: 10.1073/pnas.0906823106 . [61] C.P. Petersen, P.W. Reddien, Polarized notum activation at wounds inhibits Wnt function to promote planarian head regeneration., Science 332 (6031) (2011) 852–855, doi: 10.1126/science.1202143 . [62] G. Pezzulo, M. Levin, Re-membering the body: applications of computational neuroscience to the top-down control of regeneration of limbs and other complex organs, Integr. Biol. 7 (2015) 1487–1517, doi: 10.1039/C5IB00221D . [63] G. P ˘ aun, Computing with membranes, J. Comput. Syst. Sci. 61 (1) (20 0 0) 108–143. http://dx.doi.org/10.1006/jcss.1999.1693 . [64] G. P ˘ aun , M.J. Pérez-Jiménez , A. Riscos-Núñez , Tissue P systems with cell division, Int. J. Comput. Commun. Control 3 (3) (2008) 295–303 . [65] L.V. Rodriguez , R.A. Flickinger , Bipolar head regeneration in Planaria induced by chick embryo extracts, Biol. Bull. 140 (1) (1971) 117–124 . [66] 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 . [67] F.J. Romero-Campero, M.J. Pérez-Jiménez, Modelling gene expression control using P systems: the Lac operon, a case study, BioSystems 91 (3) (2008) 438–457. http://dx.doi.org/10.1016/j.biosystems.2007.02.011 . P-Systems Applications to Systems Biology [68] M. Rubenstein, Y. Sai, C.-M. Choung, W.-M. Shen, Regenerative patterning in swarm robots: mutual benefits of research in robotics and stem cell biology, Int. J. Dev. Biol. 53 (2009) 869–881, doi: 10.1387/ijdb.092937mr . [69] P. Ru, A. Martinez Arias, Cell dynamics and gene expression control in tissue homeostasis and development, Mol. Syst. Biol. 11 (2) (2015) n/a–n/a, doi: 10.15252/msb.20145549 . [70] E. Sal, J.F. Abril, T. Adell, F. Cebri, K. Eckelt, E. Fernndez-Taboada, M. Handberg-Thorsager, M. Iglesias, M.D. Molina, G. Rodríguez-Esteban, Planarian regeneration: achievements and future directions after 20 years of research, Int. J. Dev. Biol. 53 (2009) 1317–1327, doi: 10.1387/ijdb.072414es . [71] G. Schlosser, D. Thieffry, Modularity in development and evolution, Bioessays 22 (11) (20 0 0) 1043–1045, doi: 10.1002/1521-1878(200011)22:11  1043:: AIDBIES11  3.0.CO;2- C . [72] J. Slack, A serial threshold theory of regeneration, J. Theor. Biol. 82 (1) (1980) 105–140. http://dx.doi.org/10.1016/0 022-5193(80)90 092-2 . [73] E. Tanaka, P.W. Reddien, The cellular basis for animal regeneration, Dev. Cell 185 (1) (2011) 172–185, doi: 10.1016/j.devcel.2011.06.016 . [74] A. Tosenberger, N. Bessonov, M. Levin, N. Reinberg, V. Volpert, N. Morozova, A conceptual model of morphogenesis and regeneration, Acta Biotheor. 63 (3) (2015) 283–294, doi: 10.1007/s10441- 015- 9249- 9 . [75] A. Tseng, M. Levin, Cracking the bioelectric code: probing endogenous ionic controls of pattern formation, Commun. Integr. Biol. 6 (1) (2013) e22595, doi: 10.4161/cib.22595 . [76] Y. Umesono, J. Tasaki, Y. Nishimura, M. Hrouda, E. Kawaguchi, S. Yazawa, O. Nishimura, K. Hosoda, T. Inoue, K. Agata, The molecular logic for planarian regeneration along the anterior-posterior axis., Nature 500 (7460) (2013) 73–76, doi: 10.1038/nature12359 . [77] G. von Dassow, E. Munro, Modularity in animal development and evolution: elements of a conceptual framework for evodevo, J. Exp. Zool. 285 (4) (1999) 307–325, doi: 10.1002/(SICI)1097-010X(19991215)285:4  307::AIDJEZ2  3.0.CO;2- V . [78] G.P. Wagner, M. Pavlicev, J.M. Cheverud, The road to modularity, Nat. Rev. Genet. 8 (12) (2007) 921–931, doi: 10.1038/nrg2267 . [79] J. Wang, P. Shi, H. Peng, M.J. Pérez-Jiménez, T. Wang, Weighted {Fuzzy} {Spiking} {Neural} {P} {Systems}, IEEE Trans. Fuzzy Syst. 21 (2) (2013) 209–220, doi: 10.1109/TFUZZ.2012.2208974 . [80] S. Werner, T. Stückemann, M. Beirán Amigo, J.C. Rink, F. Jülicher, B.M. Friedrich, Scaling and regeneration of self-organized patterns, Phys. Rev. Lett. 114 (2015) 138101, doi: 10.1103/PhysRevLett.114.138101 . [81] X. Zeng, L. Xu, X. Liu, L. Pan, On languages generated by spiking neural P systems with weights, Inf. Sci. 278 (2014) 423–433, doi: 10.1016/j.ins.2014. 03.062 . [82] G. Zhang, M. Gheorghe, L. Pan, M.J. Pérez-Jiménez, Evolutionary membrane computing: a comprehensive survey and new results, Inf. Sci. 279 (2014a) 528–551, doi: 10.1016/j.ins.2014.04.007 . [83] X. Zhang, Y. Liu, B. Luo, L. Pan, Computational power of tissue p systems for generating control languages, Inf. Sci. 278 (2014b) 285–297, doi: 10.1016/ j.ins.2014.03.053 . M. García-Quismondo et al. / Information Sciences 381 (2017) 229–249 249 Manuel García-Quismondo is a post-doctoral fellow at Smithsonian Institution (USA) and member of the Research Group in Natural Computing at the University of Sevilla (Spain). His areas of research include membrane computing, population dynamics and high performance computing. He has held fellowships at universities in several countries, including the University of Sevilla (Spain), Manchester Metropolitan University (UK), Huazhong University of Science and Technology (China), Tufts University and University of Minnesota (USA). He has also published about 30 works in topics including neural networks, model checking, system dynamics, ecological modeling, analysis of gene networks and robotics. Michael Levin is the Vannevar Bush Professor at Tufts University, and director of the Tufts Center for Regenerative and Developmental Biology. His group studies information processing by cells and tissues,focusing on the molecular mechanisms by which cell behavior is integrated towards specific anatomical outcomes. Using a combination of molecular biology, developmental genetics, biophysics, and computational modeling, they work to understand endogenous bioelectrical signaling by which cellular networks store and process patterns that guide gene expression and morphogenesis. By uncovering novel controls large-scale remodeling, they hope to advance progress in regenerative medicine, birth defects, cancer biology, and synthetic bioengineering. Daniel Lobo is an Assistant Professor at the University of Maryland, Baltimore County. His research in systems biology aims to understand, control, and design the dynamic regulatory mechanisms governing complex biological processes. To this end, his group develops new computational methods, ontologies, and high-performance in silico experiments to automate the reverseengineering of quantitative models from biological data and the design of regulatory networks for specific functions. They seek to discover the mechanisms of development and regeneration, find therapies for cancer and other diseases, and streamline the application of synthetic biology. His work has received widespread media coverage including PBS, Wired, and Popular Mechanics.