Simulating FAS-induced apoptosis by using P systems

Simulating FAS-Induced Apoptosis by Using P Systems Smitha Cheruku1, Andrei P˘aun1†, Francisco J. Romero-Campero2 Mario J. P´erez-Jim´enez2, Oscar H. Ibarra3 1Department of Computer Science/IfM, Louisiana Tech University, P.O. Box 10348, Ruston, LA 71272 Email: {sch025,apaun}@latech.edu 2Department of Computer Science and Artificial Intelligence, University of Seville, Avda. Reina Mercedes s/n 41012, Sevilla, Spain Email: {fran, marper}@us.es 3Department of Computer Science, University of California - Santa Barbara, Santa Barbara, CA 93106 Email: [email protected] Abstract In contrast to differential equations, P systems are an unconventional model of computation which takes into consideration the discrete character of the quantity of components and the inherent randomness that exists in biological phenomena. The key feature of P systems is their compartmentalised structure which represents the heterogeneity of the structural organisation of the cells, and where one can take into account the role played by membranes in the functioning of the system, for example signalling at the cell surface [13], selective uptake of substances from the media [14], diffusion across different compartments [15], etc. We show here that P systems can be a reliable tool for Systems Biology and could even outperform in some cases the current simulation techniques based on differential equations. We will also use a strategy based on the well known Gillespie algorithm [4] but running on more than one compartment called Multi-compartmental Gillespie Algorithm [13]. †To whom correspondence should be addressed. E-mail: [email protected] 1 1 Introduction Understanding biosignalling pathways is essential for designing effective therapeutic approaches to several important diseases. For example the FAS- induced apoptotic signalling pathway was shown to be one of the most relevant processes for understanding and combating cancer, AIDS and neurodegenerative diseases such as Parkinson’s disease, Alzheimer, etc. With several pathways unraveled in the last years, each one with its own unique structure and complexity, there is an increasing need to model these signalling cascades due to their complex nature. Because there is usually immense data collected for only one pathway, it is almost always hard to understand the pathway without the help of computer simulators. For better understanding of the FAS-induced apoptosis we are proposing a new way (actually two different approaches) of simulating the pathway by using P systems. A typical group of biosignalling pathways is known to lead to apoptosis (also known as programmed cell death). Apoptosis is a mechanism which helps the unwanted, injured, or improperly developed cells to commit suicide playing a fundamental role in the fight of the organism against cancers. Thus, aberrations in apoptotic responses to death signals contribute to cancer development, resistance to treatment, but also the reverse problem: autoimmune diseases. One major mechanism for inducing the apoptosis is through the activation of death receptors. Among the death receptors the signalling pathways for FAS-induced apoptosis are best characterised at the moment. We believe that computational modelling can play an important role to advance understanding about the complex signalling behavior of this pathway. Modelling FAS-induced apoptosis (or any biosignalling pathway) can be done in many ways, the traditional approach being at the moment the use of differential equations. It is argued that the use of differential equations is not the best approach for simulating processes that involve low number of molecules/objects as the ordinary differential equations (ODEs) are assuming large populations of molecules and are modelling the changes in the concentration/numbers of molecules of a particular species. For low numbers of molecules it is questionable if ODEs provide an accurate modelling [2]. A model, an abstraction of the real-world onto a mathematical/computational domain, highlights some key features while ignoring others that are assumed to be not relevant. A good model should have at least four properties: relevance, computability, understandability and extensibility [18]. A model must be relevant capturing the essential properties of the phenomenon 2 investigated; and computable so it can allow the simulation of its dynamic behavior, as well as the qualitative and quantitative reasoning about its properties. An understandable model will correspond well to the informal concepts and ideas of molecular biology. Finally, a good model should be extensible to higher levels of organisations, like tissues, organs, organism, etc, in which molecular systems play a key role. We believe that P systems posses all these properties. The approach based on differential equations is usually referred to as macroscopic chemistry since they model the average evolution of the concentration of chemical substances across the whole system. The microscopic approach considers the molecular dynamics for each single molecule involved in the system taking into account their positions, momenta of atoms, etc. This approach is computationally intractable because of the number of atoms involved, the time scale and the uncertainty in many of the cellular components. Our approach is referred to as mesoscopic chemistry [22]. Like in the microscopic approach one considers individual molecules like proteins, DNA and mRNA but ignores many other molecules like water and non-regulated parts of the cellular machinery. Besides this, the position and momenta of the molecules are also not modelled, instead one deals with the statistics of which reactions occur and how often. This approach is more tractable than microscopic chemistry but it provides a finer and better understanding than the macroscopic chemistry. Another observation is that the P system paradigm focuses on the compartmental structure that is exhibited by the cells; in each compartment one has different rules and objects, and the system moves from one configuration to the next one by obeying the rules and using only the objects available in each compartment. These are the features that we want to simulate in the signalling pathways, in this respect the use of P systems appears to be a natural approach. 2 P Systems using Deterministic Waiting Times Algorithm The P systems used for modelling FAS-induced apoptosis consider the quantity of components in the cascade to be discrete. We also have that the individual chemical reactions are asynchronous and occur discretely at different intervals of time. We reach this extension of P systems by considering that the rules can take various amounts of time (not only 1 time cycle as they 3 were defined in the original variants). We also have the rules in the system obeying the Law of Mass Action: the reaction rate depends proportionally on the product of the concentrations of the reactants. An incipient model of P systems of this type was reported last year in [8]. Where we have defined a system in which molecules can associate and dissociate (the current rules are of this form) and also rules can have different “application times”, in the sense that different rules can take different amounts of time to finish. The current model extends that model further by modelling the Law of Mass Action in the sense that the same rule rcould take 5 clock cycles to finish in a particular configuration of the system, or 10 clock cycles in another configuration of the system. The difference is in the number of the reactants for the reaction; this number is influencing the actual speed with which the reaction is applied (i.e. in the second configuration we would have approximately half of the number of reactants with respect to the first configuration). In this paper we work with a variant of P systems being a tuple Π = (O, L, µ, M1, M2,...,Mn, R1,...,Rn) , where: •Ois a finite alphabet of symbols representing objects (proteins and complexes of proteins); •Lis a finite alphabet of symbols representing labels for the compartments (membranes); •µis a membrane structure containing n≥1 membranes labelled with elements from L; •Mi= (wi, li), 1 ≤i≤n, are pairs which represent the initial configuration of membrane i:li∈Lis its label, and wi∈O∗is the initial multiset. •Ri, 1 ≤i≤n, are finite sets of rules associated with the membrane iwhich are of the form u[v]li→u′[v′]li, where u, v, u′, v′∈O∗are finite multisets of objects and liis the label of membrane i. Next, we discuss in more detail the rules that we will use in this paper, to model protein–protein interactions taking place in the compartmentalised structure of the living cell. 4 Transformation, complex formation and dissociation rules: (1) [ a]l→[b]l (2) [ a+b]l→[c]l (3) [ a]l→[b+c]lwhere a, b, c ∈Σ and l∈L. These rules are used to specify chemical reactions taking place inside a compartment of type l∈L, more specifically they represent the transformation of ainto b; the formation of a complex cfrom the interaction of aand b; and the dissociation of a complex ainto band c. Diffusing in and out: (4) [ a]l→a[ ]l (5) a[ ]l→[a]lwhere a∈Σ and l∈L. When chemical substances move or diffuse freely from one compartment to another we use these types of rules, where amoves from or to a compartment of type l. Binding and unbinding rules: (6) a[b]l→[c]l (7) [ a]l→b[c]lwhere a, b, c ∈Σ and l∈L. Using rules of the first type we can specify reactions consisting in the binding of a ligand swimming in one compartment to a receptor placed on the membrane surface of another compartment. The reverse reaction, unbinding of substance from a receptor, can be described using the second rule. Recruitment and releasing rules: (8) a[b]l→c[ ]l (9) c[ ]l→a[b]lwhere a, b, c ∈Σ and l∈L. With these rules we represent the interaction between two chemicals in different compartments whereby one of them is recruited from its compartment by a chemical from the other compartment, and then the new complex remains in the latter compartment. In a releasing rule a complex, c, located in one compartment can dissociate into a and b, remaining ain the same compartment as c, and bbeing released into the other compartment. Each reaction rhas associated a rate constant crcomputed using those in [7] according to [4]. The time needed for the reaction rto finish will 5 then be computed in the following way: wtr=1 cr∗|A|∗|B|where |A|and |B| represent the number of molecules of the two reactants A, B of the reaction r. On the other hand, if the reaction simulated is a dissociation or other first order reaction in A, the reaction time is computed as wtr=1 cr∗|A|. We give in the following the algorithm: 1. Calculate the WT for all the reactions in all the membranes call it ”waiting time” – ’T’. Select the reaction, r, with the shortest WT. 2. Update the waiting time for the rest of the reactions by subtracting ’T’ and add it to the simulation time (t = t + T). 3. Apply ronly once, thereby updating the multiplicities of respective symbols. 4. Recalculate the waiting times only for those reactions which are in the compartments affected by the applied rule. For each such reaction compare the new waiting time with the existing WT and keep the shortest one. 5. Sort all reactions with respect to their waiting times and again take the reaction with the lowest waiting time. 6. Repeat steps 2 to 5 until the desired simulating time is reached. We believe that this strategy represents the most natural way of defining the evolution of the P system. It is based on the fact that the rate of reaction is directly proportional to the product of concentrations of reactants. We calculate the velocity of reactions using this fact and, likewise the waiting time. We implemented in Java a simulator for P systems with these dynamics. It accepts as input a Systems Biology Markup Language (SBML) file containing the rules to be simulated and initial number of molecules in the system. We used Cell Designer package to generate the SBML source file for the reactions. 3 Results Obtained for Modelling FAS-Induced Apoptosis We have implemented all the rules described in [7] for the Type I and Type II of FAS-induced apoptotic pathway starting with the stimulation of FASL1 1FASL is an abbreviation for FAS ligand. 6 until the activation of the effector Caspase-3. In total we simulated 99 rules in the P system working on 53 proteins and protein complexes. The aforementioned paper also provides the experimental data obtained from the Jurkat cell line, which is a type II pathway dominant cell type. We compared our results with both the experimental data and with the simulation data reported in [7]. The comparison is given in Figure 1. The simulation data from [7] was obtained by simulating the rules using ordinary differential equations (ODE). The ODEs were solved numerically by the Entelos PhysioLab Modeler implementing a Runge-Kutta method, thus calculating the concentration changes for all the molecules in the system. One of the major proteins in the pathway, Caspase-3 was compared to the experimental data. In the Entelos simulation, Caspase-3 was activated at 4 hours, and it was considered close to the experimental results where it was obtained that it activated at 6 hours, see the Figure 1. Figure 1: Comparison between experimental data (top left), previous ODE simulation data (top right) and the new DWT simulation data (center). (Top left, right) pictures used with permission. The same pathway is modelled in the membrane computing framework using the same reactions and initial conditions. The Caspase-3 activation dynamics is studied when Bcl2 is at baseline value. Caspase-3 is activated in our simulator after about 7 hours which is a very good approximation of the experimental data and also the simulated pathway in [7]. 7 We pass now to performe all the simulation runs reported in [7]: several changes in the initial concentrations were considered such as increases in the initial concentration for Bcl2 of ten fold and one hundred fold, also increase in Caspase-8 of twenty fold. These simulations were run to see the sensitivity of the overall pathway to changes in the initial concentrations of the aforementioned proteins. Another change considered in the paper [7] was to check whether Bcl2 should be defined as binding to Bax, Bid, tBid2or both Bax and tBid. Several simulations with the respective changes in the set of reactions to model each assumption were considered in the original paper and we matched the same simulations. Modelling the Type I Pathway Behavior by Modifying the Model There are cells which are not sensitive to Bcl2 over expression as described in [20]. In these cells Caspase-3 is activated through the type I pathway, independent of the type II pathway. Scaffidi et al. [20] has suggested that the type of pathway is chosen based on the concentration of Caspase-8 generated in active form following FASL binding. If the Caspase-8 concentration is high, then the Caspase-3 is activated directly, on the other hand, if the concentration of Caspase-8 is low, the type II pathway is chosen so that the system is amplifying the death signal through the mitochondria to be able to induce the cell death. To test this hypothesis, the active Caspase-8 formation is increased by having the initial concentration of Caspase-8 set to a value 20 times greater than its baseline value while everything else was kept the same in the system. We performed the same simulation with the increase in Caspase-8 initial concentration, this resulted in faster Caspase-3 activation also in our simulation and is in good agreement with the previous paper. The Bcl2 concentration is also increased 100 times to test the sensitivity of Caspase-3 activation to Bcl2. The Figure 3 shows that the Caspase-3 activation is not sensitive to the increase in Bcl2 concentration, which is the hallmark for type I pathway dominant behavior. The consistency between the framework and the experimental results in the paper [7] validate the new simulation model. 2tBid is an abbreviation for truncated Bid. 8 Figure 2: Left the DWT simulation, right the ODE simulation for the change in Caspase-8 initial concentration. Bcl2’s Effect on Type II Pathway Bcl2 is an inhibitor of the type II pathway. We now analyse the Caspase-3 activation kinetics by considering different mechanisms through which Bcl2 can block the type II pathway. In [23], [12], or [1] the authors suggest that Bcl2 might bind with (a) Bax, (b) Bid, (c) tBid, or (d) bind to both Bax and tBid to block the mitochondrial pathway. We created four different models of P systems having the rules r1,...,r96, r97 for the case (a), r1,...,r95, r′ 96, r′ 97 for the case (b), r1,...,r95, r”′′96, r′′ 97 for the case (c) or the rules r1,...,r97, r98, r99 for the case (d). All the other rules remain the same for all the cases. We refer the interested reader to [24] for the details of the model. The dynamics of Caspase-3 activation is studied by varying the Bcl2 concentration 10 times or 100 times the baseline value. It was concluded that Bcl2 binding to both Bax and tBid is the most efficient mechanism for the pathway in comparison with the results obtained for the cases (a), (b) or (c). The same conclusions were obtained also after using our simulator for all the previous changes in the pathway. Due to space limitations, we do not include here all the results of the simulations, but only the case (d) as a comparison between the ODE simulator and our simulator. We will now describe the same simulation of the FAS-induced apoptosis by using a strategy that is based on the well known Gillespie’s algorithm but running on more than one compartment and so it will be called Multicompartmental Gillespie Algorithm [13]. 9