Rainer W Kühne
email: kuehne70@gmx.de

Derivation of the Local Lorentz Gauge Transformation of a Dirac Spinor Field in Quantum Einstein-Cartan Theory

Keywords: Einstein-Cartan theory, quantum gravity, energy-momentum tensor, Lorentz transformations, gauge invariance, quantum field theories

I examine the groups which underly classical mechanics, non-relativistic quantum mechanics, special relativity, relativistic quantum mechanics, quantum electrodynamics, quantum flavourdynamics, quantum chromodynamics, and general relativity. This examination includes the rotations SO(2) and SO(3), the Pauli algebra, the Lorentz transformations, the Dirac algebra, and the U (1), SU (2), and SU (3) gauge transformations. I argue that general relativity must be generalized to Einstein-Cartan theory, so that Dirac spinors can be described within the framework of gravitation theory.

1 Introduction I will present my argumentation in a didactic way. So I will start at a very elementary level. If not denoted otherwise then Latin indices run from 1 to 3, Greek indices run from 0 to 3. I will use the Einstein summation convention, where it is summed over all indices which appear twice. Moreover I will use the natural units

h = c = ε 0 = 1 (1)
where h = h/2π denotes the reduced Planck constant, c the speed of light, and ε 0 the electric field constant. Inner indices of matrices will be dropped.

Readers who would like to learn more about the theories which will be mentioned in this paper are recommended to the comprehensive works of Goldstein [START_REF] Goldstein | Classical Mechanics[END_REF] (classical mechanics), Jackson [START_REF] Jackson | Classical Electrodynamics[END_REF] (classical electrodynamics), Kilmister [START_REF] Kilmister | Special Theory of Relativity[END_REF] (special relativity), Weinberg [START_REF] Weinberg | Gravitation and Cosmology[END_REF] and Misner, Thorne and Wheeler [START_REF] Misner | Gravitation[END_REF] (general relativity), Hehl et al. [START_REF] Hehl | General Relativity with Spin and Torsion: Foundations and Prospects[END_REF] (Einstein-Cartan theory), Messiah [START_REF] Messiah | Quantum Mechanics[END_REF] (non-relativistic quantum mechanics), Bjorken and Drell [START_REF] Bjorken | Relativistic Quantum Mechanics[END_REF] (relativistic quantum mechanics), Bjorken and Drell [START_REF] Bjorken | Relativistic Quantum Fields[END_REF] (quantum electrodynamics), Taylor [START_REF] Taylor | Gauge Theories of Weak Interactions[END_REF] (quantum flavourdynamics), Politzer [START_REF] Politzer | Asymptotic Freedom: An Approach to Strong Interactions[END_REF] and Marciano and Pagels [START_REF] Marciano | Quantum Chromodynamics[END_REF] (quantum chromodynamics), and Kühne [START_REF] Kühne | Quantum Field Theory with Electric-Magnetic Duality and Spin-Mass Duality but Without Grand Unification and Supersymmetry[END_REF] (quantum electromagnetodynamics).

Classical Mechanics

The space of classical mechanics is described by the three-dimensional Euclidian space. The scalar product of the three-vectors a i and b i is given by a

• b = δ ij a i b j (2)
where

δ ij = diag (1, 1, 1) (3)
denotes the Kronecker symbol. The vector product is given by

(a × b) k = ε ijk a i b j (4)
where ε ijk denotes the totally anti-symmetric Levi-Civita symbol. The square of the infinitesimal line element is given by

ds 2 = δ ij dx i dx j (5)
where dx i denotes the infinitesimal coordinate difference of the two space-points x i and y i

dx i = lim yi→xi (y i -x i ) (6)
A rotation R around the rotation angle ϕ in the two-dimensional subspace is described by the orthogonal rotation group SO(2)

R = exp(-iϕD) = cos ϕ -sin ϕ sin ϕ cos ϕ ∈ SO(2) (7)
where, because of the Euler equation

e -iϕ = cos ϕ -i sin ϕ (8)
the generator of the rotation is

D = 0 -i i 0 (9)
A rotation R around the rotation angle ϕ i in the three-dimensional space is described by the rotation group SO(3)

R = exp(-iϕ i D i ) ∈ SO(3) (10)
where the generators D i of the rotation satisfy the commutator relation

[D i , D j ] = iε ijk D k (11)
A representation of the generators of the SO(3) group is

D 1 =   0 0 0 0 0 -i 0 i 0   (12)
D 2 =   0 0 i 0 0 0 -i 0 0   (13)
D 3 =   0 -i 0 i 0 0 0 0 0   (14)
Note that the number of generators of the orthogonal groups is given by

dim SO(n) = n(n -1)/2 (15)
and that the generators of SO(n) are Hermitean.

Non-Relativistic Quantum Mechanics

Classical angular momentum is continuous. In quantum physics orbital angular momentum is quantized in units of h and intrinsic spin is quantized in units of h/2. Intrinsic spin is described by the unitary group SU (2). Its generators are the three Pauli matrices σ i . An often used representation of the Pauli matrices is

σ 1 = 0 1 1 0 (16) σ 2 = 0 -i i 0 (17) σ 3 = 1 0 0 -1 (18)
The commutator of the group SU (2) is given by

[σ i , σ j ] = 2iε ijk σ k (19)
The groups SU (2) and SO(3) are local isomorphic for angles 0 ≤ ϕ < 2π where the group SU (2) is covered by the group SO(3). Note also that the commutators eq. ( 11) and eq. ( 19) differ by a factor of two. The anti-commutator is

{σ i , σ j } = 2δ ij (20)
The multiplication of two three-vectors a i and b i is given by

(σ • a)(σ • b) = (σ i • a i )(σ j • b j ) = (δ ij + iε ijk σ k )a i b j = a • b + iσ • (a × b) (21)
According to non-relativistic quantum mechanics both spin and isospin are invariant under global transformations of the group SU (2). If Ψ denotes a two-component Pauli-spinor, ϕ i the three-component rotation angle vector, σ i the Pauli matrices, and x a space-time point, then

Ψ (x) = exp(-iϕ i σ i /2)Ψ(x) (22)
If Ψ denotes a two-component iso-spinor, ϕ i the three-component phase vector, and τ i the Pauli matrices, then

Ψ (x) = exp(-iϕ i τ i /2)Ψ(x) (23)
4 Special Relativity

The special theory of relativity is invariant under the semi-simple Poincare group. The parameters of its translational part are time and three-position. The parameters of its rotational part are the rotation angle three-vector and the three-component Lorentz boost.

The scalar product of the two four-vectors a µ and b µ is given by

a • b = g µν a µ b ν (24)
where g µν denotes the metric tensor. It is

g µν = g µν (25) g µ ν = g ν µ = δ µ ν = diag (1, 1, 1, 1) (26)
In Minkowski coordinates the metric tensor is represented by

g µν = diag (1, -1, -1, -1) (27)
The square of the infinitesimal line element is given by

ds 2 = g µν dx µ dx ν (28)
where dx µ denotes the infinitesimal coordinate difference of the two space-time points x µ and

y µ dx µ = lim y µ →x µ (y µ -x µ ) (29)
A rotation around the z-axis by the rotation angle ϕ is given by

a µ ν =     1 0 0 0 0 cos ϕ -sin ϕ 0 0 sin ϕ cos ϕ 0 0 0 0 1     (30)
A Lorentz boost along the x-axis by the speed v is given by

a µ ν =     γ -βγ 0 0 -βγ γ 0 0 0 0 1 0 0 0 0 1     (31)
where

β = v/c (32) γ = 1/ 1 -v 2 /c 2 (33)
In general, a Lorentz transformation a µ ν of a four-position x µ is given by

x µ = a µ ν x ν (34)
The Lorentz transformation of a four-derivative

∂ µ = ∂/∂x µ is given by ∂ µ = a µ ν ∂ ν (35)
Finally, a Lorentz transformation around the parameter ω αβ is given by

a µ ν = exp 1 2 ω αβ I αβ µ ν ∈ O(1, 3) (36)
where a representation of the six generators of the Lorentz group SO(1, 3) is

I 10 =     0 -1 0 0 1 0 0 0 0 0 0 0 0 0 0 0     (37)
I 20 =     0 0 -1 0 0 0 0 0 1 0 0 0 0 0 0 0     (38)
I 30 =     0 0 0 -1 0 0 0 0 0 0 0 0 1 0 0 0     (39)
I 13 =     0 0 0 0 0 0 0 -1 0 0 0 0 0 1 0 0     (40)
I 23 =     0 0 0 0 0 0 0 0 0 0 0 -1 0 0 1 0     (41)
I 12 =     0 0 0 0 0 0 -1 0 0 1 0 0 0 0 0 0     (42)

Relativistic Quantum Mechanics

The group which underlies the kinematics of relativistic quantum mechanics and relativistic quantum field theory is the Poincare group.

Isospin is invariant under global transformations of the group SU (2). The generators are the three Pauli matrices τ i . The three Pauli matrices σ i of spin are generalized by the four Dirac matrices γ µ . An often used representation of the Dirac matrices is

γ 0 = 1 0 0 -1 (43)
γ i = 0 σ i -σ i 0 (44)
It is

γ 0 = γ 0 ( 45
)
γ i = -γ i ( 46
)
The commutator is

[γ µ , γ ν ] = -2iσ µν (47)
which is the definition of the generalized Dirac matrices σ µν . The anti-commutator gives

{γ µ , γ ν } = 2g µν (48)
By using the representation eqs. ( 43) and (44) of the Dirac matrices and the representation eqs. ( 16), ( 17) and (18) of the Pauli matrices, the anti-commutator gives the representation eq. ( 27) of the metric tensor in Minkowski coordinates. Moreover it is

γ 5 = γ 5 = -iγ 0 γ 1 γ 2 γ 3 (49)
The multiplication of two four-vectors a µ and b µ is given by

(γ • a)(γ • b) = (γ µ a µ )(γ ν b ν ) = (g µν -iσ µν )a µ b ν (50)
The Lorentz transformation of a four-component spin 1/2 Dirac spinor field Ψ(x) is given by

Ψ (x ) = exp - i 4 σ µν 1 2 ω αβ I αβ µν Ψ(x)
(51) which is a generalization of eq. (36) which considers both the Lorentz transformation of the elementary particle and its intrinsic spin.

Relativistic Quantum Field Theory

Quantum electrodynamics is invariant under local transformations of the gauge group U (1).

If Ψ describes a four-component Dirac spinor, x a space-time point, e the elementary electric charge, ϕ the gauge phase, D µ the covariant derivative, ∂ µ the partial four-derivative, and A µ the electromagnetic four-potential, then the gauge transformation is

Ψ (x) = exp(-ieϕ(x))Ψ(x) (52)
iD µ = i∂ µ -eA µ (x) (53)
A µ (x) = A µ (x) -∂ µ ϕ(x) (54)
Note that U (1) and SO(2) are isomorphic.

Quantum flavourdynamics is invariant under local transformations of the gauge group SU (2) × U (1). If Ψ denotes an eight-component iso-spinor Dirac spinor, x a space-time point, g the weak coupling constant, ϕ i the threecomponent gauge phase isovector, τ i the Pauli matrices, and W i µ the weak isovector fourpotentials, then the gauge transformations of the SU (2) part are

