A Hybrid Approach to Modeling Biological Systems

A Hybrid Approach to Modeling Biological Systems Francesco Bernardini1, Marian Gheorghe2, Francisco Jos´e Romero-Campero3, and Neil Walkinshaw2 1Leiden Institute of Advanced Computer Science, University of Leiden Niels Bohrweg 1, 2333 CA Leiden, The Netherlands [email protected] 2Department of Computer Science, University of Sheffield Regent Court, Portobello Street, Sheffield S1 4DP, UK [email protected] 3Research Group on Natural Computing Department of Computer Science and Artificial Intelligence, University of Sevilla Avda. Reina Mercedes s/n, 41012, Sevilla, Spain [email protected] Abstract. This paper investigates a hybrid approach to modeling molecular interactions in biology. P systems, π-calculus, and Petri nets models, and two tools, Daikon, used in software reverse-engineering, and PRISM, a probabilistic model checker, are investigated for their expressiveness and complementary roles in describing and analyzing biological systems. A simple case study illustrates this approach. 1 Introduction In the last decade there has been a great interest in using theoretical computer science models in biology, based on different paradigms (process algebras, cellular automata, Lindenmayer systems, Petri nets, Boolean functions, P systems, etc.) with the aim of providing an understandable, extensible and computable modeling framework while keeping the needed formalization to perform mathematical analysis. Every such model covers certain aspects of a system and combining two or more leads to obtaining a better and more comprehensive modeling approach. In order to include quantitative and qualitative aspects, there have been suggested various variants of certain models with new features like: Petri nets [10,22], stochastic π-calculus [28], and stochastic P systems [19]. In this paper we investigate the concerted use of different methods that will reveal a new vision on modeling biological systems by combining different complementary approaches. This is quite different from the hybrid approach discussed by [1] where it is shown how to switch between deterministic and stochastic behavior. Section 2 introduces the three modeling approaches used in the paper: P systems, pi-calculus and Petri nets, as well as Daikon tool and a simple example involving a regulatory network that will be modeled within each approach. Section 3 presents Daikon’s findings and the analysis of the invariants provided. The following two sections show how PRISM and a Petri net tool, PIPE, are used in order to confirm some of the properties suggested by Daikon. The final section summarizes our findings. 2 Modeling Paradigms In this section we present three modeling approaches, namely P systems, πcalculus and Petri nets and a simple case study to illustrate the approach. This example is written directly into these three modeling paradigms and will be executed with a P systems simulator. A system of differential equations is associated to this example and the results obtained are compared to the stochastic behavior exhibited by the P systems simulator and PRISM. In the next three sections Daikon is used to reveal certain properties of our models as they appear through data sets generated by simulators, and PIPE, a Petri nets tool, to ascertain some invariants of the system. and two tools, namely PRISM and PIPE, that are used to analyze and verify properties identified by Daikon. The aim of this investigation is not to study the relationships between the results produced by using differential equations and those generated by P systems and/or PRISM. The relationships between a special class of P systems working in a deterministic manner according to a metabolic algorithm and differential equations models has been already considered [7]. In this study we are using differential equations only as a substitute for real data in order to illustrate our approach that allows us to “guess” certain properties of the model and then to verify whether they hold or not as general properties or just only happens to be true for the instances generated by simulation. Nowadays ordinary differential equations (ODE) constitute the most widely used approach in modeling molecular interaction networks in cell systems. They have been used successfully to model kinetics of conventional macroscopic chemical reactions. Nevertheless the realization of a reaction network as a system of ODEs is based on two assumptions. First, cells are assumed to be well stirred and homogeneous volumes so that concentrations do not change with respect to space. Whether or not this is a good approximation depends on the time and space scales involved. In bacteria it has been shown that molecular diffusion is sufficiently fast to mix proteins. This is not the case in eukaryotic cells where the volume is considerably bigger and it is structured in different compartments like nucleus, mitochondria, Golgi body, etc. The second basic assumption is that chemical concentrations vary continuously over time. This assumption is valid if the number of molecules of each species in the reaction volume (the cell or the subcellular compartment)are sufficiently largeand the reactionsare fast. A sufficiently large number of molecules is considered to be at least thousands of molecules; for hundreds or fewer molecules the continuous approach is questionable. Writing and solving numerically a system of ODE describing a chemical reaction network can be largely automated. Each species is assigned a single variable X(t) which represents the concentration of the species at time t. Then, for each molecular species, a differential equation is written to describe its concentration change over time due to interactions with other species in the system. The rate of each reaction is represented using a kinetic rate law, which commonly depends on one or more rate constants. Exponential decay law, mass action law and Michaelis-Menten dynamic are the most widely used kinetic mechanisms. Finally in order to solve the system of ODEs we must impose a set of initial condition representing the initial concentration of each species involved. Due to the limitations of ODEs to handle cellular systems with low number of molecules and spatial heterogeneity, some computational approaches have been recently proposed. In what follows we discuss three different approaches, P systems, π-calculus, and Petri nets. 2.1 P Systems Membrane computing is an emergent branch of natural computing introduced by Gh. P˘aun in [18]. The models defined in this context are called P systems. In the sequel we will use membrane computing and P systems with the same meaning. Roughly speaking, a P system consists of a cell-like membrane structure, in the compartments of which one places multisets of objects and strings which evolve according to given rules. Recently P systems have been used to model biological phenomena within the framework of computational systems biology presenting models of oscillatory systems [6], signal transduction [19], gene regulation control [20], quorum sensing [27] and metapopulations [21]. In this respect, P systems present a formal framework for the specification and simulation of cellular systems which integrates structural and dynamic aspects in a comprehensive and relevant way while providing the required formalization to perform mathematical and computational analysis. In the original approach of P systems the rules are applied in a maximally parallel way. This produces two inaccuracies: the reactions represented by the rules do not take place at the correct rate, and all time steps are equal and do not represent the time evolution of the real system. In order to solve these two problems stochastic P systems have been introduced in [19]. Definition 1. A stochastic P system is a construct Π=(O,L,μ, M1,M 2,...,M n,R 1,...,R n), where: –Ois a finite alphabet of symbols representing objects; –Lis a finite alphabet of symbols representing labels for compartments; –μis a membrane structure containing n≥1membranes labeled with elements from L; –Mi=(li,w i,s i),foreach1≤i≤n, is the initial configuration of membrane iwith li∈L, the label of this membrane, wi∈O∗a finite multiset of objects and sia finite set of strings over O; –Ri,foreach1≤i≤n, is a finite set of rewriting rules associated with membrane i, of one of the following two forms: •Multiset rewriting rules: obj1[obj2]l k −→ obj 1[obj 2]l(1) with obj1,obj 2,obj 1,obj 2∈O∗some finite multisets of objects and lalabel from L. A multiset of objects, obj is represented as obj =o1+o2+...+om with o1,...,o m∈O. These rules are multiset rewriting rules that operate on both sides of membranes, that is, a multiset obj1placed outside a membrane labeled by land amultisetobj2placed inside the same membrane can be simultaneously replaced with a multiset obj 1andamultisetobj 2, respectively. •String rewriting rules: [obj1+str1;...;objp+strp]l k −→ (2) [obj 1+str 1,1+...str  1,i1;...;obj p+str p,1+...str  p,ip]l Astringstr is represented as follows str =s1.s2.···.siwhere s1,...,s i∈ O. In this case each multiset of objects objjand string strj,1≤j≤p,are replaced by a multiset of objects obj jand strings str j,1...str j,ij. The stochastic constant kis used to compute the propensity of the rule by multiplying it by the number of distinct possible combinations of the objects and substrings present on the left-side of the rule with respect to the current contents of membranes involved in the rule. The propensity associated with each rule is used to compute the probability and time needed to apply it. Cellular systems consisting of molecular interactions taking place in different locations of living cells are specified using stochastic P systems as follows. Different regionsand compartments are specified using membranes. Each molecular species is represented by an object in the multiset associated with the region or compartment where the molecule is located. The multiplicity of each object represents the number of molecules of the molecular species represented by the object. Strings are used to specify the genetic information encoded in DNA and RNA. Molecular interactions, compartment translocation and gene expression are specified using rewriting rules on multisets of objects and strings - see Table 1. Table 1. Modeling principles in P systems Biochemistry PSystem Compartment Region defined by a membrane Molecule Object Molecular Population Multiset of objects Biochemical Transformation Rewriting rule Compartment Translocation Boundary rule In stochastic P systems [19] constants are associated with rules in order to compute their probabilities and time needed to be applied according to Gillespie algorithm. This approach is based on a Monte Carlo algorithm for stochastic simulation of molecular interactions taking place inside a single volume [8]. In contrast to this, in P systems we have a membrane structure delimiting different compartments (volumes), each one with its own set of rules (molecular interactions) and multiset of objects and strings (molecules). In this respect, a scheduling algorithm called the Multicompartmental Gillespie algorithm [19] is used so that each compartment evolves according to a different Gillespie algorithm. In this point our approach differs from other computational approaches which run a single Gillespie algorithm across the whole system without taking into account the compartmentalized cellular structure [10,28]. We illustrate our approach with a biomolecular system consisting in positive, negative and constitutive expression of a gene. Our model includes the specification of a gene, its transcribed RNA, the corresponding translated protein and activator and repressor molecules which bind to the gene producing an increase in Π=({gene,rna,protein,act,rep,act-gene,rep-gene},{b},[] b,(b, Mi,∅), {r1,...,r 9}) Initial multisets: M0,1=gene;M0,2=gene +act... +act and M0,3=gene +rep... +rep where act and rep occur 10 times each. Rules: r1:[gene ]b c1 −→ [gene +rna ]bc1=0.347 min−1 r2:[rna ]b c2 −→ [rna +protein ]bc2=0.174 min−1 r3:[rna ]b c3 −→ [] bc3=0.347 min−1 r4:[protein ]b c4 −→ [] bc4=0.0116 min−1 r5:[act +gene ]b c5 −→ [act-gene ]bc5=6.6412087 molec−1min−1 r6:[act-gene ]b c6 −→ [act +gene ]bc6=0.6s−1 r7:[act-gene ]b c7 −→ [act-gene +rna ]bc7=3.47 min−1 r8:[rep +gene ]b c8 −→ [rep-gene ]bc8=6.6412087 molec−1min−1 r9:[rep-gene ]b c9 −→ [rep +gene ]bc9=0.6min−1 Fig. 1. P system model of gene expression rna protein 0 100 200 300 400 500 600 700 800 900 1000 0 5 10 15 20 25 30 35 time (min) Number of molecules rna protein 0 100 200 300 400 500 600 700 800 900 1000 0 20 40 60 80 100 120 140 160 180 200 time (min) Number of molecules Fig. 2. Constitutive expression and positive regulation transcription rate or prevent the gene from being transcribed, respectively. The bacterium where the system is located is represented using a membrane. The stochastic constants used in our model are taken from the gene control system in the lac operon in E. coli [2,13,14]. In this case transcription and translation have been represented using rewriting rules on multisets of objects, a more detailed description of the concurrent processes of transcription and translation using rewriting rules on strings is presented in [20]. The P systems model is formally defined in Figure 1. It consists of one single compartment labeled b, with no strings, and consequently using only multiset rewriting rules. The model refers to three distinct initial conditions, denoted by multisets M0,i, and corresponding to constitutive expression, positive and negative regulations, respectively. Simulations of constitutive expression and positive regulation case studies are presented in Figure 2 using a tool available at [30]. A set of ordinary differential equations and their associated graphs, modeling the same examples, are provided in Figure 3. The ODE model is not used here to showits relationship to the previousP systems approach, but to providea set of data that normally is taken through biologicalexperiments. This will only be used to provide data measurements that will help identifying and validating properties of the P systems model. dr dt =c1−c3r+c7act act +Kwhere Kis the Michaelis-Menten constant dp dt =c2r−c4p rna protein 0 50 100 150 200 250 300 350 400 0 2 4 6 8 10 12 14 16 18 20 time (min) Number of molecules rna protein 0 100 200 300 400 500 600 700 800 900 1000 0 20 40 60 80 100 120 140 160 180 time (min) Number of molecules Fig. 3. Constitutive and positive expression using ODE model 2.2 π-Calculus The π-calculus approach was introduced as a formal language to describe mobile concurrent processes [17]. It is now a widely accepted model for interacting systems with dynamically evolving communication topology. The π-calculus has a simple semantics and a tractable algebraic theory. Starting with atomic actions and simpler processes, complex processes can be then constructed. The process expressions are defined by guarded processes, parallel composition P|Q, nondeterministic choice P+Q, replication !P, and a restriction operator (νx)P creating a local fresh channel xfor a process P. Different variants have been used to model molecular interactions [28]. A πcalculus specification of our system is provided by Figure 4. As usual for this type of modeling approach, each chemical element will be represented as a process and its definition will refer to all possible interactions of it. The initial process may be any of S0,i. The process called gene defines all possible interactions of a constitutive reaction, producing messenger RNA, a positive regulation, leading to a complex denoted by act-gene, or negative regulation, that gets the complex rep-gene. This process definition corresponds in a P systems model to rules r1,r5 and r8. In this way we can see, at least syntactically, similarities and differences between the two modeling approaches for expressing chemical interactions. More about the use of both P systems and π-calculus to model chemical interactions is provided by [26]. Table 2. Modeling Principles in π-calculus Biochemistry π-calculus Compartment Private communication channel Molecule Process Molecular Population Systems of communicating processes Biochemical Transformation Communication channel Compartment Translocation Extrusion of a private channel’s scope Initial processes: S0,1=gene;S0,2=gene |act |... |act and S0,3=gene |rep |... |rep Processes: gene := τc1.(gene|rna )+ac5?.act-gene +rc8?.rep-gene rna := τc2.(rna |protein )+τc3.0 protein := τc4.0 act := ac5!.0 act-gene := τc6.(act |gene)+τc7.(act-gene|rna ) rep := rc8!.0 rep-gene := τc9.(rep|gene) Fig. 4. π-calculus model of gene expression 2.3 Petri Nets Petri nets are a mathematical and computational tool for modeling and analysis of discrete event systems typically with a concurrent behavior. Petri nets offer a formal way to represent the structure of a discrete event system, simulate its behavior, and prove certain properties of the system. Petri nets have applications in many fields of system engineering and computer science. Here we only recall some basic concepts of Petri nets and refer to the current literature [9,23,24,25] for details regarding the theory and applications of Petri nets. In particular, we focus only on a specific class of Petri nets called place-transition nets or PT-nets, for short. Informally, a PT-net is a directed graph formed by two kinds of nodes called places and transitions respectively. Directed edges, called arcs, connect places to transitions, and transitions to places; each arc has associated a weight. Thus, for each transition, one identifies a set of input places, the places which have at least one arc directed to that transition, and a set of output places, the places which the outgoing arcs of that transitions are directed to. Then, a nonnegative integer number of tokens is assigned to each place; these numbers of tokens define the state of the PT-net also called the marking of the PT-net. In a PT-net, a transition is enabled when the number of tokens in each input place is greater than or equal to the weight of the arc connecting that place to the transition. An enabled transition can fire by consuming tokens from its input places and producing tokens in its output places; the number of tokens produced and consumed are determined by the weights of the arcs involved. The firing of a transition can be understood as the movement of tokens from some input places to some output places. More precisely, we give the following definition. Definition 2. APT-net is a construct N=(P,T,W, M0)where: Pis a finite set of places,Tis a finite set of transitions,withP∩T=∅,W:(P×T)∪ (T×P)→Nis the weight function,M0is a multiset over Pcalled the initial marking,andLis a location mapping. PT-nets are usually represented by diagrams where places are drawn as circles, transitions are drawn as squares, and an arc (x, y) is added between xand y if W(x, y)≥1. These arcs are then annotated with their weight if this is 2 or more. Given a PT-net N,thepre- and post-multiset of a transition tare respectively the multiset preN(t) and the multiset postN(t) such that, for all p∈P, |p|preN(t)=W(p, t)and|p|postN(t)=W(t, p). A configuration of N,whichis called a marking, is any multiset over P; in particular, for every p∈P,|p|Mrepresents the number of tokens present inside place p. A transition tis enabled at amarkingMif the multiset preN(t) is contained in the multiset M. An enabled transition tat marking Mcan fire and produce a new marking Msuch that M=M−preN(t)+postN(t) (i.e., for every place p∈P, the firing transition t consumes |p|preN(t)tokens and produces |p|postN(t)tokens). In order to reason about some basic properties, it is convenient to introduce a matrix-based representation for PT-nets. Specifically, let N=(P, T,W, M0)be a PT-net and let π:P→|P|and τ:T→|T|be two bijective functions. We call place jthe place pwith π(p)=j, and we call transition ithe transition twith τ(t)=i. Then, a marking Mis represented as a |P|vector which contains in each position jthe number of tokens currently present inside place j.Theincidence matrix of Nis the |T|×|P|matrix Asuch that, for every element aij of A, aij =|π−1(j)|postN(τ−1(i))|−|π−1(j)|preN(τ−1(i))|(i.e., aij denotes the change in the number of tokens in place jdue to the firing of transition i). A control vector uis a |T|vector containing 1 in position ito denote the firing of transition i,0 otherwise. Thus, if a particular marking Mnis reached from the initial marking M0through a firing sequence u1,u 2,...,u nof enabled transitions, we obtain Mn=M0+AT· n  k=1 uk which represents the reachable-marking equation. The aforementioned representation of a PT-net Nallows us to introduce the notions of P-invariants and T-invariants. P-invariants are the positive solutions of the equation A·y= 0; the non-zero entries of a solution yrepresents the set of places whose total number of tokens does not change with any firing sequence from M0. T-invariants instead are the positive solutions of the equation AT·x=0; a solution vector xrepresents the set of transitions which have to fire from some marking Mto return to the same marking M. Then, a PT-net is said to be bounded if there exists a |P|vector Bsuch that, for all marking Mreachable from M0,we have M≤B; a PT-net is said to be alive if, for all marking Mreachable from M0, there exists at least one transition enabled at marking M. As pointed out in [10,22], a PT-net model for a system of molecular interactions can be obtained by representing each molecular species as a different place and each biochemical transformation as a different transition. Tokens inside a place can then be used to indicate the presence of a molecule in certain proportions. This modeling approach is summarised in Table 3. Thus, a biochemical system is represented as a discrete event system whose structural properties are useful for drawing conclusions about the behavior and structure of the original biochemical system [22]. For instance, P-invariants determine the set of molecules whose total net concentrations remain unchanged during the application of certain biochemical transformations; T-invariants instead indicate the presence of cyclic reactions which lead to a condition where some reactions are in a state of continuous operation. Also, the property of liveness is useful to determine the absence of metabolic blocks which may hinder the progress of the biochemical system. Table 3. Modeling principles in PT-nets Biochemistry PT-net Molecule Place Molecular Population Marking Biochemical Transformation Transition Reactant Input Place Product Output Place Finally, we recall that it was shown in [4,15,16] how to transform a P system into a corresponding PT-net. This is done by considering a transition for each rule in the P system that has the left-hand side of the rule as pre-multiset and the right-hand side of the rule as post-multiset. In particular, in order to model the localization of rules and objects inside the membranes, one considers in the corresponding PT-net a distinct place for each object possibly present inside a membrane. Thus, the transformation of objects inside the membranes and the –For positive regulation, the number of activator molecules cannot be increased but can be decreased only by 1 - similar invariant is found by Daikon. –For negative regulation, the number of repressor molecules cannot be increased but can be decreased only by 1 - similar invariant is found by Daikon. PIPE also shows that the network is not bounded but it is alive. In fact, gene is always present and we can keep firing transition r1 to increase indefinitely the amount of rna. The liveness here comes from the above invariant and shows that the system will be working forever. The boundedness instead produces a result that apparently contradicts PRISM findings, where the probability that the number of rna’s is greater than 7 is almost 0! This comes from the fact that PIPE uses a non-deterministic system instead of a probabilistic one considered by PRISM and P system simulator. It will be interesting to check this property with a probabilistic Petri net tool. 6 Conclusions In this paper we have investigated the concerted use of different methods, and shown how these can provide complementary insights into different facets of biological system behavior. Individual modeling techniques have their own respective benefits and usually excel at reasoning about a system from a particular perspective. This paper shows how these benefits can be leveraged by using different modeling techniques in concert. As a case study, we have constructed a P system model of a small gene expression system and produced equivalent specifications using Petri net and π-calculus approaches.Simulations of the P system model were analyzed by Daikon (to identify potential rules that govern model output), and some of the most interesting suggested rules were checked using the PRISM probabilistic model checker. The Petri net model was analyzed with PIPE, a general Petri net analysis tool. The results show how analysis results from different models of the same system are useful for the purposes of both validating and improving each other. The gene expressionmodel was chosenbecause it is manageable,and thus forms a useful basis for a case study to compare different modeling techniques. Our future work will apply the techniques shown in this paper to a larger and more realistic case study. This should provide further insights into the benefits that arise in modeling increasingly complex systems when the modeler is increasingly reliant upon the use of various automated tools to study the model behavior. Acknowledgements The research of Francesco Bernardini is supported by NWO, Organization for Scientific Research of The Netherlands, project 635.100.006 “VIEWS”. The authors would like to thank the anonymous reviewers for their comments and suggestions that have helped to improve the paper. References 1. Bianco, L., Fontana, F.: Towards a Hybrid Metabolic Algorithm. In: Hoogeboom, H.J., P˘aun, G., Rozenberg, G., Salomaa, A. (eds.) WMC 2006. LNCS, vol. 4361, pp. 183–196. Springer, Heidelberg (2006) 2. Blundell, M., Kennell, D.: Evidence for Endonucleolytic Attack in Decay of Lac Messenger RNA in Escherichia Coli. J. Mol. Biol. 83, 143–161 (1974) 3. Calder, M., Vyshemirsky, V., Gilbert, D., Orton, R.: Analysis of Signalling Pathways Using Continuous Time Markov Chains. In: Priami, C., Plotkin, G. (eds.) Transactions on Computational Systems Biology VI. LNCS (LNBI), vol. 4220, pp. 44–67. Springer, Heidelberg (2006) 4. Dal Zilio, S., Formenti, E.: On the Dynamics of PB Systems: A Petri Net View. In: Mart´ın-Vide, C., Mauri, G., P˘aun, G., Rozenberg, G., Salomaa, A. (eds.) Membrane Computing. LNCS, vol. 2933, pp. 153–167. Springer, Heidelberg (2004) 5. Ernst, M., Cockrell, J., Griswold, W., Notkin, D.: Dynamically Discovering Likely Program Invariants to Support Program Evolution. IEEE Transactions on Software Engineering 27(2), 99–123 (2001) 6. Fontana, F., Bianco, L., Manca, V.: P Systems and the modeling of Biochemical Oscillations. In: Freund, R., P˘aun, G., Rozenberg, G., Salomaa, A. (eds.) WMC 2005. LNCS, vol. 3850, pp. 199–208. Springer, Heidelberg (2006) 7. Fontana, F., Manca, V.: Discrete Solutions of Differential Equations by Metabolic P Systems. Theoretical Computer Science 372(2-3), 165–182 (2007) 8. Gillespie, D.T.: Exact Stochastic Simulation of Coupled Chemical Reactions. The Journal of Physical Chemistry 81(25), 2340–2361 (1977) 9. Girault, C., Valk, R.: Petri Nets for Systems Engineering. Springer, Heidelberg (2003) 10. Goss, P.J.E., Peccoud, J.: Quantitative Modeling of Stochastic Systems in Molecular Biology using Stochastic Petri Nets. Proc. Natl. Acad. Sci. USA 95, 6750–6755 (1999) 11. Heath, J., Kwiatkowska, M., Norman, G., Parker, D., Tymchyshyn: Probabilistic Model Checking of Complex Biological Pathways. In: Priami, C. (ed.) CMSB 2006. LNCS (LNBI), vol. 4210, pp. 32–47. Springer, Heidelberg (2006) 12. Hinton, A., Kwiatkowska, M., Norman, G., Parker, D.: PRISM: A Tool for Automatic Verification of Probabilistic Systems. In: Hermanns, H., Palsberg, J. (eds.) TACAS 2006 and ETAPS 2006. LNCS, vol. 3920, pp. 441–444. Springer, Heidelberg (2006) 13. Hlavacek, S., Savageau, M.A.: Subunit Structure of Regulator Proteins Influences the Design of Gene Circuitry Analysis of Perfectly Coupled and Uncoupled Circuits. J. Mol. Biol. 248, 739–755 (1995) 14. Kennell, D., Riezman, H.: Transcription and Translation Initiation Frequencies of the Escherichia Coli Lac Operon. J. Mol. Biol. 114, 1–21 (1977) 15. Kleijn, K., Koutny, K.: Synchrony and Asynchrony in Membrane Systems. In: Hoogeboom, H.J., P˘aun, G., Rozenberg, G., Salomaa, A. (eds.) WMC 2006. LNCS, vol. 4361, pp. 66–85. Springer, Heidelberg (2006) 16. Kleijn, K., Koutny, M., Rozenberg, G.: Towards a Petri Net Semantics for Membrane Systems. In: Freund, R., P˘aun, G., Rozenberg, G., Salomaa, A. (eds.) WMC 2005. LNCS, vol. 3850, pp. 292–309. Springer, Heidelberg (2006) 17. Milner, R.: Communication and Mobile Systems: The π-calculus. Cambridge University Press, Cambridge (1999) 18. P˘aun, G.: Computing with Membranes. Journal of Computer and System Sciences 61(1), 108–143 (2000) 19. P´erez-Jim´enez, M.J., Romero-Campero, F.J.: P Systems, a New Computationl modeling Tool for Systems Biology. In: Priami, C., Plotkin, G. (eds.) Transactions on Computational Systems Biology VI. LNCS (LNBI), vol. 4220, pp. 176–197. Springer, Heidelberg (2006) 20. P´erez-Jim´enez, M.J., Romero-Campero, F.J.: Modeling Gene Expression Control Using P Systems: The Lac Operon, A Case Study (submitted, 2007) 21. Pescini, D., Besozzi, D., Mauri, C., Zandron, C.: Dynamical Probabilistic P systems. International Journal of Foundations of Computer Science 17(1), 183–204 (2007) 22. Reddy, V.N., Liebman, M.N., Mavrouniotis, M.L.: Qualitative Analysis of Biochemical Reaction Systems. Computers in Biology & Medicine 26(1), 9–24 (1996) 23. Reisig, W.: Elements of Distributed Algorithms, modeling and Analysis with Petri Nets. Springer, Heidelberg (1998) 24. Reisig, W., Rozenberg, G. (eds.): Lectures on Petri Nets I: Basic Models. LNCS, vol. 1491. Springer, Heidelberg (1998) 25. Reisig, W., Rozenberg, G. (eds.): Lectures on Petri Nets II: Applications. LNCS, vol. 1492. Springer, Heidelberg (1998) 26. Romero-Campero, F.J., Gheorghe, M., Ciobanu, G., Auld, J.M., P´erez-Jim´enez, M.J.: Cellular modeling Using P Systems and Process Algebra. Progress in Natural Science 17(4), 375–383 (2007) 27. Romero-Campero, F.J., P´erez-Jim´enez, M.J.: A Model of the Quorum Sensing System in Vibrio Fischeri Using P Systems (submitted, 2007) 28. Regev, A., Shapiro, E.: The π-calculus as an Abstraction for Biomolecular Systems. In: Ciobanu, G., Rozenberg, G. (eds.) Modeling in Molecular Biology, Springer, Heidelberg (2007) 29. Platform Independent Petri Net Editor: http://pipe2.sourceforge.net 30. P System Simulator: http://www.dcs.shef.ac.uk/∼marian/PSimulatorWeb/PSystemMF.htm APPENDIX 1 Transpose of the incidence matrix for PT-net of Figure 5: r1 r2 r4 r3 r5 r7 r6 r8 r9 gene 0 0 0 0 -1 0 1 -1 1 rna 100-101000 protein 0 1 -1 000000 act 0000-10100 act-gene 0 0 0 0 1 0 -1 0 0 rep 0000000-11 rep-gene 0 0 0 00001-1 which showsthe variations on the number of tokens determined by eachtransition. T-invariants obtained in PIPE: r110000 r201000 r401000 r310010 r500100 r700010 r600100 r800001 r900001 P-invariants computed by PIPE for some initial marking contains one token in gene n in act and min rep,withn, m ≥0: gene 1 0 0 rna 0 0 0 M(gene) + M(act-gene) + M(rep-gene) = 1 protein 0 0 0 M(act) + M(act-gene) = n act 0 1 0 M(rep) + M(rep-gene) = m act-gene 1 1 0 rep 0 0 1 rep-gene 1 0 1 APPENDIX 2 // Gene expression control // Model is stochastic stochastic // Bounds to the number of molecules const int rna_bound; const int protein_bound; const int number_activators; const int number_repressors; const int initact; const int initrep; // Stochastics constants associated with each command/rule/molecular interaction const double c1 = 0.347; // [ gene ]_b -c1-> [ gene + rna ]_b const double c2 = 0.174; // [ rna ]_b -c2-> [ rna + protein ]_b \\ const double c3 = 0.347; // [ rna ]_b -c3-> [ ]_b const double c4 = 0.0116; // [ protein ]_b -c4-> [ ]_b const double c5 = 6.6412087; // [ act + gene ]_b -c5-> [ actgene ]_b const double c6 = 0.6; // [ actgene ]_b -c6-> [ act + gene ]_b const double c7 = 3.47; // [ actgene ]_b -c7-> [ actgene + rna ]_b const double c8 = 6.6412087; // [ rep + gene ]_b -c8-> [ repgene ]_b const double c9 = 0.6; // [ repgene ]\_b -c9-> [ rep + gene ]_b // Module representing a bacterium module bacterium gene : [ 0 .. 1 ] init 1; actgene : [ 0 .. 1 ] init 0; repgene : [ 0 .. 1 ] init 0; act : [ 0 .. 1 ] init initact; rep : [ 0 .. 1 ] init initrep; rna : [ 0 .. rna_bound ] init 0; protein : [ 0 .. protein_bound ] init 0; // [ gene ]_b -c1-> [ gene + rna ]_b [ ] gene = 1 & rna < rna_bound -> c1 : (rna’ = rna + 1); // [ rna ]_b -c2-> [ rna + protein ]_b [ ] rna > 0 & protein < protein_bound -> c2*rna : (protein’ = protein + 1); // [ rna ]_b -c3-> [ ]_b [ ] rna > 0 -> c3*rna : (rna’ = rna - 1); // [ protein ]_b -c4-> [ ]_b [ ] protein > 0 -> c4*protein : (protein’ = protein - 1); // [ act + gene ]_b -c5-> [ actgene ]_b [ ] act = 1 & gene = 1 -> c5*number_activators : (gene’ = 0) & (act’ = 0) & (actgene’ = 1); // [ actgene ]_b -c6-> [ act + gene ]_b [ ] actgene = 1 & act = 0 -> c6 : (actgene’ = 0) & (act’ = 1) & (gene’ = 1); // [ actgene ]_b -c7-> [ actgene + rna ]_b [ ] actgene = 1 & rna < rna_bound -> c7 : (rna’ = rna + 1); // [ rep + gene ]_b -c8-> [ repgene ]_b [ ] rep = 1 & gene = 1 -> c8*number_repressors : (gene’ = 0) & (rep’ = 0) & (repgene’ = 1); // [ repgene ]_b -c9-> [ rep + gene ]_b [ ] repgene = 1 & rep = 0 -> c9 : (repgene’ = 0) & (rep’ = 1) & (gene’ = 1); endmodule APPENDIX 3 Ranges of molecules P=?[true U <= T rna > bound ] P=?[true U <= T protein > bound ] Constitutive regulation rna <= 7 rna >= 0 prot <= 32 prot >= 0 Positive regulation rna <= 24 rna >= 0 prot <= 205 prot >= 0 Negative regulation rna one of { 0, 1 } prot one of { 0, 1, 2, 3 } APPENDIX 4 Relationship between the number of repressors and rna and protein molecules. rna < rep rep > prot P=?[true U<=T rna > rep ] P=?[true U<=T protein > rep ] Expected number of molecules R=?[I=T] Other invariants P=?[true U gene = actgene ] ⇒Result: 0.0 P=?[true U gene = repgene ] ⇒Result: 0.0