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