Ψ (x) = exp(-igϕ i (x)τ i /2)Ψ(x) (55) iD µ = i∂ µ - g 2 W i µ (x)τ i ( 56
)
W µ i (x) = W µ i (x) -∂ µ ϕ i (x) -gε ijk ϕ j (x)W µ k (x) ( 57
)
Quantum chromodynamics is invariant under local transformations of the gauge group SU (3). If Ψ denotes a twelve-component colourvector Dirac spinor, x a space-time point, g the strong coupling constant, ϕ i the eightcomponent gauge phase vector, λ i the eight Gell-Mann matrices, and G i µ the eight gluon fourpotentials, then the gauge transformations are

Ψ (x) = exp(-igϕ i (x)λ i /2)Ψ(x) (58) iD µ = i∂ µ - g 2 G i µ (x)λ i (59) G µ i (x) = G µ i (x) -∂ µ ϕ i (x) -gf ijk ϕ j (x)G µ k (x) (60)
where the indices i, j, k run from 1 to 8.

A representation of the Gell-Mann matrices is

λ 1 =   0 1 0 1 0 0 0 0 0   (61) λ 2 =   0 -i 0 i 0 0 0 0 0   ( 62
)
λ 3 =   1 0 0 0 -1 0 0 0 0   (63)
λ 4 =   0 0 1 0 0 0 1 0 0   (64)
λ 5 =   0 0 -i 0 0 0 i 0 0   (65)
λ 6 =   0 0 0 0 0 1 0 1 0   (66)
λ 7 =   0 0 0 0 0 -i 0 i 0   (67)
λ 8 = 1 √ 3   1 0 0 0 1 0 0 0 -2   (68)
The commutator of the Gell-Mann matrices is

[λ i , λ j ] = 2if ijk λ k (69)
and the anti-commutator is

{λ i , λ j } = 4 3 δ ij + 2d ijk λ k (70)
Note that the number of the generators of the unitary groups is given by

dim SU (n) = n 2 -1 (71)
and that the generators of the SU (n) groups are Hermitean.

According to relativistic quantum mechanics and relativistic quantum field theory, the energy-momentum tensor Σ µν of a Dirac spinor Ψ is asymmetric

Σ µν (x) = - 1 2 D µ Ψ(x) γ ν Ψ(x) + 1 2 Ψ(x)γ ν D µ Ψ(x) (72)
where the covariant derivative D µ is given by the equations ( 53), ( 56) and (59).

General Relativity

The general theory of relativity is invariant under arbitrary curvilinear transformations. It is invariant under local transformations of the Lorentz group SO(1, 3). The Lorentz boosts depend on the four-position x µ . Note that the Poincare group which underlies special relativity has ten generators, whereas the Lorentz group which underlies general relativity has only six generators. Therefore general relativity is not a generalization of special relativity. The scalar product of the four-vectors a µ and b µ is given by

a • b = g µν (x)a µ b ν (73)
where the metric tensor g µν depends on the space-time point x. The square of the infinitesimal line element is

ds 2 = g µν (x)dx µ dx ν (74)
The local Lorentz transformation is

a µ ν (x) = exp 1 2 ω αβ (x)I αβ µ ν (75)
which is a generalization of eq. ( 36). When a four-vector C α is parallely displaced from the four-position x µ to the four-position x µ + dx µ , then it changes according to the prescription

dC α = -Γ α µν (x)C ν dx µ (76)
This is the definition of the four-positiondependent affine connection Γ α µν . According to general relativity it has only a symmetric part

{} α µν (x) = 1 2 Γ α µν (x) + Γ α νµ (x) (77)
which is named Christoffel symbol. The Riemann curvature tensor is given by

R λ µνκ (x) = ∂ κ Γ λ µν (x) -∂ ν Γ λ µκ (x) +Γ α µν (x)Γ λ κα (x) -Γ α µκ (x)Γ λ να (x) (78)
By contraction one gets the Ricci tensor

R µν (x) = R λ µλν (x) (79)
and the Ricci scalar

R(x) = R µ µ (x) (80)
The Einstein field equations are

R µν (x) - 1 2 g µν (x)R(x) = κΣ µν (x) (81)
where κ = -8πG

denotes the Einstein field constant and Σ µν the energy-momentum tensor.

Since the affine connection (Christoffel symbol) is symmetric it follows that the energymomentum tensor of general relativity is symmetric. This is in contrast to the asymmetric energy-momentum tensor of a Dirac spinor of relativistic quantum mechanics. This means that a Dirac spinor cannot be described by the geometry that underlies general relativity.

Einstein-Cartan Theory

Einstein-Cartan theory is a generalization of general relativity, because within the framework of this theory the anti-symmetric part of the affine connection which is named Cartan's torsion tensor

T α µν (x) = 1 2 Γ α µν (x) -Γ α νµ (x) (83)
is nonzero. In contrast to the Christoffel symbol the torsion tensor transforms as a tensor under arbitrary curvilinear transformations. Cartan's torsion tensor is related to the spin tensor by

T αµν (x) = κτ αµν (x) (84)
where κ = -8πG denotes the Einstein field constant known from general relativity. The Dirac spinor can be described by the geometry that underlies Einstein-Cartan theory, because this theory describes an asymmetric affine connection and therefore an asymmetric energy-momentum tensor.

Einstein-Cartan theory is invariant under local transformations of the Poincare group. So Einstein-Cartan theory is a generalization of special relativity.

In quantum Einstein-Cartan theory the local Lorentz gauge transformation of a Dirac spinor field Ψ is given by

Ψ (x ) = exp - i 4

Conclusion

I argued that the general theory of relativity must be generalized by the Einstein-Cartan theory, because the asymmetric energy-momentum tensor of relativistic quantum mechanics requires the existence of an asymmetric affine connection and therefore the existence of a nonzero torsion tensor. Moreover I pointed out that spin is the source of a gauge field which is associated with Cartan's torsion. This is analogous to the situation of isospin which is a pure quantum number in quantum mechanics, but the source of the weak gauge field in quantum flavourdynamics.

  This equation is a combination of eqs. (51) and (75).


Madalena Chaves
email: madalena.chaves@inria.fr

Daniel Figueiredo
email: daniel.figueiredo@ua.pt

Manuel António Martins
email: martins@ua.pt

Boolean dynamics revisited through feedback interconnections

Keywords: Boolean models, feedback interconnections, attractor computation, asynchronous vs. synchronous updates

Boolean models of physical or biological systems describe the global dynamics of the system and their attractors typically represent asymptotic behaviors. In the case of large networks composed of several modules, it may be difficult to identify all the attractors. To explore Boolean dynamics from a novel viewpoint, we will analyse the dynamics emerging from the composition of two known Boolean modules. The state transition graphs and attractors for each of the modules can be combined to construct a new asymptotic graph which will (1) provide a reliable method for attractor computation with partial information; (2) illustrate the differences in dynamical behavior induced by the updating strategy (asynchronous, synchronous, or mixed); and (3) show the inherited organization/structure of the original network's state transition graph.

Introduction

The effect of network topology in a system's dynamics and stability has long been a problem of interest [START_REF] Kauffman | Genetic networks with canalyzing boolean rules are always stable[END_REF][START_REF] Demongeot | Robustness in regulatory networks: A multidisciplinary approach[END_REF][START_REF] Demongeot | Attraction basins as gauges of robustness against boundary conditions in biological complex systems[END_REF]. In the presence of very large networks of interaction with uncertain knowledge or data involving large maps of activation/inactivation effects, Boolean models with their qualitative On/Off description, provide a practical framework to analyse and model such data [START_REF] Zhang | Network model of survival signaling in large granular lymphocyte leukemia[END_REF][START_REF] Saez-Rodriguez | A logical model provides insights into T cell receptor signaling[END_REF]. Therefore, it is very useful to develop new techniques to analyse and render Boolean models more realistic (software [START_REF] Klamt | A methodology for the structural and functional analysis of signaling and regulatory networks[END_REF][START_REF] Abou-Jaoudé | Logical modeling and dynamical analysis of cellular networks[END_REF], asynchronous analysis [START_REF] Chaves | Robustness and fragility of boolean models for genetic regulatory networks[END_REF][START_REF] Fauré | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF][START_REF] Tournier | Interconnection of asynchronous Boolean networks, asymptotic and transient dynamics[END_REF], probabilistic approaches [START_REF] Shmulevich | Probabilistic boolean networks: a rule-based uncertainty model for gene regulatory networks[END_REF][START_REF] Stoll | Maboss 2.0: an environment for stochastic boolean modeling[END_REF], or stochastic approaches [START_REF] Mori | Stochastic simulation of boolean rxncon models: towards quantitative analysis of large signaling networks[END_REF]).

A Boolean network is defined by a set of logical rules and a strategy for generating the dynamics, such as synchronous or asynchronous updates, or intermediate combinations [START_REF] Wang | Boolean modeling in systems biology: an overview of methodology and applications[END_REF][START_REF] Garg | Synchronous versus asynchronous modeling of gene regulatory networks[END_REF]. The state transititon graph and its dynamics depend not only on the network topology but also on the updating strategy. In fact, different modes of update encompass fundamentally distinct hypotheses regarding the dynamics of the network, and lead to distinct dynamical behavior and attractors [START_REF] Demongeot | Robustness in regulatory networks: A multidisciplinary approach[END_REF][START_REF] Demongeot | Attraction basins as gauges of robustness against boundary conditions in biological complex systems[END_REF]. Synchronous updates are commonly used but they assume simultaneous timescales for all phenomena, and typically generate cyclic attractors which have no biological meaning. Asynchronous updates are arguably more realistic, consider different timescales for the variables in the system, but they have the disadvantage of involving high calculation costs to obtain the state transition graph for large networks (for instance, those containing above 20 vertices).

The study of large networks, in particular the computation of its steady states and attractors, remains a hard problem. Several algorithms have been proposed to determine the singletons in a Boolean network, using approaches such as constrained programming [START_REF] Devloo | Identification of all steady states in large networks by logical analysis[END_REF], SAT formulas for special classes of functions [START_REF] Akutsu | Determining a singleton attractor of a boolean network with nested canalyzing functions[END_REF], or partitioning the system into smaller networks [START_REF] Hong | An efficient steady-state analysis method for large boolean networks with high maximum node connectivity[END_REF]. Using computational algebra, [START_REF] Veliz-Cuba | Steady state analysis of boolean molecular network models via model reduction and computational algebra[END_REF] construct a very cost-efficient algorithm to compute all singletons in a Boolean network. Finding cyclic attractors is an even more complex task and generally needs computation of the successors of all states (typically unfeasible for networks above 20 vertices), or at least the evaluation of certain partial states [START_REF] Irons | Improving the efficiency of attractor cycle identification in Boolean networks[END_REF]. The notion of seeds, for representing dynamically closed subspaces, is used in [START_REF] Klarner | Computing symbolic steady states of boolean networks[END_REF][START_REF] Klarner | Approximating attractors of boolean networks by iterative ctl model checking[END_REF] to aproximately compute attractors of a network. A new network reduction method in [START_REF] Zañudo | An effective network reduction approach to find the dynamical repertoire of discrete dynamic networks[END_REF] was shown to be very effective for computation of all types of attractors. Another approach, also based on network partitioning, was proposed by Chaves & Tournier [START_REF] Chaves | Predicting the asymptotic dynamics of large biological networks by interconnections of Boolean modules[END_REF][START_REF] Chaves | Analysis tools for interconnected boolean networks with biological applications[END_REF][START_REF] Tournier | Interconnection of asynchronous Boolean networks, asymptotic and transient dynamics[END_REF]: this uses a control theory framework to represent a large network as the feedback interconnection between two or more modules and is guaranteed to find not only the singletons but actually all the attractors of the large network, under asynchronous updates.

