Membrane Computing as a Modelling Tool: Looking Back and Forward from Sevilla

Membrane Computing as a Modelling Tool: Looking Back and Forward from Sevilla Manuel Garc´ıa-Quismondo, Carmen Graciani and Agust´ın Riscos-N´u˜nez Research Group on Natural Computing, Department of Computer Science and Artificial Intelligence, Universidad de Sevilla, Avda. Reina Mercedes s/n, 41012 Sevilla, Spain {mgarciaquismondo,cgdiaz,ariscosn}@us.es Abstract. This paper is a tribute to Prof. Mario de Jes´us P´erez- Jim´enez. An overview of modelling applications in membrane computing has been compiled, trying to narrate it from a historical perspective and including numerous bibliographical references. Since being exhaustive was obviously out of scope, this quick tour on almost two decades of applications is biased, paying special attention to the contributions in which Prof. P´erez-Jim´enez and members of his research group were involved. 1 Introduction Membrane computing (MC) is a computational discipline that takes inspiration from molecular and cellular biology. Therefore, it seems natural to consider membrane computing as a valuable modelling framework for biological processes. In this sense, several authors considered the possibility to “return some meaningful information to Biology” (see Chap. 1 in [36]) from the very initial stages of this discipline. The majority of the initial research lines in membrane computing were devoted to theoretical results investigating the computational power of different classes of P systems: either proving their Turing completeness or exploring their equivalence with elements of classical formal language hierarchies. However, from the onset the curiosity to work towards applications has been also present. The first book fully devoted to applications in membrane computing appeared in 2006 [19], although applications were carried out since 2000, as we will mention later. For a more up-to-date overview of applications of membrane computing, please visit the P systems web page [3] and the Bulletin of the International Membrane Computing Society [1]. Interested readers are also advised to refer to two of the most recent volumes: Applications of Membrane Computing in Systems and Synthetic Biology [26]andReal-life Applications with Membrane Computing [50]. There are also several other overviews of modelling applications, e.g. [30,37,38]. This paper presents a historical (and obviously non-exhaustive) overview of computational modelling approaches within membrane computing, with a special focus on the contributions in which RGNC1members were involved. 1.1 Computational Modelling Computational modelling and simulation are nowadays a cornerstone of the scientific method. Everything starts in Nature, when we identify some physical or biological phenomenon that attracts our interest for some reason. Usually, such intriguing phenomena can be seen as complex systems, in the sense that we can describe them as a collection of elements (let us call them “players”) that interact following relatively simple rules, exhibiting a complex behaviour of the system as a whole. This behaviour is commonly referred to as emergent, and means that the evolution of the system displays some special properties that are not trivially deduced from the local dynamics of its components. The first subjective decision to take when designing a model is to determine which are the relevant ingredients (players, features, variables) that will compose this model. Judging which ingredients are significant enough to be included is a decision linked to the reason why the phenomenon is of interest. More precisely, Regev and Shapiro explain in [43]thatagood model should combine the following key desirable properties: relevant, readable, extensible and computationally tractable. Thus, we have to find the balance between two conflicting goals. On the one hand, we should capture as many ingredients as possible to yield a relevant and realistic model. On the other hand, we should try to keep the model as simple as possible, so that it is easy to interpret and work with. A fundamental stage in the design of a model is validation. To this aim, we need software simulation tools that allow us to run virtual experiments, in order to carry out a reliable analysis of its dynamics under various initial conditions. This is where the concept of practical feasibility comes into play. 2 Historical Overview of Modelling Works in MC We would like to initiate the narrative of the time line on a very significant event: the Workshop on Multiset Processing that was held in Curtea de Arge¸s in 2000. This is commonly accepted as the origin of the series of workshops on membrane computing (which later on evolved and gave rise in 2010 to the International Conference on Membrane Computing). At that time, the limits of the field of membrane computing were still undefined, and the attendants to the meeting came from a variety of backgrounds. 1RGNC stands for Research Group on Natural Computing from Universidad de Sevilla, also known in the MC community as “Sevilla team”. 2.1 In the Beginning, There Was . . . During the above-mentioned workshop, what could be considered as the precursors of future works in the field of modelling in membrane computing were presented. Banˆatre et al. presented a contribution entitled Gamma and the Chemical Reaction Model: Fifteen Years After [7]. Gamma was originally proposed in 1986 as a formalism for the definition of concurrent programs. The basic idea underlying the formalism is to describe computation as a form of chemical reaction on a collection of elementary pieces of data. Indeed, the chemical abstract machine paradigm (CHAM – CHemical Abstract Machine –, introduced by Berry and Boudol in early 1990s) already includes the notion of “membrane” as a container of elements that react among them. Nishida presented a contribution entitled Multiset and K-Subset Transforming Systems [31], where he provided an example illustrating how to model a chaotic discrete dynamical system by means of a K-subset transforming system. He also provided theoretical results proving the expressive power of this kind of systems. Moreover, one year later, he participated again in the Workshop on Membrane Computing, presenting another application case study: a model for the light reactions of the photosynthesis [32]. There are a couple of observations on this paper: on the one hand, he mentioned a computer simulator as a natural addition to the theoretical model. On the other hand, he suggested the idea to go beyond standard multisets, using multiplicities on any semiring K. Last, but not least, Suzuki et al. presented a contribution entitled Artificial Life Applications of a Class of P Systems: Abstract Rewriting Systems on Multisets [45]. We will refer to this, together with some related works, in the next section. 2.2 The Dawn of Brainstorming Era The second significant event that we would like to highlight is the Brainstorming Week on Membrane Computing (BWMC). Its first edition was held in Tarragona (Spain) in February 2003. The Brainstorming Week is a meeting where participants can freely exchange ideas and open problems. It is much more dynamic than a standard workshop or conference, in the sense that the program is proactively set during the gathering and is made up of provocative talks about works in progress that will later on become papers, after fruitful discussions and joint work sessions in a friendly atmosphere. Except for this first edition, the BWMC has been held annually since 2004 in Sevilla (Spain), organised by the RGNC. In the proceedings volume that was produced after the first meeting, we can find several applications in the form of computational models. Suzuki and Tanaka revisit in [47]theirproposaltouseAbstract rewriting systems on multisets as a versatile modelling framework able to capture not only dynamics of chemical reactions (e.g. Brusselator model) and population dynamics in a tri-trophic ecological system (both of them already hinted in [45]), but also other applications in medicine such as inflammatory response. In particular, Suzuki and his collaborators have presented during the Fourth Workshop on Membrane Computing, held later the same year, a model of the p532signalling network [46]. The authors presented a very simple model using a few multiset rewriting rules, allowing symbols to move along two regions: nucleus and cytoplasm. The goal of the paper is to show that the use of P systems as an alternative to traditional rate equation models is useful and practically feasible. The proceedings are available for download via [3], however interested readers may also refer to the chapter “Modeling p53 Signaling Pathways by Using Multiset Processing”, in [19]. One can also find in the Brainstorming proceedings an interesting joint work by Ardelean and Besozzi [4], where the authors propose to pay attention to mechano-sensitive channels from a membrane computing perspective. The idea is further developed in [11,12], where one can find several considerations about the ingredients or features that should be integrated in a P system variant to be considered as a competitive modelling framework. They already mention the idea of adding variable parameters associated with a membrane (in their case related to membrane tension), whose changes are determined by environmental conditions. They point out that this concept is similar to the role of electrical charges in P systems with active membranes [34]. They also acknowledge the advantage of considering rules able to interact with objects at both sides of the membrane, similar to the approach in [9,10]. In particular, for mechano-sensitive channels it is straightforward that the concentration of reactants should be taken into account in the rules, and the authors cite how Suzuki and Tanaka defined their rules as pairs of the type [condition,action]. Another aspect that should be taken into account in order to capture the dynamics of biological processes is the replacement of maximal parallelism by some alternative semantics that express a probabilistic behaviour, and consequently the dependence on software tools to run multiple simulations [5,6,13]. 2.3 Reaching Maturity In 2003, the Institute for Scientific Information (ISI) listed the seminal paper of membrane computing [35]asfast breaking paper. In addition, membrane computing was selected as a fast emerging area in computer science. Let us try to keep track of the lines of work related to the modelling of phenomena in membrane computing, which were gradually gaining strength in those first years. Indeed, we can notice an important milestone in 2005, when the collective volume [19] was elaborated. There were several teams working on different approaches, but sharing similar interests. 2The p53 protein is closely related to the control of the apoptosis process (also known as “programmed cell death”), and therefore is quite relevant in cancer research. For example, the teams from Universities of Verona and Milano-Bicocca, leaded by Prof. Vincenzo Manca and Prof. Giancarlo Mauri, respectively, have a long tradition in this direction, both in DNA and membrane computing. Their work has produced not only a large number of publications, but also the development of the corresponding software tools and environments for modelling, as well as many PhD Theses: G. Franco Biomolecular Computing Combinatorial Algorithms and Laboratory Experiments (2006), L. Bianco Membrane Models of Biological Systems (2007), and A. Castellini Algorithms and Software for Biological MP Modeling by Statistical and Optimization Techniques (2010). In 2005, Mario laid the foundation stone of the Sevilla team modelling contributions. More precisely, he designed in collaboration with Romero-Campero (who was at that time one of his PhD students) a model of the Epidermal growth factor receptor (EGFR) signalling network, using Continuous P systems [40]. In this framework, objects multiplicity is represented by positive real numbers, and at each instant rules are considered to be applied a positive real number of times. The concept of computation step is formally replaced by an evolution function working over matrices of real numbers, although in the software simulator it is obviously necessary to work with an approximation. The robustness of this model was studied in [41]. Continuous P Systems In [40] a continuous variant of P systems was introduced to model the epidermal growth factor receptor (EGFR) signalling cascade. In contrast to the models developed so far, these systems could evolve in every instant by applying a maximal set of rules a positive real number of times each. Another significant difference was that they worked with continuous multisets, a mapping from an alphabet Σto R+, the set of non-negative real numbers. The rules that are used are of the form Rl={r:u[v]l Kr −→ u[v]l} where u, v, u,vare standard multisets (i.e. with natural multiplicities) and Kr is the function that determines the number of times that the rule will be applied, depending on the multiplicities of objects in the current configuration. This variant was inspired by the fact that in vivo chemical reactions evolve in a continuous way following a rate that depends on the concentration of the reactants (mass action law). The rules, used to model protein-protein interactions taking place in the compartmentalised structure of the living cell, are usually classified as follows: – Transformation, complex formation and dissociation rules: (1) [ a]l→[b]l (2) [ ab]l→[c]l (3) [ a]l→[bc]l – Diffusing in and out: (4) [ a]l→a[] l (5) a[] l→[a]l – Binding and unbinding rules: (6) a[b]l→[c]l (7) [ a]l→b[c]l – Recruitment and releasing rules: (8) a[b]l→c[] l (9) c[] l→a[b]l where a, b, c are objects from the working alphabet and lis a membrane label. An instantaneous configuration of a continuous P system Πis a matrix of An,m(R+) where ai,j (the element in row iand column j) represents the multiplicity of the object cjin membrane i. An evolution of a continuous P system is a mapping, E:R+−→ A n,m(R+), which associates each instant t∈R+with an instantaneous configuration of the system. The rules are applied during the evolution of the system in a continuous manner according to their rate of application function. Observe that the effects of the application of the rules are twofold: the multiplicity of objects appearing in the right-hand side of the rules (products) is increased, while at the same time the multiplicity of objects appearing in the left-hand sides (reactants) is decreased. In order to approximate the evolution of a continuous P systems in a finite set of instants t0,...,t qthe rectangle rule numerical method to approximate integrals is used. It is supposed that tj+1 −tjfor j=0,...,q−1 is fixed, and small enough to assume that all Krremains constant in any interval [tj,t j+1]. With this assumption in mind, the number of times that a rule ris applied during one of those intervals is approximately (tj+1 −tj)·K r(E(tj)). Therefore, by doing this approximation the evolution of a continuous P system during a time interval [t0,t q] is approximated by the computation of a discrete P system that performs qsteps working in a bounded parallel manner. 3 At the Crossroads of Cell Biology and Computation The title of this section was actually the motto of the Seventh WMC, which was celebrated in 2006 in Leiden (The Netherlands). After previously described work, in [18,33,42], an improved extension of those P systems was presented. This extension considers that the application of rules is not instantaneous, but takes a predefined amount of time. In literature, each chemical reaction r:A+Bcr −→ Chas an associated mesoscopic rate constant cr. The following type of P systems captures reaction times as waiting times that determine the order in which reactions place. This time will be computed in a deterministic way, for each reaction, in the following way: τr=1 cr|A||B|where |A|and |B|represent the number of molecules of the two reactants A,B.ρr=cr|A||B|is considered the probability of the rule to be applied in the next step of evolution. 3.1 Deterministic Waiting Times Algorithm The following algorithm represents the most natural way of defining the evolution of such a P system. 1. Set tas 0. 2. Calculate WT ={(τr,r,i)}for all membranes iin the structure μand for every rule r∈R li. 3. Until the time of the simulation treaches or exceeds a prefixed time T: (a) Sort WT according to their waiting time τ. (b) Select from WT the tuple with minimal τ,(τM,r M,i M) (if there are several τvalues, then select all). (c) Update WT by subtracting τMto the τof all its elements. (d) Update tby adding τM. (e) Apply the selected rules, rM, in their corresponding membranes, iM, only once. (f) Recalculate the τonly for those rules which are associated with those compartments affected by the applied rules. (g) For each of such rules, compare the new τwith its existing τand update WT as min(τ,WT). This algorithm simulated and translated the signalling cascade of the epidermal growth factor receptor (EGFR). In a similar manner, it was also used over the set of rules describing the Type I and Type II of FAS-induced apoptotic pathway starting with the stimulation of FAS ligand until the activation of the effector Caspase-3 (see [18]). 3.2 Multienvironment P Systems and Multi-compartmental Gillespie Algorithm Quorum sensing systems in bacteria are fundamental to the control and regulation of cell behaviour. In particular, in order to capture in detail the activation system of a gene regulation system which depends on cell density, the geographical information is very relevant, and this is why it is worth integrating it within the model [8,44]. Amultienvironment P system is a construct, ME=(H,Γ,G,ME,Π,k,RE), where: 1. His a finite set of labels for the environments. 2. Γ={o1,...,o mE}is a finite alphabet of objects (also for chemical substances). 3. G=(V,S) is a graph with nEnodes, V, that represent the environments (labelled with elements from L) and whose edges, Sdefine how the environments are linked. 4. ME={(MEi,h i)}1≤i≤nEis the initial configuration of ME.Sassociates each environment jin the graph Gwith a label, hi∈Hand a continuous multiset, MEi. 5. Πis a continuous P system with Γ⊆Σ. 6. k∈Nis the number of copies of the P system Πthat are non-deterministically distributed across the different environments in the initial configuration of the system. 7. RE={Rh}h∈Hare finite sets of rules of one the following forms: –r:(uρr −→ u)hwhere u, u∈Γ∗. –c:(v)h ρc −→ (v)hwhere v∈Γ∗for the case of movement of different substances from one environment to one of its neighbouring environments. –m:(Π)h ρm −→ (Π)h. In addition to the multisets of objects that represent chemical substances, a certain number of copies of P systems are placed inside the environments. These P systems, and all their contents, can move from one environment to another. Gillespie’s algorithm provides an exact method for the stochastic simulation of systems of bio-chemical reactions. The validity of this method has been rigorously proved and successfully used to simulate various biochemical processes. An extension of this algorithm, called Multi-compartmental Gillespie’s Algorithm, was introduced in [42]. Unlike the original version, this method considers the existence of multiple disjoint compartments that represent different regions where chemical reactions occur. Classical Gillespie’s Algorithm Let us consider an enumeration Rof all the rules for one of the membranes of the kP systems, including also the environments. Each one of them is considered to be a separate compartment enclosing a volume. 1. Calculate a0= ri∈R ρri 2. Generate two random numbers n1and n2over the unit interval (0,1) 3. Calculate the waiting time τ=1 a0 ln(1 n1 ) 4. Take the index jof the rule such that  1≤k≤j−1 ρrk<n 2·a0≤ 1≤k≤j ρrk 5. Return (τ,j) Multi-compartmental Gillespie’s Algorithm Let us consider an enumeration Cof all such compartments. 1. Set tas 0. 2. Calculate WT ={(τ,j,i):i∈C}using Classical Gillespie’s algorithm (as described above) to calculate (τ,j) for each compartment. 3. Until the time of the simulation treaches or exceeds a prefixed time: (a) Sort WT according to τ (b) Select from WT the tuple with minimal τ,(τM,j M,i M) (c) Update WT subtracting τMto each τ (d) Apply rjMin compartment iMonly once. (e) For those compartments iaffected by the applied rule Recalculate their corresponding (τ,j), using again Classical Gillespie’s algorithm, and update WT 4 Sevilla’s Ark: Giant Pandas, Bearded Vultures and Zebra Mussels The origin of the fascinating journey of RGNC across the ocean of computational modelling of ecosystems began in 2008. The first case study focused on an ecosystem related to the Bearded vulture (Gypaetus barbatus) in the Pyrenees. In this line, Probabilistic P systems were presented in [16,17]. Although the results qualitatively agreed with experimental data, this model was intended to be a preliminary proof of concept. Shortly after, the model was extended and improved in [15], by adding more species and features in order to improve the model’s accuracy. The next upgrade on the modelling framework enabled the modelling of geographical information, yielding the so-called Multienvironment probabilistic functional extended P systems [14,22]. This is useful to capture, for instance, scavengers moving along different areas looking for food, featuring different environmental conditions for each area, or the expansion of a disease among a population [20]. In parallel with such refinements of the technical details about the syntax and semantics of the type of P systems that were used, there was another evolution going on the software part. In particular, several simulation algorithms were engineered, implementing in different ways the hybridisation between probabilistic rules and their maximally parallel mode of application. An abstract virtual ecosystem having three trophic levels (grass, herbivorous and carnivorous) was designed, to be used as a scalable case study to perform virtual experiments and compare the characteristics of each algorithm [24,25,29]. On the other hand, after several ad-hoc software developments for simulation tools, including enduser GUIs, the RGNC started to work on a general purpose solution: MeCoSim (Membrane Computing Simulator) (see [22,39]andvisit[2]). In line with this general purpose long-term approach, a step-by-step protocol for building computational models was presented in [21], using Population Dynamics P systems (PDP systems). This is a generalisation of the syntax of previous frameworks, whose integration in pLinguaCore v3.0 was announced together with a new simulation algorithm in [28]. A noticeable case study is gene networks (dynamics of logical networks [49], or Arabidopsis thaliana’s regulating its circadian rhythms [48]). 4.1 Probabilistic Systems An extension of multienvironment P systems, encompassing them, was proposed in [16,17] with the following changes: – Each environment contains exactly one P system. In the initial configuration, the multisets associated with each P system is empty. – Environments and membranes have no associated labels. Previous models usually used exactly one label for environments and distinguished them as different elements in the set of nodes V={e1,...,e nE}of the graph Gthat defines their interconnections. Since the association between membranes and 42. P´erez-Jim´enez, M.J., Romero-Campero, F.J.: P systems, a new computational modelling tool for systems biology. In: Priami, C., Plotkin, G. (eds.) Transactions on Computational Systems Biology VI. LNCS, vol. 4220, pp. 176–197. Springer, Heidelberg (2006). https://doi.org/10.1007/11880646 8 43. Regev, A., Shapiro, E.: Cellular abstractions: cells as computations. Nature 419(6905), 343 (2002) 44. Romero-Campero, F.J., P´erez-Jim´enez, M.J.: A model of the quorum sensing system in Vibrio fischeri using P systems. Artif. Life 14, 95–109 (2008) 45. Suzuki, Y., Fujiwara, Y., Takabayashi, J., Tanaka, H.: Artificial life applications of a class of P systems: abstract rewriting systems on multisets. In: Calude, C.S., P˘aun, G., Rozenberg, G., Salomaa, A. (eds.) WMC 2000. LNCS, vol. 2235, pp. 299–346. Springer, Heidelberg (2001). https://doi.org/10.1007/3-540-45523-X 16 46. Suzuki, Y., Ogishima, S., Tanaka, H.: Modeling the p53 signaling network by using P systems. In: Alhazov, A., Martin-Vide, C., P˘aun, G. (eds.) Workshop on Membrane Computing, WMC-2003, pp. 449–454, 2003. Rovira i Virgili University, Technical report 28 47. Suzuki, Y., Tanaka, H.: Abstract rewriting systems on multisets and their application for modelling complex behaviours. In: Cavaliere, M., Martin-Vide, C., P˘aun, G. (eds.) Brainstorming Week on Membrane Computing, pp. 313–331 (2003). Rovira i Virgili University, Technical report 26 48. Valencia-Cabrera, L., Garc´ıa-Quismondo, M., P´erez-Jim´enez, M.J., Su, Y., Yu, H., Pan, L.: Analysing gene networks with PDP systems. Arabidopsis thaliana, a case study. In: Eleventh Brainstorming Week on Membrane Computing (11 BWMC), pp. 257–272, August 2013 49. Valencia-Cabrera, L., Garc´ıa-Quismondo, M., P´erez-Jim´enez, M.J., Su, Y., Yu, H., Pan, L.: Modeling logic gene networks by means of probabilistic dynamic P systems. Int. J. Unconv. Comput. 9(5–6), 445–464 (2013) 50. Zhang, G., P´erez-Jim´enez, M.J., Gheorghe, M.: Real-Life Applications with Membrane Computing. ECC, vol. 25. Springer, Cham (2017). https://doi.org/10.1007/ 978-3-319-55989-6