{'fraction_non_alphanumeric': 0.03104424545693969, 'fraction_numerical': 0.01435343692388728, 'mean_word_length': 4.644184318097362, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
From differential equations to Boolean networks: A Case Study in modeling regulatory networks

Methods of modeling cellular regulatory networks as diverse as differential equations and Boolean networks co-exist, however, without any closer correspondence to each other. With the example system of the fission yeast cell cycle control network, we here set the two approaches in relation to each other. We find that the Boolean network can be formulated as a specific coarse-grained limit of the more detailed differential network model for this system. This lays the mathematical foundation on which Boolean networks can be applied to biological regulatory networks in a controlled way.


I. INTRODUCTION
The task of molecular cell biology is to comprehend the control of cellular processes of living cells encoded in the genome of the cell. These cellular processes are guided by sophisticated networks of interactions between the macromolecules of the cell as proteins, nucleic acids and polysaccharides. Their complexes and structures define the unique properties that enable them to perform the functions of the cell as, for example, catalysis of chemical transformations, production of movement, and heredity. The complexity of these processes demands not only advanced experimental techniques, but also adequate mathematical and computational models for understanding them (Gunsalus et al., 2005;Riel, 2006).
Today, there are different methods for modeling the complex networks of biochemical interactions, ranging from master equations based on first principles (Monte-Carlo method) (Gillespie, 1976;Gillespie, 1977), ordinary differential equations (ODE) ( Thomas et al., 1995).
Among these methods, a most popular approach to modeling biochemical networks is via differential equations, based on the known chemical kinetics and successfully applied to describing numerous processes in living organisms (Chen et al., 2000;Novak and Tyson, 1993;Novak and Tyson, 1997;Novak et al., 2001;Novak and Tyson, 2004;Tyson et al., 2002). To build an ODE model, one starts with a schematic diagram representing the known interactions between components. Then this diagram is converted into a set of differential and algebraic equations. The full ODE model then consists of this set of rate equations, plus a set of parameter values as well as a set of initial conditions. The solutions of the ODEs give the time-dependence of each component of the system. In practice these solutions depend on rather detailed knowledge about all reactions and kinetic parameters.
In studies where prediction of exact reaction times is not an issue, simpler models than ODE and less parameters may be necessary for predicting the course of events in a regulatory network. For example, relevant features of cell commitment, cell cycle progression, and cell differentiation are already described in terms of a sequence of regulatory events (Albert and Othmer, 2003;Braunewell and Bornholdt, 2006;Davidich and Bornholdt, 2008;Li et al., 2004;Sanchez et al., 2001). In such cases, the much simpler modeling framework of Boolean networks may be a suitable method (Kaufmann, 1969;Thomas, 1973). Constructing a Boolean network model starts from a wiring diagram of interactions between biochemical elements as well, but no kinetic details are needed. Interactions are classified into just two classes, activation or inhibition, as well as the concentration levels being reduced to just an ON or OFF state.
Despite their extreme simplicity, such Boolean models are able to reproduce regulatory sequences, for example, models of genetic networks of A. thaliana, (Espinosa-Soto et al., 2004; Mendoza et al., 1999;Thum et al., 2003), the cell-cycle network of S. cerevisiae (Li et al., 2004), the mammalian cell-cycle (Faure et al., 2006), and the segment polarity gene network in D. melanogaster (Albert and Othmer, 2003;Sanchez et al., 2001). These examples show that the Boolean network approach provides reliable results for different organisms.
The two diverse methods just summarized (ODE and Boolean networks) are both based on the same "wiring" diagram of interactions between the components, however, use much different amounts of information about these interactions. This poses the interesting question how these two methods are related to each other. A first correspondence between Boolean network and ordinary differential equations has been drawn by Glass and collaborators (Glass and Kauffman, 1973;Glass and Hill, 1998) who explored the relationship between a class of non-linear equations representing biochemical control networks and homologous switching networks. They argued that such a correspondence can be achieved with the following requirements: 1) Rates of reactions are described by monotonic sigmoidal functions having distinct upper and lower asymptotes. 2) The parameters must be defined to match the upper or lower asymptote. 3) The target control function must correspond to the maximal or basal rate of biochemical processes. They subsequently demonstrated that there is a large variety of such functions as, e.g., the Heaviside function, the error function, or the Hill function defined for positive arguments.
This leads to a mapping between asymptotical solutions of the ODE system and the Boolean system, while omitting the exact way of transitions between dynamical states.
In this paper, we further explore the correspondence between ODE and Boolean network models considering a specific biological system and demonstrate how a working Boolean model can be derived in terms of a mathematically well defined coarse-grained limit of an underlying ODE model. As our working example we choose the fission yeast cell-cycle control network (Schizosaccharomyces Pombe). The division cell-cycle consists of four phases G1-S-G2-M, during which DNA is replicated and the cell divides itself into two cells. The main role is played by "Cyclin-Dependent-Kinases" (CDKs) and cyclins that bind to CDKs to form complexes. CDKs, while being present at all times, can only be active in complexes with cyclins. Cyclins are synthesized or degraded depending on other regulatory activities. Another important participant of the process is an enzyme complex called the "Anaphase-Promoting Complex" (APC), which targets cyclins for degradation. In summary, "to understand the molecular control of cell reproduction is to understand the regulation of CDK and APC activities" (Tyson et al., 2002). These processes, while being complicated, have been well studied for fission yeast S. Pombe and successful ODE models  exist. We choose the most widespread version of the model , as often cited in textbooks, which will serve us as a starting point for our study.
The paper is organized as follows. In Section 2 we show the passage from the ODE system of algebraic differential equations for the fission yeast cell cycle, to the limit of the corresponding Boolean model which we construct. Here also the difficulties that one can meet when works with Boolean approaches are discussed. Section 3 explores the dynamics of the derived Boolean model of the fission yeast cell-cycle.
Finally, in the discussion section the properties of the obtained system are recapitulated and the Boolean and ODE approaches are compared.