Some of these methods can be very efficient for analysis of large networks [START_REF] Zañudo | An effective network reduction approach to find the dynamical repertoire of discrete dynamic networks[END_REF][START_REF] Veliz-Cuba | Steady state analysis of boolean molecular network models via model reduction and computational algebra[END_REF], but our aim is to study an algorithm that allows to link together two known models and observe how the new dynamics emerges from each of the two network's attractors. Feedback interconnections are a widely used concept in control theory [START_REF] Sontag | Mathematical Control Theory[END_REF] and essentially mean that two systems with outputs and inputs may be composed into a large system by connecting each system's output to the input of the other system. Interconnection algorithms are especially relevant in physical and biological sciences, when the aim is to understand the behavior of a large system in terms of its different parts, such as the interactions between metabolic and genetic networks [START_REF] Baldazzi | The carbon assimilation network in Escherichia coli is densely connected and largely signdetermined by directions of metabolic fluxes[END_REF], or between two fundamental biological modules such as the cell cycle and circadian clock [START_REF] Gerard | Temporal self-organization of the cyclin/Cdk network driving the mammalian cell cycle[END_REF]). Simultaneously, new disciplines like synthetic biology aim at the construction of synthetic "biological bricks" that perform a well defined function [START_REF] Purnick | The second wave of synthetic biology: from modules to systems[END_REF] and can then be assembled together. In this context, the behavior of a system must be analysed in the light of its interactions with other systems (the question of retroactivity [START_REF] Vecchio | Modular cell biology: Retroactivity and insulation[END_REF] may arise), and the study of large networks as the interconnection of one or more smaller modules gains full significance.

The goal of this paper is to revisit the dynamical analysis of Boolean networks in the light of modular feedback interconnections. In particular, the following perspectives will be analysed: (i) develop an algorithm allowing to compute the attractors and asymptotic dynamics for the composition of (already known) networks (see Sections 3, 4, and 5); (ii) understand the dynamical behavior of a network as the feedback interconnection of two modules; gain insight on the structure and organization of the whole state space (see Section 2) and its dependence on the updating strategy. The method proposed in [START_REF] Chaves | Predicting the asymptotic dynamics of large biological networks by interconnections of Boolean modules[END_REF][START_REF] Tournier | Interconnection of asynchronous Boolean networks, asymptotic and transient dynamics[END_REF] is based on the computation of a new object, the asymptotic graph, a directed graph whose vertices are the products of attractors of the two (or more) modules that constitute the Boolean network and the edges are computed from the state transition graphs of these modules. This method is guaranteed to recover all the attractors of the large network (Section 3) but, for very particular cases, new spurious attractors may also be generated. Here, we will first show how to extend this method to avoid spurious attractors (Section 3.2). Second, we will show how to compute an asymptotic graph in the case of synchronous (Section 4) or mixed updating strategies (Section 5). Finally, several examples of interconnections of biological systems will be analyzed to illustrate our methods (Section 7).

Boolean module feedback interconnections and useful notation

Our objective is to study the dynamics of Boolean networks and their asymptotic behavior, as emergent properties of the interactions between two smaller separate networks. Throughout this paper, we would like to answer the following questions: How are the pathways constructed, in a Boolean network composed of different modules? How to compute the attractors of the large network, depending on updating strategies? Here, the focus is on synchronous versus asynchronous strategies, the comparison of their respective dynamical properties, and implications for the computation of the network's pathways and attractors.

Our approach is based on the representation of a Boolean network as the feedback interconnection of two separate Boolean modules. In fact, in physical and biological systems, a larger network is often constructed by putting together smaller, already known modules. To connect the modules, entry nodes and outgoing variables are needed, so we will start by introducing input/output Boolean modules (a notion which is common in control theory [START_REF] Sontag | Mathematical Control Theory[END_REF], for instance).

An input/output Boolean module Σ X = (Ω X , U X , H X , f X , h X ), with n X state variables, p X input variables, and q X output variables is characterized by:

• state space, Ω X = {0, 1} n X , • input space U X = {0, 1} p X , • output space H X = {0, 1} q X , • set of Boolean rules f X : Ω X × U X → Ω X , • output function h X : Ω X → H X .
For each u ∈ U X , the synchronous successor of state x ∈ Ω X is given by x + = f X (x; u) and the output associated to this state is y = h X (x).

Given two such modules, Σ A and Σ B , we will use the objects above with subscript A or B and the notation a ∈ Ω A , b ∈ Ω B for the states of each module. Assuming without loss of generality that p A = q B and p B = q A , a larger network can be formed by connecting the input of one system to the output of the other system, i.e., u A = h B (b) and u B = h A (a). The feedback interconnection of the two modules is defined as the system without outputs or inputs Σ I = (Ω I , F I ) with Ω I = {0, 1} n A +n B , states x = (a, b) and:

F I := Ω I → Ω I , F I (x) = F I (a, b) = (f A (a; h B (b)), f B (b; h A (a))), (1)
in other words, the synchronous sucessor of

x = (a, b) is x + = (a + , b + ) with a + = f A (a; h B (b)) and b + = f B (b; h A (a)
). The dynamical behavior of Boolean model Σ I (either synchronous or asynchronous) can be fully described by a state transition graph, G I , which can in turn be represented by a 2 n I × 2 n I matrix. The asymptotic behavior of the network is given by the attractors of G I , which are sets of states from which trajectories cannot escape.

In the case of large networks (eg., n I > 20), the exact computation of the graph G I becomes very costly and not feasible in practice. Hence, we propose to characterize the dynamical behavior of Σ I from the behavior of the two individual modules Σ A and Σ B , by constructing a new object called the asymptotic graph, G a . This new graph is based on a much smaller number of vertices and is less costly to compute. To do this we will first compute their state transition graphs denoted by G A,u , one for each input u ∈ U A and G B,v , for each v ∈ U B . These graphs can be constructed either for the asynchronous or synchronous strategies (see Sections 3 and 4, respectively).

Basic definitions and notation

Some notions used through the text are now recalled (see, for instance, [START_REF] Diestel | Graph Theory[END_REF]). A directed graph, G = (V, E) is defined by a set of vertices, V and a set of directed edges, E, where a directed edge is a pair of ordered vertices,

V 1 → V 2 .
A directed graph is strongly connected if any pair of vertices in V is connected by a path (where a path is a sequence of distinct edges

V 1 → V 2 → • • • V i ).
Any directed graph G can be organized into strongly connected components (SCCs), where an SCC is a maximal strongly connected subgraph of G. An SCC can contain either a single vertex or a subset of vertices (it can also contain the full graph). An SCC C can have incoming or outgoing pathways, that is a pathway from (or towards) a vertex V / ∈ C. An SCC without any outgoing pathways is called an attractor. In other words, when a trajectory reaches an attractor it remains inside for all subsequent updating times and, therefore, distinct attractors represent distinct asymptotic behaviors of the system.

The set of attractors of a given graph G is denoted A(G). The reachable set of a vertex a in G, R(a; G), is the set of all vertices ã such that there exists a path connecting a to ã in G, a G ; ã. The reachable set of a subset S ⊂ V in G is the union R(S; G) = ∪ a∈S R(a; G). Conversely, the basin of attraction of A ∈ A(G) in G, B(A; G), is the set of all vertices a such that there exists a path connecting it to A in G, a G ; A. For a system Σ X with input ν, the state transition graph is a directed graph, denoted by G X,ν = (Ω X , E X ) where the set of edges is computed from the Boolean rules f X and the updating strategy. To facilitate notation, ν takes values in lexicographic order, i.e., ν ∈ {1, 2, . . . , 2 p X } where "1 = 0 . . . 00","2 = 0 . . . 01", and "2 p X = 11 • • • 1".

For each module Σ X , the state transition graphs G X,ν will have L X,ν ≥ 1 attractors denoted X i ν . To avoid introducing heavy notation, we will use a ∈ A i u (or b ∈ B j v ) to name the elements of G A,u (or G B,v ) contained in each attractor.

Table 1 summarizes the notation used for the state transition graphs and attractors for different modules and updating strategies.

A working example

To illustrate our methodology and the differences generated by the three updating strategies, consider the interconnection between two modules, one negative loop and one double-negative loop:

Σ A : a + 1 = ¬a 2 , a + 2 = ¬u ∧ ¬a 1 , h A (a) = (a 1 , a 2 ) , Σ B : b + 1 = v 1 ∧ ¬b 2 , b + 2 = ¬v 2 ∧ b 1 , h B (b) = b 2 , (2)
Table 1: Summary of notation used throughout the paper to name state transition graphs. Three updating strategies will be compared for the analysis of the interconnection of two Boolean modules. These yield different transition graphs for each module and for the full network and asymptotic graphs. In the mixed strategy, the modules are synchronously updated.

In general, a subscript * will be used for asynchronous updates, mx will be used for a mixed strategy, and no subscript for synchronous updates. To avoid an overcharged notation, the semi-attractors are similarly labeled for both asynchronous and synchronous updates, but this is usually clear from the context.

Updates Modules

Full network Asymptotic graph graphs semi-attractors graph attractors graph attractors

Synchronous G A,u ,G B,v A i uα ,B j vβ G I Q k G a Q a,k Mixed G A,u ,G B,v A i uα ,B j vβ G I mx Q k mx G a mx Q a,k mx Asynchronous G A,u * ,G B,v * A i uα ,B j vβ G I * Q k * G a * Q a,k *
where

N A = N B = 2, p A = q B = 2, p B = q A = 1. The interconnected system Σ I is obtained by setting u = b 2 , v 1 = a 1 , v 2 = a 2 . (3)
The state transition graphs under an asynchronous or a synchronous strategy are shown in Fig. 2. The graphs G a or G I will be constructed throughout the following sections. 3 The asynchronous asymptotic graph method for attractor computation

Once the attractors of all the transition graphs are known, a new object will be constructed, called the asymptotic graph, G a * , a graph where the vertices are basically all products of attractors of the form A i u × B j v and the edges are computed based on the transition graphs G A,u

The main results in this paper establish that there is a correspondence between A(G a * ) and A(G I * ) which depends on the updating strategy. The asymptotic graph for asynchronous updates was first proposed in [START_REF] Chaves | Predicting the asymptotic dynamics of large biological networks by interconnections of Boolean modules[END_REF] and then further studied in [START_REF] Tournier | Interconnection of asynchronous Boolean networks, asymptotic and transient dynamics[END_REF][START_REF] Chaves | Attractor computation using interconnected boolean networks: testing growth rate models in e. coli[END_REF]. This Section briefly reviews the main results.

Step-by-step construction

The asynchronous strategy consists in changing exactly one coordinate at each time step. More precisely, fix u, for each state compute x + = f X (x; u) and define the set of indices whose coordinates change:

I(x; u) = {i : x + i -x i = 0}.
Then, the possible asynchronous successors of x in G X,u * are:

y(x; u) ∈ {x : xi = 1 -x i , for some i ∈ I(x; u), xj = x j , j = i}.
Step 1: Construct the set of vertices V * The first step of the asympotic graph method is to divide each attractor of the modules into subsets of states that have the same output. These are called semi-attractors of the modules and are defined as follows. Given A ī u ∈ A(G A,u * ) and

B  v ∈ B(G B,v * )
, define:

A i uα = {a ∈ A ī u : h A (a) = bool(α)} B j vβ = {b ∈ B  v : h B (b) = bool(β)},
where the superscripts i, j range over the total number of semi-attractors for the module (over all inputs), i ∈ {1, . . . , L A * }, j ∈ {1, . . . , L B * }. (Note that i and ī are not related but we have u = u(i) and α = α(i); we drop these dependences to simplify notation. A similar observation holds for j.) The subscripts range over v, α ∈ {1, . . . , 2 p B } and u, β ∈ {1, . . . , 2 p A }, and bool(ν) denotes the Boolean representation of integer ν. Observe that some of the semi-attractors may be empty and are neglected. As an example, the counting of semi-attractors for system (2) is given in [START_REF] Chaves | Attractor computation using interconnected boolean networks: testing growth rate models in e. coli[END_REF] below.

The vertices of G a * are the cross-products of non-empty semi-attractors:

V * = {A i uα × B j vβ : i = 1, . . . , L A * , j = 1, . . . , L B * }, (4)
where L A * and L B * are the total number of semi-attractors for modules A and B (respectively), under asynchronous updates (and the dependences u = u(i),α = α(i),v = v(j),β = β(j) are omitted for simplicity).

Step 2: Construct the set of edges E * To satisfy an asynchronous updating strategy, each edge corresponds to the update of a single module, i.e., only the a (or b) variables are allowed to change at each edge. Accordingly, only edges of the following form are considered

A i uα × B j vβ → A i uα × B  α β or A i uα × B j vβ → A ĩ β α × B j vβ (5)
such that one of the modules is fixed at a semi-attractor, while the dynamics of the other module evolves: in the first case, the semi-attractor of G A,u * is fixed while the dynamics of Σ B follows the path(s) of B j vβ along the graph with input v = α, G B,α * . Therefore, there is an edge between the two vertices if there is a path in the graph G B,α * that leads from some state in B j vβ to some state in B  α β . Similarly for an edge A i uα × B j vβ → A ĩ β α × B j vβ . The set of edges can be written:

E * = {(V, W ) ∈ V * : V = A i uα × B j vβ and (W = A i uα × B  α β or W = A ĩ β α × B j vβ ) },
where

A ĩ β α is in the reachable set of A i uα in the graph G A,β
Therefore, all the attractors of the asynchronous G I * can be obtained from the attractors of G a * . However, it may happen that some of the attractors of G a * are not true attractors of G I * , also called "spurious" attractors. The formation of such attractors is illustrated with the mixed updating strategy (Section 5).

A(G A,1 * ) : A 1 12 = {01}, A 2 13 = {10} A(G A,2 * ) : A 3 23 = {10} A(G B,1 * ) : B 1 11 = {00} (6)
A(G B,2 * ) :

B 2 21 = {00} A(G B,3 * ) : B 3 31 = {00, 10}, B 4 32 = {01, 11} A(G B,4 * ) : B 5 41 = {10}.
All semi-attractors are also true attractors, except for B 3 31 and B 4 32 each containing the states with same output of the attractor {00, 01, 10, 11} of G B,3 * . There are thus 3 × 5 = 15 vertices in V * , corresponding to all possible products of semi-attractors A i uα and B j vβ .

To compute the edges according to (5) consider the vertex A 3 23 × B 5 41 , for instance, which has three outgoing edges. First, fix the part A 3 23 , which is "forcing" system Σ B with input 3; hence we need to compute the successor of B 5 41 = {10} in the state transition graph G B,3 * . This graph contains a cyclic attractor, divided into two semi-attractors (see Fig. 2(a)) so trajectories of system Σ B cycle between these two, B 3 31 and B 4 33 . This gives the two edges

A 3 23 × B 5 41 → A 3 23 × B 3 31 and A 3 23 × B 5 41 → A 3 23 × B 4 32 .
Second, fix the part B 5 41 which is forcing system Σ A with input 1; the only successor of state

A 3 23 = {10} in G A,1 * is A 2 13 = {10}.vThis gives an edge A 3 23 × B 5 41 → A 2 13 × B 5 41 .
In this way, obtain the asymptotic graph shown in Fig. 3(b), which has two attractors:

Q a,1 * = {A 1 12 × B 2 21 } = {0100} Q a,2 * = {A 2 13 × B 3 31 , A 2 13 × B 4 32 , A 3 23 × B 3 31 , A 3 23 × B 4 32 } = {1000, 1010, 1001, 1011}.
Comparison with G I * (Fig. 3(a)) shows that the asymptotic graph exactly recovers the two attractors of the full interconnected system, as expected.

Extending the asynchronous asymptotic graph, G ext *

In the case of the asynchronous asymptotic graph, to improve its construction and avoid the generation of spurious attractors, the idea is to extend the vertex set V * to include further components. This new set will be called V * and the new asymptotic graph is denoted G ext * . An analysis of asynchronous asymptotic graphs with failures shows that spurious attractors may appear due to output changes along the pathways followed by the states of an attractor A i u in a state transition graph G A,ũ * with ũ = u (see the Example in Section 5). These pathways and the corresponding sequence of outputs "forgotten" in the construction of G a * (which keeps only the final state of the pathway, by definition).

To identify the appropriate new components, the idea is to analyze the pathways of each state of each attractor X i ν along the other partial graphs G X,ν * , with ν = ν, until a change in the output is detected, as follows (for clarity, the case of module X = A will be described):

(i) Compute the strongly connected components (SCCs) of each state transition graph G A,u * and its terminal SCCs (or attractors). Collect its attractors in A u = {A i u : i = 1, . . . , L A,u }.

(ii) Pick an attractor A i u , suppose it contains the states {r 1 , . . . , r k }, and look at all the possible forward pathways for each r in the other graphs G A,ũ * with ũ = u, for instance,

P 1 : r → s 1 → s 2 → • • • , in G A,ũ * ,
Along each path, each state has its own output h A (s i ).

(iii) For each pathway P k of each state r , pick the first state s î that has a different output from r , that is

h A (s î) = α = h A (r ) = α.
(iv) This state s î belongs to some strongly connected component S u of the graph G A,u * . Add this SCC to the set of attractors:

A u = A u ∪ {S u }.
Repeat steps (i)-(iv) for all u ∈ {1, . . . , 2 p A } and then for module B, for all v ∈ {1, . . . , 2 p B }, to produce the extended sets A u and B v . The next step is to divide each component of these sets into semi-SCCs with the same output: Āi

uα = {a ∈ Āī u : h A (a) = bool(α)} and Bj vβ = {b ∈ B v : h B (b) = bool(β)}.
Collecting all the semi-SCCs into:

S A = { Āi uα : i = 1, . . . , LA }, S B = { Bj vβ : j = 1, . . . , LB },
the set of vertices for the extended asynchronous asymptotic graph is thus:

V * = { Āi uα × Bj vβ : i = 1, . . . , LA , j = 1, . . . , LB }.
The new set of edges Ē * is computed as before [START_REF] Chaves | Robustness and fragility of boolean models for genetic regulatory networks[END_REF].

To test this method, we applied it to 750 randomly generated pairs of modules. The modules have N A = N B ∈ {2, 3, 4, 5, 6} and inputs p A = p B ∈ {1, 2, 3}. Each module was obtained by randomly generating the Boolean truth tables for G X,ν * , ν ∈ {1, . . . , 2 p X }. The results obtained with this extended method (Fig. 4) indicate that the spurious attractor problem is solved for about 99.8% of the cases. However, the new G ext * is still not exact and it is much larger to calculate, since |G ext * | > |G I * |, except for modules with dimension above N X = 6. These results also suggest that the addition of other SCCs to the sets S A and S B will further improve the asymptotic graph. For instance, by successively adding the SCCs corresponding to a second or third change of output along the pathways P k .

Indeed, this joins the result shown in [START_REF] Tournier | Interconnection of asynchronous Boolean networks, asymptotic and transient dynamics[END_REF] for another graph construction, called the crossgraph, G cr : this graph is contructed from all pairs of all SCCs of G X,ν * . A theorem shows that there is a bijection between the attractors of G cr and G I * (i.e. G cr recovers exactly all attractors of G I * ). This is not surprising, since G cr virtually contains all states of G I * . The difficulty is that G cr (and possibly other extensions of G ext * ) is as costly to compute as the graph G I * itself, so it is not an efficient method for large networks. Nevertheless, the successive constructions G a * , G ext * ,..., G cr provide a clear illustration of the organization of asynchronous dynamical behavior.

The synchronous asymptotic graph, G a

To construct the asymptotic graph in the synchronous case, one could follow a similar reasoning and define the edges of the asymptotic graph by asking that both the a and b states change simultaneously:

(A i uα , B j vβ ) → (A ĩ β α, B  α β ).
However, a few examples showed that this strategy may not only generate spurious attractors (as in the asynchronous case), but also miss finding some of the attractors of the synchronous G I . This is therefore not a favorable strategy for constructing a synchronous asymptotic graph.

To improve the construction and guarantee some result, the asymptotic graph can be constructed according to an iterative procedure, as follows.

Step 1: Construct the set of vertices W The initial set of vertices is V, constructed as in the asynchronous case [START_REF] Baldazzi | The carbon assimilation network in Escherichia coli is densely connected and largely signdetermined by directions of metabolic fluxes[END_REF]. Then, given a set of nodes W k , compute their direct successors according to (1) to form W k+1 :

W 0 := {A i uα × B j vβ : i = 1, . . . , L A , j = 1, . . . , L B }, W k+1 := {W : W = (a + , b + ), for (a, b) ∈ V, all V ∈ W k } \ (∪ k =0 W ) (7)
where L A , L B are the total number of semi-attractors of modules A and B (respectively) under synchronous updates, and

a + = f A (a; h B (b)) and b + = f B (b; h A (a)
). (Since the state space is finite there will eventually be some k such that W k = ∅.) The set of nodes for the synchronous asymptotic graph G a is then:

W = ∪ k-1 k=0 W k . ( 8
)
Step 2: Construct the set of edges E The set of edges is based on synchronous updates:

E = {(V, W ) ∈ W : V = (a, b) and W = (a + , b + ), V, W ∈ W}. (9)
As a last remark, in theory, the graph G a may be as long to compute as the full graph G I . This is however not the case in practice (see examples below).

Step 3: Attractor computation In the synchronous case, since G a is a subgraph of G I , all attractors of G a are guaranteed to be true attractors of G I . It may, however, happen that not all attractors of the full network are recovered, in contrast to the asynchronous case.

Theorem 2 In the case of synchronous networks, if Q a is an attractor of G a , then Q a contains an attractor of G I .

Proof: Let Q a be an attractor of G a and let B(Q a ) represent its basin of attraction. Then, by definition of the synchronous asymptotic graph, there is at least one element of Ω I of the form V = (A i uα , B j vβ ) ∈ B(Q a ) (since the set of vertices of G a contains W 0 , ( 7)). Again, by definition of the synchronous asymptotic graph, all reachable states from V in G I , R(V ; G I ), are also contained in W and, namely, Q a contains the attractor corresponding to the path in R(V ; G I ). In other words, Q a contains an attractor of G I .

