A bioinspired computing approach to model complex systems

A Bioinspired Computing Approach to Model Complex Systems Mario J. P´erez-Jim´enez(B) Research Group on Natural Computing, Department of Computer Science and Artificial Intelligence, University of Sevilla, Seville, Spain [email protected] Abstract. The use of models is intrinsic to any scientific activity. In particular, formal/mathematical models provide a relevant tool for scientific investigation. This paper presents a new Membrane Computing based computational paradigm as a framework for modelling processes and real-life phenomena. P systems, devices in Membrane Computing, are not used as a computing paradigm, but rather as a formalism for describing the behaviour of the system to be modelled. They offer an approach to the development of models for biological systems that meets the requirements of a good modelling framework: relevance, understandability, extensibility and computability. Keywords: Membrane Computing ·Multienvironment P systems · Multicompartimental P systems ·Population Dynamics P systems 1 Introduction Scientists regularly use abstractions with the aim to describe and understand the reality they are examining. Computational modelling is the process of representing real world problems in mathematical terms in an attempt to find solutions to their associated complex systems. A formal model is an abstraction of the real-world onto a mathematical/computational domain that highlights some key features while ignoring others that are assumed to be secondary. A formal model should not be considered as representation of the truth, but instead as a statement of our current knowledge of the phenomenon under research. It is desirable for a model to fulfill four properties: relevance,understandability, extensibility and computability [19]. A formal model must be relevant capturing the key featureswhileignoringothers assumedto besecondary.The abstractformalism used should adecquately match the informal concepts and ideas from the investigated phenomenon. Mathematical models should also be extensible to higher level of organizations, like tissues, organs, organisms, etc, in the case of cellular systems. Finally, a formal model should be able to be implemented in a computer so that we can run simulations to study the dynamics of the system in different scenarios, as well as the qualitative and quantitative reasoning about its properties. c Springer International Publishing Switzerland 2014 M. Gheorghe et al. (Eds.): CMC 2014, LNCS 8961, pp. 20–34, 2014. DOI: 10.1007/978-3-319-14370-5 2 A Bioinspired Computing Approach to Model Complex Systems 21 One of the main objectives of any model is to provide a predictive capability, that is, the possibility to make guesses in terms of plausible hypotheses related to the dynamics of the observed phenomenon in different scenarios that are of interest to experts. Cellular systems and population biology often depend on many parameters related to the observed behaviours. Since they define the dynamics of the system, parameters must satisfy some conditions, which can be referred to as the invariants of the associated behaviour. Some of these invariants can be expressed by rules and can be obtained by carrying out experiments, while others cannot be measured or they are very expensive to estimate. Therefore, before simulations can be performed in order to make predictions, we need to calibrate our model. Several parameters values are tested by calibration and the results corresponding to the state parameters are compared with the observed/expected behaviour of the system for the same state parameters. In some cases, the design of the model has to be reconsidered [15]. Nowadays ordinary/partial differential equations (ODEs/PDEs) constitute the most widely used approach in modelling complex systems. Nevertheless, in some cases such as molecular interaction networks in cellular systems, any model described by means of a system of ODEs/PDEs is based on two assumptions: (a) cells are assumed to be well stirred and homogeneous volumes so that concentrations do not change with respect to space; and (b) chemical concentrations vary continuously over time in a deterministic way. This assumption is valid if the number of molecules specified in the reaction volume are sufficiently large and reactions are fast. Membrane Computing is an emergent branch of Natural Computing introduced by G. Paun at the end of 1998. This new computing paradigm starts from the assumption that processes taking place in the compartmental structure of a living cell can be interpreted as computations. In contrast to differential equations, P systems explicitly correspond to the discrete character of the components of a complex system and use rewriting/evolution rules on multisets of objects which represent the variables of the system. The inherent stochasticity, external noise and uncertainty in cellular systems is captured by using stochastic or probabilistic strategies. A general bioinspired computing modelling framework, called multienvironment P systems is introduced. The paper is structured as follows. First, the framework of multienvironment P systems is defined in a formal way. Section 3 is devoted to multicompartmental P systems, the stochastic approach. Besides, four case studies at cellular level are presented in this Section. Population dynamics P systems, the probabilistic approach, are studied in Section 4 and three case studies related to real ecosystems are described. Finally, some conclusions are drawn. 2 Multienvironment P System Amultienvironment P system of degree (m, n, q) taking Ttime units is a tuple (G, Γ,Σ,µ,T,Π1,...,Πn,R,E 1,...,E m,RE) 22 M.J. P´erez-Jim´enez where: –G=(V,S) is a directed graph. Let V={e1,...,e m}whose elements are called environments; – Γ is the working alphabet and Σ Γ. –µis a rooted tree with qnodes. –Tis a natural number that represents the simulation time of the system; –Π k=(Γ,µ,M1,k,...,Mq,k,R), 1 ≤k≤n, is a basic P system of degree q,and Ris a finite set of rules of the type u[v]α i−→ u[v]β i. –Ej∈Mf(Σ), 1 ≤j≤m. –REis a finite set of communication rules among environments of the following forms (x)ej−→ (y1)ej1...(yh)ejhand (Πk)ej−→ (Πk)ej. – No rules from Rand REcompete for objects. – Each rule of the system has associated a computable function whose domain is {0,...,T}. A multienvironment P system can be viewed as a finite set of environments and a finite set of P systems, such that: (a) the links between environments are given by the arcs taken from a directed graph; (b) each P system has the same working alphabet, the same membrane structure and the same evolution rules; and (c) each environment contains several P systems, where each evolution rule has associated a computable function and each one of them has an initial multiset which depends on the environment; and (d) there is a finite set of rules among the environments. Furthermore, inside the environments, only objects from a distinguished alphabet can exist. It is worth pointing out that this bioinspired computational approach has some qualitative advantages with respect to ordinary/partial differential equations approach: – They use a language closer to experts than differential equations. – They are not affected by the usual constraints present when defining differential equations based models. – They are modular, that is, once an initial version of the model is designed, adding modifications is relatively easy. On the one hand, small changes in the system entails small changes in the model. On the other hand, when using differential equations most times we have to start from scratch. 3 Stochastic Approach: Multicompartmental P Systems A multienvironment P system of degree (m, n, q) taking Ttime units, is said to be stochastic if: (a) the computable functions associated with the rules of the P systems are propensities: they are computed from stochastic constants by applying the law of mass-action law (the reaction rate depends proportionally on the product of the concentrations of the reactants), and the stochastic constants are obtained from the kinetic constants in an easy way [18]; these rules depend on time but not on the environment; A Bioinspired Computing Approach to Model Complex Systems 23   Fig. 1. A multienvironment extended P system (b) at the initial moment, the P systems Π1,...,Πnare ramdomness distributed between the menvironments; for instance, as in the picture. These kind of P systems are called multicompartmental P systems. The dynamics of these systems is captured by either the multicompartmental Gillespie’s algorithm [18]orthedeterministic waiting time algorithm [2]. Gillespie’s algorithm [10–13] provides an exact method for the stochastic simulation of systems of bio-chemical reactions; the validity of the method is rigorously proved and it has been already successfully used to simulate biochemical processes and it is based on the inversion method of Monte Carlo theory. The deterministic waiting time algorithm is based on the fact that in vivo chemical reactions take place in parallel in a asynchronous manner. The time taken for the formation of each molecule, called waiting time, is calculated and the rule (reaction) with the least waiting time is applied, changing the concentration in their respective compartments. Each time there is a change in the concentration of a molecule in any compartment, the waiting time for reactions “using” that molecule needs to be recalculated for the compartment. Infobiotics workbench is a computational framework implementing a synergy among multicompartmental stochastic simulations, formal model analysis and structural/parametermodel optimisationforcomputationalsystems andsynthetic biology (http://www.infobiotics.org/index.html) Inthissection somepractical examples ofmulticompartmentalP systems applications for modelling cellular systems are presented. 3.1 Apoptosis Mediated by FAS Protein The FAS-induced apoptotic signalling pathway was shown to be one of the most relevant processes for understanding and combating cancer, AIDS and neurodegenrative diseases such as Parkinson’s disease, Alkzheimer, etc. 24 M.J. P´erez-Jim´enez Two pathways activated by FAS have been identified [23], and are referred to as type I (death receptor pathway) and type II (mitochondrial pathway), where caspases play a crucial role for both the initiation and execution of apoptosis (programmed cell death). The pathways diverge after activation of initiator caspasesand convergeattheend by activating executorcaspases.Inthetype Ipathway, initiator caspases activate executor caspases directly. In the type II pathway, a more complicated cascade is activated involving the disruption of mitochondrial membrane potential. We have designed a multicompartmental P system with only one environment thatconsists of 53objectsand 99 evolution rules(see [2] fordetails) in orderto study when a cell chooses the mitochondrial pathway or the death receptor pathway to produce apoptosis. AJava simulator has been implemented and it accepts as input a Systems Biology Markup Language (SBML) file containing the rules to be simulated and initial concentrations for the molecules in the system. We used the Cell Designer package to generate the SBML source file for the reactions (Cell Designer is a structured diagram editor for drawing gene-regulatory and biochemical networks). The simulator engine mimics the biological cell and it is designed in a modular way so that it can use different strategies for different pathways if needed. The specific strategy, based on the deterministic waiting time algorithm, will be executed depending on the initial concentrations of various objects present in the system. A simulator designed in Scilab, a scientific software package for numerical computations providing a powerful open computing environment for engineering and scientific applications [23], using the multicompartmental Gillespie algorithm has been also considered [2]. The consistency between the framework and the experimental results in the paper [14] validates our model. We have stated that our discrete methods handle low levels of molecules in a different way that ODE/PDE techniques. To further investigate the differences between discrete and ODE/PDE methods, we have chosen to focus on one rule from the FAS-mediated pathway (a transformation) Multicompartmental P systems constitute an alternative to ordinary/partial differential equations methods. We have argued that the discrete nature of our techniquemightbebetter forsimulatingtheevolution ofsystemsinvolvinglownumbers of molecules. 3.2 Gene Regulation Systems in Lac Operon in E. coli In most bacteria, gene expression is highly regulated in order to produce the necessary proteinic machinery to respond to environmental changes. Therefore, at a given time, a bacterial cell synthesises only those proteins necessary for its survival under the particular conditions of that time. Many of the genes in Escherichia coli (E. coli) are expressed constitutively; that is, they are always turned on. Others, however, are active only when their products are needed by the cell, so their expression must be regulated. The most direct way to control the expression of a gene is to regulate its rate of transcription. Adding a new substrate to the culture medium may induce the formation of new enzymes capable of metabolising that substrate. An example of this phenomenon happens A Bioinspired Computing Approach to Model Complex Systems 25 when we take a culture of E. coli that is feeding on glucose and transfer some of the cells to a medium containing lactose instead. In this case a revealing sequence of events takes place. A multicompartmental P system modelling the gene expression control in the Lac Operon has been designed (see [20] for details). Specifically, the system has only one environment, the total number of symbols in the working alphabet is 51 and there are 55 evolution rules. The novelty of this design is that the objects can be symbols or strings over the alphabet. In this context, finite multisets of strings within membranes represent the genetic information encoded in DNA and RNA. Thecentraldogma ofmolecularcellbiology statesthatgenetic informationisstored in the DNA. This information is transcribed into mRNA which in turn is translated into proteins. It is worth pointing out that transcription and translation have been modelled as rewriting and concurrent processes on strings. Using the Multicompartmental Gillepie’s Algorithm and a simulator developed in Scilab, we have studied the behaviour of the system for different environmental conditions to see how the system is able to sense the presence of different substrate (glucose and lactose). The delay between the sensing of the signal and the expression of different genes is not explicitly modelled but emerges as a consequence of the formulation of our approach. Our results agree well with experimental observations and results obtained by using other approaches. 3.3 Quorum Sensing in Vibrio Fischeri Bacteria are generally considered to be independent unicellular organisms. Nevertheless,in somecircumstances bacteriaexhibit coordinatedbehaviour which allows an entire population of bacteria to regulate the expression of specific genes depending on the size of the population.This phenomenon is called quorum sensing, that is, a cell density dependent gene regulation system. It was first investigated in the marine bacterium Vibrio Fischeri. The bacteria colonise specialised light organs in the squid which cause it to become luminescent. Vibrio Fischeri only causes luminescence when colonising the light organs and do not emit light when in the planktonic free-living state. Luminescence in the squid is involved in the attraction of prey, camouflage and communication between different individuals. Bacteria colonies behave like multicellular organisms. Each bacterium must be able to sense and communicate with other units in the colony to express some specific genes in a coordinated way. The cooperative activities carried out by members of the colony generate a social intelligence [17]. In this case, we have designed a multicompartmental P system of degree (25,n,1), with 25 environments containing each of them an ordinary P system only having the skin membrane (see [21] for details) This model has been represented in the Systems Biology Markup Language using Cell Designer [9]. The emergent behaviour of the system has been studied for three colonies of different size (10, 100 and 3000 bacteria) to examine how bacteria can sense the number of individuals in the colony and produce light only when such number is big enough. 26 M.J. P´erez-Jim´enez Our simulations show that Vibrio fischeri has a quorum sensing system where a single bacterium can guess that the size of the population is big enough and start to produce light. This bacterium starts to massively produce signals, but if the signal does not accumulate in the environment it means that the guess was wrong and it switches off the system. In contrast, if the signal does accumulate in the environment meaningthat thenumberof bacteriain thecolonyisbigenough, arecruitment process takes place that causes the entire population of bacteria to become luminescent. Let us stress that this emergent behaviour is a result of local interactions in the environments between different simple agents, the bacteria, which are only able to produce and receive molecular signals. These results agree well with in vitro observations. 4 Probabilistic Approach: Population Dynamics P Systems (PDP Systems) A multienvironment functional P system with active membranes of degree (q,m,n) taking Ttime units, is said to be probabilistic if: (a) the computable functions associated with the rules of the systems are probability functions verifying some conditions (these rules depend on the environment but not on the time, they are constant functions); (b) the total number of P systems Πkis equal to the number of environments: each environment contains one P system; and (c) the rules among environments are only of the form: (x)ej−→ (y1)ej1... (yh)ejh   Fig. 2. A Population Dynamics P system A Bioinspired Computing Approach to Model Complex Systems 27 This kind of P systems are called Population Dynamics P systems [5]. The dynamics/semantics of these systems is captured by ad hoc algorithms such as Binomial block based simulation algorithm (BBB)[1], Direct Non Deterministic distribution with Probabilities algorithm (DNDP)[7]andDirect distribution based on Consistent Blocks algorithm (DCBA)[16], among others. The BBB algorithm follows a strategy based on the binomial distribution and blocks of rules with the same left-hand side. In general, each simulation step is divided into two main stages: selection and execution. In the first one, the algorithm decides which rules will be applied, and the number of applications for each one (taking into account their left-hand sides and the available objects in the current configuration). In the second stage, the selected rules are applied, consuming themultisetsofthe rules’left-hand sidesandadding themultisetsof therules’ righthand sides the selected number of times, and possibly changing the polarization of membranes. The DNDP algorithm performs a non-deterministic distribution of objects along the rules, but considering the probabilities. The algorithm is split into two phases, selectionand execution.This time,selection phaseis dividedintotwomicro-phases: selection phase 1 (consistency) and selection phase 2 (maximality). Together with an initialization phase, it has a total of four phases. The first selection phase calculates a multiset of consistent applicable rules. This is performed by looping the rules in a random order, and applying each one (if consistent with the already selected rules) using the binomial distribution according to the probabilities. The second selection phase eventually increases the multiplicity of some of the rules in the previous multiset to assure maximal application, obtaining a multiset of maximally consistent applicable rules. Again, there is a loop over the remaining rules, checking the maximality condition. Although the DNDP algorithm achieves better results than its predecessor (BBB), the behaviour still produces some distortion in many situations (it is biased towards the rules with the highest probabilities). The DCBA algorithm is based on the idea of proportionally distributing the amount of objects along the rule blocks. A proportional calculus is made in such a way that rules requesting for more objects are penalized. However, this calculation can be adapted to the biological semantics to be captured by the model. Probabilities are applied to rule blocks locally. The simulation algorithm consists on two phases, selection and execution. But this time, selection is split into three microphases: phase 1 (distribution), phase 2 (maximality), and phase 3 (probabilities). Selection phase 1 uses a distribution table, where rows represent objects inside regions, and columns are rule blocks. A normalized distribution of the objects is performed over the rows. Phase 2 iterates the remaining rule blocks assuring maximality, and phase 3, once rule blocks have been selected, calculates multinomial distributions for each one (according to the selected number for it, and the probabilities of the corresponding rules). DCBA is able to reproduce the desired semantics for the model of PDP systems. However, its efficient implementation is a challenge (the distribution table can be very large). P-Lingua (http://www.p-lingua.org/wiki/index.php/Main Page)isaprogramming language for Membrane Computing which aims to be a standard to 28 M.J. P´erez-Jim´enez define P systems, in particular, population dynamics P systems. MeCoSim (http:// www.p-lingua.org/mecosim/) is a visual environment to model, simulate, analyse and verify solutions based on P systems, by defining custom apps for virtual experimentation under different scenarios. In what follows some practical examples of PDP systems applications for modelling real ecosystems are presented. 4.1 Bearded Vulture The Bearded Vulture is a cliff-nesting and territorial large scavenger distributed in mountains ranges in Eurasia and Africa. This is one of the rarest raptors in Europe (150 breeding pairs in 2007). This endangered species feeds almost exclusively on bone remains of wild and domestic ungulates. Its main food source is bone remains of dead small and medium-sized animals. The ecosystem to be modelled is in the Pyrenean and Prepyrenean mountains of Catalonia (NE Spain) and it is composed of 13 species: (a) three avian scavengers (Bearded vulture,the Egyptian vulture and the Griffon vulture as predator species); (b) six wild ungulates (Pyrenean Chamois, Red deer, Fallow deer, Roe deer, Mouflon and Wild boar); and (c) four domestic ungulates (sheep, cow, goat andhorse) thatarefound inanextensiveorsemi-extensiveregime providingcarrion for the avian scavengers and considered as prey species. Prey species are herbivores and their remains form the primary food resource for the avian scavengers in the study area. In order to model this real ecosystem, a population dynamics P system with two environments containing each of them an ordinary P system of degree 2 has been considered (see [3] for details). The model addresses: (a) the population dynamics previously mentioned; (b) the interactions among the 13 species; (c) the presence of two zones in the study area; (d) the communication protocol between the two areas; and (e) the ecosystem maximum load capacity for each of the areas. The algorithmic scheme of the proposed model is shown in Figure 3and it is structured following a series of modules which are run sequentially corresponding to the passing of 1 year in the ecosystem. We have studied the dynamics of the ecosystem modifying the initial conditions in order to analyse how the ecosystem would evolve if different biological factors were modified either by nature or through human intervention. We have designed a population dynamics P system with only one environment and 49 types of rules. For each type of animal the number of biological parameters are related to reproduction, mortality, feeding and other general processes of the species itself. We have shown the robustness of the model with respect to a modified order of application of the different processes modules. MeCoSim software has been used for the execution of the model. The population trend of the three scavenger species and the six wild ungulates has been obtained by