Journal Pre-proof
Using activity time windows and logical representation to reduce the complexity of
biological network models: G1/S checkpoint pathway with DNA damage
Ali Abroudi, Sandhya Samarasinghe, Don Kulasiri
PII:
S0303-2647(20)30034-4
DOI:
https://doi.org/10.1016/j.biosystems.2020.104128
Reference:
BIO 104128
To appear in:
Bio Systems
Received Date: 4 November 2019
Revised Date:
25 February 2020
Accepted Date: 25 February 2020
Please cite this article as: A. Abroudi, S. Samarasinghe, D. Kulasiri, Using activity time windows and
logical representation to reduce the complexity of biological network models: G1/S checkpoint pathway
with DNA damage, Bio Systems (2020), doi: https://doi.org/10.1016/j.biosystems.2020.104128.
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition
of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of
record. This version will undergo additional copyediting, typesetting and review before it is published
in its final form, but we are providing this version to give early visibility of the article. Please note that,
during the production process, errors may be discovered which could affect the content, and all legal
disclaimers that apply to the journal pertain.
© 2020 Published by Elsevier B.V.
Using Activity Time Windows and Logical Representation to Reduce
the Complexity of biological network models: G1/S Checkpoint
Pathway with DNA Damage
Abstract
Biological systems are difficult to understand complex systems. Scientists continue to create models
to simulate biological systems but these models are complex too; for this reason, new reduction
methods to simplify complex biological models into simpler ones are increasingly needed. In this
paper, we present a way of reducing complex quantitative (continuous) models into logical models
based on time windows of system activity and logical (Boolean) models. Time windows were used to
define slow and fast activity areas. We use the proposed approach to reduce a continuous ODE
model into a logical model describing the G1/S checkpoint with and without DNA damage as a case
study. We show that the temporal unfolding of this signaling system can be broken down into three
time windows where only two display high level of activity and the other shows little or no activity.
The two active windows represent a cell committing to cell cycle and making the G1/S transition,
respectively, the two most important high level functions of cell cycle in the G1 phase. Therefore, we
developed two models to represent these time windows to reduce time complexity and used
Boolean approach to reduce interaction complexity in the ODE model. The developed reduced
models correctly produced the commitment to cell cycle and G1/S transfer through the expected
behavior of signaling molecules involved in these processes. As most biological models have a large
number of fast reactions and a relatively smaller number of slow reactions, we believe that the
proposed approach could be suitable for representing many, if not all biological signaling networks.
The approach presented in this study greatly helps in simplifying complex continuous models (ODE
models) into simpler models. Moreover, it will also assist scientists build models concentrating on
understanding and representing system behavior rather than setting values for a large number of
kinetic parameters.
Keywords: Model Complexity; Model reduction; Time windows; ODEs models; Boolean models; G1/S
checkpoint; DNA damage.
1 Introduction
Biological networks such as molecular networks are complex to understand and model due to the large
number of entities interacting in complex ways to perform various functions in a largely intractable
1
manner. Experiments currently produce a large amount of a quantitative data and many quantitative
models have been developed for biological networks; but most of these models are very complex. They
are largely Ordinary Differential Equations (ODE) that also require the knowledge of a large number of
kinetic parameters of reactions. Therefore, an important current issue is how to simplify complex
quantitative models that can then be used to predict and explain results of biological experiments. This
paper simplifies and combines two complexity dimensions of biological networks - time and
interactions - to reduce complex mathematical models to simpler ones.
A particular form of reduction of interaction complexity that has already been exploited arises from
the fact that most of our understanding of regulatory and signaling networks remains qualitative in
nature. This makes it sensible to represent these networks by logical models, such as Boolean
networks, where the state of a molecule is either 0 or 1 (inactive or active) respectively. These models
have many benefits; they are simpler in nature and easier to develp, run and simulate perturbations
(Le Novere, 2015). Further, they do not need values for kinetic parameters and are capable of
capturing the vital behaviour of a system qualitatively. Therefore, they are extremely valuable for
addressing problems where the system is complex and data are lacking as most continous models
require greater level of system knowledge and data. However, there is much scope for reduction of
time complexity in these models as well as in other models. Reduction of time complexity is afforded
by the fact that in a biological network, most interactions are relatively fast involving phosphorylation
and activation relations and a smaller number of interactions are much slower involving gene
expression and protein synthesis. The interplay of such interactions can make networks produce
bursts of activity in a particular period followed by periods where activity levels are low or minimum.
The slow-fast nature along the time axis can allow reduction of time complexity of models. For
example, an ODE model needs to be solved continuously for the whole period of simulation.
However, if some stable or low activity regions exist where not much new information is produced,
these regions can be excluded in order to focus on active regions without loss of important
information while reducing the solution time. Further, as these active regions are in the temporal
order they appear in the process, it indirectly introduces time into Boolean models which is a
chanllenge at the moment for complex biological systems. Therefore, combination of the two
dimensions of complexity reduction can lead to greater simplification of complex models.
In this paper, we present an approach to simplify complex continuous models into logical models, by
simplifying time complexity and interaction complexity (by removing the kinetic parameters from a
system) through activity time windows and logical representation. The approach is standardized and
can easily be applied to other complex continuous models. In a case study, a quantitative ODE model
is used to produce a reduced logical model that describes the G1/S checkpoint with and without DNA
damage. We address how this model can explain and predict the behavior of G1/S checkpoint with and
2
without DNA damage, including some molecular oscillations and cell fate. This demonstrates that the
reduced model is still helpful to gain insights into biology and is simpler to operate and analyze.
This paper is organised as follows: the problem is introduced in the Introduction (Section 1), which is
followed by a discussion on relevant past work in the Literature Review (Section 2). The proposed
framework for model reduction is described in the Research Methodology (Section 3). The reduced
approach is analysed in the Results and Discussion (Section 4). The Conclusions (Section 5) provide a
further discussion on the highlights of the reduction method and directions for future work.
3 Background
There are several methods in use to reduce the complexity of biochemical reaction network models.
Quasi-steady state approximation (QSSA) is the most popular approach of model reduction used with
biochemical kinetics (Bodenstein, 1913). This idea was later used to produce the classical Michaelis
Menten formula (Michaelis and Menten, 1913). More improvements were made on this method that
has become the basic tool of analysis of chemical reaction mechanisms and kinetics (Segel and
Slemrod, 1989). Another approach to model reduction is the geometric singular perturbation method
(GSPM). This is useful when the models have slow and fast variables. This approach works based on
the assumption that when the model contains fast and slow variables, and the slow variables control
the fast variables, we can remove the fast variables from the model (Tikhonov, 1952; Fenichel, 1979;
Jones, 1995).
Wei and Kuo (1969) proposed a lumping technique that works based on combining reagents into quasicomponents for reducing the system size in monomolecular and pseudo monomolecular systems.
Another reduction method was proposed by Radulescu and his colleagues (2012). They suggested
some essential ideas such as dominance and limitation to impove the understanding of dynamical
and biological systems (Radulescu et al., 2012). The methods of averaging and invariant manifolds
proved their effectiveness in reducing the asymptotic dynamics of chemical reaction networks with
clear time-scale separation. Castiglione and his colleagues (2014) stressed the importance of multiscale modelling of biological systems to resolve the challenge of large amount of measurements
required in modern biology. The importance of this becomes clear when developing comprehensive
computational and mathematical models for a system with different spatial and time scales. For
instance, they used multi-scale modelling on the immune system.
Several other methods have also been used for reducing chemical reactions networks. For instance,
Noel et al. (2012) reduced and hybridised networks of biochemical reactions using the LitvinovMaslov correspondence principle. They applied this method on a cell cycle oscillator model. Lam and
Goussis (1994) implemented a singular computational perturbation technique for reduced kinetic
3
modeling focusing on its relative merits versus standard methodologies. Maas and Pope (1992; 1994)
used an intrinsic low dimensional manifold to exploit the interval between timescales of various
processes and variables.
Another model reduction method for computational biology was presented by Radulescu et al. (2015).
This method was inspired by tropical geometry and analysis that combines graphical approaches, semiquantitative reasoning and symbolic manipulation. However. it does not explain how to combine onescale approximations to obtain a multi-scale approximation that is valid for both fast and slow time
scales. O t h e r recent studies of model reduction with examples of application from systems
biology can be found in Anderson et al. (2011); Karadeniz et al. (2012); Kutumova et al. (2013); Ishizaki
et al. (2014); Kooshkbaghi et al. (2014) and West et al. (2015). There are many reduction methods but
a fully formal method that exploits the hierarchical orders in large biological network models is
missing. Known reduction methods often face difficulties when applied to complex systems.
We introduce an example biological network involving molecular interactions that control cell cycle
that implements cell division as the basis for conceptual model development and implementation. In
this review section, we briefly describe cell cycle and use models developed for cell cycle as the basis
for a discussion on model complexity and model reduction. Cell cycle progression is tightly regulated
by complex proteins in regulatory networks to achieve correct cell division (Behl & Ziegler, 2014) in
development, regeneration and wound healing of an organism. It has several checkpoints; for
example, the G1/S and G2/M checkpoints (Saltman, 2005) that ensure the integrity of the genome.
Moreover, dysfunction in cell cycle can lead to changes in the DNA, resulting in the development of
diseases such as cancer (King et al., 2003; Azimi et al., 2017; Farr et al., 2017). A better understanding
of protein regulatory networks will, therefore, not only advance our understanding of fundamental
cell cycle regulation, such as the G1/S checkpoint, but also provide a greater understanding of the
disease processes, thus increasing the efficiency of treatment of diseases. A factor that has greatly
hampered the progress in this direction is model complexity.
Biologists have used many ways to understand the mechanics of protein regulatory networks. One way
is to construct models that simulate the mechanisms of the interactions between proteins. Several
types of model have been developed to represent cell cycle that involves protein to protein
interactions. These types of models include: mathematical ordinary differential equations (ODE)
(Aguda, 1999; Novak& Tyson, 2004; Tashima et al., 2007; Iwamoto et al., 2008; Ling et al., 2010;
Iwamoto et al., 2011; Zhao et al., 2012), Boolean (Fauré et al., 2006), petri net (Kotani et al., 2002;
Herajy et al., 2013), recurrent neural networks (Ling et al., 2013) and hybrid (Singhania et al., 2011).
ODE models are the most common type used for modelling biological networks. This is because they
can provide continuous systems dynamics in the form of change in concentration of proteins ove
4
time. They are also useful for analysing the nature of systems dynamics, such as stability. As such
they offer the advantage of representing detailed system behaviour. However, over time, the
amount of data obtained from laboratory experiments has increased, and so these models have
accordingly become more complex. In particular, if they contain hundreds or thousands of variables,
ODE models become ineffective (Danos et al., 2007). These challenges with the ODE models offers
new opportunities for researchers to develop new methods to simplify biochemical reaction network
models. Many approaches to model reduction have been proposed that include the other types of
models presented above (Clarke, 1992; Maas & Pope, 1992; Lam & Goussis, 1994; Maas & Pope, 1994;
Clarke et al., 1996; Gorban & Karlin, 2005; Feret et al., 2009; Gorban et al., 2010; Noel et al., 2011;
Radulescu et al., 2012; Rao et al., 2013; Sonday et al., 2013; Rao et al., 2014; Radulescu et al., 2015).
In this work, we focus on reducing complex ODE mathematical models of large molecular signalling
networks by using time (activity) windows and simpler logical models. As a case study, a continuous
ODE model (Iwamoto et al., 2011) was reduced into a time-tagged logical model describing the G1/S
checkpoint with and without DNA damage.
The regulatory and signalling networks have been represented by logical models since the time of
Kauffman (1969) who pioneered this field. A logical model is one of the simplest models that can be
used to describe the dynamics of regulatory and signalling networks without the need for many
parameters. A logical model contains a series of interconnected elements Y1… Yn. Each element has
only two possible states of activation, 1 (active) and 0 (inactive). The regulatory and signalling network
is completely described by a set of logical equations as shown in Eq. 1(a, b, c) (Mendoza and Xenarios,
2006):
( + 1) = ( 1 ( )
2 ( )…
2 ( )…
( + 1) = 1 ( )
( 1( )
( + 1) =
2 ( )…
( ))
()
( 1( )
( ))
2 ( )…
( ))
(1 )
(1. )
(1. )
where Yi∈ {0, 1}, {Yna} is the set of activators of Yi, {Ymi} is the set of inhibitors of Yi. If Yi has activators
and inhibitors use Eq. 1.a, if Yi has only activators use Eq. 1.b and if Yi has only inhibitors use Eq. 1. c.
Figure 1 shows an example of building a logical model. There are three steps to the modelling process:
(i) find the regulatory graph of the system (Fig. 1a); (ii) define the logical conditions (Fig. 1b); (iii)
specify the logical equations (for updating (Fig. 1c)).
5
Figure 1 Interaction graph, logical model and logical equations of an example protein network (a)
Interaction graph: The nodes (Y1, Y2 … Y6) in the graph represent proteins, and the edges represent
interactions. Blue arrows represent activations and red bar arrows inhibitions. (b) Logic-based
representation of a logical model for the interaction graph given in (a). (c) Equations of the logical
model for the interaction graph in (a).
As seen in Figure 1, a logical model is represented by a logic-based graph and every node in the logical
model can be derived from the state of other nodes through a logic equation. Changing the state of
nodes from ‘on’ to ‘off’ or vice-versa through a sequence of changes in state represents a biological
process. (Chen et al., 2016). The dynamics of a logical model are defined by synchronous (all nodes are
updated at the same time) or asynchronous (nodes are updated individually in random or some
predefined order) updating. The system dynamics usually depend on the choice of the updating
scheme (Aracena et al., 2009).
Several software tools exist for building and analysing the dynamics of logical models, such as GINsim
(Gonzalez et al., 2006), SQUAD (Di Cara et al., 2007), Boolean Net (Albert et al., 2008), Chem Chains
(Helikar & Rogers, 2009), Odefy (Krumsiek et al., 2010), Bool Net (Müssel et al., 2010), Cell NOpt
(Terfve et al., 2012), MaBoSS (Stolletal., 2012) and Cell Collective (Helikar et al., 2013).
There are few logical models that have been applied to the mammalian cell cycle. For more details
about these models, refer to (Huang & Ingber, 2000; Fauré et al., 2006; Sahin et al., 2009; Singhania
et al., 2011; Mombach et al., 2014, Tanaka et al., 2017). Many logical models have been built to
model cell cycle of budding yeast (Li et al., 2004; Faur´e et al., 2009; Irons, 2009; Todd & Helikar,
2012; Alcasabas et al., 2013; Rubinstein et al., 2013; Chasapi et al. 2015; Linke et al., 2017). They
have shown the efficacy of Boolean models in producing representative qualitative systems
dynamics. Most of these models represent all interactions in Boolean form but lately hybrid logical
and ODE methods have been used for representing fast and slow reactions, respectively (Singhania et
al., 2011). These models still attempt to capture the full-time spectrum.
6
The purpose of this research is to aid further reduction of model complexity through a combined
logical and time slicing (window) method that still yields the expected behaviour of a biological
network. The model concept is rooted in the separation of the time spectrum into periods of
increased activity and low activity. Specifically, we focus on simplifying the time dimension by
dividing a signaling system into time windows (active time windows, and steady or frozen time
windows). The active time windows are then represented by logical models to aid the model
reduction process. We demonstrate the method on the G1/S checkpoint of cell cycle with and
without DNA damage.
4 Research Methodology
The proposed approach to reduce a model of a biological network into a simpler one using time
windows and logical models involves the following tasks: (i) Understanding the protein interactions
involved in the system; (ii) Dividing the dynamics of the system into time windows (active time
windows, steady or frozen time windows); (iii) Determining the key elements in each active time
window; and (iv) Building a logical model for each active time window. Last task is divided into three
sub-tasks: first, finding a regulatory graph for the system; secondly, defining the logical conditions of
the system; and thirdly, specifying the logical equations (for updating). These tasks are described in
detail in the following sections.
4.1 Understanding the Protein Interactions in the System
After reviewing previous works and models of Lev Bar-Or et al. (2000), Tashima et al. (2004), Tashima
et al.(2006), Iwamoto et al. (2008) and Iwamoto et al. (2011), we extracted the model of Iwamoto et
al. (2011) for the whole cell cycle regulation incorporating DNA damage signalling pathway. We use
this full model for extracting the G1/S checkpoint pathway integrated with DNA damage signalling
shown in Fig. 2 for our research. The model consisted of 35 dependent variables and 92 kinetic
parameters (see Appendix A). Both initial conditions and kinetic parameters in our model were from
the original Iwamoto et al. (2011) model as found in the previous works of Lev Bar-Or et al. (2000),
Tashima et al. (2006) and Iwamoto et al. (2008).
Cell cycle takes place in 4 phases of G1, S, G2 and M. In G1, a cell grows in size and prepares for DNA
replication that happens in S phase. The cell grows further in G2 phase and DNA segregation happens
in M phase followed by cytokinesis that separates the cell into two daughter cells. The transition from
G1 to S and G2 to M are regulated by checkpoints (G1/S and G2/M) to ensure the integrity of the DNA.
We consider G1 to S transition phase of cell cycle shown in Fig. 2 to demonstrate the model concept.
In a normal cell, Cyclin D (CycD) is synthesised and G1 progression is initiated in the presence of growth
factors. The activated binary complex CycD/Cdk4, forms when CycD binds to Cdk4. The CycD/Cdk4
7
complex phosphorylates Rb bound to E2F, which generates the hypophosphorylated form, RbPP/E2F. Furthermore, the hyperphosphorylated form (Rb-PPP) is generated when Rb-PP/E2F is
phosphorylated by both CycE/Cdk2 and CycA/Cdk2. Rb-PPP disassociates from Rb-PP/E2F and E2F is
activated. The transcription of CycE, CycA and Cdc25A (phosphatase that activates these 2 cyclins) are
activated by transcription factor E2F, and this is required for DNA replication. CycE and CycA bind to
inactive Cdk2 molecules to form the inactive complexes, CycE/Cdk2 and CycA/Cdk2, respectively
(Iwamoto et al., 2011; Alberts et al., 2015). Activated binary complexes, CycE/Cdk2-P and CycA/Cdk2P, form after Cdc25A dephosphorylates and activates the corresponding inactive complexes. The
further dissociation of Rb-PPP and E2F and the concomitant activation of E2F happen as a result of
phosphorylation of Rb-PP/E2F by both CycE/Cdk2-p and CycA/Cdk2-p. This results in further freeing
of E2F, thus establishing a positive feedback loop between E2F and CycE; increased concentration of
CycE (and E2F) moves the cell from G1 to S phase (Satyanarayana & Kaldis, 2009).
8
Figure 2. G1/S transitions and DNA damage signalling pathway abstracted from Iwamoto et al.
(2011) showing the proteins (green nodes) and their various interactions (red and blue arrows)
The positive feedback loop between the two types of Cyc/Cdk complexes and the Rb/E2F complex
plays an essential role in G1 progression. Driving the progression to the S phase and initiating DNA
replication happen through sufficient expression of both E2F and CycE/Cdk2-P. After completion of
DNA replication during S phase, CycA/Cdk2-P promotes the degradation of E2F that reduces the
synthesis of CycE and the formation of CycE/Cdk2-P, which causes progression to the G2 phase.
Regulation of G1/S progression depends on three Cdk inhibitors (CKIs), p16, p21, and p27. p16 inhibits
the activation of CycD/Cdk4, and both p21 and p27 repress CycE/Cdk2-P and CycA/Cdk2-P (Iwamoto
et al., 2011; Alberts et al., 2015) as shown in Figure 2.
Different protein kinases are recruited to the damage site when DNA is damaged and a signalling
pathway is initiated that triggers cell cycle arrest. Depending on the type of damage, the first kinase
9
on the damage site is ATM/ATR. Additional protein kinases, Chk1 and Chk2, are then recruited and
activated, leading to the activation (phosphorylation) of gene regulatory protein (transcription factor)
p53. Mdm2 usually binds to p53 and encourages its ubiquitylation and proteasome destruction. p53
phosphorylation blocks its binding to Mdm2; as a consequence, p53 accumulates at elevated
concentrations and stimulates the transcription of the CKI protein encoding gene, p21. p21 binds and
inactivates complexes CycE/Cdk2 and CycA/Cdk2, arresting the cell in G1 (Iwamoto et al., 2011; Alberts
et al., 2015).
We edited the ODEs in the original Iwamoto et al. (2011) whole cell cycle model to make the model
suitable to represent the DNA damage signalling pathway integrated only with G1/S checkpoint
pathway. After this base model was extracted, we estimated the time steps required for G1/S, which
was approximately 12.32 hours (Ohtsubo et al., 1995) in real time, where every 189 model time steps
corresponds to an hour, meaning that approximately 2330-time steps are needed to simulate G1/S.
To evaluate the base model, we simulated it with and without DNA damage, and compared the results
with the Iwamoto full model to ensure that the base model was correct before using it to answer the
specific questions we raised in this study. We used the base model in Fig. 2 as a complex system case
study to apply our new reduction approach.
4.2 Divide the Dynamics of the System into Time Windows (Active Time Windows, Steady
or Frozen Time Windows)
Biological processes can take place over a long timescale. For example, our G1/S network base model
generates a circadian rhythm as shown in Figure 3 where protein behaviour can be seen over three
distinct time-windows. In time window A, some proteins undergo significant changes in concentration
while other proteins see very little or no change. In time window B, there are no changes in most
proteins. In time window C, many proteins interact and undergo changes in concentration and no
change occurs in others. A system where exist such distinct regions of activity makes it amenable to
model reduction while preserving the general system behavior. To model such systems, primary
activity regions (time- windows) must be first chosen (time window A and C); and the other timewindows are then treated as relatively frozen in time (time window B).
10
Figure 3. Behaviour of proteins in multiple time windows for G1/S checkpoint with
DNA damage signalling
Before building reduced models, these separate time windows are decided from the response of the
base model. They often motivate choices on what species and processes should be included in the
model and what can be overlooked. It is possible to simplify current models that display distinct time
windows. With this, our method of model reduction aims to approximate the initial model with a
model of reduced complexity.
Model reduction by separation of time windows should lead to similar results to an original model,
where a differential equation describing a state variable is replaced by a new simpler equation. The
main idea behind model reduction technique by time windows is to treat slow variables in low or
frozen activity windows as fixed entities (constant value), rather than as state variables. In the
reduction method using time windows, the number of elements in the system depends on the time
windows and the number of elements changing with time within them.
4.3 Determine the key elements in each active time window
After studying the G1/S checkpoint protein behaviour with and without DNA damage in multiple time
windows, the key elements in each active time window (A and C) are determined. Table 1 shows the
G1/S checkpoint protein behaviour with and without DNA damage over multiple time windows.
Accordingly, time windows A and C contain a significant number of fast reactions and window B
contains a large number of slow reactions that either already had their peak activity in time window A
or are going to reach peak activity in time window C. Time window A contains 14 elements that play
11
main roles as determined by their higher level of activity. These elements are CycD, Cdk4, CycD/Cdk4,
p27, p27/CycD/Cdk4, p21, p53, Mdm2 and ATM/ATR (latter four when DNA damage occurs),
p21/CycD/Cdk4, p16, Rb/E2F, Rb-PP/E2F and Rb. Also, time window C also contains 14 elements that
play main roles: CycE, Cdk2, CycE/Cdk2, CycE/Cdk2-P, p27, p27/CycE/Cdk2-P, p21, p21/CycE/Cdk4-P,
Rb- PP/E2F, E2F, Rb-PPP, p53, Mdm2 and ATM/ATR.
Table 1 G1/S checkpoint protein behaviour with and without DNA damage in multiple time
windows
No
Molecule/complex
A
B
C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
CycD
CycE
CycA
Cdk4
Cdk2
CycD/Cdk4
CycE/Cdk2
CycE/Cdk2-P
CycA/Cdk2
CycA/Cdk2-P
p27
p27/CycD/Cdk4
p27/CycE/Cdk2-P
p27/CycA/Cdk2-P
Fast
Slow
Slow
Fast
Slow
Fast
Slow
Slow
Slow
Slow
Fast
Fast
Slow
Slow
Slow
Slow
Slow
Slow
Linear
Slow
Slow
Slow
Slow
Slow
Slow
Slow
Slow
Slow
Slow
Fast
Slow
Slow
Slow
Slow
Fast
Fast
Slow
Slow
Fast
Slow
Fast
Slow
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
p21
p21/CycD/Cdk4
p21/CycE/Cdk2-P
p21/CycA/Cdk2-P
p16
Rb/E2F
Rb-PP/E2F
E2F
Rb-PPP
Rb
p53
Mdm2
ATM/ATR
Cdc25A
Cdc25A-P
Chk1
Chk1-P
Fast *
Fast
Slow
Slow
Linear
Fast
Fast
Slow
Slow
Fast
Fast *
Fast *
Fast *
Slow
Slow
Linear
Slow
Slow
Linear
Slow
Slow
Linear
Slow
Slow
Slow
Slow
Slow
Slow
Slow
Slow
Linear
Slow
Slow
Slow
Fast *
Slow
Fast
Slow
Linear
Slow
Fast
Fast
Fast
Slow
Fast *
Fast *
Fast *
Slow
Slow
Slow
Slow
32
33
34
35
NF-Y
B-Myb
B-Myb-P
Im
Slow
Linear
Slow
Fast
Slow
Linear
Slow
Slow
Slow
Slow
Slow
Slow
* When DNA damage occurs
12
4.4 Building Logical Models for each Active Time Window
According to Table 1, we need to build two logical models, R1 and R2, representing reduced models
for time windows A and C, respectively. The first step to construct any logical model is to find a
regulatory graph for the system, which in our case depends on the ODE model. See Table B1 in
Appendix B for the ODEs of the reduced model (R1) and Table C1 in Appendix C for the ODEs of the
reduced model (R2). The new ODEs contain only the relations of the 14 elements that play main
roles. Figures 4 and 5 show the regulatory graphs for reduced models R1 and R2, respectively.
Figure 4. Regulatory graph for the reduced model (R1) for time window A of G1/S checkpoint with
and without DNA damage signaling (yellow nodes on the top right represent DNA
damage part)
The second step is to convert each regulatory graph into a logic-based graph. Figure 6 and 7 show the
logical models for the interactions in regulatory graphs in Figures 4 and 5, respectively. The last step is
to specify the logical equations (for updating) and the initial values of the proteins. The reduced model
(R1) represents 14 logical equations (see Table B3 in Appendix B). Table B2 shows the initial conditions
of the reduced model (R1). Meanwhile, the reduced model (R2) represents 14 logical equations (see
Table C3 in Appendix C). Table C2 shows the initial conditions of the reduced model (R2). Synchronous
method is used to update the status of nodes where nodes are updated simultaneously. See
Appendix D for an example detailing the steps to convert an ODE equation to Boolean equation.
13
Figure 5. Regulatory graph for the reduced model (R2) for time window C for G1/S checkpoint with
and without DNA damage signaling (yellow nodes on the top right represent DNA
damage part)
Figure 6. Logic-based graph for the reduced model (R1)
14
Figure 7. Logic-based graph for the reduced model (R2)
5 Results and Discussion
The simulation results of the reduced models R1 and R2 with and without DNA damage agreed with
the biological knowledge from previous studies as described below. CycD and CycE play the main roles
in mammalian G1-to-S transition. As shown in Figure 8, for the reduced model (R1) without DNA
damage, the state of CycD becomes active in early G1/S. CycD binds to Cdk4 forming the binary
complex, CycD/Cdk4, as shown (at time step 2), and CycD/Cdk4 becomes active. Furthermore,
CycD/Cdk4 complex phosphorylates E2F bound Rb, which becomes hypophosphorylated Rb-PP/E2F
(Helin, 1998). We note that Rb/E2F becomes inactive and Rb-PP/E2F becomes active (at time step 3)
and Rb becomes inactive (at time step 4). This behvaiour that initiates the first phase of release of
E2F by hypophosphorylating Rb/E2F to form Rb-PP/E2F represents the early part of the cell cycle up
to the commitment point where a cell makes an irreversible decision to commit to the completion of
cell cycle. p53 and p21 stay inactive for all iterations in the simulation of the reduced model (R1) as
there was no DNA damage.
15
Figure 8. Simulation results of the reduced model (R1) without DNA damage (green: active, red:
inactive)
After a damage occurs in DNA, the simulation results of the reduced model (R1) change as shown in
Figure 9. The DNA damage signal activates ATM/ATR (at time step 1). ATM/ATR activates p53 (at time
step 2), which regulates transcription of a large number of genes required for various purposes; cell
cycle arrest, cell death (apoptosis) or damage recovery (Ciliberto et al., 2005; Harris & Levine, 2005;
Geva-Zatorsky et al., 2006).
In the initial part of the cell cycle described by this model, p53 activation leads to the activation of p21
(Yu et al., 1999) (step 3). The function of p21 is to inhibit CDK activity by inhibiting Rb phosphorylation
in order to maintain E2F inactive (Campisi & Fabrizio, 2007). This explains the activation of
p21/CycD/Cdk4 (at time step 4). Furthermore, as seen in Figure 9, there are oscillations in the
behaviour of p53, p21 and Mdm2. These results agree with the results as described in the Iwamoto et
al. (2011) model. Specifically, activated p53 activates Mdm2 and p21. However, in the next step MdM2
deactivates p53 which results in reduction or loss of p21 in the next step. But ATM/ATR reactivates
p53 and suppresses Mdm2. This oscillatory behaviour of the 3 proteins repeats as clearly shown in
the figure.
16
Figure 9. Simulation results of the reduced model (R1) with DNA damage (green: active, red:
inactive)
As shown in Figure 10, the simulation results of the reduced model (R2) without DNA damage are for
mid to late G1 phase of cell cycle. Activation of CycE/CDK2 is the primary reason for further
hyperphosphorylation of Rb-pp/E2F leading to the dissociation of Rb-PPP from E2F, thereby releasing
E2F. As shown (at time step 1) in Figure 10, E2F and Rb-PPP become active. The enhanced
concentration of E2F promotes CycE synthesis in the mid-to-late G1 phase (time step 2). Produced
CycE then binds to CycE to CDK2 to create more CycE/CDK2 (at time step 3) so CycE/CDK2-P becomes
active. This leads to a further release of E2F and thus to a positive feedback loop between E2F and
CycE; increased concentration of CycE (and E2F) moves the cell from the G1/S checkpoint to the S
phase (Satyanarayana & Kaldis, 2009). These results agree with the reported experimental results
(Kohn, 1999; Hochegger et al., 2008).
Now, there is an interesting interplay between cyclins and p27 in relation to the above processes.
One of CycD/CDK4's roles is to maintain CycE/CDK2-P active by (i) competing with CycE/CDK2-P for
free p27 to form p27/CycD/CDK4 complex as shown (at time step 3) in Figure 9 and by (ii)
sequestering p27 from p27/CycE/CDK2-P. The resulting active CycE/CDK2-P causes additional Rbpp/E2F hypophosphorylation, that encourages Rb-PPP and E2F dissociation leading to the release of
E2F (Obaya and Sedivy, 2002). Increasing the amount of free E2F promotes CycE synthesis in the midto-late G1 phase, facilitating the CycE and CDK2 to create more of the CycE/CDK2 complex. This
results in further increased E2F activity contributing to the establishment of the aforementioned
positive feedback loop between CycE and E2F; increased E2F and CycE concentrations allow the cell to
pass the G1/S transition and initiate the S phase (Hiebert et al., 1992). CycD degradation happens in
the mid- G1 stage promoting the release of p27 bound to the CycD/CDK4 complex. Free p27 is also
redistributed to the new CycE/CDK2-P complex. Although p27 can inhibit CycE/CDK2-P, when p27
binds to CycE/CDK2-P, the large amount of CycE/CDK2-P can initiate p27 degradation by
17
phosphorylation at
18
threonin (Coqueret, 2003). This explains a significant degradation of p27 at the end of G1 due to the
accumulation of CycE/CDK2-P in the cell as shown (at time step 1 onwards) in Figure 10.
Figure 10. Simulation results of the reduced model (R2) without DNA damage (green: active, red:
inactive)
As shown in Figure 11 for the simulation results of the reduced model (R2) with DNA damage, p21
becomes active. The function of p21 is to inhibit CDK activity by inhibiting Rb phosphorylation to
prevent E2F from being active to effect cell cycle (Campisi & Fabrizio, 2007). Here, Rb-PP/E2F is the
focus which is further phosphorylated by CycE/Cdk2-p. Therefore, p21 is bound to CycE/Cdk2-P
forming the complex, p21/CycE/Cdk2-P. As shown (at time step 1), p21/CycE/Cdk2-P becomes active.
This delays the release of E2F from Rb-PP/E2F for one iteration more, which means a delay in the
synthesis of CycE. As shown, E2F becomes active (at time step 2) and CycE becomes active (at time
step 3). This delay in activation represents cell cycle arrest. These results in terms of the steps agree
with experimental results (Lev Bar-Or et al., 2000; Batchelor et al., 2008). The negative feedback loop
of p53 and Mdm2 is completely restored with the removal of DNA damage and p53 returns to a low
level (step 4). The decrease in p53 reduces the p21 level, releasing the CycE/CDK2-P complexes (and
CycA/CDK2-P, not shown in the figure) and returning the cell cycle to normal condition.
Figure 11. Simulation results of the reduced model (R2) with DNA damage (green: active, red:
inactive)
19
Results show that the reduced model R1 has well captured the main processes in the early G1 phase
up to the commitment point where Rb/E2F is hypophosphorylated. The reduced model R2 has well
captured the processes beyond the commitment point where E2F is released and E2F and CycE
promote each other in a feedback loop resulting in increased levels of CycE/Cdk2-P which denotes the
cell cycle transition into S phase. CycE/Cdk2-P has the role of activating DNA synthesis in the S phase
and therefore reaching the required level of its activation is important for this transition (Dulic et al.,
1994). The results from the simulation indicates that the reduced models together can depict the most
important processes converging into important aspects such as commitment point and phase
transition of cell cycle in G1 phase with reduced complexity in each model. Further, reduced
complexity means that each model can be developed more quickly, understood more clearly and
interpreted more easily.
6 Summary and Conclusions
There has been much interest in the reduction of complex ODE regulatory network models for many
reasons: (i) Most of the models that have been built to understand regulatory networks were ODEbased mathematical models; (ii) Most of the ODE-based mathematical models were complex and
needed kinetic information for reactions that was not easily gathered; (iii) Simplifying complex ODEbased mathematical models can lead to better understanding and control of these systems.
We presented a way of reducing complex quantitative (ODE) models into logical models where the use
of time windows allowed the reduction of time complexity that facilitated the reduction of model
complexity. Specifically, it exploits the fact that many interactions of a signaling network are of
activation/inhibition nature occurring relatively fast and fewer interactions involve protein synthesis
that is a slow process. Therefore, the operation of a network can pass through windows of time with
high, low or minimum activity and this could reveal boundaries for modularising or reducing networks.
We showed this in our example G1/S network that revealed windows of time with high or little to no
activity that enabled the separation of time into active and inactive time windows in a qualitative
sense. This allowed extraction of segments of the network (reduced models) that are active in windows
characterised by high activity. This has significantly reduced the number of time steps needed to run
the model. For example, for the part in Time window A, the ODE model needs 670 time steps while
our R1 model needs only 10 time steps. Similarly, for the part in Time window C, ODE model needs
290 time steps and our R2 model takes only 4 time steps. Further model reduction was enabled by
logical representation of the reduced models to do away with the kinetic parameters in a system. The
method is general and can be readily applied to any complex quantitative model. In our case study, a
continuous ODE model was transformed into a reduced logical model to describe the G1/S
checkpoint with and without DNA damage. The results revealed that the two reduced models
represent the major events in cell cycle: the first reduced model represents events up to the
20
commitment point in cell cycle where a cell makes an irreversible decision to complete cell cycle. The
second reduced model properly represents the activation of CycE/Cdk2 that denotes transition into S
phase as it is the major controller of G1/S transition. The reduced model also explains and predicts
the G1/S checkpoint behaviour with DNA damage, including correct representation of known
oscillations of the state of molecules in the DNA damage pathway and cell fate after damage, such as
cell cycle arrest. This showed that the reduced model was able to capture the main events in the
G1/S checkpoint and therefore it is helpful to gain insights into biology. It was much easier and
simpler to model, run and analyse than the original ODE model.
The approach presented here could greatly help to simplify complex quantitative models into simpler
models. Further, the interaction between modelling and experiments could be better facilitated due
to its simplicity and interpretability. Moreover, researchers and those who construct these models
will be helped to concentrate on understanding and representing system behavior rather than
determining kinetic parameter values until they are accessible. Both ODE and Boolean network
modelling approaches have advantages and disadvantages. For large complex biological problems
such as cell cycle, using simpler Boolean models are more appropriate for gaining a system level
qualitative understanding than complex ODE models. This advantage can be further exploited by
intergating Boolean with time windows, as presented in this study, which implicitly brings time into
the picture which is a challenge for Boolean models in general.
However, Boolean models being
discrete and binary, they cannot produce intermeidate states of protein concentrations as ODE
models do.
Therefore, in future, the proposed reduced method could be improved by using time
windows and multi-valued logical models (semi-quantitative models). Further, in our model, there is
no loss of essential information as the reactions in the frozen time window reached full activity (half
expression of proteins were considered as full activity in Boolean sense) either in time window A or
C. Therefore, all important activities are caaptured by the model while reducing the time needed for
simulation. Appropriate selection of time windows is important for getting good results from the
proposed method.
7 References
Aguda, B. D. (1999). A quantitative analysis of the kinetics of the G2 DNA damage checkpoint system.
Proceedings of the National Academy of Sciences, 96(20), 11352-11357.
Albert, I., Thakar, J., Li, S., Zhang, R., & Albert, R. (2008). Boolean network simulations for life scientists.
Source code for biology and medicine, 3(1), 16.
Alcasabas, A. A., de Clare, M., Pir, P., & Oliver, S. G. (2013). Control analysis of the eukaryotic cell cycle
using gene copy-number series in yeast tetraploids. BMC genomics, 14(1), 744.
Anderson, J., Chang, Y. C., & Papachristodoulou, A. (2011). Model decomposition and reduction tools
for large-scale networks in systems biology. Automatica, 47(6), 1165-1174.
21
Aracena, J., Goles, E., Moreira, A., & Salinas, L. (2009). On the robustness of update schedules in
Boolean networks. Biosystems, 97(1), 1-8.
Azimi, I., Milevskiy, M. J., Kaemmerer, E., Turner, D., Yapa, K. T., Brown, M. A., ... & Monteith, G. R.
(2017). TRPC1 is a differential regulator of hypoxia-mediated events and Akt signaling in PTENdeficient breast cancer cells. J Cell Sci, jcs-196659.
Bodenstein, M. (1913). Eine theorie der photochemischen reaktionsgeschwindigkeiten. Zeitschrift für
physikalische Chemie, 85(1), 329-397.
Campisi, J., Fabrizio, d.A.d.F., (2007). Cellular senescence: when bad things happen to good cells. Nat.
Rev. Mol. Cell Biol. 8 (9), 729–740.
Castiglione, F., Pappalardo, F., Bianca, C., Russo, G., & Motta, S. (2014). Modeling biology spanning
different scales: an open challenge. BioMed Research International.
Chasapi, A., Wachowicz, P., Niknejad, A., Collin, P., Krapp, A., Cano, E., ... & Xenarios, I. (2015). An
extended, Boolean model of the Septation Initiation Network in S. pombe provides insights
into its regulation. PloS one, 10(8), e0134214.
Chen, H., Wang, G., Simha, R., Du, C., & Zeng, C. (2016). Boolean models of biological processes explain
cascade-like behavior. Scientific reports, 7.
Ciliberto, A., Novak, B., Tyson, J., (2005). Steady states and oscillations in the p53/Mdm2 network. Cell
Cycle 4 (3), 488–493.
Coqueret O., (2003). New roles for p21 and p27 cell-cycle inhibitors: a function for each cell
compartment? Trends in Cell Biology, 13(2), 65-70.
Di Cara, A., Garg, A., De Micheli, G., Xenarios, I., & Mendoza, L. (2007). Dynamic simulation of
regulatory networks using SQUAD. BMC bioinformatics, 8(1), 462.
Dulic V., Kaufmann W.K., Wilson S.J., Tisty T.D., Lees E., Harper J.W., Elledge S.J., & Reed S.I., (1994).
p53-dependent inhibition of cyclin-dependent kinase activities in human fibroblasts during
radiation-induced G1 arrest. Cell, 76(6), 1013-1023.
Farr, K. J., Gallaway, P. J., & Hongu, N. (2017). Breast Cancer Prevention: Exercise and Healthy Diet. The
University of Arizona Cooperative Extension.
Fauré, A., Naldi, A., Chaouiya, C., & Thieffry, D. (2006). Dynamical analysis of a generic Boolean model
for the control of the mammalian cell cycle. Bioinformatics, 22(14), e124-e131.
Fauré, A., Naldi, A., Lopez, F., Chaouiya, C., Ciliberto, A., & Thieffry, D. (2009). Modular logical modelling
of the budding yeast cell cycle. Molecular BioSystems, 5(12), 1787-1796.
Fenichel, N. (1979). Geometric singular perturbation theory for ordinary differential equations. Journal
of Differential Equations, 31(1), 53-98.
Geva-Zatorsky, N., Rosenfeld, N., Itzkovitz, S., Milo, R., Sigal, A., Dekel, E., Yarnitzky, T., Liron, Y., Polak,
P., Lahav, G., Alon, U., (2006). Oscillations and variability in the p53 system. Mol. Syst. Biol., 2.
Gonzalez, A. G., Naldi, A., Sanchez, L., Thieffry, D., & Chaouiya, C. (2006). GINsim: a software suite for
the qualitative modelling, simulation and analysis of regulatory networks. Biosystems, 84(2),
91-100.
Harris, S.L., Levine, A.J., (2005). The p53 pathway: positive and negative feedback loops. Oncogene 24,
2899–2908.
Helikar, T., & Rogers, J. A. (2009). ChemChains: a platform for simulation and analysis of biochemical
networks aimed to laboratory scientists. BMC systems biology, 3(1), 58.
Helikar,T., Kowal,B., & Rogers,J.A. (2013). A cell simulator platform: the cell collective. Clin. Pharmacol.
22
Ther. 93, 393–395. doi:10.1038/clpt.2013.41.
Helin, K., 1998. Regulation of cell proliferation by the E2F transcription factors. Curr. Opin. Genet. Dev.
8 (1), 28–35.
Hiebert, S.W., Chellappan, S.P., Horowitz, J.M., & Nevins, J.R., (1992). The interaction of RB with E2F
coincides with an inhibition of the transcriptional activity of E2F. Genes & Development, 6(2),
177-185.
Huang, S., & Ingber, D. E. (2000). Shape-dependent control of cell growth, differentiation, and
apoptosis: switching between attractors in cell regulatory networks. Experimental cell
research, 261(1), 91-103.
Irons, D. J. (2009). Logical analysis of the budding yeast cell cycle. Journal of theoretical biology, 257(4),
543-559.
Ishizaki, T., Kashima, K., Imura, J. I., & Aihara, K. (2014). Model reduction and clusterization of largescale bidirectional networks. IEEE Transactions on Automatic Control, 59(1), 48-63.
Iwamoto, K., Hamada, H., Eguchi, Y., & Okamoto, M. (2011). Mathematical modeling of cell cycle
regulation in response to DNA damage: exploring mechanisms of cell-fate determination.
Biosystems, 103(3), 384-391.
Iwamoto, K., Tashima, Y., Hamada, H., Eguchi, Y., & Okamoto, M. (2008). Mathematical modeling and
sensitivity analysis of G1/S phase in the cell cycle including the DNA-damage signal
transduction pathway. Biosystems, 94(1), 109-117.
Jones, C. (1995). Geometric singular perturbation theory. Dynamical Systems, 44-118.
Karadeniz, H., Soyhan, H. S., & Sorusbay, C. (2012). Reduction of large kinetic mechanisms with a new
approach to the necessity analysis method. Combustion and Flame, 159(4), 1467-1480.
Kauffman, S. A. (1969). Metabolic stability and epigenesis in randomly constructed genetic nets.
Journal of theoretical biology, 22(3), 437-467.
King, M. C., Marks, J. H., & Mandell, J. B. (2003). Breast and ovarian cancer risks due to inherited
mutations in BRCA1 and BRCA2. Science, 302(5645), 643-646.
Kooshkbaghi, M., Frouzakis, C. E., Boulouchos, K., & Karlin, I. V. (2014). Entropy production analysis for
mechanism reduction. Combustion and Flame, 161(6), 1507-1515.
Krumsiek, J., Pölsterl, S., Wittmann, D. M., & Theis, F. J. (2010). Odefy-from discrete to continuous
models. BMC bioinformatics, 11(1), 233.
Kutumova, E., Zinovyev, A., Sharipov, R., & Kolpakov, F. (2013). Model composition through model
reduction: a combined model of CD95 and NF-κB signaling pathways. BMC Systems Biology,
7(1), 1.
Lam, S. H., & Goussis, D. A. (1994). The CSP method for simplifying kinetics. International Journal of
Chemical Kinetics, 26(4), 461-486.
Le Novere, N. (2015). Quantitative and logic modelling of gene and molecular networks. Nature
Reviews. Genetics, 16(3), 146.
Li, F., Long, T., Lu, Y., Ouyang, Q., & Tang, C. (2004). The yeast cell-cycle network is robustly designed.
Proceedings of the National Academy of Sciences of the United States of America, 101(14),
4781-4786.
Ling, H., Kulasiri, D., & Samarasinghe, S. (2010). Robustness of G1/S checkpoint pathways in cell cycle
regulation based on probability of DNA-damaged cells passing as healthy cells. Biosystems,
101(3), 213-221.
23
Ling, H., Samarasinghe, S., & Kulasiri, D. (2013). Novel recurrent neural network for modelling biological
networks: oscillatory p53 interaction dynamics. Biosystems, 114(3), 191-205.
Linke, C., Chasapi, A., González-Novo, A., Al Sawad, I., Tognetti, S., Klipp, E., ... & Barberis, M. (2017). A
Clb/Cdk1-mediated regulation of Fkh2 synchronizes CLB expression in the budding yeast cell
cycle. npj Systems Biology and Applications, 3, 1.
Maas, U., & Pope, S. B. (1992). Simplifying chemical kinetics: intrinsic low-dimensional manifolds in
composition space. Combustion and Flame, 88(3), 239-264.
Maas, U., & Pope, S. B. (1994). Laminar flame calculations using simplified chemical kinetics based on
intrinsic low-dimensional manifolds. In Symposium (International) on Combustion (Vol. 25, No.
1, pp. 1349-1356). Elsevier.
Mendoza, L., & Xenarios, I. (2006). A method for the generation of standardized qualitative
dynamical systems of regulatory networks. Theoretical Biology and Medical Modelling, 3(1),
13.
Michaelis, L., & Menten, M. L. (1913). Die kinetik der invertinwirkung. Biochem. z, 49(333-369), 352.
Mombach, J. C., Bugs, C. A., & Chaouiya, C. (2014). Modelling the onset of senescence at the G1/S cell
cycle checkpoint. BMC genomics, 15(7), S7.
Müssel, C., Hopfensitz, M., & Kestler, H. A. (2010). BoolNet—an R package for generation,
reconstruction and analysis of Boolean networks. Bioinformatics, 26(10), 1378-1380.
Noel, V., Grigoriev, D., Vakulenko, S., & Radulescu, O. (2012). Tropical geometries and dynamics of
biochemical networks application to hybrid cell cycle models. Electronic Notes in Theoretical
Computer Science, 284, 75-91.
Novak, B., & Tyson, J. J. (2004). A model for restriction point control of the mammalian cell cycle.
Journal of theoretical biology, 230(4), 563-579.
Obaya A.J., & Sedivy J.M., (2002). Regulation of cyclin-Cdk activity in mammalian cells. Cellular and
Molecular Life Sciences, 59, 126-142.
Radulescu, O., Gorban, A. N., Zinovyev, A., & Noel, V. (2012). Reduction of dynamical biochemical
reaction networks in computational biology. arXiv preprint arXiv:1205.2851.
Radulescu, O., Vakulenko, S., & Grigoriev, D. (2015). Model reduction of biochemical reactions
networks by tropical analysis methods. Mathematical Modelling of Natural Phenomena, 10(3),
124-138.
Rubinstein, A., Hazan, O., Chor, B., Pinter, R. Y., & Kassir, Y. (2013). The effective application of a
discrete transition model to explore cell-cycle regulation in yeast. BMC research notes, 6(1),
311.
Sahin, Ö., Fröhlich, H., Löbke, C., Korf, U., Burmester, S., Majety, M., . . . Thieffry, D. (2009). Modeling
ERBB receptorregulated G1/S transition to find novel targets for de novo trastuzumab
resistance. BMC systems biology, 3(1), 1.
Satyanarayana, A., Kaldis, P., (2009). Mammalian cell-cycle regulation: several Cdks, numerous cyclins
and diverse compensatory mechanisms. Oncogene 28 (33), 2925–2939.
Segel, L. A., & Slemrod, M. (1989). The quasi-steady-state assumption: a case study in perturbation.
SIAM Review, 31(3), 446-477.
Singhania, R., Sramkoski, R. M., Jacobberger, J. W., & Tyson, J. J. (2011). A hybrid model of mammalian
cell cycle regulation. PLoS computational biology, 7(2), e1001077.
Stoll, G., Viara, E., Barillot, E., & Calzone, L. (2012). Continuous time Boolean modeling for biological
signaling: application of Gillespie algorithm. BMC Syst. Biol. 6: 116. doi: 10.1186/1752-0509-624
116.
Tanaka, H., Fauré, A., & Matsuno, H. (2017). Boolean modeling of mammalian cell cycle and cancer
pathways. The 2017 International Conference on Artificial Life and Robotics (ICAROB 2017),
Seagaia Convention Center, Miyazaki, Japan, (pp. 507-510).
Tashima, Y., Hanai, T., Hamada, H., Eguchi, Y., & Okamoto, M. (2007). Mathematical modeling of G2/M
phase in the cell cycle with involving the p53/Mdm2 oscillation system. In World Congress on
Medical Physics and Biomedical Engineering 2006 (pp. 197-200). Springer Berlin Heidelberg.
Terfve, C., Cokelaer, T., Henriques, D., MacNamara, A., Goncalves, E., Morris, M. K., ... & SaezRodriguez,
J.
(2012). CellNOptR:
a
flexible toolkit to using multiple logic
formalisms. BMC systems biology, 6(1), 133.
Tikhonov, A. N. (1952). Systems of differential equations containing small parameters in the
derivatives. Matematicheskii sbornik, 73(3), 575-586.
Todd, R. G., & Helikar, T. (2012). Ergodic sets as cell phenotype of budding yeast cell cycle. PloS one,
7(10), e45780.
Wei, J., & Kuo, J. C. (1969). Lumping analysis in monomolecular reaction systems. Analysis of the exactly
lumpable system. Industrial & Engineering Chemistry Fundamentals, 8(1), 114-123.
West, S., Bridge, L. J., White, M. R., Paszek, P., & Biktashev, V. N. (2015). A method of ‘speed
coefficients’ for biochemical model reduction applied to the NF-κB system. Journal of
Mathematical Biology, 70(3), 591-620.
Yu, J., Zhang, L., Hwang, P. M., Rago, C., Kinzler, K. W., & Vogelstein, B. (1999). Identification and
classification of p53-regulated genes. Proceedings of the National Academy of Sciences,
96(25), 14517-14522.
Zhao, Y., Lou, I. C., & Conolly, R. B. (2012). Computational modeling of signaling pathways mediating
cell cycle checkpoint control and apoptotic responses to ionizing radiation-induced DNA
damage. Dose-Response, 10(2), dose-response.
25
Appendix A
Initial conditions, kinetic parameters and mass balance equations of
the ODE mathematical model for the DNA damage signalling
pathway and G1/S checkpoint used as the basis for the current
study
A.1 Initial conditions
Table A.1 Initial conditions of the mathematical model used as the basis for the current study
Chemical
species
Initial value
Chemical
species
Initial value
Y1
Y2
Y3
Y4
Y5
Y6
Y7
Y8
Y9
Y10
Y11
Y12
3.00e-02
1.00e-03
4.00e-05
5.00e+00
1.50e+01
7.50e+00
1.00e-03
1.00e-03
4.00e-04
1.00e-04
1.40e+01
1.00e-03
Y13
Y14
Y15
Y16
Y17
Y18
Y19
Y20
Y21
Y22
Y23
Y24
1.00e+00
1.00e-04
0
0
0
0
1.00e-03
1.95e+00
1.00e-03
0
1.00e-02
5.00e-02
Chemical
species
Initial value
Y25
Y26
Y27
Y28
Y29
Y30
Y31
Y32
Y33
Y34
Im
2.65e-02
2.35e-04
0
1.00e-03
1.00e-04
9.90e-01
1.00e-02
0
0
0
0
Abbreviations are as follows: Y1: CycD, Y2: CycE, Y3: CycA, Y4: Cdk4, Y5: Cdk2, Y6: CycD/Cdk4, Y7:
CycE/Cdk2, Y8: CycE/Cdk2-P, Y9: CycA/Cdk2, Y10: CycA/Cdk2-P, Y11: p27, Y12: p27/CycD/Cdk4, Y13:
p27/CycE/Cdk2-P, Y14: p27/CycA/Cdk2-P, Y15: p21, Y16: p21/CycD/Cdk4, Y17: p21/CycE/Cdk2-P, Y18:
p21/CycA/Cdk2-P, Y19: p16, Y20: Rb/E2F, Y21: Rb-PP/E2F, Y22: E2F, Y23: Rb-PPP, Y24: Rb, Y25: p53,
Y26: Mdm2, Y27: ATM/ATR, Y28: Cdc25A, Y29: Cdc25A-P, Y30: Chk1, Y31: Chk1-P, Y32: NF-Y, Y33: BMyb,
Y34:
B-Myb-P,
Im:
Intermediate,
DDS:
DNA
damage
26
signal.
A.2 Kinetic parameters
Table A.2 Kinetic parameters of the ODE mathematical model used as the basis for the current
study
Kinetic
parameter
Value
Kinetic
parameter
Value
Kinetic
parameter
Value
k1
k2
k3
k4
k5
k6
k7
k8
k9
k10
k11
k12
k13
k14
k15
k16
k17
k18
k19
k20
k21
k22
k23
k24
5.00e-04
5.00e-04
5.00e-03
2.50e-03
1.00e-01
2.50e-03
2.50e-03
2.50e-05
3.00e-04
5.00e-04
5.00e-04
2.00e-04
5.00e-04
7.50e-03
5.00e-03
5.00e-03
5.00e-02
5.00e-04
5.00e-03
5.00e-04
5.00e-05
6.00e-03
1.75e-03
2.25e-02
k25
k26
k27
k28
k29
k30
k31
k32
k33
k34
k35
k36
k37
k38
k39
k40
k41
k42
k43
k44
k45
k46
k47
k48
1.75e-04
2.25e-02
1.75e-04
9.00e-04
5.00e-05
2.50e-03
1.75e-04
2.50e-03
1.75e-04
5.00e-08
5.00e-02
1.50e-03
5.00e-05
1.00e-03
5.00e-03
2.00e-03
5.00e-05
1.00e-04
5.00e-04
5.00e-04
5.00e-05
2.50e-03
2.50e-03
2.50e-03
k49
k50
k51
k52
k53
k54
k55
k56
k57
k58
k59
k60
k61
k62
k63
k64
k65
k66
k67
k68
k69
k70
k71
k72
4.00e-02
2.50e-03
5.00e-08
5.00e-07
5.00e-05
1.00e-02
5.00e-08
5.00e-05
5.00e-03
5.00e-05
5.00e-04
1.00e-04
7.00e-02
1.00e-03
9.40e-04
2.00e-02
9.50e+00
1.00e+01
5.00e-03
5.00e-02
5.00e-02
6.00e+00
4.00e-03
1.00e-08
Kinetic
parameter
Value
k73
k74
k75
k76
k77
k78
k79
k80
k81
k82
k83
k84
k85
k86
k87
k88
k89
k90
k91
k92
DDS
2.00e-03
7.72e-01
1.00e-05
5.56e-02
2.00e-02
2.00e-01
1.00e-02
4.00e-02
1.00e-03
5.00e-02
5.00e-03
1.00e-03
5.00e-03
5.00e-04
1.00e+00
1.00e+00
1.00e-03
5.00e-04
5.00e-03
3.00e-06
*
* The values of DDS were as follows: 0 (no-damage), 0.002 (low-damage), 0.004 (mediumdamage),
0.008 (high-damage), and 0.016 (excess-damage).
A.3 Biochemical meaning of the kinetic parameters
Table A.3 Biochemical meaning of the kinetic parameters of the G1/S base model
Parameter
Biochemical Meaning
K1
K2
K3
K4
K5
K6
K7
K8
K9
K10
Synthesis rate of CycD
Degradation rate of CycD
Association rate of CycD/CDK4
Dissociation rate of CycD/CDK4
Synthesis rate of CycE through E2F
Degradation rate of CycE
Association rate of CycE/CDK2
Dissociation rate of CycE/CDK2
Synthesis rate of CycA through B-Myb-P
Degradation rate of CycA
27
K11
K12
K13
K14
K15
K16
K17
K18
K19
K20
K21
K22
K23
K24
K25
K26
K27
K28
K29
K30
K31
K32
K33
K34
K35
K36
K37
K38
K39
K40
K41
K42
K43
K44
K45
K46
K47
K48
K49
K50
K51
K52
K53
K54
K55
K56
K57
K58
K59
K60
K61
Association rate of CycA/CDK2
Dissociation rate of CycA/CDK2
Rate of CDK4 production through CycD/CDK4
Rate of CDK2 production through CycA/CDK2-P
Rate of CDK2 production through CycA/CDK2
Rate of CDK2 production through CycE/CDK2
Rate of CDK2 production through CycE/CDK2-P
Association rate of p21/CycD/CDK4
Disassociation rate of p21/CycD/CDK4
Association rate of p27/CycD/CDK4
Disassociation rate of p27/CycD/CDK4
Phosphorylation rate of CycE/CDK2 to form CycE/CDK2-P
De-phosphorylation rate of CycE/CDK2-P to form CycE/CDK2
Association rate of p27/CycE/CDK2-P
Disassociation rate of p27/CycE/CDK2-P
Association rate of p21/CycE/CDK2-P
Disassociation rate of p21/CycE/CDK2-P
Phosphorylation rate of CycA/CDK2 to form CycA/CDK2-P
De-phosphorylation rate of CycA/CDK2-P to form CycA/CDK2
Association rate of p27/CycA/CDK2-P
Disassociation rate of p27/CycA/CDK2-P
Association rate of p21/CycA/CDK2-P
Disassociation rate of p21/CycA/CDK2-P
Synthesis rate of p27
Association rate of p27/CycE/CDK2-P
Association rate of p27/CycA/CDK2-P
Synthesis rate of p21
Rate of synthesis of p21 through p53
Degradation rate of p21
Synthesis rate of p16
Constant as influx or precursor
Rate of inhibition of synthesis p16 by Rb
Degradation rate of p16
Degradation rate of CycD/CDK4 by p16
Association rate of Rb/E2F
Phosphorylation rate of Rb/E2F to form Rb-PP/E2F through CycD/CDK4
Phosphorylation rate of Rb/E2F to form Rb- PP/E2F through p27/CycD/CDK4
Phosphorylation rate of Rb/E2F to form Rb-PP/E2F through p21/CycD/CDK4
Rate of activation of E2F by CycE/CDK2-P
Rate of activation of E2F by CycA/CDK2-P
Rate of synthesis of E2F promoted by itself
Synthesis rate of E2F
Degradation rate of E2F
Rate of E2F degradation by CycA/CDK2-P
De-phosphorylation rate of Rb-PPP to Rb
Synthesis rate of Rb
Degradation rate of Rb
Constant as influxes or precursor
Rate of inhibition of synthesis Rb by p16
Synthesis rate of p53
Rate of synthesis of p53 through ATM/ATR
28
K62
K63
K64
K65
K66
K67
K68
K69
K70
K71
K72
K73
K74
K75
K76
K77
K78
K79
K80
K81
K82
K83
K84
K85
K86
K87
K88
K89
K90
K91
K92
Degradation rate of p53
Synthesis rate of Mdm2
Degradation rate of Mdm2
Dissociation constant in the Hill function
Rate of synthesis of Mdm2 through Im
Degradation rate of Im
Synthesis rate of B-Myb through E2F
Phosphorylation rate of B-Myb to form B-Myb-P through CycA/Cdk2-P
Rate of p53’s sequence-specific DNA binding activity by DNA-damage signal
Association rate of p53 and Mdm2
Rate of DNA-damage repair
Degradation rate of B-Myb-P
Rate of inhibition of degradation of p53 and/or Mdm2 by DNA-damage signal
Synthesis rate of CycA through NF-Y
Strength rate of Mdm2’s ability to promote p53 degradation
Rate of inhibition of Mdm2-mediated p53 degradation under the initial damage signal
Rate of synthesis of ATM/ATR through DNA-damage signal
Degradation rate of ATM/ATR
Rate of Cdc25A production through E2F
Degradation rate of Cdc25A through Chk1-P
Phosphorylation rate of Cdc25A to form Cdc25A-P through CycA/Cdk2-P and
CycA/Cdk2-P
Degradation rate of Cdc25A
Degradation rate of Cdc25A-P through Chk1-P
De-phosphorylation rate of Cdc25A-P to form Cdc25A through CycE/Cdk2to be
CycE/Cdk2-P
De-phosphorylation rate of Cdc25A-P to form Cdc25A through CycA/Cdk2to be
CycA/Cdk2-P
De-phosphorylation rate of Chk1-P to form Chk1
Phosphorylation rate of Chk1 to form Chk1-P
Synthesis rate of NF-Y through CycA/Cdk2-P
Degradation rate of NF-Y
Degradation rate of CycA
Synthesis rate of CycA through E2F
29
A.4 Mass balance equations of the ODE mathematical model
Table A.4 Mass balance equations of the ODE mathematical model used as the basis for the current
study
CycD dY1/dt = k1 + k4Y6 − (k2 + k3Y4)Y1
CycE dY2/dt = k5Y22 + k8Y7 − (k6 + k7Y5)Y2
CycA dY3/dt = k9Y34 + k92Y22 + k12Y9 + k75Y32 − (k10 + k11Y5 + k91)Y3
Cdk4 dY4/dt = k4Y6 + k13Y6 − k3Y1Y4
Cdk2 dY5/dt = k8Y7 + k12Y9 + k14Y10 + k15Y9 + k16Y7 + k17Y8Y8 − (k7Y2 + k11Y3)Y5
CycD/Cdk4 dY6/dt = k3Y1Y4 + k19Y16 + k21Y12 − (k4 + k13 + k18Y15 + k20Y11 + k44Y19)Y6
CycE/Cdk2 dY7/dt = k7Y2Y5 + k23Y8 − (k8 + k22Y29 + k16)Y7
CycE/Cdk2-P dY8/dt = k22Y7Y29 + k25Y13 + k27Y17 − (k23 + k24Y11 + k26Y15 + k17Y8)Y8
CycA/Cdk2 dY9/dt = k11Y3Y5 + k29Y10 − (k12 + k28Y29 + k15)Y9
CycA/Cdk2-P dY10/dt = k28Y9Y29 + k31Y14 + k33Y18 − (k29 + k30Y11 + k32Y15 + k14)Y10
p27 dY11/dt = k34 + k21Y12 + k25Y13 + k31Y14 − (k35Y8 + k36Y10 + k20Y6 + k24Y8 + k30Y10)Y11
p27/CycD/Cdk4 dY12/dt = k20Y6Y11 − k21Y12
p27/CycE/Cdk2-P dY13/dt = k24Y8Y11 − k25Y13
p27/CycA/Cdk2-P dY14/dt = k30Y10Y11 − k31Y14
p21 dY15/dt = k37 + k38Y25 + k19Y16 + k27Y17 + k33Y18 − (k39 + k18Y6 + k26Y8 + k32Y10)Y15
p21/CycD/Cdk4 dY16/dt = k18Y6Y15 − k19Y16
p21/CycE/Cdk2-P dY17/dt = k26Y8Y15 − k27Y17
p21/CycA/Cdk2-P dY18/dt = k32Y10Y15 − k33Y18
p16 dY19/dt = k40 + k41/(1 + k42Y24) − (k43 + k44Y6)Y19
Rb/E2F dY20/dt = k45Y22Y24 − (k46Y6 + k47Y12 + k48Y16)Y20
Rb-PP/E2F dY21/dt = k46Y6Y20 + k47Y12Y20 + k48Y16Y20 − (k49Y8 + k50Y10)Y21
E2F dY22/dt = k49Y8Y21 + k50Y10Y21 + k51Y22 + k52 − (k45Y24 + k53 + k54Y10)Y22
Rb-PPP dY23/dt = k49Y8Y21 + k50Y10Y21 − k55Y23
Rb dY24/dt = k56 + k58/(1 + k59Y19) + k55Y23 − (k57 + k45Y22)Y24
p53 dY25/dt = k60 + k61Y27 − (deg (t)Y26 + k62)Y25
Mdm2 dY26/dt
= k + (k Im50)/(k 50 + Im50) − k Y
63
66
65
64 26
ATM/ATR dY27/dt =
Cdc25A dY28/dt =
k78sig (t) − k79Y27
k80Y22 + k85Y29 − (k81Y31 + k82(Y8 + Y10) + k83)Y28
Cdc25A-P dY29/dt = k82(Y8 + Y10)Y28 − (k84Y31 + k85 + k86)Y29
Chk1 dY30/dt = k87Y31 − k88Y27Y30
Chk1-P dY31/dt = k88Y27Y30 − k87Y31
NF-Y dY32/dt = k89Y10 − k90Y32
B-Myb dY33/dt = k68Y22 – k69Y10Y33
B-Myb-P dY34/dt = k69Y10Y33 – k73Y34
Im dIm/dt = k70Y25sig(t)/(1 + k71Y25Y26) − k67Im
Sig sig(t) = DDS × exp(−k72×time)
Deg deg(t) = k76 − k74×(sig(t) – DDS × exp(−k77×DDS×time))
Abbreviations are as shown in Table A1.
30
Appendix B
ODEs, initial conditions and logical equations of the reduced model
(R1) (used for time slicing and logical modelling)
B.1 ODEs of the reduced model (R1)
Table B.1 ODEs of the reduced model (R1)
1
CycD
dY1/dt = k1 + k4Y6 − (k2 + k3Y4)Y1
2
Cdk4
dY4/dt = k4Y6 + k13Y6 − k3Y1Y4
3
CycD/Cdk4
4
p27
5
p27/CycD/Cdk4
6
p21
7
p21/CycD/Cdk4
8
p16
9
Rb/E2F
10
Rb-PP/E2F
dY21/dt = k46Y6Y20 + k47Y12Y20 + k48Y16Y20
11
Rb
dY24/dt = k56 + k58/(1 + k59Y19) − (k57)Y24
12
p53
13
Mdm2
14
ATM/ATR
15
Im
16
Sig
17
Deg
dY6/dt = k3Y1Y4 + k19Y16 + k21Y12 − (k4 + k13 + k18Y15 + k20Y11 + k44Y19)Y6
dY11/dt = k34 + k21Y12 − (k20Y6)Y11
dY12/dt = k20Y6Y11 − k21Y12
dY15/dt = k37 + k38Y25 + k19Y16 − (k39 + k18Y6)Y15
dY16/dt = k18Y6Y15 − k19Y16
dY19/dt = k40 + k41/(1 + k42Y24) − (k43 + k44Y6)Y19
dY20/dt = − (k46Y6 + k47Y12 + k48Y16)Y20
dY25/dt = k60 + k61Y27 − (deg (t)Y26 + k62)Y25
dY26/dt = k63 + (k66Im50)/(k6550 + Im50) − k64Y26
dY27/dt = k78sig (t) − k79Y27
dIm/dt = k70Y25sig(t)/(1 + k71Y25Y26) − k67Im
sig(t) = DDS × exp(−k72×time)
deg(t) = k76 − k74×(sig(t) – DDS × exp(−k77×DDS×time))
Abbreviations are as follows: Y1: CycD, Y4: Cdk4, Y6: CycD/Cdk4, Y11: p27, Y12: p27/CycD/Cdk4, Y15:
p21, Y16: p21/CycD/Cdk4, Y19: p16, Y20: Rb/E2F, Y21: Rb-PP/E2F, Y24: Rb, Y25: p53, Y26: Mdm2, Y27:
ATM/ATR, Im: Intermediate, DDS: DNA damage signal.
31
B.2 Initial conditions of the reduced model (R1)
Table B.2 Initial conditions of the reduced model (R1)
No
Chemical species
Initial value
No
Chemical species
Initial value
1
2
3
4
5
6
7
CycD
Cdk4
CycD/Cdk4
p27
p27/CycD/Cdk4
p21
p21/CycD/Cdk4
0
1
0
1
0
0
0
8
9
10
11
12
13
14
p16
Rb/E2F
Rb-PP/E2F
Rb
p53
Mdm2
ATM/ATR
0
1
0
1
0
0
0
B.3 Logical equations of the reduced model (R1)
Table B.3 Logical equations of the reduced model (R1)
1
CycD
2
3
4
5
6
7
8
9
10
11
12
13
14
Cdk4
CycD/Cdk4
p27
p27/CycD/Cdk4
p21
p21/CycD/Cdk4
p16
Rb/E2F
Rb-PP/E2F
Rb
p53
Mdm2
ATM/ATR
= GFs And Not (CycD And Cdk4)
= Cdk4 And Not (CycD And Cdk4)
= (CycD/Cdk4 Or (CycD And Cdk4)) And Not p16
= p27
= (p27 And CycD/Cdk4) Or p27/CycD/Cdk4
= p53
= (p21 And Cyc/DCdk4) Or p21/CycD/Cdk4
= p16 Or (Rb And CycD/Cdk4)
= Not (CycD/Cdk4 Or p27/CycD/Cdk4 Or p21/CycD/Cdk4)
= Rb-PP/E2F Or (CycD/Cdk4 And Rb/E2F)
= Not p16
= ATM/ATR And Not Mdm2
= p53 And Not Mdm2
= Sig
*** GFs: Growth factors equals 1 at t=1, Sig: DNA damage signal (0 or 1).
32
Appendix C
ODEs, initial conditions and logical equations of the reduced model
(R2) (used for time slicing and logical modelling)
C.1 ODEs of the reduced model (R2)
Table C.1 ODEs of the reduced model (R2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
CycE dY2/dt = k5Y22 + k8Y7 − (k6 + k7Y5)Y2
Cdk2 dY5/dt = k8Y7 + k16Y7 + k17Y8Y8 − (k7Y2)Y5
CycE/Cdk2 dY7/dt = k7Y2Y5 + k23Y8 − (k8 + k16)Y7
CycE/Cdk2-P dY8/dt = k25Y13 + k27Y17 − (k23 + k24Y11 + k26Y15 + k17Y8)Y8
p27 dY11/dt = k34 + k25Y13 − (k35Y8 + k24Y8)Y11
p27/CycE/Cdk2-P dY13/dt = k24Y8Y11 − k25Y13
p21 dY15/dt = k37 + k38Y25 + k27Y17 − (k39 + k26Y8)Y15
p21/CycE/Cdk2-P dY17/dt = k26Y8Y15 − k27Y17
Rb-PP/E2F dY21/dt = − (k49Y8)Y21
E2F dY22/dt = k49Y8Y21 + k51Y22 + k52
Rb-PPP dY23/dt = k49Y8Y21 + k55Y23
p53 dY25/dt = k60 + k61Y27 − (deg (t)Y26 + k62)Y25
Mdm2 dY26/dt = k63 + (k66Im50)/(k6550 + Im50) − k64Y26
ATM/ATR dY27/dt = k78sig (t) − k79Y27
Im dIm/dt = k70Y25sig(t)/(1 + k71Y25Y26) − k67Im
Sig sig(t) = DDS × exp(−k72×time)
Deg deg(t) = k76 − k74×(sig(t) – DDS × exp(−k77×DDS×time))
Abbreviations are as follows: Y2: CycE, Y5: Cdk2, Y7: CycE/Cdk2, Y8: CycE/Cdk2-P, Y11: p27, Y13:
p27/CycE/Cdk2-P, Y15: p21, Y17: p21/CycE/Cdk2-P, Y21: Rb-PP/E2F, Y22: E2F, Y23: Rb-PPP, Y25: p53,
Y26: Mdm2, Y27: ATM/ATR, Im: Intermediate, DDS: DNA damage signal.
33
C.2 Initial conditions of the reduced model (R2)
Table C.2 Initial conditions of the reduced model (R2)
No
Chemical species
Initial value
No
Chemical species
Initial value
1
CycE
0
8
p21/CycE/Cdk2-P
0
2
Cdk2
1
9
Rb-PP/E2F
1
3
CycE/Cdk2
0
10
E2F
0
4
CycE/Cdk2-P
0
11
Rb-PPP
0
5
p27
1
12
p53
0
6
p27/CycE/Cdk2-P
0
13
Mdm2
0
7
p21
0
14
ATM/ATR
0
*** Initial values of P53, p21, Mdm2 and ATM/ATR are 1 if there is DNA damage.
C.3 Logical equations of the reduced model (R2)
Table C.3 Logical equations of the reduced model (R2)
1
CycE
= E2F And Not (CycE And Cdk2)
2
Cdk2
= Cdk2 And Not (CycE And Cdk2)
3
CycE/Cdk2
= CycE/Cdk2 Or (CycE And Cdk2)
4
CycE/Cdk2-P
5
p27
6
p27/CycE/Cdk2-P
7
p21
8
p21/CycE/Cdk2-P
9
Rb-PP/E2F
10
E2F
11
Rb-PPP
= (p21/CycE/Cdk2-P Or p27/CycE/Cdk2-P Or CycE/Cdk2)
= p27 And Not (p27 And CycE/Cdk2-P)
= p27/CycE/Cdk2-P Or (p27 And CycE/Cdk2-P)
= p53
= p21/CycE/Cdk2-P Or (p21 And CycE/Cdk2-P)
= Not (CycE/Cdk2-P And Rb-PP/E2F)
= (CycE/Cdk2-P And Rb-PP/E2F) Or E2F
= (CycE/Cdk2-P And Rb-PP/E2F) Or Rb-PPP
(Without DNA Damage)
12
p53
13
Mdm2
14
ATM/ATR
= ATM/ATR And Not Mdm2
= p53 And Not Mdm2
= Sig
*** Sig: DNA damage signal (0 or 1).
*** Rb-PPP = ( (CycE/Cdk2-P And Rb-PP/E2F) Or Rb-PPP ) And Not(p21) (With DNA Damage)
34
Appendix D
Steps to convert mathematical equations (ODEs) to logical
equations
1. Take the ODE equation for the protein from the base model.
Ex: ODE for CycD (in R1 model)
dY1/dt = k1 + k4Y6 − (k2 + k3Y4)Y1
2. Remove the proteins that are not included in the reduced model (R1) and and remove
terms that have a very small effect.
Equation for CycD will become
dY1/dt = k1 + k4Y6 − k3Y4Y1
We remove − k2Y1 because degradaZon of CycD in G1 is very small and no other terms are removed
because all proteins Y6: CycD/Cdk4, Y1: CycD and Y4: Cdk4 are included in the reduced model (R1)
3. Draw the regulatory graph for this equation
4. Convert the ODE equ. to Boolean equ. based on the following rules (Mendoza, L., &
Xenarios, I. (2006):
( + 1) = ( 1 ( )
2 ( )…
( ))
( 1()
2 ( )…
( ))
(1. )
( + 1) =
( + 1) =
1 ()
( 1()
2 ( )…
2 ( )…
()
(1. )
( ))
(1. )
35
Explain for 1.a, 1.b and 1.c rules:
A. (+) in ODE equ. it will be (or) in Boolean equ.
Ex: ODE equ. x= y + z
Boolean equ. x = y or z
B. (.) in ODE equ. it will be (and) in Boolean equ.
Ex: ODE equ. x= yz + k
Boolean equ. x = (y and z) or k
C. (-) in ODE equ. it will be (and not) in Boolean equ.
Ex: ODE equ. x= y - z
Boolean equ. x = y and (not z)
------------------------------------------------------------------------------------------------------For CycD ODE equation:
dY1/dt = k1 + k4Y6 − k3Y4Y1
CycD Boolean equ:
CycD
= (GFs or CycD/Cdk4) And Not (CycD And Cdk4)
5. Apply Boolean Algebra laws to reduce the Boolean equation:
Ex: CycD = (GFs or CycD/Cdk4) And Not (CycD And Cdk4)
The initial value of GFs= 1 and stay 1 for all time steps in the reduced model (R1). When applying
Boolean Algebra laws:
(A + 1) = 1
(GFs or CycD/Cdk4)
(1 or CycD/Cdk4)
(1)
GFs
we reduce the equ. to be:
CycD
= GFs And Not (CycD And Cdk4)
6. Draw the regulatory graph for this equ.
36