Example For the network (2), the synchronous transition graphs G A,u and G B,u are shown in Fig. 2(b) and the sets of semi-attractors are:

A(G A,1 ) : A 1 11 = {00}, A 2 14 = {11}, A 3 12 = {01}, A 4 13 = {10} A(G A,2 ) : A 5 23 = {10} A(G B,1 ) : B 1 11 = {00} (10)
A(G B,2 ) : All semi-attractors are also true attractors, except for two pairs A 1 11 ∪ A 2 14 = {00, 11} and B 3 31 ∪B 4 32 = {00, 01, 10, 11} organized into same-output semi-attractors. There are thus 5×5 = 25 vertices in V. According to the iterative procedure [START_REF] Chaves | Hierarchy of models: From qualitative to quantitative analysis of circadian rhythms in cyanobacteria[END_REF], the following four vertices need to be added: 00 × 00, 00 × B 2 21 , A 2 14 × 01, A 5 23 × 01. Fig. 5(a), represents depicts the synchronous asymptotic graph which has three attractors:

B 2 21 = {00} A(G B,3 ) : B 3 31 = {00, 10}, B 4 32 = {01, 11} A(G B,4 ) : B 5 41 = {10}.
Q a,1 = A 3 12 × B 2 21 = {0100} Q a,2 = {A 1 11 × B 5 41 , A 2 14 × 01, A 2 14 × B 1 11 , 00 × 00} = {0000, 0010, 1100, 1101}, Q a,3 = {A 4 13 × B 4 32 , A 5 23 × B 3 31 } = {1000, 1010, 1001, 1011}.
Comparison with G I (Fig. 5(b)) shows that the synchronous asymptotic graph exactly recovers the attractors of the full interconnected system.

5 The asymptotic graph, for mixed strategy networks, G a mx

The results in Sections 3 and 4 suggest a natural application of the asymptotic graph method for mixed-strategy networks. Assume that a network Σ M with graph G I mx is synchronous in blocks, for instance composed of two groups of variables, x = (a, b), and within each group a 1 , . . . , a k or b 1 , . . . , b the variables are synchronously updated. The choice to update either group a or group b is according to an asynchronous strategy.

To apply the asymptotic graph method, it is intuitive to partition the network into two modules, each containing one of the synchronous blocks of variables. For each module compute the synchronous graphs G X,u and their corresponding attractors A(G X,u ). Using these attractors, construct the asymptotic graph as in the asynchronous case, G a mx = (V mx , E mx ). Then, the attractors of G a mx will contain (a representative of) all attractors of G I mx : A(G a mx ) ⊃ A(G I mx ). Intuitively, this is because the attractors of G X,u represent each basin of attraction of each of the synchronous blocks.

To facilitate the analysis, we now introduce a projection function, π : A(G X,u ) → Ω X which returns all the states contained in a given attractor of G X,u and a grouping function, ψ : A(G I mx ) → V mx , which returns all the vertices of G a mx "contained" in an attractor of G I mx :

ψ(Q) = {(A uα , B v,β ) ∈ V mx : (a, b) ∈ Q where a ∈ π(A uα ), b ∈ π(B v,β )}.
The equivalent of Theorem 1 for a mixed block-synchronous strategy can now be stated. The proof follows the same argument as in [START_REF] Chaves | Predicting the asymptotic dynamics of large biological networks by interconnections of Boolean modules[END_REF]. Example Suppose that the two modules Σ A and Σ B given in (2) are both synchronously updated, but their interconnection is asynchronous. To construct the asymptotic graph G a mx , use the synchronous modules transition graphs G X,u shown in Fig. 2(b) and then construct the asymptotic graph using an asynchronous strategy. In this case, the sets of semi-attractors are the same as in the synchronous case [START_REF] Comet | Simplified models for the mammalian circadian clock[END_REF], so that V has 5 × 5 = 25 vertices in G a mx . Fig. 6(a), The mixed strategy asymptotic graph has three attractors: As before, this result can be extended to the interconnection between 3 or more synchronous blocks. Furthermore, if the strategy among blocks is asynchronous, it is clear that each block could follow any updating strategy, since only its state transition graphs G X,u , u ∈ {1, 2, . . . , 2 p X }, are needed to construct the asymptotic graph.

Q a,1 mx = {A 3 12 × B 2 21 } = {0100}, Q a,2 mx = {A 1 11 × B 1 11 , A

The synchronous vs. asynchronous dynamics dichotomy

The different dynamics induced by synchronous or asynchronous updating strategies is well illustrated through the construction of the asymptotic graphs.

In practice, we can compute all or most of the attractors of the interconnected network using only the knowlegde of the two individual modules, through the computation of a new graph of (much) reduced dimension, the asymptotic graph. Each strategy -either updating all variables simultaneouly or only one at a time -leads to a different construction for the asymptotic graph and to two apparently opposite results:

A(G I ) ⊇ A(G a ), A(G I mx ) ⊆ A(G a mx ), A(G I * ) ⊆ A(G a * )
. In asynchronously interconnected networks, it appears that the information on the asymptotic behavior is completely included in the state transition graphs of the two separate modules, and this for any possible pair of modules. The asymptotic graph scheme may sometimes introduce a false attractor. Roughly speaking, the asynchronous nature of the dynamics is reminiscent of a slow/fast timescale separation, since each part of an attractor (either the A or B part) is fixed while the other part evolves rapidly inside its own module: the dynamics within each module is "fast" when compared to the timescale of the interconnection. Thus, the asynchronous asymptotic graph contains one representative of each family of pathways leading to any given attractor. In some particular cases, the families of pathways are not adequately represented in G a * , thereby leading to the appearance of a spurious attractor as illustrated by the example in Section 5.

In contrast, in synchronously interconnected networks, by construction the asymptotic graph contains actual pathways of G I (those starting from an initial choice of vertices, W 0 ), so it cannot introduce any spurious attractors. In synchronous networks, each state has a unique successor and converges to a single attractor. The basins of attraction are pairwise disjoint and deterministically associated with a unique attractor. So, in general, the attractors associated to a pair of modules Σ A , Σ B for Σ I may not sample appropriately all the basins of attraction, hence some attractors of G I will be "forgotten" in G a . This is the case even in straigthforward networks such as the positive loop of size n ≥ 4.

This difference in the structure of asynchronous and synchronous networks joins the work of Demongeot and colleagues [START_REF] Demongeot | Attraction basins as gauges of robustness against boundary conditions in biological complex systems[END_REF] who have analysed the basins of attraction of synchronous networks, their sizes and the distances between them, as indicators of robustness. A small distance between two basins of attraction means that perturbations to an initial configuration are more likely to induce a switch from one basin to another, and hence to change the asymptotic behavior of the network.

The main disadvantage of the asymptotic graph construction is the possibility of generating spurious attractors in the case of asynchronous updates, or the failure to recover all attractors under synchronous updates. However, our tests with randomly generated networks show that these errors are very small and decrease with size for networks with some structure (below 2%, see Figure 7) and applications to biological networks modules has never shown any of these two "errors".

A quick comparison of the three stategies is shown in Figure 7. Two classes of Boolean functions were considered: unate cascade functions as characterized in [START_REF] Jarrah | Nested canalyzing, unate cascade, and polynomial functions[END_REF], and Horn functions as characterized in [START_REF] Crama | Boolean functions. Theory, algorithms, and applications[END_REF], for n variables: Unate cascades:

f (x) = xσ(1) 1 xσ(2) 2 . . . xσ(n-1) n-1 xσ(n) . . . ,
Horn functions:

f (x) = m i=1   j∈P i x j k∈Q i ¬x k  
where σ is a permutation of {1, . . . , n}, x is either x or ¬x, and i represents either the logical AND or the logical OR; for Horn functions 1 ≤ m ≤ 6, 0 ≤ |P i | ≤ n and 0 ≤ |Q i | ≤ 1, that is, each term in a Horn function has at most one negated variable. To generate Horn functions, the integers m, |P i |, and |Q i | were randomly chosen out of a uniform distribution in the given intervals. To obtain unate cascade functions, a permutation and a sequence of n -1 logical operations were randomly generated from uniform distributions. For the computation of each of the three asymptotic graphs for each class, A total of 500 × 5 networks were randomly generated as pairs of modules with N A = N B ∈ {3, 4, 5, 6, 7} and N in = 2, for each class (each symbol corresponds to an average over 500 networks). For each of these networks, we computed the three asymptotic graphs and the full interconnected graph and compared the attractors. It is immediate to observe that the asynchronous strategy gives the most accurate results, with the lowest percentage of errors.

For the analysis of large networks, it is clear that synchronous updates have a lower computational cost. Thus, in view of Section 5, the most advantageous solution is to use a mixed strategy for analysis of interconnections: consider synchronous individual modules and then construct an asynchronous asymptotic graph. In cases where each individual module has its own time-scale and their communication happens at a slower pace, mixed strategies are a good option.

Composition of biological networks by feedback interconnection

A recurrent problem in biology is to study the emergent dynamics arising from the coupling of two or more different phenomena [START_REF] Gerard | Temporal self-organization of the cyclin/Cdk network driving the mammalian cell cycle[END_REF][START_REF] Baldazzi | The carbon assimilation network in Escherichia coli is densely connected and largely signdetermined by directions of metabolic fluxes[END_REF]. Although many biological phenomena are currently represented by mathematical models, there are no specific tools to systematically explore the interconnections between two existing models. Our asymptotic graph method addresses this problem by providing a framework to analyse the joint dynamics of two known networks. Given two or more previously developed models, our method requires only an investigation of the input-output connections between the two networks (in Fig. 1, finding u = h B (b) and v = h A (a)).

To illustrate our method, four pairs of Boolean networks that model distinct phenomena in several organisms have been coupled to analyze the corresponding joint dynamics. These pairs are listed in Table 2 by order of size/complexity. (For a complete description and construction of the models, please refer to the bibliography cited in the Table .)

Following the discussion on Section 6, we will consider synchronous modules and then compare the synchronous and mixed asymptotic graphs. The notation used throughout these examples is the same as before, using the lexicographical order for Boolean numbers: A 3 25 means that this semi-attractor belongs to the graph G A,2 , it has input "2" (Boolean 001) and output "5" (Boolean 100). The superscript "3" means that it is semi-attractor number 3 in the total count of semi-attractors of module A.

The cyanobacteria example

Circadian rhythms in cyanobacteria are widely studied and the corresponding mechanism is very well known: it is characterized by an ordered sequence of phosphorylations in the protein Kai C under the regulation of another protein, Kai A. This mechanism can be faithfully reproduced in vitro (see [START_REF] Rust | Ordered phosphorylation governs oscillation of a three-protein circadian clock[END_REF] and references therein). In living cyanobacteria, the clock mechanism interacts with the transcription/translation cycle (TTC) [START_REF] Dong | Elevated atpase activity of kaic applies a circadian checkpoint on cell division in synechococcus elongatus[END_REF].

The two cyanobacteria modules The Boolean clock mechanism (Σ A ) was developed in [START_REF] Chaves | Hierarchy of models: From qualitative to quantitative analysis of circadian rhythms in cyanobacteria[END_REF] and a simple TTC circuit is composed by unphosphorylated protein Kaic C (C U ) and its mRNA (R), and a "lumped" protein P closes the circuit:

Σ A : A + 1 = ¬C S ∨ A 2 A + 2 = ¬C S ∧ A 2 C + T = u ∧ A 1 C + T S = A 1 ∧ C T C + S = ¬A 2 ∧ C T S Σ B : R + = ¬v ∧ ¬P C + U = R P + = C U
Based on [START_REF] Dong | Elevated atpase activity of kaic applies a circadian checkpoint on cell division in synechococcus elongatus[END_REF], the interactions are as follows:

u = h B (b) = C U , v = h A (a) = C T ∨ C T S .
Computation of the state transition graphs for each module shows that there are five semiattractors for each module, as follows: where the Kai C protein is sequencially phosphorylated on its site T (10100), then T and S (10110), then only on site S (10111), which inhibits Kai A (00111), leading to repeat the cycle. Similarly, the attractor B 1 11 ∪ B 2 12 represents the full TTC cycle with Kai C mRNA production (100) followed by protein production (110), which induces the full pathway (111) and leads to cycle repetition.

G A,1 : A 1 11 = {10000}, A 2 11 = {11000}, G A,2 : A 3 21 = {00000,
The cyanobacteria synchronous asymptotic graph Computation of the synchronous asymptotic graph leads to the transition graph shown in Fig. 11. This graph contains 52 states (instead of the 2 8 = 256 of the full interconnected system), represented by products of semi-attractors or states a × b (i.e., Boolean states a of Σ A and b of Σ B ) and two cyclic attractors:

Q a,1 : 00000 × 100 → 10000 × 110 → 10100 × 111 → 10110 × 011 → 10111 × 001 → → 00011 × 000 → 00001 × 000 → • • • Q a,2 : 11000 × 100 → 11000 × 110 → 11100 × 111 → 11110 × 011 → 11110 × 001 → → 11010 × 000 → 11000 × 000 → • • • .
The first attractor oscillates simultaneaously through the clock and TTC orbits (except that 00111 is replaced by 00011, an equivalent state at this step), showing that the interconnected system has synchronized clock and transcription/translation oscillations, as expected for wild type cells.

The second attractor still cycles through TTC but the clock evolves among states which contain a high concentration of protein Kai A (A 2 11 and A 5 22 ), thus preventing the formation of C S and subsequent re-initialization of the Kai proteins phosphorylation cycle. Such dynamics can indeed be observed experimentally [START_REF] Rust | Ordered phosphorylation governs oscillation of a three-protein circadian clock[END_REF]: a high expression of Kai A prevents circadian clock oscillations.

The cyanobacteria mixed mode asymptotic graph To compare updating strategies, the mixed mode asymptotic graph is shown in Fig. 12. It is composed of 25 states (only products of semi-attractors) and also contains two attractors: From equations ( 11) and ( 12), observe that attractors Q a,i mx correspond to attractors Q a,i . Merging semi-attractors that belong to the same attractor, Q a,1 mx has the following compact dynamics:

Q a,1 mx = {A 1 11 × B 1 11 , A 1 11 × B 2 12 , A 1 11 × B 5 21 , A 3 21 × B 2 12 , A 3 21 × B 1 11 , A 4 22 × B 2 12 , A 4 22 × B 5 21 } Q a,2 mx = {A 5 22 × B 2 12 , A 5 22 × B 5 21 , A 2 11 × B 2 12 , A 2 11 × B 5 21 , A 2 11 × B 1 11 }.
A 4 22 × B 5 21 ← A 3 2• × B 1 1• ↓ ↑↓ A 1 11 × B 5 21 → A 1 11 × B 1 1•
meaning that there is (as in the synchronous asymptotic graph) a simultaneous clock and TTC oscillation (A 3 2• × B 1 1• ), but also the possibility that the clock oscillator still follows a clear sequence with Kai C successively phosphorylated in its three possible combinations (A 4 22 ) and then all unphosphorylated A 1 11 , while the TTC is arrested (fixed at B 5 21 ). This last case may represent a situation protein production is slower or has a higher period than the clock.

Mammalian example

The interactions between the mammalian cell cycle and circadian clock are currently generating a large amount of interest in biology, as the progress in single cell measurements allows answering specific questions [START_REF] Feillet | Phase locking and multiple oscillating attractors for the coupled mammalian clock and cell cycle[END_REF]. These new experiments are uncovering unexpected bi-directional links between the two modules.

The two mammalian modules

The cell cycle Boolean model (module Σ A ) is taken from [START_REF] Fauré | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF] and has the following variables:

CycD + = CycD Rb + = (¬CycD ∧ ¬CycE ∧ ¬CycA ∧ ¬CycB) ∨ (p27 ∧ ¬CycD ∧ ¬CycB) E2F + = (¬Rb ∧ ¬CycA ∧ ¬CycB) ∨ (p27 ∧ ¬Rb ∧ ¬CycB) CycE + = ¬u (E2F ∧ ¬Rb) (13)
CycA + = (E2F ∧ ¬Rb ∧ ¬Cdc20 ∧ ¬(Cdh1 ∧ U bc)) ∨ (CycA ∧ ¬Rb ∧ ¬Cdc20 ∧ ¬(Cdh1 ∧ U bc)) p27 + = ¬CycD ∧ ¬CycE ∧ ¬CycA ∧ ¬CycB Cdc20 + = CycB Cdh1 + = (¬CycA ∧ ¬CycB) ∨ Cdc20 ∨ (p27 ∧ ¬CycB) U bc + = ¬Cdh1 ∨ (Cdh1 ∧ U bc ∧ (Cdc20 ∨ CycA ∨ CycB))
CycB + = ¬Cdc20 ∧ ¬Cdh1, as described in [START_REF] Fauré | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF], where CycX represent the four cyclins involved in the cell cycle. The clock model (module Σ B ) is based on [START_REF] Comet | Simplified models for the mammalian circadian clock[END_REF]:

BM AL + = ¬P Cnuc mP ER + = ¬v ∧ BM AL mCRY + = ¬v ∧ BM AL pP ER + = mP ER pCRY + = mCRY P C + = pP ER ∧ pCRY P Cnuc + = P C.
In the clock model, mX and pX denotes mRNA and protein corresponding to gene X, while P C denotes the complex formed by the proteins PER and CRY, and P Cnuc denotes this complex in the nucleus.

Testing interactions between modules Although there are currently many questions on the form of interaction between mammalian cell cycle and circadian clock [START_REF] Feillet | Coupling between the circadian clock and cell cycle oscillators: Implication for healthy cells and malignant growth[END_REF], it seems clear that clock's BM AL acts on the cell cycle, possibly at different stages. Conversely, no conclusive biological knowledge exists on how the cell cycle may affect the clock. To include a basic link, we have considered that during mitosis all gene expression is arrested (in this model, mitosis can be modeled as Cdc20 ∧ CycB). The interconnection between modules is thus given by (see Fig. 9):

u = h B (b) = BM AL, v = h A (a) = Cdc20 ∧ CycB.
BM AL acts negatively during G1, S or M-phases of the cell cycle. We analyzed three different interconnection models, one for each of these possibilities, and compared the results (see Table 4). The three cases give similar results in terms of joint behavior. For a closer analysis consider, for instance, BM AL acting during G1 (here represented by a negative efect on cyclin

The mixed mode graph has also five attractors, which correspond to each of Q a,i . In fact, we can see from Fig. 14 that there is set equality Q a,i = Q a,i mx for i = 1, 3, 4, 5 but strictly inequality for Q a,2 ⊂ Q a,2 mx . The attractors Q a,i i = 1, 3, 4, 5 represent quiescent cells with ongoing circadian oscillations (as observed in hair cells [START_REF] Plikus | Local circadian clock gates cell cycle progression of transient amplifying cells during regenerative hair cycling[END_REF]), while Q a,2 represents the joint cell division cycle and clock oscillatory behavior. The mixed mode asymptotic graph shows how this joint behavior emerges from the combination of the module's periodic orbits (see left attractor in Fig. 14). Merging semi-attractors that belong to the same attractor, Q a,2 mx can be represented in the compact form:

A 3 12 × B 9 22 ← A 2 1• × B 1 1• ↓ ↑↓ A 5 2• × B 9 22 A 5 2• × B 1 1•
where

A 2 1• = A 2 11 ∪ A 3 12 , A 5 2• = A 5 21 ∪ A 6 22 , and B 1 1• = B 1 11 ∪ B 2 12 .
The cell cycle and clock may jointly oscillate and alternate states with a regular cycle of cyclin E (present mostly through S phase and mitosis) or eventually switch to a joint cycle with absence of cyclin E. However, at mitosis (A 3 12 ), the clock may switch to its steady state (B 9 22 ), which leads directly to a full degradation of cyclin E in the cell cycle (A 5 2• ).

T helper cells regulatory network

The third example analyzes the regulatory network developed in [START_REF] Mendoza | A method for the generation of standardized qualitative dynamical systems of regulatory networks[END_REF] to describe the differentiation of T helper cells, which are part of the vertebrate immune system, into one of three classes: Th0 precursor cells and Th1, Th2 effector cells. Network reconstruction is discussed in [START_REF] Mendoza | A method for the generation of standardized qualitative dynamical systems of regulatory networks[END_REF], it involves 19 variables and four external input variables, X e ∈ {IF N β, IL12, IL18, T CR}. For the analysis in [START_REF] Mendoza | A method for the generation of standardized qualitative dynamical systems of regulatory networks[END_REF], these four inputs are all set to 0 but here, these external inputs are treated as constant variables, setting X + e = X e to allow for greater generality in the results. Thus, the 23 variable network was first split into two modules, of 10 and 13 variables, respectively:

Σ A : GAT A3 + = (GAT A3 ∨ ST AT 6) ∧ ¬T bet IF N β + = IF N β IF N βR + = IF N β IL4 + = GAT A3 ∧ ¬ST AT 1 IL4R + = IL4 ∧ ¬SOCS1 JAK1 + = u ∧ ¬SOCS1 SOCS1 + = ST AT 1 ∨ T bet ST AT 1 + = IF N βR ∨ JAK1 ST AT 6 + = IL4R T bet + = (ST AT 1 ∨ T bet) ∧ ¬GAT A3 Σ B : IF N γ + = (IRAK ∨ N F AT ∨ ST AT 4 ∨ v 3 ) ∧ ¬ST AT 3 IF N γR + = IF N γ IL10 + = v 1 IL10R + = IL10 IL12 + = IL12 IL12R + = IL12 IL18 + = IL18 IL18R + = IL18 ∧ ¬v 2 IRAK + = IL18R N F AT + = T CR ST AT 3 + = IL10R ST AT 4 + = IL12R ∧ ¬v 1 T CR + = T CR
where the interconnection between modules is given by:

u = h B (b) = IF N γR, v = h A (a) = (GAT A3, ST AT 6, T bet) .
Module A is single input and has the following 12 semi-attractors: Module B has 3 inputs, and therefore eight fixed input graphs G B,v , with altogether 64 single state attractors (not listed here).

G A,1 : A 1 11 = {0000000000}, A 2 12 = {0000001001}, A 3 12 = {0110001101}, A 4 17 = {1001100010}, A 5 15 = {1110001100}, G A,2 : A 6 22 = {0000001001}, A 7 22 = {0110001101}, A 8
For the interconnection of the two systems, the synchronous and mixed asymptotic graphs give similar results, both with 33 single state attractors, which can be divided into the following five groups: where the symbol * indicates that the corresponding digit may equal either 0 or 1. In the particular case where all the external input variables are set to 0, X e ≡ 0, there are only 3 attractors for the full T helper cell regulatory network, as expected. Indeed, this was the only case considered in [START_REF] Mendoza | A method for the generation of standardized qualitative dynamical systems of regulatory networks[END_REF], which lists the three attractors corresponding to the three T helper cell phenotypes, Th0, Th1, and Th2. Our analysis was more general, since the four external input variables (X e ) were allowed to take all possible combinations of 0,1 values, leading to a total of 33 distinct attractors for the T helper cell network. However, notice that these attractors follow some specific patterns, and five distinct groups can be identified, according to the cytokines expressed or not in each attractor. Furthermore, each group of attractors can then be classified as a representative of a distinct T helper cell phenotype, as a function of 10 cytokines and proteins that typically discriminate between the phenotypes.

A 1 11 × B 1 11 = {0000000000} × {0000000000000}, A 4
As summarized in [START_REF] Mendoza | A method for the generation of standardized qualitative dynamical systems of regulatory networks[END_REF], precursor T helper cells, or Th0, do not produce any of the cytokines represented in this model. The interconnected network (Σ A , Σ B ) describes the differentiation into the two types of effector cells. Th1 type cells are characterized by high expression of cytokine IF N γ, protein SOCS1, and transcription factor T bet. In contrast, the dominant characteristic of Th2 cells is a high expression of cytokine IL4 and its receptor IL4R. It is known that SOCS1 inhibits the IL4 pathway, hence SOCS1 is hardly present in Th2 cells. Moreover, it appears that GATA3 forms a positive feedback loop with IL4, hence it is also expressed in Th2 cells. Another cytokine, IL10, is also observed in Th2 cells: its expression is enhanced by IL4 and depends on the presence of ST AT 3.

A 4 17 × A 5 15 × A 6 22 × A 7 22 × {0011 • • • 00 • 10•} {0011 • • • • • •10•} {1100 • • • • • •0 • •} {1100 • • • • • •0 • •} IF N γ 0 0 0 1 1 IF N γR 0 0 0 1 1 IL10 0 1 1 0 0 IL10R 0 1 1 0 0 IL4 0 1 0 0 0 IL4R 0 1 0 0 0 ST AT 3 0 1 1 0 0 GAT A3 0 1 1 0 0 SOCS1 0 0 1 1 1 T bet 0 0 0 1 1 IF N β 0 0 1 0 1 Cell type Th0 Th2 Th1p Th1 Th1
Indeed, all the cytokines and proteins that distinguish between effector cell types are constant in each of the five attractor families, as seen in Table 3. Concerning the external inputs, only IF N β has constant expression in each group.

It is clear that A 1 11 ×B 1 11 represents precursor cells, Th0, with no expressed variables. There are two attractor groups (columns five and six) that correspond to Th1 phenotype, with all the required proteins expressed: IF N γ, IF N γ, SOCS1, and T bet. There is one group (column three) which exactly corresponds to phenotype Th2, as it expresses both IL4, IL10 and their receptors, as well as GAT A3 and ST AT 3.

A fifth attractor group (column four) is identified, representing an intermediate type of cell which we classified as "pre-T h1" cell phenotype. While this group expresses most of Th2 characteristic proteins, it also expresses the Th1 protein SOCS1, which inhibits the dominant protein in Th2, IL4. Therefore, this group of attractors cannot be said to represent Th2. One of the differences between Th2 and Th1p is in the fact that Th1p is stimulated by external input IF N β. Experimental evidence [START_REF] Baccala | Interferons as pathogenic effectors in autoimmunity[END_REF] shows that IF N β promotes Th1 cell maturation, by increasing IF N γ production.

The T helper cell network does reproduce this fact, as shown in Fig. 10: since both GAT A3 and T bet have self-loops, it suffices to apply a perturbation (to T bet for instance) to invert the expression levels of GAT A3 and T bet, which will in turn lead to the expression of the other proteins in a Th1 phenotype. Conversely, it can also be remarked that suppression of external input IF N β will lead to loss of SOCS1 expression and subsequent de-repression of IL4R, which would determine a Th2 phenotype.

Arabidopsis thaliana example

In the last two years, two fundamental Arabidopsis thaliana signaling networks have been developed using Boolean models, one comprehensively representing the cell cycle [START_REF] Ortiz-Gutiérrez | A dynamic gene regulatory network model that recovers the cyclic behavior of Arabidopsis thaliana cell cycle[END_REF] and another studying the formation of the root apical meristem (RAM) through a genetic-hormonal regulatory network (GHRN) [START_REF] García-Gomez | A dynamic genetic-hormonal regulatory network model explains multiple cellular behaviors of the root apical meristem of Arabidopsis thaliana[END_REF]. Both have about 15 Boolean variables (we will not reproduce the Boolean networks here, as they are easily accessible through these references) the first has only one attractor while the later has 17. Several points of possible interactions between these two networks are discussed in [START_REF] García-Gomez | A dynamic genetic-hormonal regulatory network model explains multiple cellular behaviors of the root apical meristem of Arabidopsis thaliana[END_REF], notably two main actions of components of GHRN on the cell cycle: the first through auxin which modulates the transition from mitotic to endocycle in Arabidopsis, and the second through cytokinin which promotes the expression of an APC activator in the root meristem. In the GHRN1 model, auxin is represented by WOX5 (transcription factor Wuschel-related-homeobox-5) and cytokinin by ARR2 (Arabidopsis Response Regulator 2) . The action of the cell cycle on the GHRN was again represented by arresting all gene expression during mitosis, here modeled by AP C ∧ CycB.

The interconnection between modules is thus given by:

u = h B (b) = [ARR2, W OX5] , v = h A (a) = AP C ∧ CycB
where u inhibits expression of the five genes in GHRN1: SHORT-HYPOCOTYL2 (SHY2), SHORTROOT (SHR), SCARECROW (SCR), JACKDAW (JKD), and MAGPIE (MGP). Conversely, v 1 positively affects APC and v 2 negatively affects cyclin D. As expected, the transition graphs G A,1 and G B,1 have exactly the same attractors obtained in [START_REF] Ortiz-Gutiérrez | A dynamic gene regulatory network model that recovers the cyclic behavior of Arabidopsis thaliana cell cycle[END_REF][START_REF] García-Gomez | A dynamic genetic-hormonal regulatory network model explains multiple cellular behaviors of the root apical meristem of Arabidopsis thaliana[END_REF]. The mixed mode asymptotic graph G a mx has seven attractors, of the form:

Q a,k mx = A 1 11 × B k 11 ← A 1 11 × B 21 ↑↓ ↑↓ A 2 12 × B k 11 → A 2 12 × B 21 , Q a,QC mx = A 3 21 × B 6 12 (14)
where A 1 11 ∪ A 2 12 is the cell cycle attractor [START_REF] Ortiz-Gutiérrez | A dynamic gene regulatory network model that recovers the cyclic behavior of Arabidopsis thaliana cell cycle[END_REF] and A 3 21 is an attractor representing quiescent cells, and satisfying E2F a = 1, SCF = CycB = CycA = 0. The attractors B vβ represent some well known hormonal activity profiles, as listed in [START_REF] García-Gomez | A dynamic genetic-hormonal regulatory network model explains multiple cellular behaviors of the root apical meristem of Arabidopsis thaliana[END_REF] The mixed mode asymptotic graph thus predicts that the interconnected cell cycle / GHRN1 network will have seven possible types of asymptotic behavior, characterized by [START_REF] Devloo | Identification of all steady states in large networks by logical analysis[END_REF], where the pairs (k, ) ∈ {(13, 20), [START_REF] Devloo | Identification of all steady states in large networks by logical analysis[END_REF][START_REF] Garg | Synchronous versus asynchronous modeling of gene regulatory networks[END_REF], [START_REF] Diestel | Graph Theory[END_REF][START_REF] Gerard | Temporal self-organization of the cyclin/Cdk network driving the mammalian cell cycle[END_REF], [START_REF] Dong | Elevated atpase activity of kaic applies a circadian checkpoint on cell division in synechococcus elongatus[END_REF][START_REF] Irons | Improving the efficiency of attractor cycle identification in Boolean networks[END_REF], [START_REF] Fauré | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF][START_REF] Jarrah | Nested canalyzing, unate cascade, and polynomial functions[END_REF], [START_REF] Feillet | Phase locking and multiple oscillating attractors for the coupled mammalian clock and cell cycle[END_REF][START_REF] Hong | An efficient steady-state analysis method for large boolean networks with high maximum node connectivity[END_REF]}. In summary, there is a quiescent cell steady state A 3 21 ×B 6 12 , where cell cycle is arrested. All other dynamical behaviors consist of a normal cell cycle together with a hormonal profile coresponding either to Root Cap formation or to central pro-vascular tissues, in the transition domain (TD) of the root apical meristem. The synchronous asymptotic graph obtains all attractors in [START_REF] Devloo | Identification of all steady states in large networks by logical analysis[END_REF] as well as some variations around the quiescent cell state Q a,QC mx (not shown). To explore whether other dynamical behaviors may be obtained, a deep study of the interactions between the cell cycle and GHRN is needed, but this is out of the scope of our paper. Nevertheless, it stands out that our method constitutes a very practical tool to attack this question.

