Document Type : Original Research

Authors

1 PhD Candidate, Department of Biophysics & Biomedical Engineering, School of Medicine, Tehran University of Medical Sciences, Tehran, Iran

2 PhD Candidate, Research center for Science and Technology in Medicine (RCSTIM), Tehran University of Medical Sciences, Tehran, Iran

3 PhD Candidate, Department of Medical Genetics, School of Medicine, Mashhad University of Medical Science, Mashhad, Iran

4 PhD, Department of Biophysics & Biomedical Engineering, School of Medicine, Tehran University of Medical Sciences, Tehran, Iran

5 PhD, Research center for Science and Technology in Medicine (RCSTIM), Tehran University of Medical Sciences, Tehran, Iran

Abstract

Background: Interactions of many key proteins or genes in signalling pathway have been studied qualitatively in the literature, but only little quantitative information is available.
Objective: Although much has been done to clarify the biochemistry of transcriptional dynamics in signalling pathway, it remains difficult to find out and predict quantitative responses. The aim of this study is to construct a computational model of epidermal growth factor receptor (EGFR) signalling pathway as one of hallmarks of cancer so as to predict quantitative responses.
Material and Methods: In this analytical study, we presented a computational model to investigate EGFR signalling pathway. Interaction of Arsenic trioxide (ATO) with EGFR signalling pathway factors has been elicited by systematic search in data bases, as ATO is one of the mysterious chemotherapy agents that control EGFR expression in cancer. ATO has dichotomous manner in vivo, dependent on its concentration. According to fuzzy rules based upon qualitative knowledge and Petri Net, we can construct a quantitative model to describe ATO mechanism in EGFR signalling pathway.
Results: By Fuzzy Logic models that have the potential to trade with the loss of quantitative information on how different species interact, along with Petri net quantitatively describe the dynamics of EGFR signalling pathway. By this model the dynamic of different factors in EGFR signalling pathway is achieved.
Conclusion: The use of Fuzzy Logic and PNs in biological network modelling causes a deeper understanding and comprehensive analysis of the biological networks.

Keywords

Introduction

The major challenge in biology is to understand how interactions among genes, proteins or molecules in signal transduction pathways and gene regulatory networks regulate different functions of a cell and how a cell specifies its various fates, such as apoptosis (death), proliferation, differentiation and quiescence [ 1 ]. The answer is that interactions are inherently stochastic. Actually two kinds of stochasticity in a biological system exist: 1) inherent stochasticity in the biochemical procedures of gene expressions (intrinsic noise) and 2) variations in other cellular components (extrinsic noise). Genetic factors, regulatory dynamics and transcription rates control the amplitude of noises. Stochastic effects in gene expression play a crucial roles in biological procedures and determine different cell fates. There is much experimental evidence that proves stochasticity in gene expression s in biochemical procedures [ 2 , 3 ]. Although previous studies demonstrate the biochemistry of transcriptional dynamics in signaling pathways, it remains difficult to find out and predict quantitative responses. Integration of biological experiments with computational models can identify predictive models of signal transduction pathways and gene regulatory networks [ 4 , 5 ].

There are many tools and methods used for modeling and analysis of biological networks such as ordinary differential equation (ODE) [ 6 ], stochastic differential equation (SDE) [ 7 , 8 ], partial differential equation (PDE) [ 9 ], Markov chain [ 10 - 12 ], queuing networks [ 12 ], Boolean network (BN) [ 13 - 15 ], Probabilistic Boolean network [ 16 - 18 ], agent-based modeling (ABM) [ 19 - 23 ], Cellular automata (CA) [ 24 , 25 ], Fuzzy logic based models [ 26 ] and Petri net (PN) [ 27 - 29 ].