II. BOOLEAN VARIABLES AS STATIONARY STATES
The passage from a differential equations model to a Boolean model requires the mapping of continuous solutions  into the ON/OFF states of a Boolean network's nodes. In order to achieve this, the time evolution of a function, determined by the rate functions and kinetic constants, has to be replaced with a discrete mapping of the node set into itself. Moreover, the rules of this self-mapping have to be governed by logical functions, connecting the binary states of interacting nodes. The dynamics of the resulting Boolean network is an ordered sequence of states of the network nodes, instead of the continuous time output of the ODE model. Having this in mind, let us find the conditions, which allow to perform the transition from a differential equations model to a Boolean system. In the following, we will first describe the passage of the continuous variables to discrete states and, in a second step, construct the logical functions representing the dynamics.

A. Stationary states of ODE system
The model of the fission yeast cell-cycle ) is based on the antagonist interaction of CDKcyclin complex with AP C (via proteins Ste9 and Slp1) and CDK inhibitor Rum1. The CDK-cyclin complex (Cdc2/Cdc13) is represented by two variables -preM P F and M P F ("Maturation Promoting Factor"). Furthermore, helper molecules (Start Kinase SK, transcription factor T F , kinase k wee , phosphatase k 25 , and time-delay enzyme IE) participate in the process. Regulatory interaction between these macromolecules are described by two types of equations -differential equations for Ste9, Slp1, IEP , M , SK, Rum1, preM P F and algebraic equations for k wee , k 25 , T F and M P F variables. Among the first ones, some are of Michaelis-Menten type (Ste9, Slp1, IEP ) and some are exponential growth (M , SK) equations. The mass variable M plays a special role as it parametrizes the time evolution in the system. The ODE model ) uses arbitrary units for concentrations in all equations, since there are few data of actual protein concentrations. The kinetic constants deretmine the right timing of the processes. Solutions of the system show that the concentrations of the major proteins in general rise or decrease steeply.
To make the transition to a Boolean system, we first need to rescale the differential equations such that their solutions assume values between 0 (inactive) and 1 (maximum activity). This is a first step towards mapping these variables onto Boolean OFF/ON variables with values 0 and 1. The rescaling does not change the form of equations, it only affects the values of kinetic constants. To do this, let us divide all functions by their respective maximum value. For example, for Slp1 we introduce the new rescaled function Slp1 1 =Slp1Ampl, where Ampl = 2.1 is the amplitude of the original solution. Rescaling all variables except M we obtain: where the Goldbeter-Koshland (GK) (Goldbeter and Koshland, 1981;Novak et al., 2001) function has the following general form: . (16) Square brackets denote the concentrations of their elements. The subscript 1 marks the rescaled variables with maximum 1. The new values of parameters are shown in Table 1.
Next we map the continuous solution of the ODE model ) into the discrete states of a Boolean network's nodes. Since in a Boolean model there is no continuous time, but rather a sequence of switching events between two stationary states of the nodes, one needs to reduce (whereever possible) the initial dynamics of the system to a sequence of the evolution of stationary states. For this, one can use the results of a bifurcation analysis  of the transitions during the cell cycle. Thus, there are some variables (Ste9, Slp1, IEP ) that are described by Goldbeter-Koshland (GK) functions  in the stationary state: The characteristic properties of the GK function imply that its variable mainly resides in two limiting states: High and Low. The transition time between them is short (as the variables J are small). Therefore, we approximate them as Boolean (binary) variables.
It is easy to see that Slp1 T 1 determines only the amplitude and the smoothness of the transition in (19), therefore, we neglect Slp1 T 1 as a first step towards a Boolean model and write The states of the remaining variables can be described by Golbeter-Koshland functions, as well. For SK, while the equation for this variable is exponential, one can evaluate the right-hand part through a GK function with the stationary solution : Furthermore, there are three algebraic equations for T F (10), k wee (11), and k 25 (12) that contain Golbeter-Koshland functions. Here, again, the two limiting states of the GK function will be related to the binary ON/OFF version of the corresponding variables in the Boolean limit.
Finally, we will simplify the functional behavior for the Cdc13 T , preM P F , and Rum1 T , as well. Again we want to neglect the exact path of their transitions, keeping the limiting stationary states, eventually enabling us to take the limit of Boolean functions as a simplified description of the dynamics. Equations with the following requirements will allow us to take this limit: a) In a small neighborhood of the switching point, the functional behavior can be approximated by an exponential rise. b) On the larger interval, it has a stationary solution with the steep transition between two limiting stationary states.
c) This function converges to the Heaviside function in the limit of steep transition. In the following we will see that the Michaelis-Menten dynamics fulfils these requirements, resulting in exact conformity of initial and final states and permitting a well controlled passage to a Boolean function. Let us start from the Michaelis-Menten equation and first check the condition a). Expanding (22) in the neighborhood of switching points where X << J 1,2 << 1 and keeping leading order terms yields This is a common equation of exponential growth/decrease. This allow us to take equations of exponential growth as an expansion of (23) in the neighborhood of switching ON/OFF points. For Cdc13 T , for example, the equation is the expansion of the equation For illustration, in Fig. 1, this function is compared to the initial Cdc13 T 1 . Condition b) is satisfied as well, since the stationary states of (22) are described by a Goldbeter-Koshland function. Validity of condition c) is shown in the next section.
Thus, for the above equations we have a system of GK-functions which are responsible for the transitions between stationary states: One can see in Fig. 1 that the new and initial functions start to grow and start to decrease at the same times, respectively. Note that the obtained substituted equations (26-28) play only a helper role and cannot be directly applied to the initial system of differential equations. The complete correspondence of the obtained system to the initial one is achieved by the limit transition shown in section II B. Summarizing all information above we obtain the system of equations (10)(11)(12), (17)(18), (20)(21), (26)(27)(28), describing the stationary states of the corresponding variables. The next step is to perform a transition from the continuous functions to discrete functions. In functions with continuous time, the stepwise change corresponds to the Heaviside function.

B. Passage to Boolean variables
We now show that there is an exact passage from the function of Goldbeter-Koshland (Goldbeter and Koshland, 1981) to the indicator function (Heaviside step function). Let us remark that in (16)  3. Very small parameters c and d mean that the enzymesubstrate complex is tightly bound and hardly dissociates. Thereby in  an assumption is made that the enzyme-substrate complexes involved are very stable (Aguda, 2006). For this reason, let us consider the behavior of the corresponding GK-functions in the limiting case c → 0, d → 0 while a and b take finite values.
First note that, both, numerator and denominator depend on d. Moreover, a is a numerator's factor. At the same time, c appears in a sum with the finite terms in the denominator, only (except at the point b = a). Therefore, we can assume, without loss of generality, that c = 0 and consider below .
Using the Taylor expansion of the square root in the denominator and neglecting higher powers of d, we obtain for all points a = b There are two possible cases: • a) if b − a > 0 then |b − a| = b − a and (29) takes a form ) and (29) is simply G(a,b,0,d)=1.
This implies that there is a passage to the limit from Goldbeter-Koshland function G(a, b, c, d) to the Heaviside function θ(a − b): Thus, in this limit, the variable a plays the role of an activator input and the variable b is an inhibitor input. The output is active/inactive (its Boolean value is equal to one/zero) if the total value of activator inputs is larger/smaller than the total value of inhibitor inputs. Thus, the Goldbeter-Koshland-function converges to the Heaviside function in the limit of steep transitions.

C. Logical Boolean functions
Let us now rewrite the system of equations in the limit of the parameters J → 0 (J 3 , J 5 , J 7 , J 8 , J 1 0, J 16 , J awee , J iwee , J i25 , J a25 and the corresponding parameters J in equations (20), (26)(27)(28) equations) as Let us add two simplifications to the equations (32)  3) If Cdc13 = 1, preM P F = 0, and Rum1 = 1, then M P F = 1, with its value slightly decreasing to M P F = 0.93 if Rum1 = 0. This corresponds to a high level of activation, when M P F is activated by Cdc13 and this activation is not reduced by preM P F .
Let us reformulate these rules in the following. Assume there are two variables -M P F and M P F 2. The first one is activated by Cdc13. For activation of the second variable M P F 2, one assumes that M P F has to be present as a low-level and preM P F should be inactive. Thereby, one needs to rewrite the system of equations (30), (36-37) taking into account which level of excitation of M P F is crucial for each particular variable.
Second, as in the model based on differential equations the cell mass M takes a special role in the present model. The solution  treats it during a cell growth as an independent variable, which is described by an exponential growth function (9). Thus, the variable M corresponds to a time in this system, which drives the evolution between stationary states ). In the system, M directly influences Cdc13 and T F . As soon as M reaches a threshold value, it activates Cdc13 and induces the sequence of consecutive transitions between stationary states. For T F it plays a role of constantly positive input, T F is always active unless M P F has a high activity.
As a criterium for the end of the cycle, Novak et al. ) determine when the cell divides by monitoring the values of the other variables. When these chosen variables have certain values that indicate the end of the cell cycle, the current value of M is divided by two manually, as at the end of mitosis the cell divides into two daughter cells of approximately equal masses. Subsequently, M continues its exponential growth, again.
Following this strategy, one needs to distinguish M between two principal different values -M and 2M in the Boolean model. Here M works at the beginning of the cell cycle as a trigger of switching events, whereas 2M play a role of an indicator for the end of the cell cycle. Correspondingly, M becomes 2M at the end of mitosis, when Slp1, Ste9 and IEP all have high concentrations.
Thus, one can add the following Boolean rule: (53) Thus, we have a system of equations (43-53), where each variable can take values 0 or 1, only. It is easy to simplify this system, reducing the kinetic coefficients to 0 or 1 and adding thresholds. Consider, for example, Cdc13 T . In Table 2, based on equation (36), we show all possible cases for Cdc13. In a more compact form, where the kinetic constants are reduced to 1, these rules become (54) Repeating the same procedure for all variables, we obtain the system of equations:

III. BOOLEAN MODEL
We now have a system of algebraic equations (54-63), which describe the switch-like transitions between stationary states, plus Boolean equations (43-44),(64-65) for M and M P F . Note that in this discrete system, no information about continuous time is present any more, except the sequence of events. To obtain this discrete dynamical sequence, we iteratively solve the system. We start from the known initial conditions : 2M = 1, Slp1 1 = 1, IEP 1 =1, Ste9 1 = 1, k wee1 = 1, with all other variables being 0. Following the terminology of Boolean models, each variable is represented by one node. The network of nodes is shown in Fig.2. Each node i has only two states, S i (n) = 1 (active) and S i (n) = 0 (inactive). The index n is the number of iterations. The iterative solution of the system has the following general form: with the Heaviside function All nodes are updated synchronously, which corresponds to the iteration of the full dynamical system. The interaction matrix T ik and the state vector Q i determine the transition rule between states.
The specific values of the interactions T ik are determined as follows. All interactions between a pair of nodes are defined by an interaction strength given as the elements T ik and an activation threshold Q i . Positive (negative) arguments of the θ-functions have T ik = 1 (T ik = −1). T ii = 1 for M , this rule is true only for M node, since it is described by a growing exponential function.
The resulting matrix T ik and vector Q i have the following form:

IV. RESULTS OF BOOLEAN SIMULATION OF THE FISSION YEAST CELL CYCLE
First we run the model described in the previous section with the initial conditions . The temporal evolution of the protein states is presented in Table 3. One can see that the iterative solution of the system (43-44) (54-65) is the switching between unstable stationary states, which coincide with the corresponding evolution in the ODE model. The final state is a stable stationary state of the system. One notices that the initial and end states are identical except for the activation of the nodes M and 2M . The update of nodes M and 2M , keeping all other nodes in the same states, starts the new cycle. This cycling of the model is similar to the realization of cycling in the original model of differential equations.
Let us briefly summarize our coarse-graining strategy that we followed in this article. In Fig. 3, we show consecutive abstractions of the model for the Ste9 and IEP variables as an example. For this we first plot the dynamics of the differential equations model, then the evolution between the stationary states solutions of the ODE model, and finally the sequence of states obtained from the iterated Boolean network model.
In the next step we run the model starting from each of the 2 15 possible initial states. We find that from all initial states 67% flows into one big attractor. This attractor is the same stable stationary state that one obtains starting with biological conditions described above.

A. Mutations
Let us also compare the behavior of mutants. We model two mutations -W ee − and W ee − Cdc25 described in . In Boolean models one cannot distinguish between reduced activity and no activity. This is why we model W ee − as W ee in both cases. We run the model starting with wild type initial conditions, but with deleted nodes k wee , and in the second case k wee and k 25 . In both mutations, the number of steps is reduced to 12, as compared to 14 in the wild type cell cycle described in the previous section. This suggests that the cell can divide at a smaller size than the wild type, where both mutations are viable. Our results are in accordance with the predictions of the earlier differential equation model .

B. Comparison with an existing Boolean model for the fission yeast cell cycle
It is interesting to compare this model with another Boolean model for the fission yeast cell cycle (Davidich and Bornholdt, 2008) that was built on known biochemical interactions between proteins, only. Both models are quite similar and have the same connections between homologue nodes. Their dynamics matches the wild type signaling sequence during the cell cycle. The difference is that in the current model the nodes Cdc13, preM P F , M P F 1, M P F 2 correspond to only two nodes Cdc2/Cdc13 and Cdc2/Cdc13* in (Davidich and Bornholdt, 2008). So in the model (Davidich and Bornholdt, 2008), the complex Cdc2/Cdc13 can have two levels of activation -medium and high. The intermediate level corresponds to sole activation of the Cdc2/Cdc13 node, whereas a high level of activation is represented by activation of both, Cdc2/Cdc13 and Cdc2/Cdc13* .The last one plays the role of a dephosphorylated Cdc2/Cdc13, that closely corresponds to M P F in . In the current model, the node M P F was separated into two nodes M P F 1 and M P F 2 as well, to distinguish different levels of M P F activation.
The formulas responsible for evolution of proteins are similar in a sense that in both a threshold activation rule is used. In the model (Davidich and Bornholdt, 2008), proteins remain active if the corresponding node was not switched-off by other incoming inhibiting signals. This rule means that if the protein was activated it requires some other signals to change its state. Whereas in the current model one needs to have always a positive incoming signal in order to keep the protein in its active form.

V. DISCUSSION AND CONCLUSION
In this paper, our aim is to make a connection between two successfully used methods -the ODE and the Boolean models for predicting properties of a real biological process. In particular, we show a possible limit transition between the ODE and the Boolean model for the fission yeast cell cycle. For this purpose the known ODE model  of fission yeast has been chosen.
In order to do the transition, first the differential equations are rewritten in a such way that the solutions of equations reach 1 in their maximum values. Then the obtained equations are transformed in a limiting procedure to Boolean functions. Firstly, Michaelis-Menten equations and equations with switch-like behavior can be directly reduced to Boolean functions. Secondly, a set of equations with sigmoidal transfer functions can be replaced with Michaelis-Menten equations without changing the sequence of states through which the system evolves. Thirdly, there are also some cases that cannot be reduced to the two previous ones. It happens when the variable is described by a constantly growing function or a function which has distinctly different levels. In this case we propose in the Boolean model to substitute those variables with two labeling intermediate and high activity of it. Finally, all continuous solutions of equations are mapped into ON/OFF states of Boolean network and the transition between states are described by Boolean functions.
This Boolean model reproduces the results of the initial ODE model . In particular, starting with initial conditions as in , the system evolves through the same sequence of states. The second evidence of similar behavior of the ODE and the Boolean model is the robustness to the initial conditions. The Boolean model has a dominant attractor (67%), attracting most of the trajectories, starting from different initial conditions. The dominant attractor coincides with initial biological conditions of the system. The ability to model mutations in the Boolean approach additionally confirms a good correspondence between the ODE and the Boolean model.
We find that the transition to a Boolean model is possible for differential equations, which have monotonic sigmoidal functions with distinct upper and lower asymptotes. In particular, firstly, in our case Michaelis-Menten equations are reduced to S-shaped GK-functions which have the necessary asymptotes (Glass and Kauffman, 1973;Glass and Hill, 1998). This function works as a switch between the cases when parameters are defined as the upper or lower asymptote and the target control function corresponds to the maximal or basal rate of biochemical processes. Secondly, here substituting some equations that have monotonic sigmoidal functions on the right-hand-side with Michaelis-Menten functions, we also find, that the exact form of the sigmoidal function does not strongly influence the behavior of the system. The comparison of the current model with a previous Boolean model for fission yeast reveals that they both have a similar set of variables (proteins) and similar Boolean functions responsible for update.
Our results also confirm the idea that some molecular control networks are so robustly designed that timing is not a critical factor (Braunewell and Bornholdt, 2006). In our case it is possible to reproduce the main results of  without including time, but reproducing the right sequence of events. It supports the idea that the Boolean approach could contain sufficient information. Thereby one needs less information about the system, the knowledge about reactions on the level of activation/inhibition is sufficient, which eliminates the problem of finding the right kinetic constants. Another advantage is the low computational cost of Boolean networks. The problems one meets working with the Boolean approach are that it is sometimes difficult to reduce the concentration level of some proteins only to ON/OFF states. Sometimes there are intermediate states of concentration which need to be separated from high concentration. In this case two methods are possible. One is, as we implemented it here, to divide this variable into two and to perform as two different nodes in a system. Doing this, one needs to take into account the differences in influences of this protein when it has intermediate and high concentration. Another solution for a such situation could be the introduction of two discrete levels of concentration that the protein can have, for example 1 for intermediate and 2 for high concentration, which was used in modeling gene regulatory network model for Arabidopsisthaliana flower development (Espinosa-Soto et al., 2004).
We would like to note that the ODE and the Boolean approach are both useful methods. The advantage of the ODE approach is that it provides detailed information about the system at any given time in contrast to the Boolean method, which reproduces only the right sequence of events. But the costs for this information are the following. One needs to have exhaustive information about the reactions, where the most difficult part is to find the right kinetic constants. Also it demands more computational costs to find the solution of the system. One could say that the ODE approach is appropriate when the system is well studied and it is necessary to make a detailed study of all reactions that take place. On the other hand if the task is to understand the main principles of some process and one has less information, the Boolean approach is very suitable to use.