Discussion and conclusions

The state transition graph of a large network is a computationaly expensive object, especially for asynchronous networks. Likewise, computing the attractors of large Boolean networks is a NP-hard problem [START_REF] Akutsu | Determining a singleton attractor of a boolean network with nested canalyzing functions[END_REF][START_REF] Hong | An efficient steady-state analysis method for large boolean networks with high maximum node connectivity[END_REF]. Here, we proposed a new method that computes the attractors of large Boolean networks which are constructed by assembling together sets of smaller, already known, modules, as is often the case in physics or biology, for instance.

Theoretical advantages The asymptotic graph for the feedback interconnection of two modules is a new object computed using only the state transition graphs and attractors of each of the modules. This is a much less demanding state transition graph to compute, as attested by our studies.

Other methods are very fast in computing all singletons for very large networks (up to ≈ 1000 variables) [START_REF] Veliz-Cuba | Steady state analysis of boolean molecular network models via model reduction and computational algebra[END_REF], or very efficient in computing most types of attractors (up to ≈ 100 under variables) [START_REF] Zañudo | An effective network reduction approach to find the dynamical repertoire of discrete dynamic networks[END_REF] but our method is guaranteed to recover all attractors both singletons and cyclic attractors in the case of asynchronous updates. In addition, it sheds light on Boolean organization and the role of synchronous, asynchronous, or mixed updates.