The complex interactions in gene regulatory networks or in signaling pathways exhibit characteristics which are hard to describe and implement mathematically using traditional tools like differential equations, while by computational methods such as CA, ABM, BN, fuzzy logic based models and PN, more details of interactions among agents can be modelled, without needing many kinetic parameters [ 30 , 31 ]. PN has been extended to handle more complexity in modeling of biological networks. PN is weighted, directed, bipartite multi graphs which includes places, transitions and arcs. In PN, enabled transitions (a reaction to be possible) fire (a reaction occur) immediately, while timed petri net (TPN) [ 32 , 33 ] regards a deterministic delay for firing of enabled transitions, and stochastic petri net (SPN) [ 29 , 34 ] applies a stochastic delay with the exponential distribution for firing enabled transitions. One of the ways for modeling stochasticity in biological networks is regarding stochastic delay for the occurrence of reactions [ 35 ]. Fuzzy petri net (FPN) [ 36 ] considers a linguistic value for a token such as low, medium and high one and defines the membership degree of token to each of the predefined fuzzy subsets, with a certain membership function. In Fuzzy stochastic petri net (FSPN) [ 29 ], the kinetic parameter related to exponential distribution (stochastic delay) is described by a fuzzy number and through regarding the different conditions of modeled system, the membership degree of kinetic parameters to each of the fuzzy subsets is being determined.

In this study, we briefly reviewed Fuzzy logic and Petri net respectively for biological networks modeling.

Fuzzy logic is capable to effectively capture systems, incomplete quantitative or qualitative information by characterizing the dynamic behavior of a system by a set of fuzzy rules. It has been used to simulate and build a model of a number of biological networks [ 4 , 26 ]. Fuzzy logic can model gene interactions at different levels of detail in the gene regulatory network across various conditions. Fuzzy model is constructed according to expert knowledge to quantify uncertainty, which is an inherent feature of biological networks.

Analysis and design of biological networks such as signaling pathways often include two kinds of uncertainty: fuzziness and randomness. Fuzziness models the measurement of imprecision due to incomplete information or linguistic structure. On the other hand, randomness models stochastic variability in interactions (or in rate of interaction) among agents (protein or gene). Fuzzy logic can encode probabilistic and dynamic transitions among network states so as to construct simple and realistic depiction of cell signaling pathways. Large body of evidences, turned the attention on fuzzy logic and fuzzy set theory which have been applied successfully in designing and modeling of many biological, and gene regulatory networks with regarding uncertainty and imprecision [ 4 , 37 ].

Petri net was introduced by Adam petri (1962), as a mathematical tool for modeling and analyzing of complex systems that can be characterized as simultaneous, parallel, synchronous, distributed, nondeterministic or stochastic.

Petri net is weighted, directed, bipartite multi graphs, which includes places, transitions and arcs. In models of gene regulatory networks or in signaling pathways, places represent chemical species, e.g., protein or protein complexes, gene, while transitions describe chemical reactions, e.g., transcription and translocation. The precursors and products of a chemical reaction correspond to the prior transition and the following one that they are input and output of transitions respectively. The interaction between two nodes (places) is mediated by one transition, that arcs connect transitions to places and places to transitions. Actually straight connection of two places or two transition is not allowed. The weight of an arc depicts its multiplicity, reflecting e.g., stoichiometric coefficient of a chemical reaction. Tokens describe the value of each place and are natural numbers or real positive numbers, reflecting e.g., the concentration of proteins or number of specific cells. A distribution of tokens over all places is called marking of PN, which describes a state of PN. In most PN, these states are finite [ 27 ]. Petri net provides state equations describing system behavior and model static and dynamic system characteristics. The static part of PN is its structure which is constructed from places, transitions and arcs, while the dynamic part models flow of tokens among places [ 38 ].

Material and Methods

In this analytical study, we presented a computational model to investigate EGFR signalling pathway. In this section, we introduce a pathway of key proteins in EGFR signalling and effect of ATO as medication on it.

1. ATO mechanisms and EGFR signaling

ATO is one of the mysterious chemotherapy agents with dichotomous manner in vivo depending on its concentration. Low concentration levels of ATO lead to cell proliferation and angiogenesis whereas high concentration induce cellular apoptosis [ 39 - 41 ].

1. 1. High concentrations of ATO

One of the crucial roles of ATO which is currently unraveled, is ATO-induced C-Src activation triggered EGFR-Y845/ERK phosphorylation and p21 expression. The major effect of ATO in cells is being modulated by NADPH oxidase and its products such as hydrogen peroxide, superoxide, and peroxynitrite. Activation of C-Src in response to ATO exerts this function in the cells. C-Src activity might results in phosphorylation and up regulation of NADPH oxidase component such as p67phox and p47phox that it elevating NADPH oxidase activity and generates ROS (reactive oxygen species). Furthermore C-Src can effectively transact by EGFR and effectively phosphorylates tyrosine residues such as Y845 which activating ERK1/2 [ 42 , 43 ] (Figure 1).

