A Multiscale Modeling Framework Based on P Systems

A Multiscale Modeling Framework Based on P Systems Francisco Jos´e Romero-Campero1, Jamie Twycross1,2, Hongqing Cao1, Jonathan Blakes1, and Natalio Krasnogor1 1Automated Scheduling, Optimisation and Planning Research Group School of Computer Science, Jubilee Campus, University of Nottingham Nottingham NG8 1BB, United Kingdom 2Centre for Plant Integrative Biology Sutton Bonington Campus, University of Nottingham Nottingham LE12 5RD, United Kingdom {fxc,jpt,hxc,jvb,nxk}@cs.nott.ac.uk Abstract. Cellular systems present a highly complex organization at different scales including the molecular, cellular and colony levels. The complexity at each one of these levels is tightly interrelated. Integrative systems biology aims to obtain a deeper understanding of cellular systems by focusing on the systemic and systematic integration of the different levels of organization in cellular systems. The different approaches in cellular modeling within systems biology have been classified into mathematical and computational frameworks. Specifically, the methodology to develop computational models has been recently called executable biology since it produces executable algorithms whose computations resemble the evolution of cellular systems. In this work we present P systems as a multiscale modeling framework within executable biology. P system models explicitly specify the molecular, cellular and colony levels in cellular systems in a relevant and understandable manner. Molecular species and their structure are represented by objects or strings, compartmentalization is described using membrane structures and finally cellular colonies and tissues are modeled as a collection of interacting individual P systems. The interactions between the components of cellular systems are described using rewriting rules. These rules can in turn be grouped together into modules to characterize specific cellular processes. One of our current research lines focuses on the design of cell systems biology models exhibiting a prefixed behavior through the automatic assembly of these cellular modules. Our approach is equally applicable to synthetic as well as systems biology. 1 Introduction Models in systems biology has been recently classified according to their semantics into denotational and operational models [6]. Models with denotational semantics are the classical approach in modeling cellular systems which uses a set of equations to describe how the quantities of the different molecular species are related to each other over time. The classical example are ordinary and partial differential equations. In this case the behavior of the system is obtained by approximating numerically these equations. On the other hand, computational models have operational semantics which describe the behavior of the system using an algorithm or list of instructions that can be executed by an abstract machine. The models developed within this last framework has been termed recently executable biology [6]. In this case a more detailed description of the processes producing the behavior of the system is provided. Several formal computational approaches have been proposed to model cellular systems like Petri nets [10] and process algebra [19]. They mainly focus on system specification at the molecular level: membranes, compartmentalization and cellular colonies are seldom described. This fact makes it difficult to study multicellular systems whose function is determined by molecular interactions. Membrane computing is a branch of natural computing inspired directly from the structure and functioning of the living cell [14]. It has been applied to cellular modeling as one of the few computational frameworks which presents an integrative approach to multiscale systems ranging from the molecular to the multicellular level. Specifically, it represents the molecular interaction level of living cells using objects or strings and rewriting rules; the compartmental/cellular level using membranes; and the colony level using collections of membranes called membrane structures. The devices of this computational paradigm are referred to as Psystems. Although most research in P systems focuses on the study of the computational power of the different proposed variants, recently their application as a modeling formalism to cellular systems is emerging [3,4,7,11,17,20,21,15]. In this paper we discuss through a running example the use of P systems as a multiscale modeling framework for cell systems biology models. The paper is organized as follows. Stochastic P systems for cellular modeling are introduced in Section 2. Section 3 presents the running example used throughout this paper. The modeling principles in P systems are described in Section 4. Modularization in P systems is briefly discussed in Section 5. Finally, conclusions and future work are discussed in Section 6. 2 Stochastic P Systems The original strategy for the application of the rewriting rules in P systems was based on maximal parallelism and non-determinism [13]. This strategy does not represent the rate at which molecular interactions take place as every object that can evolve according to any rule must evolve in a single computation step, without taking into account that some molecular interactions are more frequent than others. Moreover, the real time evolution of cellular systems is not captured as all the computation steps are assumed to be of the same time length, neglecting the fact that some molecular interactions are faster than others. Different strategies for the application of the rewriting rules in P systems have been studied [5,8]. Specifically, a sequential stochastic strategy based on Gillespie’s theory of stochastic kinetics [9] was introduced in order to overcome the two previous problems when developing a modeling framework for cellular systems biology based on P systems [16]. Here we refer to this variant as stochastic P systems. Definition 1 (Stochastic P Systems). A Stochastic P system is a construct: Π=((Σobj,Σ str),L,μ,M l1,...,M lm,(Robj l1,R str l1),...,(Robj lm,R str lm)), where: •Σobj is a finite alphabet of objects representing molecular species whose internal structure is not relevant in the functioning of the system under study. •Σstr is a finite alphabet of objects representing relevant parts of some molecular species in the system. These objects are arranged into strings describing the structure of molecular species. •L={l1,...,l m}is a finite alphabet of symbols representing compartment labels used to identify compartment classes. Compartments with the same label share the same class, i.e., set of rewriting rules and initial multisets. •μis a membrane structure consisting of n≥1membranes defining compartments identified in a one to one manner with values from {1,...,n}and labeled with elements from L. •Mlt=(wt,s t),foreach1≤t≤m, is the initial state of the compartments from the class identified by label lt,wherewt∈Σ∗ obj is a finite multiset of individual objects and stis a finite set of strings over Σstr.Amultisetof objects, obj is represented as obj =o1+o2+...+opwith o1,...,o p∈Σobj. Strings are represented as follows s1·s2···sqwhere s1,...,s q∈Σstr. •Robj lt={robj,lt 1,...,r obj,lt kobj,lt},foreach1≤t≤m, is a finite multiset of rewriting rules on multisets of objects associated with compartments of the type specified by the label lt. The rewriting rules on multisets of objects are of the following form: robj,lt j:obj1[obj2]l cobj,lt j −→ obj 1[obj 2]l(1) with obj1,obj 2,obj 1,obj 2some finite multisets of objects from Σobj and lalabel from L. These rules are multiset rewriting rules that operate on both sides of membranes, that is, a multiset obj1placed outside a membrane labeled by landamultisetobj2placed inside the same membrane can be simultaneously replaced with a multiset obj 1andamultisetobj 2, respectively. Note that a constant cobj,lt jis associated specifically with each rule. This constant will be referred to as stochastic constant and is key to provide P systems with a stochastic extension as it will be used to compute the probability and time needed to apply each rule. This constant depends only on the physical properties of the molecules and compartments involved in the reaction described by the rule like temperature, pressure, pH, volume, etc. •Rstr lt={rstr,lt 1,...,r str,lt kstr,lt},foreach1≤t≤m, is a finite set of rewriting rules on multisets of strings and objects associated with compartments of the type defined by ltand of the following form: rstr,lt j:[obj +str ]l cstr,lt j −→ [obj+str;str 1+...+str s]l(2) with obj, objmultisets of objects over Σobj and str, str,str 1,...,str  sstrings over Σstr. These rules operate on both multisets of objects and strings. The objects obj are replaced by the objects obj. Simultaneously a substring str is replaced by strwhereas the strings str 1,...,str  sare produced to form part of the content of the compartment. In the same way as for rewriting rules on multisets of objects a stochastic constant cstr,lt jis associated with each rule. The previous definition is provided with a stochastic strategy for the application of the rewriting rules by extending the Gillespie algorithm to the multicompartmental structure of P systems. The resulting algorithm has been referred to as the Multicompartmental Gillespie Algorithm (MGA) [16]. The Gillespie algorithm [9] can only be applied directly in a single, fixed and well mixed volume. In our approach the first step consists of treating each compartment defined by a membrane as a fixed and well mixed volume where the rewriting rule to be applied and the elapsed time before its application is computed using the Gillespie Direct Method. Our algorithm then applies the corresponding rules following the order determined by these waiting times. After the application of each rule the algorithm recomputes the rules to be applied and the waiting times in the compartments affected by the application of the last rule using the Gillespie Direct Method. Finally, the MGA halts when a prefixed simulation time is reached or no further rules can be applied. 3 Running Example In order to illustrate our modeling framework we will use an abstract gene regulation system inspired from the functioning and structure of the lac operon in Escherichia coli (E. coli). This operon consists of three structural genes, lacZ, lacY and lacA, located sequentially on the genome and transcribed into one single mRNA. Their protein products are involved in the sensing, uptake and metabolism of lactose. The transcription of the lac operon is both positively and negatively regulated and it is considered a canonical example of gene transcription regulation in prokaryotes [18]. The linear structure of the lac operon (Figure 1) starts with a region called cap where the activator protein CRP binds and increases the rate of transcription. Following this site there is an operator sequence that we will refer to as op where the repressor protein LacI binds to stop transcription. The structural genes lacZ, lacY and lacA then follow. The first gene lacZ codifies the enzyme β-galactosidase involved in the metabolism of lactose by cleaving it into glucose and galactose; allolactose appears as a byproduct of this reaction. The protein product of the second gene lacY is a permease that associates to the cell membrane and acts as a pump transporting lactose into the cell. The function of the protein coded in the last gene lacA is not yet fully understood. The regulation of the lac operon allows E. coli to express the genes in the operon only when it is more beneficial for the cell. In the absence of lactose in the media the repressor LacI binds to the operator op preventing the structural ⇓ cap.op.lacZ.lacY.lacA Fig. 1. A schematic representation of the structure of the lactose operon (top) and its representation as a string (bottom) genes from being transcribed since they are not needed under these conditions. Nevertheless, occasionally the repressor drops from the operator producing a basal transcription of the operon. When lactose becomes available it starts to be transported inside the cell by the basal number of LacY proteins on the cell surface. Once in the cytoplasm it interacts with the basal number of β-galactosidase producing as a byproduct allolactose. Allolactose in turn binds to the repressor LacI and changes its state so it cannot bind to the operator allowing transcription of the structural genes. The resulting increase in production of LacY and β-galactosidase forms a positive feed-back loop increasing the number of allolactose molecules which interact with the repressors preventing premature termination of transcription. The lac operon is also under positive regulation by the protein CRP. This protein is activated by the glucose transport system and when active it binds to the cap site facilitating transcription. Even in the presence of lactose if glucose is present in the media CRP will not be active as the transport system will be occupied, pumping glucose into the cell. Therefore CRP will not bind to the operon to assist transcription. Only in the presence of lactose and absence of glucose will CRP be active and bound to the operon, producing the full transcription of the operon. This gene regulation system will be used in the following section as the running example illustrating our modeling principles. 4 Modeling Principles The complexity of cellular systems is organized into different levels ranging from the molecular to the cellular and colony scales. These levels of complexity are not independent instead they are tightly interrelated influencing each other directly. In this respect, stochastic P systems present an integrating multiscale modeling framework which explicitly specifies the molecular, cellular and colony levels in cellular systems in a relevant and understandable manner. One of our research lines consists of the development of integrative modeling principles within the modeling framework of stochastic P systems. More specifically we will present some ideas on how to describe molecular species, cellular regions and compartments, molecular interactions, gene expression control and cell colonies. Our running example will be used to illustrate our modeling principles. •Molecular species: These are specified as individual objects or strings of objects. Molecules with an internal structure that is relevant in the functioning of the system are specified using strings. For example, gene operons with a linear structure consisting of promoters, operators, transcription/translation starting points, etc, otherwise molecular species are described using individual objects. Table 1. Specification of the molecular species in the lac operon Molecular Species Object RNA Polymerase RNAP Ribosome Rib Repressor LacI Activator CRP CRP∗ LacZ product LacZ LacY product LacY LacA product LacA Lactose Lact Allolactose Allolac Glucose Gluc Glucose transport system Gluc Complex glucose transport system Gluc-GT S Complex lactose LacY product Lact-LacY Complex lactose LacZ product Lact-LacZ Complex lactose LacZ product Lact-LacZ Complex allolactose repressor Allolac-LacI Operon site Object Activator binding site cap Occupied activator binding site capCRP∗ Repressor binding site op Occupied repressor binding site opLacI lacZ gene lacZ lacY gene lacY lacA gene lacA lacZ mRNA mlacZ lacY mRNA mlacY lacA mRNA mlacA Running example: The different molecular species in our example will be specified according to this modeling principle. On the one hand, the proteins and complexes of proteins involved in the regulation and expression of the lac operon are specified as individual objects since we are not interested in their internal structure (Table 1). On the other hand, each component of the lac operon will be described using an object such that the lac operon structure is specified as a string containing these objects in the specific order they can be found in E. coli’s genome (Figure 1). •Membranes: Compartmentalization and membranes are fundamental in the structural organization and functioning of living cells. Membranes do not act as passive boundaries of cells and compartments; instead they play a key role in the regulation of the metabolism and information processing between the outside and the inside of compartments. P systems constitutes one of the few computational frameworks which explicitly specifies compartments and membranes. For instance, P systems have been used to study selective uptake of molecules from the environment [20], signalling at the cell surface [12] and colonies of interacting bacteria which communicate by sending and receiving diffusing signals [2,21]. In general P system membranes are used to define relevant regions in cellular systems and therefore they do not always correspond to real cell membranes although normally they do. Running example: In the lac operon gene regulation system there are two relevant regions. Namely, the bacterium surface where LacY and GT S act as pumps transporting lactose and glucose into the cell and the aqueous interior or cytoplasm where the operon is located together with the different transcription factors and proteins. These two regions are represented using two membranes embedded one inside the other to described the structure of an E. coli bacterium (Figure 2). Fig. 2. Graphical representation of the membrane structure specifying an E. coli bacterium •Molecular processes consisting of protein-protein interactions and protein translocation: Such processes are normally described in P systems using rewriting rules on multisets of objects. Our P system modeling framework aims at providing a comprehensive and relevant rule-based schema for the most common molecular interactions taking place in living cells. More specifically, our approach focuses on the transformation and degradation of molecular species, the formation and dissociation of complexes, and the basic processes of communication and transport between different compartments in cellular systems (Table 2). Running example: The protein-protein interactions in our gene regulation system are described using the rewriting rules on multisets of objects presented in Table 3. Rules r29,r 30,r 31 and r32 are examples of complex formation and dissociation rules. The degradation and dilution of different proteins Table 2. P system rule-based schemas for the most common molecular interactions Molecular Interaction PSystemRules Transformation and Degradation [a]l c −→ [b]l[a]l c −→ [] l Complex formation and dissociation [a+b]l cf −→ [c]l[c]l cd −→ [a+b]l Diffusion in and out a[] l cin −→ [a]l[a]l cout −→ a[] l Binding and debinding a[b]l clb −→ [c]l[c]l cld −→ a[b]l Recruitment and releasing a[b]l crt −→ c[] lc[] l crl −→ a[b]l is specified in rules r22,r 23 and r24. Finally, active uptake of glucose and lactose are modeled using the binding and releasing rules r27,r 28,r 33 and r34. •Gene expression control: The sensing of signals and the processing of the information they convey is performed in living cells through molecular interactions of the type presented in Table 2. The response of cells to these signals consists of the expression of appropriate proteins codified in specific genes. Gene expression control has been described in P systems using either rewriting rules on multisets of objects or rewriting rules on multisets of objects and strings according to the structural organization of the genes in the system under study. Tables 4 and 5 presents these two alternatives for the specification of the most important processes in gene expression control; transcription factor binding and debinding, transcription and translation. From a simplistic point of view the processes involved in transcription factor binding and debinding, transcription and translation can be represented by individual rewriting rules on multisets of objects (Table 4). Nevertheless, these processes are very complex and they consist of different stages like operator/promoter recognition by transcription factors and RNA polymerase, transcription/translation initiation/termination, elongation, etc. A more accurate and detailed description of all these processes is achieved by using rewriting rules on multisets of strings and objects of the form of the rules in Table 5. Running Example: The gene regulation control in the lac operon is modeled using the rewriting rules on multisets of objects and strings given in Table 3. More specifically, the binding and debinding of the activator and repressor to their corresponding binding sites is represented using rules r3,r 4,r 7and r8. Transcription initiation in the presence and absence of the promoter site occupied by the activator CRP∗is specified using rules r1,r 2,r 5and r6. The transcription of the structural genes lacZ,lacY and lacA is described by the rules r9,r 10,r 11 and r12. Finally, translation and mRNA degradation is modeled with the rules r13 -r21. •Cell colonies: The last level of organization that has been represented using P systems consists of cellular systems where cells form colonies by interacting Table 3. Lac Operon Regulation Rules Nr. Rule Stochastic Constant r1:[RNAP +cap]b c1 −→ [cap.RNAP ]bc1=5×10−3min−1 r2:[cap.RNAP ]b c2 −→ [RNAP +cap]bc2=1min−1 r3:[CRP∗+cap]b c3 −→ [capCRP∗]bc3=16.6min−1 r4:[capCRP∗]b c4 −→ [CRP∗+cap]bc4=10min−1 r5:[RNAP +capCRP∗]b c5 −→ [capCRP∗.RNAP ]bc5=0.2min−1 r6:[capCRP∗.RNAP ]b c6 −→ [RNAP +capCRP∗]bc6=1min−1 r7:[LacI +op]b c7 −→ [opLacI ]bc7= 166min−1 r8:[opLacI ]b c8 −→ [LacI +op]bc8=0.1min−1 r9:[RNAP.op]b c9 −→ [op.RNAP ]bc9=3min−1 r10 :[RNAP.lacZ]b c10 −→ [lacZ.RNAP;mlacZ]bc10 =0.78min−1 r11 :[RNAP.lacY ]b c11 −→ [lacY.RNAP;mlacY ]bc11 =1.92min−1 r12 :[RNAP.lacA]b c12 −→ [RNAP +lacA;mlacA]bc12 =4min−1 r13 :[Rib +mlacZ]b c13 −→ [Rib.mlacZ]bc13 =0.12min−1 r14 :[Rib +mlacY ]b c14 −→ [Rib.mlacY ]bc14 =0.12min−1 r15 :[Rib +mlacA]b c15 −→ [Rib.mlacA]bc15 =0.12min−1 r16 :[Rib.mlacZ]b c16 −→ [Rib +LacZ +mlacZ]bc16 =0.12min−1 r17 :[Rib.mlacY ]b c17 −→ [Rib +LacY +mlacY ]bc17 =1.73min−1 r18 :[Rib.mlacA]b c18 −→ [Rib +LacA +mlacA]bc18 =3.55min−1 r19 :[mlacZ]b c19 −→ [] bc19 =6×10−3min−1 r20 :[mlacY ]b c20 −→ [] bc20 =6×10−3min−1 r21 :[mlacA]b c21 −→ [] bc21 =6×10−3min−1 r22 :[LacZ ]b c22 −→ [] bc22 =6.9×10−2min−1 r23 :[LacY ]b c23 −→ [] bc23 =6.9×10−2min−1 r24 :[LacA ]b c24 −→ [] bc24 =6.9×10−2min−1 r25 :[LacY ]b c25 −→ LacY [] bc25 =1min−1 r26 :LacY [] b c26 −→ [LacY ]bc26 =0.7min−1 r27 :Lact [LacY ]s c27 −→ [Lact-LacY ]sc27 =10min−1 r28 :Lact-LacY [] b c28 −→ LacY [Lact ]bc28 =10min−1 r29 :[Lact +LacZ ]b c29 −→ [Lact-LacZ ]bc29 =10min−1 r30 :[Lact-LacZ ]b c30 −→ [Allolac +LacZ ]bc30 =10min−1 r31 :[Allolac +LacI ]b c31 −→ [Allolac-LacI ]bc31 =1min−1 r32 :[Allolac-LacI ]b c32 −→ [Allolac +LacI ]bc32 =10 −4min−1 r33 :Gluc [GT S ]s c33 −→ [Gluc-GT S ]sc33 =1min−1 r34 :Gluc-GT S [] b c34 −→ GT S [Gluc ]bc34 =10min−1 r35 :GT S [CRP ]b c6 −→ GT S [CRP∗]bc35 =6.9×10−3min−1 r36 :[CRP∗]b c6 −→ [] bc36 =0.069min−1 and exhibiting coordinated behavior. The specification of the environment where the colony of cell is located cannot always be described by a single membrane since it is normally too big to be considered a well mixed volume or region where the Gillespie Algorithm can be applied. In this respect,