Gene Regulatory Networks represent the interactions among genes regulating the activation of speci c cell functionalities and they have been successfully modeled using threshold Boolean networks. In this paper we propose a systematic translation of threshold Boolean networks into Reaction Systems. Our translation produces a non redundant set of rules with the minimal number of objects. This translation allows us to simulate the behavior of a Boolean network simply by executing the (closed) Reaction System we obtain. Moreover, it allows us to investigate the role of di erent genes simply by \playing" with the rules related to di erent genes. Starting from a well-known Boolean network model of the Yeast-Cell Cycle, we construct the corresponding Reaction System and start an investigation on causal dependencies among genes.
Introduction In the context of molecular biology of cells, Gene Regulatory Networks (GRNs) represent the interactions among genes regulating the activation of speci c cell functionalities. More speci cally, genes in a GRN can be either active (i.e. the corresponding protein is expressed) or not, and each active gene can either stimulate or inhibit the activation of a number of other genes. Moreover, the activation of some genes is also usually in uenced by other factors such as the availability of some substances in the environment or the reception of a signal form neighbor cells. As a consequence, gene regulatory networks can be seen as the mechanism that allow a cell to react to external stimuli. When a stimulus is received, it causes a change in the activation state of a few genes, that, in turn, in uence other genes, allowing a new con guration of active genes (corresponding to a new set of active cell functionalities) to be reached.
Several approaches have been proposed to model and analyze GRNs (see [24] for a survey on this topic). Modeling techniques can either deal only with the qualitative aspects of such networks (treating them essentially as logic circuits), or can describe also the quantitative aspects, such as the rates of the interactions. The latter approach is for sure more precise, but requires many additional details of the network dynamics to be taken into account, such as the rates of transcription and translation of genes' DNA into proteins and the rates of protein-protein and protein-DNA interactions. Qualitative models are often su cient to reason on the behavior of the regulatory networks, although with some degree of approximation.
In the qualitative modeling setting, one of the most successful modeling
frameworks for gene regulatory networks are Boolean networks. In this setting,
a particular simple form of Boolean networks, the so called threshold Boolean
networks [
Boolean networks allow dynamical properties of GRNs to be investigated. Starting from an initial con guration of active genes, the dynamics of a GRNs is expressed as a sequence of steps in which such a con guration is updated according to the in uences among the genes described by the Boolean network. Dynamical properties can be investigated either by performing simulations, or by constructing the ( nite) graph representing all possible dynamical evolutions. Example of properties that are often studied on these models are reachability and stability of con gurations, and con uence of evolutions started from di erent initial con gurations into stable con gurations (attractors).
Analysis of dynamical properties may become computationally very
expensive. In order to reduce the state space to be analyzed, minimization techniques
can be applied. The Boolean function represented by a Boolean network can
be synthesized in any framework of logic minimization. The classical approach
to logic minimization that produces sum of products two level formulas can be
used (see e.g., Espresso [
Other analysis methods for GRNs could be applied by changing the
representation of the Boolean networks describing them. In this paper we propose
a systematic translation of threshold Boolean networks into Reaction Systems
[
In this setting the dynamic run of the Reaction System simulates the
evolution of the Boolean network. This correspondence allows us to \play" with the
rules of the Reaction System related to di erent genes in order to detect
dynamic causality dependencies between genes activation/deactivation. Moreover,
we believe that this correspondence will allow us to apply to Boolean networks
well-known techniques to detect causality relationships between objects in
biological systems. The understanding of causality relationships among the events
happening in a biological (or bio-inspired) system is an issue investigated in the
context both of systems biology (see e.g. [
In [
The translation we propose in this paper produces a non redundant set of rules with the minimal number of objects.
The paper is organized as follows. Section 2 introduces the main concepts of (Closed) Reaction Systems. In Section 3 we describe how Boolean networks are de ned and how they work. Section 4 presents the systematic translation from Boolean network to Reaction Systems we propose. Finally, in Section 5 we apply our methodology to simulate and study the Yeast-Cell Cycle Boolean Network. 2
Closed Reaction Systems
In this section we recall the basic de nition of Reaction Systems [
The state of a Reaction System is described by a set of objects. Let a = (Ra; Ia; Pa) be a reaction and T a set of objects. The result resa(T ) of the application of a to T is either Pa, if T separates Ra from Ia (i.e. Ra T and Ia \ T = ;), or the empty set ; otherwise. The application of multiple reactions at the same time occurs without any competition for the used reactants (t hreshold supply assumption). Therefore, each reaction which is not inhibited can be applied, and the result of the application of multiple reactions is cumulative. Formally, given a Reaction System A = (S; A), the result of application of A to a set T S is de ned as resA(T ) = resA(T ) = Sa2A resa(T ).
An important feature of Reaction System is the assumption about the nonpermanency of objects: the objects carried over to the next step are only those produced by reactions. All the other objects vanish, even if they are not involved in any reaction.
The dynamics of a Reaction System is generally driven by the contextual objects, namely the objects supplied to the system by the external environment at each step. Closed Reaction Systems are the subset of general Reaction Systems where the external environment provides objects at the rst step only.
This allows us to simplify the dynamics of a (closed) Reaction System A = (S; A). Indeed, given the initial set D0 the semantics can be simply de ned as the result sequence, = D1; : : : ; Dn where each set Di, for i 1, is obtained from the application of reactions A to the state obtained at the previous step Di 1 ; formally Di = resA(Di 1) for all 1 i < n. For the sake of simplicity, we write Di 1 !A Di as a shorthand for Di = resA(Di 1). In this case the sequence of states of the Reaction System coincides with the result sequence = D1; : : : ; Dn. 3
Boolean Networks We present a formal de nition of threshold Boolean networks [22] considering a set M of n elements, S1; S2; : : : Sn to be nodes of a network. We assign to each element, at each time instant t, a value Si(t) 2 f0; 1g denoting if the element Si is present at that instant or not. The interactions among elements are given by the set of edges of the network called E. Each edge in E can be either activating or inhibiting. An edge from element Sj to element Si is denoted aij (where i 6= j given that a element cannot either activate or inhibit itself). An activating edge has value 1 while an inhibiting one has value 1. Elements M can be partitioned in two sets Msa and Mnsa of self-activating and non-self-activating elements, respectively, M = Msa [ Mnsa. A self-activating elements if is present at time t and it is not inhibited will be present at time t + 1, while a non-self-activating one will not. Moreover we assume that each element Si has associated a value i 2 which is called the threshold parameter of Si. The pair (M; E) is called a threshold Boolean network.
The states of the nodes in the network are updated in parallel in discrete time. The rules for updating the values of nodes are the following, for i 2 f1; : : : ; ng >8>> 1 if Pj aij Sj (t) > i Si(t + 1) = >>>>>>>< 0Si(t) iiff SPji a2ijMSjs(at)^< Pj aiij Sj (t) = i
>>>>: 0 if Si 2 Mnsa ^ Pj aij Sj (t) = i where the value i is the threshold parameter associated to the element Si.
Typically the threshold parameter i associated to Si is equal to 0 so that the switch is inactive if there is no input signal, and it switches on when signals are present. A node which needs more than one incoming signals to be activated can be represented in the model by setting i to a value greater than 0.
Starting from an initial condition, the network produces a dynamical sequence of states and it can reach a periodic attractor or a xed point.
There is a natural representation of Boolean networks by means of graph where the nodes represent the elements and the edges represent the interactions between the elements; an activating edge is indicated by ! while an inhibiting one is indicated by a. Non self-activating genes are represented by nodes with half-arrow (*) self loops.
We introduce an example to illustrate threshold Boolean networks and their dynamic evolution.
Example 1. Let us consider the threshold Boolean network (M; E) with elements M = fA; B; C; Dg such that Msa = fA; B; Cg and Mnsa = fDg and with the edges depicted in Fig. 1. Thus, the element D is non self-activating while the elements A; B and C are self-activating. We also assume that the threshold parameter for each element is 0.
We describe the temporal evolution considering an initial state in which the element D is present while the others are not. We have Initially the element D stimulate the activation of both elements A and B because the element C, their inhibitor is not present. Note that at the second step the element D is inactive because it is non self-activating. Then, at step 3, the elements C is present because it is activated by B while A and B are still present because they are self-activating and at step 2 C was not present. At step 4, the element C actives D and inhibits A which is not present anymore. By contrast, C does not inhibit the element B which is still present because it was present at step 3 but also A was present. The step 5 is similar. Finally at step 6 the element B is not present anymore because in the previous con guration A and D were absent and the inhibitor C was activated. The last one is also a stationary state of the system, since no more evolutions of the network are possible from such state.
Translation of Boolean Networks into Reaction Systems We present a translation of threshold boolean networks in closed Reaction System. Given a Boolean network (M; E) with M = fS1; S2; : : : Sng we de ne for Si 2 M , Act(Si) = fSj j j 2 [1; n] ^ aij = 1g In(Si) = fSj j j 2 [1; n] ^ aij = 1g
We recall that aij denotes an edge from element Sj to element Si. Hence, Act(Si) reports the elements Sj which activates Si and analogously In(Si) reports the elements Sj which inhibits it.
De nition 1. Let (M; E) be a threshold Boolean network with elements M = fS1; S2; : : : Sng and threshold parameters = f 1; 2; : : : ng. We de ne its translation as the closed Reaction System RS((M; E)) = (M; A), where reactions in A are constructed according to the following inference rules: 1)
Pi #Pi
Act(Si) Ii In(Si) #(In(Si) nIi) = i 1 2)
Si 2 Msa #Pi
Pi
Act(Si) Ii #(In(Si)nIi) = i
In(Si) (Pi; Ii; fSig) 2 A (Pi [ fSig; Ii; fSig) 2 A
The closed Reaction System RS((M; E)) simulates the threshold Boolean network (M; E) using reactions obtained by applying either the inference rule 1) or the inference rule 2). Rule 1) de nes the reactions which simulate the production of a element Si at time t + 1 whenever at time t the number of the elements which activate Si minus the number of the elements which inhibit Si is greater than i (according to the rule given in Section 3). This behavior is simulated by a reaction which has as product Si, as reactants Pi and as inhibitors Ii where Pi is a subset of the elements which activates Si and analogously Ii is a corresponding subset of the elements which inhibits it. Note that this reaction can be applied if in the Reaction System state none of the elements in Ii is present. As a consequence, the set of the elements which are inhibitors of Si and which may be present is given by In(Si)nIi. Therefore we guarantee that the cardinality of Pi minus that of In(Si)nIi is greater than i. More speci cally we require that this di erence is equal to i 1 in order to built a minimal set of reactions in the corresponding Reaction System.
Rule 2) is de ned for self-activating elements which remains active at time t + 1 if they are present at time t and they are not inhibited (according to the rule given in Section 3). In Reaction System, due to the non-permanency of objects, the objects carried over to the next step are only those produced by reactions. Therefore, in this case, the reaction which simulates the behavior has Si as reagent and also as product. Similarly as in the case 1), Pi is a subset of the elements which activates Si and Ii is a corresponding subset of the elements which inhibits Si. In this case, however, we require that the cardinality Pi minus the number of inhibitors which might be present (i.e. (In(Si)nIi) ) is exactly i. Example 2. We give the translation of the threshold Boolean network (M; E) presented in Example 1. Fig. 1) illustrates the interactions between the elements of the network M = fA; B; C; Dg such that Msa = fA; B; Cg and Mnsa = fDg.
By assuming again that the threshold parameter for each element is 0 we obtain the closed Reaction System RS((M; E)) = (M; A) with reactions A are de ned as follows: (fCg; ;; fDg); (fBg; ;; fCg) (fCg; ;; fCg); (fA; Dg; ;; fAg) (fDg; fCg; fAg); (fAg; fCg; fBg) (fAg; fCg; fAg); (fB; Dg; ;; fBg) (fDg; fCg; fBg); (fA; Dg; ;; fBg) (fB; Ag; ;; fBg); (fBg; fCg; fBg) The reactions on the left are obtained by applying the rule 1) while those on the right by applying the rule 2). The two columns on the left contains one or more reaction for each elements which produces the element. For the elements D and C there is exactly one reaction given that they do not have inhibitors. By contrast, the element C inhibits both A and B. In the rst case, there is just one reaction which says that D produces A if not inhibited by C. The case of B is similar but it can be activated by two element, A and D. Hence, there are three di erent reactions corresponding to the possible conditions of elements which can activate and inhibit B. Note that the requirements of rule 1) guarantee that only minimal subsets are considered. For instance a reaction such as (fA; Dg; fCg; fBg) is subsumed by (fA; Dg; ;; fBg) given that the latter reaction can be applied regardless of the presence of C.
The two columns on the right shows the reactions for self-activating elements, A, B and C. Similarly as in the previous case there is one or more reaction for each self-activating elements which has the element both as reagent and as product.
We can now prove the soundness of the translation of threshold Boolean networks in closed Reaction System. To relate the con gurations of a threshold Boolean network with the state of the associated Reaction System we introduce the following de nition.
De nition 2. Given a threshold Boolean network (M; E) with M = fS1; : : : Sng, and a state at time t, S(t) = fS1(t); S2(t); : : : Sn(t)g. The translation of the state S(t) in a corresponding Reaction System state is given by RS(S(t)), de ned as follows: RS(S(t)) = fSi j Si(t) = 1; i 2 [1; n]g.
Theorem 1. Let (M; E) be a threshold Boolean network with elements M = fS1; S2; : : : Sng and threshold parameters = f 1; 2; : : : ng. Given a state at time t, S(t) = fS1(t); S2(t); : : : Sn(t)g we have that
RS(S(t + 1)) = resA(RS(S(t))) where A = RS((M; E)) = (M; A) is the Reaction System obtained by the translation.
Due to Theorem 1 a threshold Boolean network (M; E) with a state S(0) at time 0 can be simulated by the corresponding closed Reaction System RS((M; E)) by considering the initial state RS(S(0)).
At this point, it important to count the number reactions of the closed Reaction System necessary to simulate a threshold Boolean network (M; E). Such number depends on the number of the nodes M and of the edges of the network E and also on the threshold parameters . For each Si 2 M the number of the reactions which have Si as a product depends on the cardinality of Act(Si) and In(Si) and on i. Indeed, Act(Si) and In(Si) represents the number of the incoming edges which activates and inhibits Si respectively.
Proposition 1. Given a threshold Boolean network (M; E) with elements M = fS1; S2; : : : Sng and threshold parameters = f 1; 2; : : : ng.
Moreover, let RS((M; E)) = (S; A) be the corresponding closed Reaction System. { For each i 2 f1; : : : ; ng, let N (Si) = #(f(R; I; P ) 2 A j Si 2 P g) be the number of the reactions which produce the element Si. We have that N (Si) = 8 Pmin(mi;(li+1+ i)) mi >>> k=1+ i k > > < Pmin(mi;(lj+1+ i)) mi
k=1+ i k > > > >:> Pmin(mi;(lj+ i)) mi h= i h k 1li i ; if Si 2 Mnsa;
li k 1 i + li ; k i if Si 2 Msa: where mi = #(Act(Si)) and li = #(IN (Si)). { We have that #(A) = Pn
i=1 N (Si).
The previous result is a direct consequence of De nition 1.
Example 3. Let us consider the closed Reaction System RS(M; E) = (M; A), presented in the Example 2, which is the translation of the threshold Boolean network (M; E) of Example 1. The reaction system has 12 reactions. Indeed, since A; B and C belong to Msa while D belongs to Mnsa and the threshold parameter is 0 by applying Proposition 1, we obtain:
N (A) = Pmin(1;2) 1
i=1 i N (B) = Pmin(2;2) 2
i=1 i N (C) = Pmin(1;1) 1
i=1 i N (D) = Pmin(1;2) 1 i=1 i i 11 + Pmin(1;1) 1
i=0 i i 11 + Pmin(2;1) 2
i=0 i i 01 + Pmin(1;0) 1
i=0 i 0 i 1 1 = 1 + (1 + 1) i 1 = (2 + 1) + (1 + 2) = 6 i 0 = 1 + 1 i = 3 = 2 = 1 5
Simulating the Yeast-Cell Cycle Boolean Network The cell-cycle process by which a cell goes and divides into two cells is a vital process the regulation of which is conserved among the eukaryotes [21]. The process mainly consists in four phases depicted in Figure 2.
In phase G1 the cell grows and, under appropriate conditions, commits to division, in phase S the DNA is synthesized and chromosomes replicated, G2 is the phase where the cell checks the duplicated chromosomes, and nally in the M (Mitosis) phase the cell is divided into two. After the Mitosis phase, the cell enters the G1 phase, hence completing a cycle. There are about 800 genes involved in the cell-cycle process of the budding yeast [25]. However, the number of key regulators that are responsible for the control and regulation of this complex process is much smaller. Based on extensive literature studies, the authors in [20] constructed a network of key regulators involving 11 genes. The relations between genes are described by the boolean network (MCell; ECell) depicted in Figure 3, where the threshold parameter is always 0. The boolean network was used to study the time evolution of the protein states. Starting from the 211 = 2048 possible initial states describing a con guration for gene activation, they discover that all of them ow into one of seven attractor stationary states. In particular, among the seven xed points there is one big attractor that attracts 1764 initial states. We translated the Boolean network (MCell; ECell) of Figure 3 into a Reaction System Acell = RS(MCell; ECell) = ((MCell; ACell) by applying the procedure described in De nition 1. The Reaction System Acell has 52 reactions. For simplicity, we just show the translation of reactions describing the production (activation) of a single node of the Boolean network. Consider the central node named Sic1 in the Boolean network of Figure 3. It has 2 activating incoming arcs and 3 inhibiting incoming arcs. Since Sic1 2 Msa, by Proposition 1 there will be Pim=i1n(2;4) 2i i 31 +Pim=i0n(2;3) 2i 3i = (2+3)+(1+6+3) = 15. Indeed, by applying our translation we nd the following rules for the production of Sic1: (fCdc20g; fClb5; 6; Clb1; 2; Cln1; 2g; fSic1g) (fSwi5gfClb5; 6; Clb1; 2; Cln1; 2g; fSic1g) (fCdc20; Swi5g; fClb1; 2; Cln5; 6g; fSic1g) (fCdc20; Swi5g; fClb1; 2; Cln1; 2g; fSic1g) (fCdc20; Swi5g; fCln1; 2; Clb5; 6g; fSic1g) (fSic1g; fClb5; 6; Clb1; 2; Cln1; 2g; fSic1g) (fSic1; Cdc20g; fClb5; 6; Clb1; 2g; fSic1g) (fSic1; Cdc20g; fClb5; 6; Cln1; 2g; fSic1g) (fSic1; Cdc20g; fClb5; 6; Clb1; 2g; fSic1g) (fSic1; Swi5g; fClb5; 6; Clb1; 2g; fSic1g) (fSic1; Swi5g; fClb5; 6; Cln1; 2g; fSic1g) (fSic1; Swi5g; fClb5; 6; Clb1; 2g; fSic1g) (fSic1; Cdc20; Swi5g; fCln1; 2g; fSic1g) (fSic1; Cdc20; Swi5g; fClb1; 2g; fSic1g) (fSic1; Cdc20; Swi5g; fClb5; 6g; fSic1g) Note that the behavior of the reactions producing Sic1 faithfully model the activation of gene Sic1 in the Boolean network. Consider the case where genes Cdc20; Swi5 and Clb5; 6 are all active according to the Boolean network of Figure 3 after one step Sic1 becomes active. The previous state is represented in the Reaction System as the set of activated genes D0 = fCdc20; Swi5; Clb5; 6g. Now starting from D0, we can apply rule (fCdc20; Swi5g; fClb1; 2; Cln1; 2g; fSic1g) to obtain the production(activation) of gene Sic1, therefore Sic1 2 D1 where D0 !ACell D1. Note that if Cln1; 2 was also active in the initial state then gene Sic1 could not be activated according to the Boolean network. This is modeled in the Reaction System by the fact that none of the 15 rules producing Sic1 could be applied to the set D0 = fCdc20; Swi5; Clb5; 6; Cln1; 2g. Once we have obtained the complete Reaction System Acell that simulates the entire Boolean network we can run it with any initial state D0 in order to study the behaviour of the cell when some genes are activated. As a rst experiment we executed the Reaction System with the initial state D0 that it was observed in nature triggers the cell-cycle. Indeed, usually the cell stays in a stationary state where just genes Sic1 and Cdh1 are active. When the cell grows, the external cell size signal Cell size arrives and activates Cln3. This "excites" the cell from its stationary state and triggers the cycle. We can observe the di erent states of activation/deactivation of genes during the cell cycle by executing the Reaction System with an initial state where Sic1, Cdh1 and Cln3 are active. Thus, the following evolution can be observed. At this point the evolution reaches the stationary state fSic1; Cdh1g and the cell waits for another external stimulus to arrive, that is an external new cell size signal that activates gene Cln3 and triggers a new cycle. The evolution of the Reaction System represents the evolution of the Boolean network depicted in Figure 4 that describes the entire cell cycle. The Reaction System ACell can now be used for studying the in uence that each gene has in the cell cycle. Each gene can be silenced in turn simply by deleting the rules that produces such gene. Note that this corresponds to simulate the Boolean network where we canceled the node representing the gene together with all his arcs. As a second example consider the case where gene SBF was silenced. To this aim, let the Reaction System ACell = (MCell; ACell); we consider the Reaction System ACell SBF = (MCell; ACell=f(Ra; Pa; fSBFg)j a 2 ACellg). In this case, starting from the stationary state fSic1; Cdh1; Cln3g the following evolution can be observed.
fSic1; Cdh1; Cln3g !ACell SBF fCdh1; MBF; Sic1g !ACell SBF fCdh1; MBF; Sic1g The corresponding evolution on the Boolean network is depicted in Figure 5.
It can be observed that if gene SBF was silenced the cell could not perform the cycle because after one step it reaches a new stationary state. This shows that the gene SBF was, in some way, necessary for the cycle to be performed. As a last example consider the case where gene Mcm1 was silenced. The evolution we obtain by considering the Reaction System ACell Mcm1 = (MCell; ACell=f(Ra; Pa; fMcm1g)j a 2 ACellg) starting with the stationary state fSic1; Cdh1; Cln3g is directly depicted in Figure 6.
In this case we obtain a very di erent result from the previous one. Indeed, it can be observed that even if gene Mcm1 was silenced the cell could perform most of its cycle and go back to the initial stationary state. This suggests that the gene Mcm1 was not necessary for the cycle to be performed. Indeed, the cell can recover even if, for some reasons, gene Mcm1 could not be activated. 6
Conclusions In this paper we proposed a systematic translation of boolean networks into reaction systems. This allows us to simulate the behaviour of a boolean network simply by executing the reaction system we obtain. We applied our method to model the extensively studied Cell cycle that allows a cell to split. The translation of the boolean network describing the Cell cycle into a reaction system allows us to investigate the role of di erent genes simply by "playing" with the rules related to the production of such gene. In this way, we studied the e ects of the silencing of some genes such SBF and Mcm1 on the entire cell cycle.
However we think that the main advantage of our translation will be the
application of well known techniques for detecting dynamic causalities relations
in reaction systems (see [