Figure 1. Signaling pathway of different concentrations of Arsenic trioxide (ATO). A) High concentration of ATO lead to P21 expression and induce cellular apoptosis. B) Low concentration of ATO lead to cell proliferation and angiogenesis.

Interestingly, various species of ROS can stimulate phosphorylation of EGFR. Altogether, these mechanisms lead to p21 expression. It’s astonishing that an oncogene acts as a tumor suppressor gene. This phenomena is known as soncogene-induced senescence, which controls activation of transformed oncogenes [ 44 ].

It has been demonstrated that ATO-induced sustained ERK1/2 phosphorylation can increase C-Fos protein stability to dimerize with c-Jun and form the AP-1 complex. This complex cooperates with p300/CBP. These actions result in chromatin remodeling, and consequently enhance p21 expression which leading to cellular apoptosis [ 39 - 41 ].

Moreover, it has been shown that EGFR can enhance DNA-damage repair such as non-homologous end-joining and homologous end-joining. As a matter of fact, inhibiting EGFR by erlotinib leads to down regulation of BRCA1, Rad50 and Rad 51 activity which results in enhancement of DNA damage based on ATO-induced ROS in ATO/erlotinib combination [ 45 ].

1. 2. Low concentration of ATO

In contrast, there is large body of evidence demonstrating, in low concentration of ATO, tumor growth and angiogenesis attributed to suppression of p21 expression [ 43 ]. Activation of JNK pathway by ATO phosphorylates c-Jun, and then through recruitment of TGIF/HDAC1 to sp1 binding sites of p21WAF1/CIP1 (p21) gene JNK pathway attenuates p21 expression [ 43 , 44 ]. In addition, in low concentration of ATO, FOXM1 (acting as a cell proliferation specific transcription factor) regulates the expression of CDC6, CDC25A, and Cyclin D1 through binging to its promoters and E2F/RB pathway. These genes by dephosphorylating and activation of CDK, phosphorylate FOXM1 that increased its activity. FOXM1 has an imperative role in the signal transduction network in low concentration of ATO [ 40 ].

In section 2 we suggest a method for modeling and analysis of dynamic networks which are based on Fuzzy Logic and PNs to exhibit both dimensions of uncertainty which are imprecision (fuzziness) and probabilistic (stochastic) variability.

2. EGFR model

EGFR signaling pathway model of the current paper is composed of two parts. First part, the relationship between key components of EGFR signaling pathway is modeled statically by fuzzy logic, considering ATO, P300/CBP and TGIF/HDAC (green circles of Figure 2) as normalized input of fuzzy logic model. The normalized concentration or activation levels of other components (blue circles of Figure 2) of Fuzzy logic model with fuzzy rules are resulted. More details about Fuzzy Logic model are illustrated in section 2.1.

Figure 2. Fuzzy logic with expert knowledge and fuzzy rule sets construct a quantitative model of epidermal growth factor receptor (EGFR) signaling pathway and Arsenic trioxide (ATO) as medication. This fuzzy logic model regard normalized values for each of the circles. This Fuzzy Logic model consists of three inputs: ATO, P300/CREB-binding protein (CBP) and TG-interacting factor (TGIF)/Histone deacetylase (HDAC) (green circles). The normalized value of blue circles is determined according to the value of green circles by fuzzy rule set. For each value of Forkhead box protein M1 (FOXM1) continuous petri net is executed to model inherent stochasticity in interactions between FOXM1, Cell division cycle 25A (DC25A), Cell division control protein 6 (CDC6), CyclinD1 and Cyclin-dependent kinase (CDK).

The second part of this model involves the interactions between FOXM1, CDC25A, CyclinD1, CDC6 and CDK during time into account. These interactions are inherently stochastic; therefore, continuous petri net is suitable for modeling this section. For each value of FOXM1 (as initial value of FOXM1 that is determined by Fuzzy logic model), continuous petri net is executed to model inherent stochasticity in interactions between FOXM1, CDC25A, CDC6, CyclinD1 and CDK and determine variations of concentration or activation levels of these proteins during time. More details about continuous petri net are illustrated in section 2. 2.

2. 1. Fuzzy Logic Model