An interesting observation is the complementarity of the synchronous and mixed mode asymptotic graphs, as particularly illlustrated by the mammalian example: the mixed mode asymptotic graph recovers the general structure of the atttractors in terms of the modules' attractors, while the synchronous asymptotic graph provides details of the state sequence within an attractor. The mixed asymptotic graph can be interpreted as the graph of a network, Table 4: Summary of interconnection results between each pair of networks listed in Table 2. The column "Interconnections" states the input/output connections used for each pair of models. The column "Attractors" contains the total number of attractors of each module (attractors of G A,u and G B,v , for all u, all v), and in the synchronous (G a ) and mixed (G a mx ) strategy asymptotic graphs. The column "Size" reports the size of the synchronous and mixed asymptotic graphs, to be compared with the total number of states, 2 N A ×N B , of the full interconnected systems. formed by an interconnection of two modules, where the dynamical timescale of each module is faster relative to the timescale of communication between modules. This property is often observed in the analysis of biological processes: fast metabolism / slow gene expression; fast post-translational processes/slow transcription.

Organism Interconnections Attractors Size A → B B → A G A,u G B,v G a G a mx G a G a mx G I - Cyanobacteria C T ∨ C T S R C U → C T 4
Synchronous modules with asynchronous interconnections The dichotomy between asynchronous and synchronous networks shows that, from a dynamical point of view, the asymptotic behavior of an asynchronous network can be fully predicted by that of its component modules (any pair). In contrast, a synchronous interconnection can generate new asymptotic behavior, not predicted by its component modules (some pairs). In fact, the asynchronous strategy ensures that any partition of a network contains representatives of each basins of attraction while under the synchronous strategy some network partitions fail to represent all basins of attraction, hence some attractors are "forgotten". In other words, the asymptotic graph method is guaranteed to give exact results whenever W 0 (of G a ) contains a representative of each basin of attraction in the synchronous case or G a * contains a representative of each family of pathways leading to an attractor.

The dynamical behavior of large networks is less costly to evaluate with synchronous updates. Our results suggest that an interesting coupling strategy is to use the mixed asymptotic graph method: each network module is evaluated according to a synchronous updating strategy (or any other, since only the state transition graphs of the modules are needed), and then the two networks are interconnected through an asynchronous asymptotic graph, which guarantees computation of all attractors of the full network and is much smaller than a synchronous asymptotic graph (see Table 4).

A method for analysis of emergent behaviors State-of-the-art mathematical models are presently available for a large selection of biological phenomena. One of the next steps in mathematical biology is to study the emergent dynamics arising from the coupling of two or more different models and thus predict the joints effects of two known biological processes. Our asymptotic graph method is an ideal framework to analyse emergent phenomena due to the coupling of previously studied networks. Indeed, given the state transition graphs of each module network (computed according to any updating strategy), an advantageous method is to adopt an asynchronous strategy for the coupling of the two modules. This will preserve the specific dynamics of each module and construct an asymptotic graph that recovers all the emergent behavior of the fully interconnected network. As illustrated by several examples throughout different organisms, the (mixed) asymptotic graph method is a powerful tool for studying the interaction between two known networks and predicting their joint behavior. As in the mammalian or vertebrate T helper cell examples, our method is also useful for interrogating models and testing new hypotheses, thereby checking the validity of the proposed interaction mechanisms.

Figure 1 :

 1 Figure 1: A dynamical system Σ I viewed as the feedback interconnection of two modules Σ A and Σ B .

Figure 2 :

 2 Figure 2: The state transition graphs for modules Σ A and Σ B of (2), using (a) asynchronous or (b) synchronous updates.

Figure 3 :

 3 Figure 3: (a) The asynchronous asymptotic graph and (b) the full network state transition graph for the interconnection of modules Σ A and Σ B in (2).

Figure 4 :

 4 Figure 4: Performance of the extended asymptotic graph method for asynchronous networks (darker curves, G ext ) and comparison with original method (lighter curves, G a * ), for randomly generated modules, and interconnected networks with p A = 1, 2, 3. Left: the percentage of failure (i.e., the generation of spurious). Right: the sizes of both asymptotic graphs relative to the size of the full network G I * .

Figure 5 :

 5 Figure 5: (a) The synchronous asymptotic graph and (b) the full network state transition graph for the interconnection of modules Σ A and Σ B in (2).

Theorem 3 Figure 6 :

 36 Figure 6: (a) The mixed updates asymptotic graph and (b) the full network state transition graph for the interconnection of modules Σ A and Σ B in (2).

Figure 7 :

 7 Figure 7: Performance of the asymptotic graph method for synchronous (blue curves), mixed (black curves), or asynchronous (red curves) interconnected networks. We considered two classes of networks: Horn Boolean functions (left column) and unate cascade functions (right column). For each class and each N A , 500 networks were generated, formed by a pair of modules of equal size N A = N B ∈ {3, 4, 5, 6, 7} and N in = 2. For the three strategies, the corresponding full graphs and asymptotic graphs were computed. The plots compare percentage of failures (top) and the sizes of the asymptotic graph relative to the size of the full network, |G a |/|G I | (bottom).

Table 2 :Figure 8 :

 28 Figure 8: The interconnected cyanobacteria clock (left) and transcription/translation cycle (right). Bold blue arrows denote input/output connections. The protein KaiA affects KaiC-S (C S ) at a different rate than KaiC-T (C T ) and KaiC-TS (C T S ). Therefore, KaiA is represented by two variables, A 1 (low KaiA level) and A 2 (high KaiA level). The unphosphorylated protein Kai C is denoted by C U .

  00001, 10000}, A 4 22 = {00111, 10100, 10110, 10111}, A 5 22 = {11110}, attractor of Σ A (under input 2), which corresponds to the circadian clock oscillatory orbit: 00000 → 10000 → 10100 → 10110 → 10111 → 00111 → 00001 → 00000

Figure 9 :

 9 Figure9: The interconnected mammalian cell cycle (left, adapted from[START_REF] Fauré | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF]) and clock (right, adapted from[START_REF] Comet | Simplified models for the mammalian circadian clock[END_REF]). Bold blue arrows denote input/output connections.

25 =

 25 {1000001000, 1001000000, 1001110000}, A 9 27 = {1000001110, 1000111110, 1001110110}, A 10 25 = {1000111000}, A 11 27 = {1001000110}, A 12 25 = {1110001100}.

17 ×

 17 {0011 * * * 00 * 10 * }, A 5 15 × {0011 * * * * * * 10 * }, A 6 22 × {1100 * * * * * * 0 * * }, A 7 22 × {1100 * * * * * * 0 * * },

Figure 10 :

 10 Figure 10: Partial visualization of the T helper cell network, containing only the nodes and edges involved in maturation of Th1p towards Th1 or Th2, as discussed in Section 7.3. Σ A nodes are coloured red and Σ B nodes are coloured green. Dashed arrows denote a simplified pathway, composed of several steps.

Figure 11 :

 11 Figure 11: The synchronous asymptotic graph of the Kai proteins interconnection. Dark shaded states belong to an attractor.

Figure 12 :

 12 Figure 12: The mixed mode asymptotic graph of the Kai proteins interconnection. Dark shaded states belong to an attractor.

Figure 13 :

 13 Figure 13: The five attractors of the synchronous asymptotic graph of the mammalian modules interconnection.

Figure 14 :

 14 Figure 14: The five attractors of the mixed mode asymptotic graph of the mammalian modules interconnection. Left: Q 2 mx

  B4 32 is not allowed in G a mx , by definition (see edge rule[START_REF] Chaves | Robustness and fragility of boolean models for genetic regulatory networks[END_REF]). Note also that 0010 → 0001 → 1001 is the only pathway in G I mx that allows transitions from the states in set Q a,2 mx to those in set Q a,3 mx . The definition of the edges of G a mx together with the lack of redundancy among pathways in G I mx originates a spurious attractor in G a mx .

	and note that the edge A 1 11 × B 5 41 → A 1 11 ×
		1 11 × B 5 41 , A 2 14 × B 1 11 , A 2 14 × B 5 41 } = {0000, 0010, 1100, 1110},
	Q a,3 mx = {A 4 13 × B 3 31 , A 4 13 × B 4 32 , A 5 23 × B 3 31 , A 5 23 × B 4 32 } = {1000, 1010, 1001, 1011}.
	However, Fig. 5(b) shows that the full network G I mx has only two attractors, corresponding to Q a,1 mx and Q a,3 mx . It follows that Q a,2 mx is a spurious attractor, which does not correspond to
	any true attractor of G I mx . To understand the origin of this spurious result, compare the edges linking the following equivalent vertices in G I mx and G a mx :
	0010	0000	A 1 11 × B 5 41 → A 1 11 × B 1 11
		↑	↑
		0001	A 1 11 × B 4 32
		↓	↓
		1001	A 5 23 × B 4 32

Table 3 :

 3 Classification of attractors into T helper cell phenotypes.

	Variable A 1 11 × B 1 11

  :

	B 6 12 :	Quiescent center cells, QC
	B 13 11 , B 20 21 :	Central pro-vascular tissues, in Meristem Transition Domain (C. Prov. TD3)
	B 14 11 , B 21 21 :	Root Cap (R.C. 2)
	B 15 11 , B 22 21 :	Cyclic attractor 2
	B 16 11 :	Central pro-vascular tissues, in Meristem Transition Domain (C. Prov. TD1)
	B 17 11 :	Root Cap (R.C. 1)
	B 18 11 :	Cyclic attractor 1
	B 23 11 :	C. Prov. TD1, with SHY2=0
	B 24 11 :	R.C. 1, with SHY2=0
	B 25 11 :	Cyclic attractor 1, with SHY2=0

Acknowledgements M. Chaves was partially supported by the French agency for research through project ICycle ANR-16-CE33-0016-01. D. Figueiredo and M.A. Martins are partially supported by the ERDF European Regional Development Fund through the Operational Programme for Competitiveness and Internationalisation -COMPETE 2020 and by National Funds through the Portuguese funding agency, FCT -Fundação para a Ciência e a Tecnologia within project POCI-01-0145-FEDER-016692 and project UID/MAT/ 04106/2013 at CIDMA. D. Figueiredo acknowledges the support of FCT via the PhD scholarship PD/BD/114186/2016. This work was partially supported by a France-Portugal partnership PHC PESSOA 2018 between M. Chaves (Campus France #40823SD) and M.A. Martins (POCI-01-0145-FEDER-030947).

 E). In this case, the attractors of each module are: G A,1 :

A 1 11 = {0100010100}, A 2 11 = {1000001110, 1000100011, 1001100000, 1010000110, 1011000100, 1011100100}, A 3 12 = {1000101011}, G A,2 :

A 4 As expected, the attractors of the graph G A,1 correspond exactly to those listed in [START_REF] Fauré | Dynamical analysis of a generic boolean model for the control of the mammalian cell cycle[END_REF]. Attractors A 1 11 and A 4 21 correspond to the quiescent cell state while the (full) attractor A 2 11 ∪ A 3 12 corresponds to the known cell cycle progression. A similar cycle is also recovered in G A,2 , with the difference that states 1001100000, 1011000100, and 1011100100 are replaced, respectively, by 1000100000, 1010000100, and 1010100100: the Boolean sequence of the latter three equals that of the former three except in the 4th digit, which is set to 0. This is because u = 2 implies CycE + ≡ 0 (and CycE is the 4th variable). In both cases, the semi-attractors A ) where only the first two correspond to circadian oscillations.

Mammalian synchronous and mixed mode asymptotic graphs Among the different interaction schemes tested, all allow at least one cycle that captures joint oscillations of cell cycle and clock. For the scheme where BMAL affects CycE (13), the synchronous asymptotic graph has five attractors, but only one of these recovers joint oscillations (Q 2 ), all others corresponding to a quiescent cell with oscillatory clock: