Modelling EGFR signalling cascade using continuous membrane systems

Modelling EGFR signalling cascade using continuous membrane systems M.J. P´erez-Jim´enez; F.J. Romero-Campero Research Group on Natural Computing Department of Computer Science and Artificial Intelligence University of Sevilla Avda. Reina Mercedes s/n, 41012, Sevilla, Spain {marper, fran}@us.es Abstract. The complexity of networks of biological signalling pathways is such that the development of simplifying models is essential in trying to understand the wide-ranging cellular responses they can generate. In this paper a continuous variant of membrane systems is introduced and used to model the epidermal growth factor receptor signalling network which is known to play a key role in tumour cell proliferation, angiogenesis and metastasis. Keywords: membrane computing, EGFR signalling network, signal transduction 1 Introduction Membrane Computing is an emergent branch of Natural Computing introduced by G. P˘aun in [9]. Since then it has received important attention from the scientific community. In fact, Membrane Computing has been selected by the Institute for Scientific Information, USA, as a fast Emerging Research Front in Computer Science, and [8] was mentioned in [17] as a highly cited paper in October 2003. This new non-deterministic model of computation starts from the assumption that the processes taking place in the compartmental structure of a living cell can be interpreted as computations. The devices of this model are called P systems. Roughly speaking, a P system consists of a cell-like membrane structure, in the compartments of which one places multisets of objects which evolve according to given rules in a synchronous non-deterministic maximally parallel manner. Most variants of membrane systems have been proved to be computationally complete, that is equivalent in power to Turing machines, and computationally efficient, that is being able to solve computationally hard problems in polynomial time by trading time for space. P systems as a discrete model of computation have also been used to model biological phenomena (see the volume in [1]); and as a continuous model in [7]. A first formalization of non-discrete P system and a way to approximate them was introduced in [2]. 2 In this paper we introduce a continuous variant of P systems different from that in [7] and we use it to model the epidermal growth factor receptor (EGFR) signalling cascade. Up to now the usual mathematical formalization of biochemical signalling networks has been done using differential equations which are focused on the global description of the change in concentration of chemical compounds. Here we use a new formalization of these phenomena in a computational framework which focuses on the compartmental structure (membrane structure) of the cell and on the chemical reactions (rules) that take place in different regions of the cell. Thus this new approach focuses in the local interactions between the components of the system as so it makes possible a topological and modular modelling of intracellular signalling networks. In this framework expansion of an existing model is done by adding new rules (reactions) and, if it is necessary, new membranes to represent new regions (organella) of the cell; therefore the previous model does not need to be changed. Besides modularity and easy extensibility, in favour of our approach we also mention the easy understandability and programmability. All these features are not easily achieved in models which use differential equations. The epidermal growth factor receptor (EGFR) belongs to the tyrosine kinase family of receptors. The binding of the epidermal growth factor (EGF) to the extracellular domain of EGFR induces receptor dimerization and autophosphorylation of intracellular domains. Then a multitude of proteins are recruited starting a complex signalling cascade and the receptor follows a process of internalization and degradation in endosomals. Two principal pathways lead to activation of Ras-GTP by hydrolization of Ras-GDP. One of these pathways depends on the Src homology and collagen domain protein (Shc) and the other one is Shc-independent. Ras-GTP acts like a switch that stimulates the Mitogen Activated Protein (MAP) kinase cascade by phosphorylating the proteins Raf, MEK and ERK. Subsequently phosphorylated ERK regulates several cellular proteins and nuclear transcription factors. Deregulated EGFR expression, ligand production and signalling have been proved to have a strong association with tumourgenesis. As a result of this, EGFR has been identified as a key biological target for the development of novel anticancer therapies. The paper is organised as follows. Continuous P systems are introduced in the next section. In section 3 we discuss how continuous P systems can be approximated by discrete systems in order to implement them in computers; a description of the EGF signalling network is given in section 4. In section 5 the model of the EGF signalling network is presented. Some results and discussion are exposed in the next section. Finally, conclusions and future work are given in the last section. 2 Continuous P Systems Usual variants of P systems are discrete models of computation where in every step the rules are applied in a maximal way an integer number of times, we refer to [10] for details. Here we introduce a variant whose systems can evolve 3 in every instant applying a maximal set of rules a positive real number of times determined by a certain function K. This variant is 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. Roughly speaking, a continuous P system consists in a membrane structure, a hierarchically arranged set of membranes. More formally, a membrane structure is a rooted tree, where the nodes are called membranes, the root is called skin, and the leaves are called elementary membranes. Informally we can represent a membrane structure using Venn diagrams. In the membrane structure one places multisets of objects; usual P systems deal with discrete multisets but here we work with continuous multisets. A continuous multiset over an alphabet Σis a mapping from Σto R+. Next we give a formal definition of continuous P systems. A continuous P system is a construct, Π= (Σ, µ, w1,...,wn,R,K), where: 1. n ≥1 is the degree of the system (number of membranes); 2. Σ ={c1,...,cm}is the alphabet of objects; 3. µ is a membrane structure consisting of nmembranes labelled with 1,...,n (often, we identify the membranes with labels from a finite set H). 4. w1,...,wnare continuous multisets associated with each membrane of the membrane structure µ 5.Ris a finite set of rules of the form r≡(u, v, u′, v′, i) where u, v, u′, v′∈Σ∗, and 1 ≤i≤n. We represent a rule ras follows: u[v]i→u′[v′]i Notation: (r)1=u; (r)2=v; (r)3=u′; (r)4=v′and (r)5=i. 6.Kis the rate of application function which associates with each rule and multiplicity of the objects in µthe rate of application of the rule: K:R×Mn×m(R+)→R+ where Mn×m(R+) is the set of matrixes of order n×mover R+. An instantaneous configuration of a continuous P system Πis a matrix of Mn×m(R+) where the object in row iand column j(ai,j ) represents the multiplicity of the object cjin the membrane i. We interpret the configurations as assignments of continuous multisets to the membranes of the system, that is, the association of each region with the number of molecules of chemical substances present in it. For usual P systems we talk about computations but for continuous P systems we prefer to think of evolutions. An evolution of a continuous P system is a mapping from R+to Mn×m(R+). That is, an evolution Eassociates with each instant t∈R+an instantaneous configuration E(t) of the system: E(t) = (ai,j(t)) 1≤i≤n 1≤j≤m For each t∈R+and i, 1 ≤i≤n, we denote by vi(t) the continuous multisets over Σ={c1,...,cm}defined as follows: (vi(t))(cj) = aij(t) for 1 ≤j≤m. That is, we can describe E(t) by a tuple (v1(t),...,vn(t)). 4 In order to describe how to determine the evolution of a P system we need to define the relevant rules to a membrane. Given a continuous P system, Π= (Σ, µ, w1,...,wn,R,K), and a membrane i(1 ≤i≤n) we denote: Ri={r: (r)5=i}, R∗ i={r:f((r)5) = i}, where f((r)5) is the father of the membrane (r)5in the membrane structure µ. We say that the rules in Ri∪ R∗ iare the relevant rules to the membrane i. The way a continuous P system, Π= (Σ, µ, w1,...,wn,R,K), evolves is determined by the initial multisets w1,...,wnand the rate of application function K. We define the initial configuration of Πas the tuple (w1,...,wn). The rules are applied during the evolution of the system in a continuous way according to the rate of application function K. During an infinitesimal interval of time dt, a rule r∈ R is applied exactly K(r, E(t))dt times (in this sense, we can say that the rules are applied in a K-maximal way); that is, K(r, E(t))dt units of the reactants are consumed and K(r, E(t))dt units of the products are produced. Observe that the effect of the rule rdecreases the multiplicity (number of molecules) of its reactants and increases the multiplicity of its products. More precisely, we define the effect of a rule rduring an interval of time [t, T ] as follows: Ef(r, t, T ) = ZT t K(r, E(s)) ds. Therefore, given an object cj∈Σ, 1 ≤j≤m, and a membrane i, 1 ≤i≤n, (vi(t))(cj), denoted by |cj|i(t), is determined by the next formula: |cj|i(t) = |cj|i(0) + X r∈Ri∧cj∈alph((r)4) Ef(r, 0, t) + (1) +X r∈R∗ i∧cj∈alph((r)3) Ef(r, 0, t)−(2) −X r∈Ri∧cj∈alph((r)2) Ef(r, 0, t)−(3) −X r∈R∗ i∧cj∈alph((r)1) Ef(r, 0, t) (4) Observe that on the one hand the effect of the application of the rules in (1) and (2) increases |cj|i(t) because cjappears in the right-hand side of the rules (cjis a product) but on the other hand (3) and (4) decrease |cj|i(t) because cj appears in the left-hand side of the rules (cjis a reactant). 3 Approximating Continuous P Systems In computers, real numbers are represented by a finite set of rational numbers. Therefore, like in most continuous models we need to develop approximations in order to simulate evolutions of continuous P systems in computers. 5 As shown before, in order to determine the effect of a rule on the evolution of a system during an interval of time [t, T ] we only need to compute an integral of the rate of application function K. Hence, in order to approximate the evolution of a continuous P systems in a finite set of instants t0,··· , tqwe can use any suitable known numerical method to approximate integrals. Here for simplicity we use the rectangle rule; that is, we suppose tl+1 −tl=pis small enough to assume that Kremains constant and equal to K(r, E(tl)) in the interval [tl, tl+1] for l= 0,...,q −1. With this assumption we can approximate the effect of a rule during an interval of time of length pby Ef(r, tl, tl+1)≈pK(r, E(tl)). By doing this approximation we reach an usual P systems that performs qsteps (t0,...,tq) and in each steps the rules are applied pK(r, E(tl)) times. Therefore, we have approximated the evolution of a continuous P system by the computation of an usual discrete P system working in a pKbounded parallel manner. 4 EGFR Signalling Cascade In this section we will describe briefly the EGFR signalling cascade following the cascade depicted in the next page. During the signal transduction which takes place in this cascade, the information about the concentration of the EGF in the outside of the cell is translated into kinetic information inside the cell by EGFR phosphorylation. The epidermal growth factor receptor (EGFR) belongs to the tyrosine kinase family of receptors. The binding of the epidermal growth factor (EGF) to the extracellular domain of EGFR induces receptor dimerization and autophosphorylation of intracellular domains. Then, on the one hand, a multitude of proteins are recruited starting a complex signalling cascade and, on the other hand, the receptor follows a process of internalization, ubiquitination and degradation. In our model we consider two marginal pathways and two principal pathways starting from the phosphorylated receptor. In the first marginal pathway phospholipase C-γ(PLCγ) binds to the phospholyrated receptor, then it is phosphorylated (PLC∗ γ) and released into the cytoplasm where it can be translocated to the cell membrane or dephosphorylated. In the second marginal pathway the protein PI3K binds to the phospholyrated receptor, then it is phosphorylated (PI3K∗) and released into the cytoplasm where it regulates several proteins that we do not include in our model. Both principal pathways lead to activation of Ras-GTP. The first pathway does not depend on the concentration of the Src homology and collagen domain protein (Shc). This pathway consist of a cycle where the proteins growth factor receptor-binding protein 2 (Grb2) and Son of Sevenless homolog protein (SOS) bind to the phosphorylated receptor. Later the complex Grb2-SOS is released in the cytoplasm where it dissociates into Grb2 and SOS. In the other main pathway Shc plays a key role, it binds to the receptor and it is phosphorylated. Then either Shc∗is released in the cytoplasm or the proteins Grb2 and SOS binds to the receptor yielding a four protein complex (EGFR-EGF2*-Shc*-Grb2-SOS). Subsequently this complex dissociates into the 6 complexes Shc∗-Grb2-SOS, Shc∗-Grb2 and Grb2-SOS which in turn can also dissociate to produce the proteins Shc∗, Grb2 and SOS. Finally, Ras-GTP is activated by these two pathways and in turn it stimulates the Mitogen Activated Protein (MAP) kinase cascade by phosphorylating the proteins Raf, MEK and ERK. Subsequently phosphorylated ERK regulates 7 several cellular proteins and nuclear transcription factors that we do not include in our model. Note that there exist cross-talks between different parts and cycles of the signalling cascade which suggest a strong robustness of the system. For a more detailed description of the cascade see the literature listed in the bibliography. 5 Modelling EGFR Signalling Cascade by Continuous P Systems We have developed a model of the signalling cascade described in the previous section using a continuous P system, ΠEGF = (Σ, µ, we, ws, wc,R,K). Our model consists of more that 60 proteins and complexes of proteins and 160 chemical reactions. Supplementary information and details about the model are available on the web page www.gcn.us.es/egfr.pdf. •Alphabet: In the alphabet Σwe collect all the proteins and complexes of proteins that take part in the signalling cascade. In table 1 of the Supplementary information all the objects of the alphabet and the chemical compounds that they represent are listed. •Membrane Structure: In the EGFR signalling cascade described in the previous section, there are three relevant regions, namely the environment, the cell surface and the cytoplasm. We represent them in the membrane structure as the membranes labelled with: efor the environment, sfor the cell surface and cfor the cytoplasm. A Venn diagram and a rooted tree representing the membrane structure can be seen below: •Initial Multisets: In the initial multisets we represent the initial number of molecules of the chemical substances in the environment, the cell surface and the cytoplasm. These estimations has been obtained from the references listed in the bibliography. A detailed presentation of the initial multisets and the reference from where they were taken are shown in table 2 of the Supplementary information. •Rules and Rate of application function: In the rules we model the chemical reactions which form the signalling cascade. To model the reactions we use 8 the Law of Mass Action which states that the speed of a reaction is proportional to the product of the concentrations of the reactants. That is, if we have a reaction of the form: r1+···+rk→p1+···+pk′, then the rate of this reaction is k|r1|···|rn|, where kis called kinetic constant. In tables 3-11 of the Supplementary information (www.gcn.us.es/egfr.pdf) all the rules are listed as well as the kinetic constants and the references from where they were taken. As an example of the procedure we have followed to develop our model, we next present the derivation of one of the 160 rules. Let us consider the binding of EGF to EGFR: EGF EGFR →EGF-EGFR We know from biological experiments that EGF, which is present in the environment, binds to EGFR, which is present in the cell surface at a rate of 0.003 nM−1s−1. According to this, the relevant membrane in this reaction is the cell-surface because it separates the two regions involved in this reaction. Besides following the Mass Action Law the reaction takes place at a velocity of 0.003|EGF ||EGF R|. Therefore, in our model we represent this chemical reaction by the following rule and rate of application: EGF [EGFR]s→[EGF-EGFR]sK(r, E(t)) = 0.003|EGF (t)|e|EGF R(t)|s 6 Results and discussion The model presented in the previous section has been implementing using CLIPS, a productive development and delivery expert system tool which provides a complete environment for the construction of rule and/or object based expert systems. CLIPS uses the Rete Algorithm for the firing and application of rules, for details we refer to [16]. To implement our model we have approximated the evolutions of the continuous P system ΠEGF by computations of an usual P system working in a bounded parallel manner. The parameter p, chosen for the approximation, was fixed to 10−3after testing different values until the results obtained did not change. In this section we present the evolution of a number of key proteins in the EGFR signalling cascade to study two of its characteristics; namely its robustness and redundancy. Finally, we study the role of the number of receptors in the cell surface. First, we investigate the effect of different EGF concentrations on the signalling cascade. To illustrate this effect we depict the evolution of the number of molecules (nM) of the phosphorylated receptor and the complex EGFR-EGF2*- Shc*. It can be seen that the receptor autophophorylation is clearly concentration dependent showing a high peak in the first 5 seconds to decay rapidly afterwards to very low levels. 9 010 20 30 40 50 60 0.00 0.04 0.08 0.12 0.16 0.20 0.24 0.28 0.32 Receptor Autophosphorylation Time (s) Molecules (nM) 100nM 200nM 300nM 010 20 30 40 50 60 0 2 4 6 8 10 12 14 16 EGFR−EGF2*−Shc* Time (s) Molecules (nM) 100nM 200nM 300nM According to the variance in the receptor activation it is intuitive to expect different cell responses to different EGF concentrations. Here we will show that this is not the case. Two of the most important proteins in the cascade are the mitogenic kinases ERK and MEK. Next we present the evolution of these two phosphorylated proteins over 180 seconds in our model. 020 40 60 80 100 120 140 160 180 0.0e+000 2.0e-005 4.0e-005 6.0e-005 8.0e-005 1.0e-004 1.2e-004 1.4e-004 ERK Evolution Time (s) Molecules (nM) 100nM 200nM 300nM 020 40 60 80 100 120 140 160 180 0 1 2 3 4 5 MEK Evolution Time (s) Molecules (nM) 100nM 200nM 300nM Above it is shown the surprising robustness of the signalling cascade. The signals from outside due to different EGF concentration have been either attenuated or amplified to get the same concentration of the most relevant kinases in the signalling cascade. Note that after 100 seconds, when the response gets sustained, the three lines representing the response to different external EGF concentrations are identical. Currently, the robustness of the EGFR signalling cascade is proposed to be a product of receptor internalization and cross-talk between different pathways in the cascade. According to the literature listed in the bibliography receptor internalization produces signal attenuation protecting the cell from high external EGF concentrations, meanwhile amplification of the signal due to low EGF concentrations is performed in several nodes where there exists cross-talk between different pathways of the cascade. The results obtained using our model are in accordance with this hypothesis. Another important feature of this signalling cascade is its redundancy, there exist two different pathways leading to Ras-GTP activation. To study this feature