According to expert knowledge about gene interactions in epidermal growth factor receptor (EGFR) signaling pathway and effects of ATO as medication in EGFR signaling, we constructed a fuzzy logic based model.

Fuzzy logic can model the interactions between major genes or molecules or components of EGFR signaling regarding the effect of different concentration of ATO. Therefore, without having kinetic data which are often incomplete or even unavailable, we can simulate the behavior of this pathway to achieve quantitative description of concentration with regarding uncertainty. Fuzzy theory deals with uncertainty associated with imprecision rather than with randomness. Actually many quantities such as protein concentration or its level of activation are not crisp deterministic quantities, and uncertainty due to imprecision and vagueness causes these variables to be fuzzy and could be represented by membership functions.

We use Mamdani-type fuzzy inference process that includes five stages as follow respectively:

Fuzzify input variables, apply fuzzy operator, apply implication method, apply aggregation method, defuzzification.

Stage 1: Fuzzify input variables.

The first stage is to take the input and specify the degree to which the input belongs to each of the suitable fuzzy sets via membership functions. In the following section, there is a definition of a fuzzy set θ on a universset X=μθ:X→[0,1] which assigns to each crisp point x∈X a real value μθ(X) ∈ [0,1].

μθ(X) is a membership function on fuzzy set θ. The membership function is a one-input/one-output function that determines the membership degree of each value of the variable χ to the predefined subsets of fuzzy set. There are different fuzzy membership functions such as: sigmoid, Gaussian curve, Bell curve, Pi-shaped curve, S-shaped curve, Trapezoidal, Triangular and Z-shaped curve.

Membership functions which are different in degree of freedom in the definition and implementation and cause different levels of performance and sensitivity in the system. We used triangular fuzzy membership function, denoted by θ=(a,b,c),a≤b≤c (1).

Stage 2: Apply fuzzy operator

Once the inputs have been fuzzified, the degree to which each section of the antecedent of every rule satisfied is determined. If the antecedent of each rule has more than one section, the fuzzy operator is applied to obtain one number which indicates the outcome of the antecedent for that rule. The fuzzified input variables include of two or more membership values.

Stage 3: Apply implication method

The input for the implication procedure is a single number determined by the antecedent, and its output is a fuzzy set. Implication is executed for each rule.

if <FP1> then <FP2> rule(R)

p = μFP1(x) q = μFP2(y) (2)

With regarding Mamdani implication method and Min for and method, the following formula is resulted:

μp→q(x,y)=min{p,q} (3)

Stage 4: Apply aggregation method

Aggregation is the procedure by which the fuzzy sets describing the outputs of each rule are merged into a single fuzzy set. The input of the aggregation procedure in our model is the list of truncated outcome functions returned by the implication procedure for each rule. The output of the aggregation procedure is one fuzzy set for each outcome variable. Aggregation only occurs once for each output variable and aggregation method is Max.

Stage 5: Defuzzification

Defuzzification is a mapping from fuzzy set to the crisp point. There are five built-in methods supported: bisector, centroid, middle of maximum (the average value of the maximum of the output set), smallest of maximum and largest of maximum. The defuzzification method implemented is centroid or center of gravity (Figures 3).

Figure 3. Adapted from Bree et al (46]. Fuzzy logic uses rule-based gates to quantify mechanisms that regulate network species. Fuzzy logic based models by expertized knowledge and fuzzy rule sets construct quantitative models of network species interactions.

2. 2. Continuous Petri net

The key protein FOXM1 in low concentration of ATO regulates the expression of CDC6, CDC25A, and CyclinD1 through binding to its promoters and E2F/RB pathway. These genes phosphorylate FOXM1 by dephosphorylating and activating of CDK. FOXM1 has an imperative role in the signal transduction network in low concentration of ATO. In this section for modeling interactions between: FOXM1, CDC6, CDC25A, CyclinD1 and CDK, we have constructed a continuous petri net (CPN). CPN with regarding nondeterministic firing rate for enabled transitions can model inherent stochasticity in the biochemical procedures of expression of CDC6, CDC25A, and Cyclin D1 by FOXM1 and regulation of CDK by CDC6, CDC25A and CyclinD1. Finally, dephosphorylation and activation of CDK phosphorylate FOXM1 as depicted in Figure 4. Therefore, CPN model variability in rate of interactions that this variability is made by intrinsic noise. Actually genetic factors, regulatory dynamics and transcription rates control the amplitude of noise. CPN captures randomness which is an inherent trait of biological networks.

