A circuit-preserving mapping from multilevel to Boolean dynamics
AA CIRCUIT-PRESERVING MAPPING FROM MULTILEVEL TO BOOLEANDYNAMICS
ADRIEN FAUR ´E AND SHIZUO KAJIA bstract . Many discrete models of biological networks rely exclusively on Boolean variables andmany tools and theorems are available for analysis of strictly Boolean models. However, multilevelvariables are often required to account for threshold e ff ects, in which knowledge of the Boolean casedoes not generalise straightforwardly. This motivated the development of conversion methods for mul-tilevel to Boolean models. In particular, Van Ham’s method has been shown to yield a one-to-one,neighbour and regulation preserving dynamics, making it the de facto standard approach to the prob-lem. However, Van Ham’s method has several drawbacks: most notably, it introduces vast regions of“non-admissible” states that have no counterpart in the multilevel, original model. This raises specialdi ffi culties for the analysis of interaction between variables and circuit functionality, which is believedto be central to the understanding of dynamic properties of logical models. Here, we propose a newmultilevel to Boolean conversion method, with software implementation. Contrary to Van Ham’s, ourmethod doesn’t yield a one-to-one transposition of multilevel trajectories; however, it maps each andevery Boolean state to a specific multilevel state, thus getting rid of the non-admissible regions and, atthe expense of (apparently) more complicated, “parallel” trajectories. One of the prominent features ofour method is that it preserves dynamics and interaction of variables in a certain manner. As a demon-stration of the usability of our method, we apply it to construct a new Boolean counter-example to thewell-known conjecture that a local negative circuit is necessary to generate sustained oscillations. Thisresult illustrates the general relevance of our method for the study of multilevel logical models.
1. B ackground
Boolean models have proved very useful in the analysis of various networks in biology. However, itis often convenient to introduce multilevel variables to account for multiple threshold e ff ects. We areoften faced with choices between using Boolean variables or multilevel variables. This can be crucialsince theoretical results are sometimes proved only for Boolean or multilevel networks. A particularexample of this situation is in Ren´e Thomas’ conjecture that a local negative circuit is necessary toproduce sustained (asynchronous) oscillations. This paper stems from the simple idea that a Booleancounter-example to that conjecture could be found by transposing a multilevel counter-example foundearlier by Richard and Comet. However, we believe the method developed in this paper, together witha handy script which implements it, is widely applicable to other theoretical studies which involvesdiscrete networks. We also find the notion of asymptotic evolution function defined in this paper shedslight on the understanding of relation between the state transition graph and the interaction graph.1.1. Introduction.
Introduced in the 1960s-70s to model biological regulatory networks , the log-ical (discrete) formalism has gained increasing popularity, with recent applications as diverse asdrosophila development, cell cycle control, or immunology (see [1] for a survey). While many ofthese models rely exclusively on Boolean variables, it is often useful to introduce multilevel variablesto account for more refined behaviour. However, many tools and theoretical results are restricted tothe Boolean case (see e.g. [2, 3, 4]) This situation motivated the development of methods to convertmultilevel models to Boolean ones [5, 6]. A simple idea for such a conversion was introduced byVan Ham [6], and this method has been shown to be essentially the only one that could provide a“one-to-one, neighbour and regulation preserving map” [7]. One problem with the conversion is that
Mathematics Subject Classification.
Primary 68R05; Secondary 92D99.
Key words and phrases.
Discrete dynamical system, Automata network, Boolean network, Regulatory network, Ge-netic regulation, Thomas’ conjecture.The authors contributed equally to this work. a r X i v : . [ q - b i o . M N ] M a y ADRIEN FAUR ´E AND SHIZUO KAJI the resulting Boolean model is defined only on a sub-region of the whole Boolean state space, calledthe admissible region , and how to extend the model outside that region is not trivial. This leads topotential problems with analytical tools designed to deal with the whole state space, as a property thatis true in the restricted domain may be false on the whole state space, and vice versa. The primarygoal of the present paper is to address this issue by introducing an extension of Van Ham’s method.More precisely, we introduce a new method for multilevel to Boolean model conversion which ex-tends the domain of Van Ham’s model to the whole state space while preserving edge functionalityand, therefore, local circuits. Our mapping yields a state transition graph with “parallel” trajectoriesthat contains the one obtained by Van Ham’s mapping as a sub-graph in such a way that attractors ofthe dynamics are preserved.We apply our method to investigate a particular class of theoretical results that connect the asyn-chronous behaviour of a model to the presence of regulatory circuits in the interaction graph. In theearly 1980s, R. Thomas conjectured that the presence of a positive circuit ( i.e. a circuit where eachcomponent directly or indirectly has a positive e ff ect upon itself) in the interaction graph is a necessarycondition for multi-stability, and a negative circuit (where each component has a negative e ff ect uponitself) is necessary for sustained oscillations [8]. One particular formulation of the conjecture focuseson local or “type-1” circuits [9], i.e. circuits whose arcs are all functional in the same single pointof the system’s state transition graph – as opposed to global circuits whose arcs may be functionalanywhere. While the conjecture holds for positive circuits both at the global and local levels, andfor multilevel as well as Boolean models [10, 11], in the negative case the conjecture could only beproved true at the global level [10]. At the local level, a counter-example has first been published formultilevel models [12], while the Boolean case remained open [9] until a Boolean counter-examplewas eventually discovered [13], showing that contrary to expectations, a local negative circuit wasnot necessary to generate sustained oscillations. Interestingly, the approaches taken by P. Ruet and A.Richard are rather di ff erent, and their counter-examples have little in common. Applying our methodto the Richard-Comet multilevel counter-example, we obtain a new Boolean counter-example to theconjecture that a local negative circuit is necessary to produce sustained oscillations.1.2. Definitions.
Evolution function and State transition graph.
We work within the generalised logical frame-work introduced by Ren´e Thomas and collaborators [14]; see Abou-Jaoud´e et al. [1] for a recentreview. Here, we introduce the notation we use throughout this paper. Fix positive integers n and m i (1 ≤ i ≤ n ). Consider a system consisting of mutually interacting n genes, indexed by the set I = { , , . . . , n } . Each gene a i takes expression levels in the integer interval { , , . . . , m i } . The stateof the system evolves depending on the current state. This leads to a discrete dynamical systemrepresented by a evolution function over Mf = ( f , f , . . . , f n ) : M → M , where M = { ( x , . . . , x n ) | x i ∈ { , , . . . , m i }} . As a special case when m i = i ∈ I , we denote M = B n with B = { , } and call the system Boolean . A basic question asks what we can tell aboutthe asymptotic global behaviour of the dynamics, which is encoded in the state transition graph , fromlocal data of f , which are encoded in the partial derivatives of f or the interaction graph .The evolution of the whole system can be formally modelled by a certain kind of directed graphon M . We equip M with the usual metric d ( x , x (cid:48) ) = (cid:80) ni = | x i − x (cid:48) i | for x , x (cid:48) ∈ M . Denote by e = (1 , , , . . . ) , e = (0 , , , , . . . ) , etc. the coordinate vectors of M . A grid graph Γ over M is a graphwith the vertex set M satisfying that • each directed edge connects a pair of vertices of distance one • at each vertex x there are no two parallel outward edges; that is, x − e j ← x → x + e j is notallowed.The state of the whole system is represented by the levels of genes, and corresponds to a vertex in Γ .At each time step, the state evolves to one of its neighbouring vertices connected by an arrow in the CIRCUIT-PRESERVING MAPPING FROM MULTILEVEL TO BOOLEAN DYNAMICS 3 following way. To an evolution function over M , we associate a grid graph Γ ( f ) over M called the (asynchronous) state transition graph with the edge set(1) ( x , x , . . . , x j , . . . , x n ) → ( x , x , . . . , x j + δ, . . . , x n ) , δ = − f j ( x ) < x j ) + f j ( x ) > x j ) . Note that here we follow the standard convention that transition of states is unitary (see [12, § x → x (cid:48) implies d ( x , x (cid:48) ) =
1; that is, at each step the level of a single genechanges at most by one.Asymptotic behaviour of the evolution of a system can be captured in a graph theoretical entityof the state transition graph. An attractor is a terminal strongly connected sub-graph of Γ ; that is,any two elements of it are connected by a path and there is no edge from its elements to one in thecomplement. An attractor consisting of a single vertex is called a stable state , otherwise it is called a cyclic attractor . Intuitively, attractors are domains in Γ in which the system eventually resides; thereis no way to escape once the system arrives in it, but each state in the domain can be visited afterarbitrarily many steps.1.2.2. Interaction graph and circuit functionality.
A common practice in analysing interactions amonggenes in a network is to encode it in the form of a labelled directed graph called the interaction graph,where interaction is measured by the partial derivatives of the evolution function f = ( f , f , . . . , f n ) : M → M .The forward partial derivative of f i along the j -th coordinate at x = ( x , . . . , x n ) with x j < m j isdefined by ∂ + j f i ( x ) = f i ( x , . . . , x j + , . . . , x n ) − f i ( x , . . . , x j , . . . , x n ) = f i ( x + e j ) − f i ( x ) . The backward partial derivative along the j -th coordinate at x with x j > ∂ − j f i ( x ) = f j ( x , . . . , x j , . . . , x n ) − f j ( x , . . . , x j − , . . . , x n ) = f i ( x ) − f i ( x − e j ) . Partial derivatives ∂ + j f i ( x ) and ∂ − j f i ( x ) are non-trivial when the i -th gene’s target value changes alongthe change of the j -th gene. They encode the dependence between genes locally at the state x ∈ M . Remark 1.
For a Boolean network, only one of the forward or the backward partial derivative existsat each x , so we just put them together to define the ordinary partial derivative denoted by ∂ j . On theother hand, in multilevel case, we have both the forward and the backward partial derivatives at some x . It is important to consider both of them (c.f. [12, Definition 8]). Definition 1.
The (local) interaction graph G f ( x ) of f at x is a graph over the vertex set I such thatthere exists an edge from j to i • with label “ + ” if ∂ + j f i ( x ) > ∂ − j f i ( x ) > • with label “ − “ if ∂ + j f i ( x ) < ∂ − j f i ( x ) < j to i at the same time. We define the global interaction graph G f ( M ) as the union of edges of G f ( x ) for all x ∈ M . Definition 2 ([9, Definition 1, Proposition 1]) . • A cycle C in G f ( M ) is called a positive (resp. negative) circuit if it contains an even (resp.odd) number of negative edges. • A circuit C is said to be type-1 functional if C ⊂ G f ( x ) for some x .As in the continuous case, a function f is recovered up to constant by its partial derivatives: For twoevolution functions f , g : M → M satisfying ∂ + j f i = ∂ + j g i for all i , j ∈ I (or ∂ − j f i = ∂ − j g i for all i , j ∈ I ),the di ff erence f i ( x ) − g i ( x ) is constant for any i ∈ I . In particular, two distinct Boolean evolutionfunctions have the same partial derivatives if and only if they are constant and do not coincide at anypoint. This means the partial derivatives have almost all the information of the network However, thenext example shows that the partial derivatives are not enough to determine the asymptotic behaviourof the dynamics. ADRIEN FAUR ´E AND SHIZUO KAJI
Example 1.
Consider the the Boolean evolution functions defined by f ( x , x ) = (0 ,
0) (( x , x ) = (1 , ,
0) (otherwise) g ( x , x ) = (0 ,
1) (( x , x ) = (1 , ,
1) (otherwise) . Since they di ff er by a constant, their partial derivatives agree. There exists the unique cyclic attractor(0 , ↔ (1 ,
0) in Γ ( f ), whereas there exists the unique stable state (1 ,
1) in Γ ( g ).2. M ethods Asymptotic evolution function.
The correspondence between evolution functions and statetransition graphs is not bijective. In fact, as discussed by Streck et al. [15], for a given (multilevel)grid graph Γ , there are multiple evolution functions which have Γ as their state transition graphs. Tohave a bijective correspondence between the two representations of the system, we restrict ourselvesto a certain class of evolution functions. There are two major conventions:(1) We say f is stepwise or unitary if | f i ( x , . . . , x n ) − x i | ≤ i ∈ I and x ∈ M .(2) We say f is asymptotic if f i ( x , . . . , x i , . . . , x n ) ∈ { , x i , m i } for all i ∈ I and x ∈ M .In both cases, f i encodes only the sign of f i ( x ) − x i .For any evolution function, there exists a unique asymptotic and a unique stepwise evolution func-tions having the same state transition graph. Proposition 1.
For any evolution function f , define¯ f i ( x ) = m i ( f i ( x ) > x i ) x i ( f i ( x ) = x i )0 ( f i ( x ) < x i ) ˆ f i ( x ) = x i + f i ( x ) > x i ) x i ( f i ( x ) = x i ) x i − f i ( x ) < x i ) . Then, ¯ f is asymptotic and ˆ f is stepwise with Γ ( f ) = Γ ( ¯ f ) = Γ ( ˆ f ).Similarly, for any grid graph Γ , there exists a unique asymptotic function ¯ f Γ and a unique stepwisefunction ˆ f Γ such that Γ = Γ ( ¯ f Γ ) = Γ ( ˆ f Γ ). Proof.
We see how to define ¯ f Γ . For a vertex x ∈ Γ and i ∈ I , we have only one of the threepossibilities: • x → x − e i • x → x + e i • there is no edge from x in the direction of e i .We define an asymptotic evolution function ¯ f Γ by setting f Γ i ( x ) = , m , x i accordingly. (cid:3) If we are interested in the evolution of a system, which is encoded in the state transition graph,we can restrict ourselves to either the class of stepwise evolution functions or the class of asymptoticevolution functions. Our choice in this paper is to restrict ourselves to the latter, and we identify anasymptotic evolution function with its state transition graph and vice versa. In the rest of the paper,we assume functions are asymptotic unless otherwise stated and denoted just by f without a bar overit. There is a little di ff erence in the interaction graph when we consider the stepwise case instead ofthe asymptotic case. When i (cid:44) j , ∂ + i f j ( x ) and ∂ + i ˆ f j ( x ) have the same sign, and same is true for thebackward partial derivatives. However, when i = j , we have the following di ff erence. Example 2.
Consider the asymptotic evolution function f (0) = , f (1) = , f (2) = M = { , , } . The corresponding stepwise evolution function is ˆ f (0) = , ˆ f (1) = , ˆ f (2) =
2. At x = f while there is a positive self-loop in the one of ˆ f . Onthe other hand, consider the asymptotic evolution function f (0) = , f (1) = , f (2) =
0. At x = f , whereas there is no arrow in the one of thecorresponding stepwise evolution function ˆ f (0) = , ˆ f (1) = , ˆ f (2) = CIRCUIT-PRESERVING MAPPING FROM MULTILEVEL TO BOOLEAN DYNAMICS 5
In short, the interaction graphs of an asymptotic function and its corresponding stepwise functionare the same only up to self-loops.A function which is neither asymptotic nor stepwise has in general more non-trivial partial deriva-tives than the asymptotic and the stepwise functions sharing the same state transition graph given inProposition 1. (See [15] for a detailed discussion. The stepwise function in our paper can be seen asa special case of the canonical function defined there.)
Example 3.
An asymptotic function f : { , , } → { , , } defined by ( f , f )( x , x ) = (2 ,
2) havethe same state transition graph with( g , g )( x , x ) = (1 ,
2) ( x , x ) = (0 , ,
2) (otherwise) . However, ∂ f (0 , = ∂ g (0 , =
1. This means, there is a positive arrow x → x in Gg (0 , G f (0 , A mapping from multilevel to Boolean networks.
Fix the set of states M and a natural number l . We consider mappings from the set Evo( M ) of asymptotic evolution functions on M to the setEvo( B l ) of l -dimensional Boolean evolution functions. Mappings between grid graphs are obtainedfrom them by the correspondence given in Proposition 1. Following [7], we introduce two preferableproperties of such mappings. Definition 3.
A mapping Ψ : Evo( M ) → Evo( B l ) is said to be • neighbour-preserving if there exists a map b : M → B l and ψ : B l → M such that ψ ◦ b is the identity on M and b and ψ induce graph homomorphisms ˜ b : Γ ( f ) → Γ ( Ψ ( f )) and˜ ψ : Γ ( Ψ ( f )) → Γ ( f ) for any f ∈ Evo( M ). • globally regulation-preserving if G Ψ ( f )( B l ) (cid:44) G Ψ ( f (cid:48) )( B l ) for any f , f (cid:48) ∈ Evo( M ) with G f ( M ) (cid:44) G f (cid:48) ( M ). • locally regulation-preserving if there exists a map b : M → B l such that G Ψ ( f )( b ( y )) (cid:44) G Ψ ( f (cid:48) )( b ( y )) for any y ∈ M and any f , f (cid:48) ∈ Evo( M ) with G f ( y ) (cid:44) G f (cid:48) ( y ).These properties are practically useful. For a neighbour-preserving mapping, the two maps b and ψ give correspondence between the multilevel states and the Boolean states in such a way that the statetransition graph of any multilevel model is embedded in that of a Boolean model. With a regulation-preserving mapping, one can recover the interaction graph of a multilevel network from the corre-sponding Boolean one.A naive idea to convert an evolution function f : M → M to a Boolean one is to use an embedding(one-to-one map) b : M → B l of the set of multilevel states to a higher dimensional set of Booleanstates. Then, the conjugate of f with respect to b is defined as f b ( x ) : = b ◦ f ◦ b − ( x ) , which is defined only on Im( b ) ⊂ B l , the image of b . The domain Im( b ) is called the admissible region for f b . The state transition graph Γ ( f b ) in this case is defined to be the full sub-graph on Im( b ) of theone defined by Eq. (1). Van Ham [6] proposed one particular embedding b : M → B m + m + ··· + m n which is defined as the direct product of(2) ( b ) i : { , , . . . , m i } → B m i , ( b ) i ( k ) = (1 , , . . . , (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) k , , . . . ,
0) for all i ∈ I . Didier et al. [7] showed that Van Ham’s embedding is essentially the only one satisfying niceproperties which they call neighbour preservation and regulation preservation (see Remark 2 below).On the other hand, an apparent inconvenience of this method is that the resulting evolution function isdefined only on the restricted domain Im( b ), the set of admissible states [7]. In contrast, we will givea construction which produces a Boolean network defined on the whole state space B m + ··· + m n . Theidea is to use a surjective map ψ : B m + ··· + m n → M rather than an embedding in the opposite direction. ADRIEN FAUR ´E AND SHIZUO KAJI
Remark 2.
Properties similar to the first two in Definition 3 were introduced by Didier et al. [7]but only for embeddings b : M → B l (and mappings obtained by conjugation with embeddingsEq. (2)). Recall that an embedding b : M → B l is said to be neighbour-preserving if it satisfies d ( b ( y ) , b ( y (cid:48) )) = y , y (cid:48) ∈ M with d ( y , y (cid:48) ) =
1. Also an embedding is said to be regulation-preserving if
G f b ( B l ) (cid:44) G f (cid:48) b ( B l ) when G f ( M ) (cid:44) G f (cid:48) ( M ) for any f , f (cid:48) : M → M ; that is, the globalinteraction graphs of the Boolean networks obtained by conjugation di ff er when so do those of themultilevel networks. Our definitions are modified versions of theirs which apply to any mapping.Here, we define a mapping from Evo( M ) to Evo( B l ) with l = ( m + · · · + m n ), and a mapping fromgrid graphs over M to grid graphs over B l , which possesses all three above properties.Define a map ψ : B l → M by ψ ( x , , x , , . . . , x , m , x , , . . . , x n , m n ) = ( | y | , . . . , | y n | ) , where y i = ( x i , , . . . , x i , m i ) ∈ B m i and | y i | = (cid:80) m i k = x i , k . We denote the index set of B l = B m + ··· + m n by I B = { ( i , j i ) | ≤ i ≤ n , ≤ j i ≤ m i } . Definition 4.
For an asymptotic multilevel evolution function f ∈ Evo( M ), its binarisation B ( f ) ∈ Evo( B l ) is defined by B ( f ) i , j ( x ) : = f i ( ψ ( x )) < ψ ( x ) i ) x i , j ( f i ( ψ ( x )) = ψ ( x ) i )1 ( f i ( ψ ( x )) > ψ ( x ) i ) . Conversely, we have(3) f i ( ψ ( x )) = m i (cid:88) j = B ( f ) i , j ( x ) . The binarisation of a grid graph Γ on M , denoted by B ( Γ ), is defined to be the grid graph on B l suchthat there exists a directed edge x → x (cid:48) in B ( Γ ) if and only if d ( x , x (cid:48) ) = ψ ( x ) → ψ ( x (cid:48) ) in Γ .It is trivial to see B ( f Γ ) = f B ( Γ ) and B ( Γ ( f )) = Γ ( B ( f )). We now identify the image of binarisationEvo( M ) → Evo( B l ). The symmetric group S m i acts on B m i by permuting the coordinates. We considerthe coordinate-wise action of S = S m × S m × · · · × S m n on B l = B m × B m × · · · B m n . Since the map ψ is invariant under this action, the binarisation B ( f ) has symmetry with respect to this action. Definition 5.
A Boolean network f (cid:48) : B l → B l is said to be S -symmetric if f (cid:48) i , j ( x , , , . . . , x n , m ) = f (cid:48) i ,σ i ( j ) ( x ,σ (1) , x ,σ (2) , . . . , x ,σ ( m ) , x ,σ (1) , . . . , x n ,σ n ( m n ) )for any σ = ( σ , . . . , σ n ) ∈ S . Similarly, a Boolean grid graph Γ (cid:48) is said to be S -symmetric when anedge x → x (cid:48) exists if and only if so does σ ( x ) → σ ( x (cid:48) ).For an S -symmetric Boolean evolution function f (cid:48) we obtain a well-defined evolution function ψ ◦ f (cid:48) ◦ ψ − . Similarly, for an S -symmetric Γ (cid:48) , we obtain a grid graph over M as the image under ψ . Theorem 1.
The binarisation induces a bijective mapping between the set Evo( M ) of asymptoticevolution functions on M and the set Evo( B m + ··· + m n ) of Boolean evolution functions on B m + ··· + m n which are S -symmetric.It immediately follows that Van Ham’s embedding b : M → B l and our own map ψ : B l → M induce graph homomorphisms (cid:101) b : Γ ( f ) → Γ ( B ( f )) and ˜ ψ : Γ ( B ( f )) → Γ ( f ) for any f ∈ Evo( M )such that ˜ ψ ◦ (cid:101) b is the identity. Thus, the binarisation is neighbour-preserving. We will see that itis also locally, and hence globally as well, regulation-preserving. We show that the dynamics of thesystem, namely, attractors in the state transition graph are preserved under binarisation.In what follows, we often make use of the following two obvious facts: Lemma 1.
CIRCUIT-PRESERVING MAPPING FROM MULTILEVEL TO BOOLEAN DYNAMICS 7 • When there exists an edge x → x (cid:48) in B ( Γ ), there exists an edge ψ ( x ) → ψ ( x (cid:48) ) in Γ . • When there exists an edge y → y (cid:48) in Γ , for any x ∈ ψ − ( y ) there exists an edge x → x (cid:48) in B ( Γ )for some x (cid:48) ∈ ψ − ( y (cid:48) ). Proposition 2.
The strongly connected components of B ( Γ ) map surjectively onto those of Γ via ˜ ψ .Moreover, attractors of B ( Γ ) map surjectively onto those of Γ via ˜ ψ . Proof.
Assume that x , x (cid:48) ∈ B ( Γ ) are in the same strongly connected component. This means, there isa cycle containing x , x (cid:48) and it maps to a cycle containing ψ ( x ) , ψ ( x (cid:48) ) ∈ Γ . Therefore, ψ ( x ) , ψ ( x (cid:48) ) arein the same strongly connected component. Conversely, assume that there exists a cycle containing y , y (cid:48) ∈ Γ . For any vertex x ∈ ψ − ( y ), there exists a vertex x (cid:48) ∈ ψ − ( y (cid:48) ) and a cycle containing both x and x (cid:48) . To sum up, the image of a strongly connected component of B ( Γ ) is a strongly connectedcomponent of Γ , and for any strongly connected component of Γ there exists a strongly connectedcomponent of B ( Γ ) which maps to it. Since ˜ ψ is a surjective graph homomorphism, attractors of B ( Γ )map surjectively onto those of Γ via ˜ ψ . (cid:3) Corollary 1. • A stable state exists in B ( Γ ) if and only if it does in Γ . • A cyclic attractor exists in B ( Γ ) if and only if it does in Γ . Proof.
The first statement is trivial. For any cyclic attractor in Γ , there exists an attractor in B ( Γ )which maps to it by the previous proposition. Since it contains more than one element, it is a cyclicattractor in B ( Γ ). Conversely, assume that there is a cyclic attractor in B ( Γ ). It contains at least twoelements x , x (cid:48) with d ( x , x (cid:48) ) =
1. Their images ψ ( x ) , ψ ( x (cid:48) ) should be di ff erent since any two distinctelements in a single fibre (the inverse image of a point) ψ − ( ψ ( x )) have at least distance two. Thus, theimage of the cyclic attractor contains at least two distinct elements ψ ( x ) , ψ ( x (cid:48) ) and is a cyclic attractorin Γ . (cid:3) Theorem 2.
The map I B → I defined by ( i , j i ) (cid:55)→ i induces a surjective graph homomorphism on G B ( f )( x ) → G f ( y ) for y = ψ ( x ) and any x ∈ B l . More precisely, the following two statements hold.(1) At any y ∈ M , if a positive (resp. negative) edge i → i (cid:48) exists in G f ( y ), so does a positive(resp. negative) edge ( i , j ) → ( i (cid:48) , j (cid:48) ) in G B ( f )( x ) for some j and j (cid:48) at any x ∈ ψ − ( y ).(2) At any x ∈ B l , if a positive (resp. negative) edge ( i , j ) → ( i (cid:48) , j (cid:48) ) exists in G B ( f )( x ), so does apositive (resp. negative) edge i → i (cid:48) in G f ( ψ ( x )). Proof.
We only show the statements for the case of a positive edge, as the case of a negativeedge follows by a similar argument.For the first statement, assume that there exists a positive edge i → i (cid:48) in G f ( y ). We have twocases ∂ + i f i (cid:48) ( y ) > ∂ − i f i (cid:48) ( y ) >
0. When ∂ − i f i (cid:48) ( y ) = f i (cid:48) ( y ) − f i (cid:48) ( y − e i ) >
0, for any x ∈ ψ − ( y )there exists j such that x i , j = y i >
0. By Eq. (3) there must exist some j (cid:48) such that B ( f ) i (cid:48) , j (cid:48) ( x ) − B ( f ) i (cid:48) , j (cid:48) ( x − e i , j ) =
1. This means there exists a positive edge ( i , j ) → ( i (cid:48) , j (cid:48) ) in G B ( f )( x ). A similar argument applies when ∂ + i f i (cid:48) ( y ) > i , j ) → ( i (cid:48) , j (cid:48) ) in G B ( f )( x ).When x i , j =
0, this means B ( f ) i (cid:48) , j (cid:48) ( x ) = B ( f ) i (cid:48) , j (cid:48) ( x + e i , j ) =
1. Since ψ ( x ) + e i = ψ ( x + e i , j ), we have f i (cid:48) ( ψ ( x ) + e i ) − f i (cid:48) ( ψ ( x )) = (cid:80) m i (cid:48) k = (cid:16) B ( f ) i (cid:48) , k ( x + e i , j ) − B ( f ) i (cid:48) , k ( x ) (cid:17) . Since B ( f ) i (cid:48) , k ( x + e i , j ) − B ( f ) i (cid:48) , k ( x ) ≥ k and B ( f ) i (cid:48) , j (cid:48) ( x + e i , j ) − B ( f ) i (cid:48) , j (cid:48) ( x ) =
1, we have f i (cid:48) ( ψ ( x ) + e i ) − f i (cid:48) ( ψ ( x )) >
0. This in turn means that there exists a positive edge i → i (cid:48) in G f ( ψ ( x )). A similar argument applies to the case when x i , j = (cid:3) Intuitively speaking, (1) says all the regulation in the original multilevel network is captured in theconverted Boolean network, while (2) says all the regulation in the converted network comes from theoriginal multilevel network.
Corollary 2.
An asymptotic evolution function f over M has a positive (negative) type-1 functionalcircuit if so does its binarisation B ( f ). ADRIEN FAUR ´E AND SHIZUO KAJI
Note that two negative arrows ( i , j ) → ( i , j (cid:48) ) and ( i , j (cid:48) ) → ( i , j ) in G B ( f )( x ), both of which corre-spond to a negative self-loop i → i in G B ( f )( x ), can be composed to produce a positive circuit. Thispositive functional type-1 circuit corresponds the one which is the composition of a single negativeself-circuit with itself in G f ( ψ ( x )). Example 4.
We give two characteristic examples of the binarisation.Consider the evolution function f ( y ) = − y over M = { , , } . Its binarisation is B ( f )( x ) = (1 ,
1) ( x = (0 , x ( x = (1 , , (0 , ,
0) ( x = (1 , . The corresponding state transition graphs are0 (cid:47) (cid:47) (cid:111) (cid:111) , (cid:62) (cid:62) (cid:32) (cid:32) (cid:96) (cid:96) (cid:126) (cid:126) y = x = (0 ,
0) respectively look:
G f (0) = y − (cid:5) (cid:5) , G B ( f )(0 , = x − (cid:35) (cid:35) x − (cid:102) (cid:102) Notice that the self-loop on y corresponds to each of the edges x → x and x ← x . In particular,the converse to Corollary 2 does not hold.Consider another evolution function over M = { , , } defined by f ( y ) = y ( y = , y = . Its binari-sation is B ( f )( x ) = x ( x (cid:44) (1 , ,
0) ( x = (1 , . The corresponding state transition graphs are0 1 2 (cid:111) (cid:111) , (cid:96) (cid:96) (cid:126) (cid:126) y = x = (0 ,
1) respectively look:
G f (1) = y − (cid:5) (cid:5) + (cid:87) (cid:87) , G B ( f )(0 , = x − (cid:35) (cid:35) x + (cid:120) (cid:120) Another extension method.
Recently, Tonello [16] has independently constructed a mappingwhich also extends Van Ham’s while preserving the dynamics and the local regulations in a morestringent sense than ours. Her method was used to produce a counter-example to Conjecture 1 aswell. Her mapping can be described in our context as follows: f (cid:55)→ b ◦ f ◦ ψ, where f is a stepwise function. Compared to ours, her method yields fewer arrows in the statetransition graph. Her strategy was to stipulate the converted function to take values in the admissible CIRCUIT-PRESERVING MAPPING FROM MULTILEVEL TO BOOLEAN DYNAMICS 9 region Im( b ), whereas ours was to equip the converted function with the symmetry described inTheorem 1. 3. R esults Lambda phage.
As an illustration, we first apply our method to the 2-variable lambda phagemodel proposed by Thie ff ry and Thomas [17].The lambda phage is a bacterial virus that infects E. coli . It is a temperate phage, i.e. it caneither multiply and eventually kill the host cell ( lytic phase), or integrate its DNA into the bacterialchromosome ( lysogenic phase), conferring the cell immunity against super-infection by other lambdaphage. The switch between lysis and lysogeny, as modelled by Thie ff ry and Thomas, is essentiallycontrolled by a positive feedback circuit between genes cI and cro . The two genes inhibit each other,such that cI dominates the lysogenic phase, whereas cro is active during the lytic phase. The genecro further inhibits its own activity. The system is modelled by a discrete system with the state space M = { ( cI , cro ) ∈ { , } × { , , }} . The dynamics displays a stable state with high cI and low cro activity, and a two-state cyclic attractor with low cI , and cro oscillating around its activity threshold(Fig. 1, left), which the authors describe as homeostasis.F igure State transition graphs for the lambda phage model.
Left, the original,multilevel model [17]; centre, the Boolean version obtained using Van Ham’s method;right, the Boolean version obtained using our method.Table 1 shows how the same dynamics can be encoded using a stepwise or asymptotic evolutionfunction. Notice that the stepwise function creates a positive feedback on cro ( ˆ f cro (0 , − ˆ f cro (0 , > f cro (1 , − ˆ f cro (1 , >
0) that is not visible in the asymptotic function. This di ff erence in the globalregulatory graphs is shown in Fig. 2.T able
1. Truth table of the lambda phage model cI cro ˆ f cI ˆ f cro ¯ f cI ¯ f cro cI and cro in each state, encoded using either the stepwise orasymptotic evolution function.Boolean systems are generated by Van Ham’s method and ours with the state space { ( cI , cro , cro ∈ B } . However, since Van Ham’s method yields a dynamics with as many states as the original, mul-tilevel dynamics (Fig. 1, centre), the system thus obtained includes a ”non-admissible” region (grey area in the Figure) whose states do not have any counterpart in the multilevel model. To make compar-ison, we extended the domain of the Boolean model obtained by Van Ham’s method (Fig. 1, centre)by completing the grey dashed arrows in such a way that • it does not create any extra arrows in the global interaction graph that is not already visibleelsewhere within the admissible region • and there is no outgoing arrow in the state transition graph from an admissible state to anon-admissible state.This extension is based on our understanding of Van Ham’s original publication [6]. The correspond-ing global regulatory graph is shown in Fig. 2 (centre).F igure Global interaction graphs for the lambda phage model.
Left, the origi-nal, multilevel model [17]: top, stepwise evolution function, bottom, asymptotic evo-lution function; centre, the Boolean version obtained using Van Ham’s method; right,the Boolean version obtained using our method.In contrast, the dynamics produced by our method is more complex, and it occupies the whole statespace: every state has a counterpart in the original multilevel model. However, the two-state attractorof the original model is now represented by a three-state attractor, where state 011 corresponds tostate 02 and both states 001 and 010 correspond to state 01. Another di ff erence is that, in the modelobtained with Van Ham’s method, variables cro cro cro . In the Boolean model generated with our method, cro cro ff erence appears at the local level (Fig. 3). Using Van Ham’s method, localgraphs may include edges that have no visible counterpart in the local graph of the correspondingmultilevel state. For example, while in 02 the only visible regulation is the negative loop on cro ,Van Ham’s method adds an edge from cro cI in 011, whereas the corresponding local graphobtained using our method includes only regulations between cro cro
2. Similarly, in 10, theonly visible edges occur between cI and cro , and the same is true in 100 using our method; however,using Van Ham’s method an additional edge between cro cro cro cro
2, which in multilevel model correspondsto the composition of the negative self circuit on cro with itself. Such positive circuits can only appearbetween variables that represent the same multilevel variable (see Corollary 2).3.2.
Boolean counter-example to Conjecture 1.
One of the main goals in genetic regulatory net-work analysis is to find relations between circuits in the interaction graphs and attractors in the statetransition graphs. Here we list a few known results in this direction:(1) If f has no type-1 functional circuit, Γ ( f ) has a unique stable state x ([18]).(2) If f has no type-1 functional positive circuit, Γ ( f ) has a unique attractor ([11, 10]).(3) If there exists a cyclic attractor in Γ ( f ), G f ( M ) must contain a negative circuit ([12]).Notice that the first two statements connect the asymptotic behaviour of the dynamics with the local interaction graph at some x ∈ M , while the last one does so with the global interaction graph G f ( M ).This naturally gives rise to the following conjectures: Conjecture 1 ([9, Question 1],[12, Question 1,2]) . CIRCUIT-PRESERVING MAPPING FROM MULTILEVEL TO BOOLEAN DYNAMICS 11 F igure Local interaction graphs for the lambda phage model.
Left, the original,multilevel model [17]; centre, the Boolean version obtained using Van Ham’s method;right, the Boolean version obtained using our method. Top, local graphs for low cI and high cro (states 02 or 011); bottom, local graphs for states with high cI and low cro (states 10 and 100). Local graphs in 02 and 10 are identical under the stepwise orasymptotic evolution function.(1) If f has no type-1 functional negative circuit, then Γ ( f ) has no cyclic attractors.(2) If f has no type-1 functional negative circuit, then Γ ( f ) has a stable state.Note that the second statement is weaker in the sense that it follows from the first one.Richard, together with Comet [12, Example 6], gave a counter-example to the conjectures when M = { , , , } . The grid graph Γ over M = { , , , } given in [12, Example 6] has no stable stateand there exists no type-1 negative functional circuit in the interaction graph of ¯ f Γ (see Fig. 4). Theexistence of Boolean counter-examples was left as an open problem in [12].F igure A multilevel counter-example given in [12]By Corollary 2, the Boolean grid graph B ( Γ ) yielded by our method (Fig. 5) has no type-1 negativefunctional circuit in the interaction graph of the binarisation B ( ¯ f Γ ). Furthermore, by Corollary 1 thereis no stable state in the state transition graph of B ( ¯ f Γ ). Thus, we obtain a new Boolean counter-example to Conjecture 1.Finally, we note that Ruet recently gave a systematic way to produce counter-examples to the con-jectures for the Boolean case [13]. Our method is very di ff erent from Ruet’s and while his example F igure A Boolean counter-example to Conjecture 1
The binarisation of the mul-tilevel counter-examples gives a Boolean counter-examplepossesses a special property that it has an attractive cycle , our counter-example has a smaller dimen-sion of 6. For comparison, we include an example based on Ruet’s theory [13] Fig. 6.F igure Antipodal attractive cycle with no negative circuit.
States not shown areall stable states. Adapted from [13].3.3.
Implementation.
We implemented our method in the form of a Perl script. For comparison,Tonello’s method [16] is also implemented in the same script. It is available at [19] under the MITlicense. The input is a multilevel evolution function described in the Truth table format [20] and theoutput is a Boolean evolution function described in the same format. The maximum levels m i for each CIRCUIT-PRESERVING MAPPING FROM MULTILEVEL TO BOOLEAN DYNAMICS 13 gene is automatically detected. When the input function is defined on { , , . . . , m } × { , , . . . , m } ×· · · × { , , . . . , m n } , the converted Boolean network is defined on B m + m + ··· + m n .4. D iscussion and conclusions In the discrete formulation of gene regulatory networks, a system is commonly modelled by afunction. When some genes take more than two levels, there are multiple choices for functions hav-ing the same (asynchronous) state transition graph. We single out a unique choice, which we callthe asymptotic evolution function (Proposition 1). Then, we introduced a mapping which convertsan asymptotic evolution function to a Boolean evolution function (Definition 4). This mapping pre-serves dynamical and regulatory properties (Proposition 2, Theorem 2), thus allowing us to analysemultilevel networks by methods developed for Boolean networks.Mappings from multilevel to Boolean networks have been used in the study of gene regulatory net-works. In particular, Van Ham’s mapping has been shown by Didier et al. to be essentially the onlymethod to provide a one-to-one, neighbour-preserving and regulation-preserving Boolean representa-tion of multilevel models [7]. However, although the authors did suggest that the mapping could beuseful to study the role of regulatory circuits, the question of how interaction functionality contextsare preserved had not been studied so far.One such instance is Thomas’s conjecture, which states that the existence of a cyclic attractor inthe asynchronous state transition graph requires that of local negative circuit. The conjecture hasrecently been given a Boolean counter-example by P. Ruet [13]. Until then, although A. Richard andand J-P. Comet had produced a multilevel counter-example [12], the Boolean case remained open. Itis straightforward to apply Van Ham’s method to the counter-example in order to obtain a Booleanmodel, but it is defined only on the admissible region. To extend the model to the whole Booleanstate space, while preserving its dynamics and regulatory relation, is highly non-trivial. The methodwe propose has been designed specifically to circumvent the problem. The idea was to avoid extrainteractions and circuits by extending the state transition graph obtained with Van Ham’s method insuch a way that we have “parallel trajectories” going through the whole state space. This was achievedby loosening the one-to-one criterion such that states including intermediate values in the multilevelmodel would match with several states in the Boolean version, e ff ectively creating equivalent Booleantransitions for each multilevel transition in the original model.In a sense, our method works opposite to Van Ham’s: instead of embedding multilevel states intoBoolean states, we define a Boolean model such that each Boolean state can be mapped to a multilevelstate. With our method, interaction functionality is preserved, and thus all local interaction graphs inthe Boolean model come from their counterpart in the original, multilevel model. In contrast, VanHam’s method only preserves the global interaction graph.One limitation of our method is that the synchronous dynamics of the original multilevel model cannot be directly retrieved from the Boolean model produced by the conversion. Here, by synchronousdynamics we mean the state transition graph having edges x → f ( x ) instead of (1) (see, for exam-ple, [21]). Synchronous state transition is often deemed unrealistic since it assumes all processes arerealised simultaneously with the same delay (see discussion by Abou-Jaoud´e et al. [1]). Neverthe-less, the synchronous mode is still a popular update method in simulation due to its simplicity, andoccasionally used for multilevel models (see e.g. Chifman et al. [22] for a recent example). If ourbinarisation is used for a synchronous simulation, any increase in a variable would translated into itsincrease to the maximum value. It is worth noting that the counter-examples for Conjecture 1 consid-ered in § ff erent ways to represent a multilevel system, fordi ff erent ways can represent the same model [15], which causes ambiguities in the notation. A cknowledgements The authors are grateful to Paul Ruet for explaining his result in [13], to Elisa Tonello for fruitfuldiscussion and careful reading of our draft, and to Yuki Ikawa and Sergey Tishchenko for their help inthe early stage of this work. The second named author is partially supported by JST PRESTO GrantNumber JPMJPR16E3, Japan. R eferences [1] Abou-Jaoud´e, W., Traynard, P., Monteiro, P.T., Saez-Rodriguez, J., Helikar, T., Thie ff ry, D., Chaouiya, C.: Logi-cal modeling and dynamical analysis of cellular networks. Front Genet , 94 (2016). doi: [2] Stoll, G., Viara, E., Barillot, E., Calzone, L.: Continuous time boolean modeling for biological signaling: applicationof gillespie algorithm. BMC Syst Biol , 116 (2012). doi: [3] MacNamara, A., Terfve, C., Henriques, D., Pe˜nalver Bernab´e, B., Saez-Rodriguez, J.: State–time spectrum of signaltransduction logic models. Physical Biology (4), 045003 (2012)[4] Helikar, T., Kowal, B., McClenathan, S., Bruckner, M., Rowley, T., Madrahimov, A., Wicks, B., Shrestha, M.,Limbu, K., Rogers, J.A.: The cell collective: toward an open and collaborative approach to systems biology. BMCSyst Biol , 96 (2012). doi: [5] Remy, ´E., Ruet, P., Thie ff ry, D.: Positive or negative regulatory circuit inference from multilevel dynamics. LectureNotes in Control and Information Sciences , 263–70 (2006)[6] Van Ham, P.: Kinetic logic: A boolean approach to the analysis of complex regulatory systems. Lecture notes inBiomathematics., vol. 29, pp. 326–43. Levin, S. (1979). Chap. How to deal with variables with more than two levels[7] Didier, G., Remy, E., Chaouiya, C.: Mapping multivalued onto boolean dynamics. J. Theor. Biol. (1), 177–84(2011). doi: [8] Thomas, R.: On the relation between the logical structure of systems and their ability to generate multiplesteady states or sustained oscillations. In Springer series in synergies: Numerical methods in the study of crit-ical phenomena, vol. 9, pp. 180–193. Della Dora, Jean and Demongeot, Jacques and Lacolle, Bernard (1981).doi: [9] Comet, J.-P., Noual, M., Richard, A., Aracena, J., Calzone, L., Demongeot, J., Kaufman, M., Naldi, A., Snoussi,E.H., Thie ff ry, D.: On circuit functionality in boolean networks. Bull. Math. Biol. (6), 906–19 (2013). doi: [10] Remy, ´E., Ruet, P., Thie ff ry, D.: Graphic requirements for multistability and attractive cycles in a boolean dynamicalframework. Advances in Applied Mathematics , 335–350 (2008). doi: [11] Richard, A., Comet, J.-P.: Necessary conditions for multistationarity in discrete dynamical systems. Discrete AppliedMathematics , 2403–13 (2007). doi: [12] Richard, A.: Negative circuits and sustained oscillations in asynchronous automata networks. Adv. Appl. Math. (4), 378–392 (2010). doi: [13] Ruet, P.: Negative local feedbacks in boolean networks. Discrete Applied Mathematics , 1–17 (2017). doi: [14] Thomas, R., D’Ari, R.: Biological Feedback. CRC Press, Inc., 2000 Corporate Blvd, N.W., Boca raton, Florida,33431, United States of America (1990). http://hal.archives-ouvertes.fr/hal-00087681/en/ [15] Streck, A., Lorenz, T., Siebert, H.: Minimization and equivalence in multi-valued logical models of regulatorynetworks. Natural Computing , 555–566 (2015). doi: [16] Tonello, E.: On the conversion of multivalued gene regulatory networks to Boolean dynamics. ArXiv e-prints (2017). [17] Thie ff ry, D., Thomas, R.: Dynamical behaviour of biological regulatory networks–ii. immunity control in bacterio-phage lambda. Bull. Math. Biol. (2), 277–97 (1995). doi: [18] Shih, M.-H., Dong, J.-L.: A combinatorial analogue of the jacobian problem in automata networks. Advances inApplied Mathematics (1), 30–46 (2005). doi: [19] Kaji, S.: A multilevel to Boolean gene regulatory network converter. GitHub (2017)[20] Naldi, A., Berenguier, D., Faur´e, A., Lopez, F., Thie ff ry, D., Chaouiya, C.: Logical modelling of regulatory networkswith ginsim 2.3. BioSystems (2), 134–9 (2009). doi: [21] Remy, E., Moss´e, B., Chaouiya, C., Thie ff ry, D.: A description of dynamical graphs associated to elementaryregulatory circuits. Bioinformatics
19 Suppl 2 , 172–8 (2003). doi: [22] Chifman, J., Arat, S., Deng, Z., Lemler, E., Pino, J.C., Harris, L.A., Kochen, M.A., Lopez, C.F., Akman, S.A.,Torti, F.M., Torti, S.V., Laubenbacher, R.: Activated oncogenic pathway modifies iron network in breast epithelialcells: A dynamic modeling perspective. PLoS Comput. Biol. (2), 1005352 (2017). doi: (A. Faur´e) D epartment of P hysics and I nformation S cience CIRCUIT-PRESERVING MAPPING FROM MULTILEVEL TO BOOLEAN DYNAMICS 15 Y amaguchi U niversity oshida , Y amaguchi apan E-mail address : [email protected] (S. Kaji) D epartment of M athematics Y amaguchi U niversity oshida , Y amaguchi apan/ JST PRESTO
E-mail address ::