A Study of the Robustness of the EGFR Signalling Cascade Using Continuous Membrane Systems

A Study of the Robustness of the EGFR Signalling Cascade Using Continuous Membrane Systems M.J. P´erez-Jim´enez and 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. Many approaches to anticancer treatment have had a limited success. A fundamental hurdle to cancer therapy is the robustness of the signalling networks involved in tumourgenesis. 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 model of the epidermal growth factor receptor signalling cascade is developed using continuous membrane systems. This model is used to study the robustness of this signalling cascade which is known to play a key role in tumour cell proliferation, angiogenesis and metastasis. Keywords: membrane computing, EGFR signalling network, signal transduction, robustness. 1 Introduction Membrane Computing is an emergent branch of Natural Computing introduced by Gh. P˘aun in [9]. Since then it has received important attention from the scientific community. In fact in 2003 the Institute for Scientific Information (ISI) has considered the seminal paper [9] as fast breaking and Membrane Computing has been selected as a fast emerging area in computer science. 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 trading time for space. P systems as a discrete model of computation have also been used to model biological phenomena (see the volume [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]. In this paper we use a continuous variant of P systems, different from that in [7], 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 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 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, features which are not easily achieved in models which use differential equations. The epidermal growth factor receptor (EGFR) belongs to the tyrosine kinase family of receptors. 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 concentration of 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. Disregulated 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 the EGFR signalling cascade is briefly described. We present our model using continuous P systems in section 4. Results and discussion are exposed in section 5. 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 use a variant whose systems can evolve 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 of a membrane structure, a hierarchically arranged set of membranes, where one places multisets of objects that represent the concentration of chemical substances. Usual P systems deal with discrete multisets over an alphabet Σbut here we work with continuous multisets (mappings from Σto R+, the set of non-negative real numbers). These multisets evolve according to a finite set of rules that represent chemical reactions. 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 (a rooted tree) consisting of nmembranes (nodes of the tree) 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: u[v]i→u0[v0]i, where u, v ∈Σ∗represent the reactants, u0, v0∈Σ∗represent the products, and i∈His the label of the relevant membrane of the reaction that is modelled. 6. Kis the rate of application function which associates with each rule and multiplicity of the objects in µ, a non-negative real number considered as the rate of application of the rule: K:R×Mn×m(R+)→R+, where Mn×m(R+) is the set of matrices of order n×mover R+. Aconfiguration 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 cj in 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 concentration 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)). 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. At an instant t∈R+, a rule r∈ R is applied exactly K(r, E(t)) times (in this sense, we can say that the rules are applied in a K-maximal way); that is, K(r, E(t)) units of the reactants are consumed and K(r, E(t)) units of the products are produced. Observe that the effect of the rule rdecreases the multiplicity (concentration) of its reactants and increases the multiplicity (concentration) 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. More formally, given an object cj∈Σ, 1 ≤j≤m, and a membrane i, 1≤i≤n, we denote by productioni(cj) (resp. consumptioni(cj)) the set of rules where cjis a product in membrane i(resp. a reactant). Therefore the real number (vi(t))(cj), denoted by |cj|i(t), is determined by the next formula: |cj|i(t) = |cj|i(0) + X r∈productioni(cj)Zt 0 K(r, E(s)) ds − −X r∈consumptioni(cj)Zt 0 K(r, E(s)) ds, with vi(0) = wi, that is, |cj|i(0) = wi(cj)). 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. As shown above, 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. 3 EGFR Signalling Cascade In this section we will describe briefly the EGFR signalling cascade following the network depicted below. 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 desphosphorylated. 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 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 several cellular proteins and nuclear transcription factors that we do not include in our model. 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. 4 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 c for the cytoplasm. •Initial Multisets: In the initial multisets we represent the initial concentrations of the chemical substances in the environment, the cell surface and the cytoplasm. These concentration has been obtained from the references listed in the bibliography. A detailed presentation of initial multisets is shown in table 2 of the Supplementary information. •Rules and Rate of application function: In the rules we model the chemical reactions described which form the signalling cascade. To model the reactions we use the Law of Mass Action which states that the rate 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+· · · +pk0, 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 5 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. 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. We study the effect of different EGF concentrations on the signalling cascade. To illustrate this effect we depict the evolution of the concentration of the most relevant proteins in the signalling cascade over time. 010 20 30 40 50 60 0.00 0.04 0.08 0.12 0.16 0.20 0.24 0.28 0.32 time (s) Concentration (nM) 100nM 200nM 300nM Above it can be seen that the receptor activation by autophophorylation is clearly concentration dependent showing a high peak in the first 5 seconds to decay rapidly afterwards to very low levels of concentration. 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. Next we show the evolution of the phosphorylation of Shc after binding to the receptor over 60 seconds. 010 20 30 40 50 60 0 2 4 6 8 10 12 14 16 time (s) Concentration (nM) 100nM 200nM 300nM It can be observed that the responses to EGF stimulation get more sustained as we get deeper in the cascade and that the network has almost managed to attenuate the overstimulation with EGF. Observe that the red and blue line are almost identical. The activation of the Ras protein is a key node in the EGFR signalling cascade. We depict its evolution over 180 seconds in the next graphic. 020 40 60 80 100 120 140 160 180 0 10 20 30 40 50 time (s) Concentration (nM) 100nM 200nM 300nM Note that at this point the response to the overstimulation with EGF has been completely attenuated and that an amplification of the response to the low concentration has been performed. Finally, the goal of the cascade is the activation of the mitogenic kinases MEK and ERK. These proteins regulates the transcription of several proto-oncogenes like c-fos and c-Jun. Next the evolution of phosphorylated MEK over 180 second is presented. 020 40 60 80 100 120 140 160 180 0 1 2 3 4 5 time (s) Concentration (nM) 100nM 200nM 300nM In this picture it is shown the surprising robustness of the signalling cascade. The signals from outside due to EGF concentration have been either attenuated or amplified to get the same concentration of the most relevant kinases. Observe that after 100 seconds, when the response gets sustained, the three line representing the response to different external EGF concentrations are identical. For more