Figure 4. Petri net (PN) model all interactions and dynamic of Forkhead box protein M1 (FOXM1), Cell division control protein 6 (CDC6), Cell division cycle 25A (DC25A), CyclinD1 and Cyclin-dependent kinase (CDK) with regarding initial value of FOXM1 by Fuzzy Logic model.

A CPN is a four-tuple V=<P,T,F,M0> where:

P is a finite set of places.

T is a finite set of transitions at which P∩T=∅.

F Pre ∪ Post ; Pre=P × T Post=T × P are finite set of directed arcs, which formerly describe arcs from places to transitions and later describe arcs from transitions to places.

M0: P → R≥0, represents the initial marking, which assigns a positive real number to each place, p ∈ P.

The enabling degree of transition t at marking m, enab(t,m) ∈ R≥0 ∪ ∞ is defined by:

{enab(t,m)=min{m[p]Pre[p,t]|Pt}enab(t,m)= ift= (4)

At which t represents input places of transition t.

t is enabled in m if enab(t,m) > 0.

An enabled transition t can be fired immediately by any values α ∈ R 0≤α≤ enab(t,m), and after firing of enabled transitions, marking defined by following rules all of p∈P:

m'[p]=m[p] + αC[p,t]

C=Post-Pre (5)

When a transition is fired, tokens immediately move from their pre-places to post-places according to α value and enabling degree, and updating the marking to a new reachable one. Our aim in this section is to construct a CPN for describing all interactions and dynamic of FOXM1, CDC6, CDC25A, CyclinD1 and CDK as depicted in Figure 4. This CPN consists of five places and five transitions, and Pre and Post matrices that describe the structure of interactions of all proteins as follow:

Pre=[0.30.30.30000010000100001000001]Post=[0000110000010000010000010] (6)

Actually Fuzzy Logic model according to different concentration of ATO regulate FOXM1. When Fuzzy logic model determines the concentration of FOXM1 (as initial value), then CPN models dynamic of FOXM1, CDC6, CDC25A, CyclinD1 and CDK during time.

Actually EGFR signaling pathway is composed of two part, static and dynamic part. At static part, relationships between key components of EGFR signaling pathway statically by fuzzy logic are modeled, and by regarding ATO, P300/CBP and TGIF/HDAC as inputs of fuzzy model, the concentration or activation level of all components of Figure 2 with fuzzy rules is resulted. The dynamic part of model is composed of CPN that describes the interactions between FOXM1, CDC25A, CDC6, CyclinD1 and CDK during time.

Results

In this work, we have applied Fuzzy Logic model and PN model to elucidate the interactions between some key components of EGFR signalling pathway. According to last studies low concentration of ATO through recruitment of TGIF/HDAC to sp1 binding sites of p21WAF1/CIP1 (p21) gene, JNK pathway attenuates p21 expression, and causes tumor growth and angiogenesis. Also ATO-induced sustained ERK1/2 phosphorylation can increases C-Fos protein stability to dimerize with c-Jun to form the AP-1 complex. This complex cooperates with p300/CBP. These actions result in chromatin remodeling, and consequent enhance p21 expression, which lead to cellular apoptosis. Therefore, with increment of P300/CBP and TGIF/HDAC, the value of P21 increased and decreased respectively as depicted in Figure 5. Also the dynamics of P21 with respect to dynamics of P300/CBP and ATO is computed with this model and depicted in Figure 6.

Figure 5. Normalized concentration or activation level of P21 related to normalized value of P300/CREB-binding protein (CBP) and TG-interacting factor (TGIF)/Histone deacetylase (HDAC). With increment of P300/CBP and increment of TGIF/HDAC, the value of P21 is increased and is decreased respectively.

Figure 6. Normalized concentration or activation levels of P21 related to normalized value of P300/CREB-binding protein (CBP) and Arsenic trioxide (ATO). As depicted, low concentrations of ATO suppressed p21 expression and high concentration of ATO lead to P21 expression and induced cellular apoptosis.

As mentioned previously, different concentration of ATO regulates different signalling pathways. The key protein FOXM1 in low concentration of ATO regulates the expression of CDC6, CDC25A, and Cyclin D1 and these genes by dephosphorylating and activating of CDK, phosphorylate FOXM1. Increment of ATO concentration inhibits upregulation of FOXM1. Figure 7 depicts dynamic of FOXM1, CDC25A, CDC6, CyclinD1 and CDK related to different concentration of ATO during time. As depicted in Figure 7 elevation of ATO decreases the value of FOXM1, consequently FOXM1 regulates the expression of other predefined proteins.

Figure 7. Fuzzy Logic model according different concentrations of Arsenic trioxide (ATO) determine the initial concentration or activation levels of Forehead box protein M1 (FOXM1). Average of dynamic of FOXM1, Cell division control protein 6 (CDC6), Cell division cycle 25A (DC25A), CyclinD1 and Cyclin-dependent kinase (CDK) during time is computed by Continuous petri net (CPN).

Discussion

In this study, we suggested a method for modeling and the analysis of dynamic networks which are based on Fuzzy Logic and PN to exhibit both dimensions of uncertainty which are imprecision (fuzziness) and probabilistic (stochastic) variability. This uncertainty is due to vague or missing kinetic data, incomplete data or naturally vary, e.g., between different experimental conditions, etc. Fuzzy Logic model is flexible, capable to incorporate qualitative and noisy data to describe the dynamics of signaling pathway of EGFR, and to produce quantitative predictions and novel biological insights about the function of signaling pathway. Along with Fuzzy Logic model, we have used PNs to model inherent stochasticity (noise) in the biochemical procedures of gene expressions. PNs model variability in rate of interactions that this variability is made by intrinsic noise. Actually genetic factors, regulatory dynamics and transcription rates control the amplitude of noise. Although the proposed method was applied to EGFR signaling pathway, it is not confined to modeling and analysis of this pathway. It can be utilized for modeling and the analysis of any signaling pathway or the biological network.

Model fitness can be conducted by calculating the sum of the squared difference between a model and the data. Fuzzy logic model parameters can be estimated if data on concentration of proteins during time is available. Therefore calibrated models with experimental data produce quantitative predictions and novel biological insights about the function of signaling pathways. Degree of fuzziness and membership function thresholds are parameters that could be estimated.

Conclusion

We believe that the major contribution of our survey, in addition to the offer of the use of Fuzzy Logic and PNs in biological network modeling, are a deeper understanding and the comprehensive analysis of the biological networks that can be achieved.

References

  1. Huang S, Ingber D E. Shape-dependent control of cell growth, differentiation, and apoptosis: switching between attractors in cell regulatory networks. Exp Cell Res. 2000; 261:91-103. DOI | PubMed
  2. Elowitz M B, Levine A J, Siggia E D, Swain P S. Stochastic gene expression in a single cell. Science. 2002; 297:1183-6. DOI | PubMed
  3. Blake W J, KAErn M, Cantor C R, Collins J J. Noise in eukaryotic gene expression. Nature. 2003; 422:633-7. DOI | PubMed
  4. Aldridge B B, Saez-Rodriguez J, Muhlich J L, Sorger P K, Lauffenburger D A. Fuzzy logic analysis of kinase pathway crosstalk in TNF/EGF/insulin-induced signaling. PLoS Comput Biol. 2009; 5:e1000340. Publisher Full Text | DOI | PubMed
  5. Hatakeyama M, Kimura S, Naka T, Kawasaki T, et al. A computational model on the modulation of mitogen-activated protein kinase (MAPK) and Akt pathways in heregulin-induced ErbB signalling. Biochem J. 2003; 373:451-63. Publisher Full Text | DOI | PubMed
  6. Chen T, He H L, Church GM. Modeling gene expression with differential equations. Pac Symp Biocomput. 1999;29-40. PubMed
  7. Chen K C, Wang T Y, Tseng H H, Huang C Y, Kao C Y. A stochastic differential equation model for quantifying transcriptional regulatory network in Saccharomyces cerevisiae. Bioinformatics. 2005; 21:2883-90. DOI | PubMed
  8. Climescu-Haulica A, Quirk M D. A stochastic differential equation model for transcriptional regulatory networks. BMC Bioinformatics. 2007; 8(Suppl 5):S4. Publisher Full Text | DOI | PubMed
  9. Baker R E, Gaffney E, Maini P. Partial differential equations for self-organization in cellular and developmental biology. Nonlinearity. 2008; 21:R251-90. DOI
  10. Cruz-Monteagudo M, Gonzalez-Diaz H, Aguero-Chapin G, et al. Computational chemistry development of a unified free energy Markov model for the distribution of 1300 chemicals to 38 different environmental or biological systems. J Comput Chem. 2007; 28:1909-23. DOI | PubMed
  11. Anderson D F, Kurtz T G. Continuous time Markov chain models for chemical reaction networks. Springer: Design and analysis of biomolecular circuits ; 2011.
  12. Bolch G, Greiner S, De Meer H, Trivedi K S. Queueing networks and Markov chains: modeling and performance evaluation with computer science applications. John Wiley &amp; Sons: New Jersey ; 2006.
  13. Li F, Long T, Lu Y, Ouyang Q, Tang C. The yeast cell-cycle network is robustly designed. Proc Natl Acad Sci U S A. 2004; 101:4781-6. Publisher Full Text | DOI | PubMed
  14. Davidich M I, Bornholdt S. Boolean network model predicts cell cycle sequence of fission yeast. PLoS One. 2008; 3:e1672. Publisher Full Text | DOI | PubMed
  15. Bornholdt S. Boolean network models of cellular regulation: prospects and limitations. J R Soc Interface. 2008; 5 Suppl 1:S85-94. Publisher Full Text | DOI | PubMed
  16. Shmulevich I, Dougherty E R, Kim S, Zhang W. Probabilistic Boolean Networks: a rule-based uncertainty model for gene regulatory networks. Bioinformatics. 2002; 18:261-74. PubMed
  17. Shmulevich I, Dougherty E R, Zhang W. From Boolean to probabilistic Boolean networks as models of genetic regulatory networks. Proceedings of the IEEE. 2002; 90:1778-92. DOI
  18. Liang J, Han J. Stochastic Boolean networks: an efficient approach to modeling gene regulatory networks. BMC Syst Biol. 2012; 6:113. Publisher Full Text | DOI | PubMed
  19. Figueredo G P, Siebers P O, Aickelin U. Investigating mathematical models of immuno-interactions with early-stage cancer under an agent-based modelling perspective. BMC Bioinformatics. 2013; 14 Suppl 6:S6. Publisher Full Text | DOI | PubMed
  20. Pogson M, Holcombe M, Smallwood R, Qwarnstrom E. Introducing spatial information into predictive NF-kappaB modelling--an agent-based approach. PLoS One. 2008; 3:e2367. Publisher Full Text | DOI | PubMed
  21. Gorochowski T E, Matyjaszkiewicz A, Todd T, et al. BSim: an agent-based tool for modeling bacterial populations in systems and synthetic biology. PLoS One. 2012; 7:e42790. Publisher Full Text | DOI | PubMed
  22. Figueredo G P, Siebers P O, Owen M R, Reps J, Aickelin U. Comparing stochastic differential equations and agent-based modelling and simulation for early-stage cancer. PLoS One. 2014; 9:e95150. Publisher Full Text | DOI | PubMed
  23. An G, Mi Q, Dutta-Moscato J, Vodovotz Y. Agent-based models in translational systems biology. Wiley Interdiscip Rev Syst Biol Med. 2009; 1:159-71. Publisher Full Text | DOI | PubMed
  24. Souza-e-Silva H, Savino W, Feijoo R A, Vasconcelos A T. A cellular automata-based mathematical model for thymocyte development. PLoS One. 2009; 4:e8233. Publisher Full Text | DOI | PubMed
  25. Monteagudo A, Santos J. Treatment Analysis in a Cancer Stem Cell Context Using a Tumor Growth Model Based on Cellular Automata. PLoS One. 2015; 10:e0132306. Publisher Full Text | DOI | PubMed
  26. Morris M K, Saez-Rodriguez J, Clarke D C, et al. Training signaling pathway maps to biochemical data with constrained fuzzy logic: quantitative analysis of liver cell responses to inflammatory stimuli. PLoS Comput Biol. 2011; 7:e1001099. Publisher Full Text | DOI | PubMed
  27. Hardy S, Robillard P N. Petri net-based method for the analysis of the dynamics of signal propagation in signaling pathways. Bioinformatics. 2008; 24:209-17. DOI | PubMed
  28. Ruths D, Muller M, Tseng J T, Nakhleh L, Ram P T. The signaling petri net-based simulator: a non-parametric strategy for characterizing the dynamics of cell-specific signaling networks. PLoS Comput Biol. 2008; 4:e1000005. Publisher Full Text | DOI | PubMed
  29. Liu F, Heiner M, Yang M. Fuzzy Stochastic Petri Nets for Modeling Biological Systems with Uncertain Kinetic Parameters. PLoS One. 2016; 11:e0149674. Publisher Full Text | DOI | PubMed
  30. Fisher J, Henzinger T A. Executable cell biology. Nat Biotechnol. 2007; 25:1239-49. DOI | PubMed
  31. Butcher J C. Numerical methods for ordinary differential equations. John Wiley &amp; Sons: New Jersey ; 2016.
  32. Popova-Zeugmann L. Time petri nets. Time and Petri nets. Springer: New York ; 2013.
  33. Andreychenko A, Magnin M, Inoue K. Modeling of Resilience Properties in Oscillatory Biological Systems using Parametric Time Petri Nets, Supplementary Information. ArXiv Preprint ArXiv:150606299. 2015.
  34. Heidary Z, Ghaisari J, Moein S, Naderi M, Gheisari Y. Stochastic Petri net modeling of hypoxia pathway predicts a novel incoherent feed-forward loop controlling sdf-1 expression in acute kidney injury. IEEE Transactions on Nanobioscience. 2016; 15:19-26. DOI
  35. Josic K, Lopez J M, Ott W, Shiau L, Bennett M R. Stochastic delay accelerates signaling in gene networks. PLoS Comput Biol. 2011; 7:e1002264. Publisher Full Text | DOI | PubMed
  36. Windhager L. Modeling of dynamic systems with Petri nets and fuzzy logic. LMU: Munich ; 2013.
  37. Bordon J, Moskon M, Zimic N, Mraz M. Fuzzy Logic as a Computational Tool for Quantitative Modelling of Biological Systems with Uncertain Kinetic Data. IEEE/ACM Trans Comput Biol Bioinform. 2015; 12:1199-205. DOI | PubMed
  38. Peterson J L. Petri net theory and the modeling of systems. 1981.
  39. Huang H S, Liu Z M, Hong D Y. Blockage of JNK pathway enhances arsenic trioxide-induced apoptosis in human keratinocytes. Toxicol Appl Pharmacol. 2010; 244:234-41. DOI | PubMed
  40. Liu Y, Hock J M, Van Beneden R J, Li X. Aberrant overexpression of FOXM1 transcription factor plays a critical role in lung carcinogenesis induced by low doses of arsenic. Mol Carcinog. 2014; 53:380-91. DOI | PubMed
  41. Liu Z M, Huang H S. Arsenic trioxide phosphorylates c-Fos to transactivate p21(WAF1/CIP1) expression. Toxicol Appl Pharmacol. 2008; 233:297-307. DOI | PubMed
  42. Lunghi P, Tabilio A, Lo-Coco F, et al. Arsenic trioxide (ATO) and MEK1 inhibition synergize to induce apoptosis in acute promyelocytic leukemia cells. Leukemia. 2005; 19:234-44. DOI | PubMed
  43. Sato K. Cellular functions regulated by phosphorylation of EGFR on Tyr845. Int J Mol Sci. 2013; 14:10761-90. Publisher Full Text | DOI | PubMed
  44. Tseng H Y, Liu Z M, Huang H S. NADPH oxidase-produced superoxide mediates EGFR transactivation by c-Src in arsenic trioxide-stimulated human keratinocytes. Arch Toxicol. 2012; 86:935-45. DOI | PubMed
  45. Kryeziu K, Jungwirth U, Hoda M A, Ferk F, et al. Synergistic anticancer activity of arsenic trioxide with erlotinib is based on inhibition of EGFR-mediated DNA double-strand break repair. Mol Cancer Ther. 2013; 12:1073-84. DOI | PubMed
  46. Aldridge B B, Saez-Rodriguez J, Muhlich J L, Sorger P K, Lauffenburger D A. Fuzzy logic analysis of kinase pathway crosstalk in TNF/EGF/insulin-induced signaling. PLoS Computational Biology. 2009; 5(4):e1000340.