Madalena Chaves
email: madalena.chaves@inria.fr

Luca Scardovi
email: scardovi@scg.utoronto.ca

Eleni Firippi
email: eleni.firippi@inria.fr

† L Scardovi

Coupling and synchronization of piecewise linear genetic regulatory systems

published or not. The documents may come

INTRODUCTION

The two interaction motifs known as the positive and the negative feedback loops are present at the core of a wide diversity of regulatory networks, involved in pattern formation, memory phenomena, calcium waves, circadian rhythms, or cell cycle [START_REF] Yu-Chen | Robust, tunable biological oscillations from interlinked positive and negative feedback loops[END_REF]. The positive and negative feedback loops are circuits whose topology is generally organized in a ring composed of n elements, where the i-th element activates/inhibits the (i + 1)th element. In addition, each element also has a natural degradation. In a positive (resp., negative) feedback loop, the number of inhibitions is even (resp., odd), so that the composed effect of all the interactions is "positive" or "negative".

From the work initiated by R. Thomas and followed by other authors, it is known that these two types of circuits give rise to fundamentally distinct dynamical behaviors: the topology of a positive loop leads to multistability while a negative loop leads to oscillatory behavior (see, for instance, [START_REF] Thomas | Dynamical behaviour of biological regulatory networks -I. biological role of feedback loops and practical use of the concept of feedback loop-characteristic states[END_REF]). Depending on the number of elements in the circuit, the form of inhibition functions, and the parameters, the oscillations may be damped or sustained [START_REF] Goodwin | Oscillatory behavior in enzymatic control processes[END_REF][START_REF] Michael | A synthetic oscillatory network of transcriptional regulators[END_REF][START_REF] Farcot | Periodic Solutions of Piecewise Affine Gene Network Models with Non Uniform Decay Rates: The Case of a Negative Feedback Loop[END_REF]. These theoretical predictions have been verified not only in different "natural" systems [START_REF] Yu-Chen | Robust, tunable biological oscillations from interlinked positive and negative feedback loops[END_REF] but also in synthetically constructed molecular networks, the repressilator [START_REF] Michael | A synthetic oscillatory network of transcriptional regulators[END_REF] and the toggle switch [START_REF] Gardner | Construction of a genetic toggle switch in escherichia coli[END_REF].

In living organisms, similar cells in a group within an organ, a tissue or a specific part of the organism, also express the same regulatory networks. In response to external or internal signals, cells communicate by exchanging proteins and other elements, and evolve to adapt their states to new conditions. A natural question is to ask what is the behavior of an interconnected network of cells: do the cells synchronize their dynamics? Does the coupled network generate new emergent solutions? And under which conditions?

An example is the synchronization of circadian clocks in mammalians. Many experiments have shown that individual cells throughout the various tissues of the body have identical autonomous clocks [START_REF] Saini | The mammalian circadian timing system: synchronization of peripheral clocks[END_REF]. Some experiments show that a population of identical cells in a culture may show independent, not synchronized oscillations [START_REF] Feillet | Phase locking and multiple oscillating attractors for the coupled mammalian clock and cell cycle[END_REF]. But other experiments have shown that a group of neurons in the brain, the suprachiasmatic nuclei, transmit signals that organize and synchronize the circadian rhythms throughout the body [START_REF] Silver | A diffusible coupling signal from the transplanted suprachiasmatic nucleus controlling circadian locomotor rhythms[END_REF].

Another example is pattern formation, for instance the segment polarity network in early fly developement [START_REF] Dassow | The segment polarity network is a robust developmental module[END_REF]: the embryo goes through different stages where neighbor cells communicate by exchanging molecules of specific proteins (cf. Wingless). Each cell will finally stabilize at distinct steady states, thus forming the striped pattern of the embryo.

These two examples illustrate that the result of coupling identical cells may lead either to synchronized behavior (in the case of mammalian circadian clock) or, conversely, to the formation of an expression pattern where cells are not synchronized to the same states.

In this paper, we will study the interconnection of N similar compartments, all containing the same species, coupled through a diffusive process. We will focus mainly on the interconnection of positive feedback loops and show that this compartmental system has the capacity to generate several new locally stable steady states, such as those observed at the basis of phenomena such as pattern formation.

PIECEWISE LINEAR BASIC MOTIFS

Some classical models have been studied for the positive and negative loops. For instance, [START_REF] Goodwin | Oscillatory behavior in enzymatic control processes[END_REF] studies a negative loop with one nonlinear sigmoidal inhibition term for one element and linear activation terms for all the other elements. To study their molecular networks, both [START_REF] Michael | A synthetic oscillatory network of transcriptional regulators[END_REF] and [START_REF] Gardner | Construction of a genetic toggle switch in escherichia coli[END_REF] have used increasing or decreasing Hill functions to represent their activation and inhibition terms, together with linear degradation rates, a widely used formalism in genetic regulatory networks.

Most of the existing work on synchronization of coupled compartmental models used nonlinear ODEs or nonlinear operators as mathematical models (see, e.g., [START_REF] Scardovi | Synchronization of interconnected systems with applications to biochemical networks: an Input-Output approach[END_REF] and references within). Throughout this paper, we will use piecewise linear functions to approximate the sigmoidal Hill functions: this approach has the advantage of providing both theoretical and algorithmic standpoints for analysis of the dynamics [START_REF] Glass | Prediction of limit cycles in mathematical models of biological oscillations[END_REF][START_REF] Farcot | Periodic Solutions of Piecewise Affine Gene Network Models with Non Uniform Decay Rates: The Case of a Negative Feedback Loop[END_REF][START_REF] Poignard | A stability result for periodic solutions of some nonmonotonic smooth negative feedback systems[END_REF]. In this framework, increasing and decreasing Hill functions are thus approximated by step functions:

s + (x, θ) =    0, x < θ [0, 1], x = θ 1, x > θ, s -(x, θ) = 1 -s + (x, θ).

Positive and negative feedback loops

The feedback loop circuits are then represented by:

ẋ1 = V 1 s * (x n , θ n ) -γ 1 x 1 ẋi = V i s + (x i-1 , θ i-1 ) -γ i x i , i = 2, . . . , n, (1)
where s * = s + for the positive loop and s * = s -for the negative loop. We have the following assumptions: A1. V i γ i > θ i ; where A1 garantees that each variable may express increasing or decreasing dynamics. It is clear that the

interval [0, V 1 /γ 1 ] × • • • × [0, V n /γ n ]
is forward-invariant for the system [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF]. Observe that this rectangular region can be partitioned into 2 n subregions, of the form

D p = D 1 p 1 × D 2 p 2 × • • • D n
pn , called regular domains where p i ∈ {0, 1} are defined by the threshold values θ i :

D i 0 = (0, θ i ), D i 1 = θ i , V i γ i .
It follows immediately that, in each regular domain, the step functions are constant, and the solutions to the system can be exactly computed. The boundaries between regular domains are also called switching domains, and correspond to regions of the space where at least one variable is equal to its threshold (x i = θ i , for some i). At these boundaries, the system is defined as a differential inclusion and solutions may be constructed according to [5] (for more details see, for instance [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF]).

It is well known that positive feedback loops exhibit two locally stable steady states corresponding, for instance, to ON and OFF states [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF]; while the negative feedback loop with n ≥ 3 exhibits a single and unstable steady state, at (θ 1 , . . . , θ n ) and converges to a periodic orbit [START_REF] Farcot | Periodic Solutions of Piecewise Affine Gene Network Models with Non Uniform Decay Rates: The Case of a Negative Feedback Loop[END_REF].

Coupling identical motifs

We consider a compartmental model with N compartments, where each compartment is composed by n species with dynamics (1). The compartments are linked by a diffusive process. This is a reasonable hypothesis in biological networks, for instance the segment polarity network model developed in [START_REF] Dassow | The segment polarity network is a robust developmental module[END_REF] considers diffusion of Wingless proteins across cell membranes. From their experiments on the suprachiasmatic nuclei, [START_REF] Silver | A diffusible coupling signal from the transplanted suprachiasmatic nucleus controlling circadian locomotor rhythms[END_REF] conclude that a diffusible signal is sufficient to control circadian rhythms.

For simplicity we will assume that:

A2. The systems are coupled by diffusion terms through their first variables.

We will denote the state of species i in compartment j with x i,j and with

x i = [x i,1 , x i,2 , . . . , x i,N
] the N -dimensional vector of the states of species i across the N compartments. We will study the following model of N -coupled identical systems

ẋ1,k = V 1 s * (x n,k , θ n ) -γ 1 x 1,k + N j=1 a k,j (x 1,j -x 1,k ) ẋi,k = V i s + (x i-1,k , θ i-1 ) -γ i x i,k , i = 2, . . . , n, (2)
where k = 1 . . . , N . The nonnegative coefficients a k,j define the coupling topology between compartments: if a k,j = 0 there is no diffusion between compartment k and j, if a k,j > 0 the diffusion coefficient is proportional to a k,j . It is possible to associate a weighted graph with N nodes and a link between node k and node j whenever a k,j > 0. The resulting graph is strongly connected if any two nodes are connected by a path. In this paper we will make the following assumption on the coupling topology.

A3. The graph encoding the diffusive interconnection is symmetric and strongly connected.

Let's define the N -dimensional Boolean vector function

q(x i ) = [s * (x i,1 , θ i ), . . . , s * (x i,N , θ i )] , i = 1, . . . n,
where for notational convenience we dropped the dependency from the parameter θ i . Note that s * = s + for i = 1, . . . n -1. Then, (2) can be rewritten as

ẋ1 = V 1 q(x n ) -(Γ 1 + L) x 1 (3) ẋi = V i q(x i-1 ) -Γ i x i , i = 2 . . . , n, (4)
where L is the symmetric N × N Laplacian matrix defined as

l kj =      N i=1 a ki , j = k -a kj , j = k,
and Γ i are diagonal matrices with elements γ i > 0.

The following properties hold. Under assumption A3, the matrix (Γ 1 + L) is invertible and strictly positive and -(Γ 1 + L) is a Hurwitz matrix [START_REF] Dörfler | Electrical Networks and Algebraic Graph Theory: Models, Properties, and Applications[END_REF]. Similarly to the previous section we can define the regular domains for the coupled system as

D P = D p 1 × • • • D p N ,
where D p i is the regular domain of compartment i and P is a n × N matrix with columns p i , i = 1, . . . , N . For example, p i,j = 1 corresponds to regions where θ i < x i,j < V i γ i .

STABLE STEADY STATES

Proposition 1 Consider system (2) under assumptions A1-A3. Then:

(a) for negative loops, there are no steady states in regular domains;

(b) positive feedback loops have at least two locally stable steady states corresponding to synchronized ON and synchronized OFF states, respectively.

Proof:

Let x = {x i,j } be an equilibrium of system (2), contained in some regular domain, D P . Then it must satisfy the following equalities

xi = V i Γ i -1 q(x i-1 ), i = 2, . . . , n.
As mentioned in the previous section (Γ 1 + L) is an invertible matrix, which implies

x1 = V 1 (Γ 1 + L) -1 q(x n ). ( 5
)
Suppose that x1,j > θ 1 (which implies P1,j = 1). Then s + (x 1,j , θ 1 ) = 1 which implies x2,j = V 2 /γ 2 and, by assumption A1, also xi,j = V i /γ i for all i = 2, . . . , n (and, consequently, Pi,j = 1 for all i = 2, . . . , n). This implies q j (x i ) = 1 for every i = 1, . . . n. Conversely, if x1,j < θ 1 , it follows that q j (x i ) = 0 for every i = 1, . . . n. Therefore, defining a non-negative sign function sign(r) := 0, r ≤ 0 1, r > 0, the steady states must also satisfy, for negative and positive loops, respectively:

(a) 1 -q j (x n ) = sign(V 1 (Γ 1 + L) -1 j q(x n ) -θ 1 ) (b) q j (x n ) = sign(V 1 (Γ 1 + L) -1 j q(x n ) -θ 1 ),
for every j = 1, . . . , N . Now, define two vectors that differ only on the j entry as q(x n ) if q j (x n ) = 0 and q(

x n ) if q j (x n ) = 1, with qk (x n ) = qk (x n ) for k = j. Define likewise x1 = V 1 (Γ 1 + L) -1 q(x n ) and x1 = V 1 (Γ 1 + L) -1 q(x n ). Since (Γ 1 + L) -1
is strictly positive, we have x1,j < x1,j for any j.

Let's consider case (a). If P1,j = 1, then q j (x 1 ) = 1, which implies that q j (x i ) = 1 for i = 1, . . . , n -1 and q j (x n ) = 0. This implies that x1,j = V 1 (Γ 1 + L) -1 j q(x n ) = x1,j > θ 1 . But now also x1,j > θ 1 with q j (x n ) = 1 yields a steady state in the same domain P , contradicting the fact that q j (x n ) = 0. If P1,j = 0, then q j (x n ) = 1 and an analogous argument applies, implying that x1,j < x1,j < θ 1 and reaching the contradiction q j (x n ) = 0. We conclude that for the negative loop, there are no equilibria. In contrast, for positive loops, it is possible to find at least two locally stable steady states, corresponding to the synchronized ON and OFF configurations xi = (0, 0, . . . , 0) and xi =

V i γ i , . . . , V i γ i , (6)
for all species i = 1, 2 . . . , n. Since -(Γ 1 + L) and -Γ i , i = 2 . . . , n are Hurwitz matrices, both the equilibria are locally asymptotically stable.

As an example, consider a compartmental system with N = 4 positive loops, of dimension n = 3, connected in a "ring-like" topology

L =     -2 1 0 1 1 -2 1 0 0 1 -2 1 1 0 1 -2     ,
with k = 0.5. The parameters are as follows, with randomly chosen V and θ: γ = [0.4, 0.4, 0.4] , V = [0.9670, 2.4922, 0.3681] ,

θ = [0.9376, 0.6552, 0.7529] .

With these parameters, the coupled system has 10 distinct locally stable steady states, shown in Table 1. Fig. 1 shows the first coordinates converging to the steady state with q(x) = 0010. This means that compartment j = 1 converges to the point (0.3597, 0, 0) , compartments j = 2, 4 both converge to the point (0.5036, 0, 0) and compartment j = 3 converges to (1.0505, V 2 /γ 2 , V 3 /γ 3 ) .

SYNCHRONIZATION OF POSITIVE LOOPS

In this section, set s * ≡ s + . The first observation is that the first coordinates are either all zero or all positive. Proposition 2 Consider system (2) under assumptions A1-A3. Suppose x1,j = 0 for some compartment j, then all compartments are synchronized at the origin.

Figure 1: The trajectories of the four compartments, each converging to a distinct steady state. Parameters as in [START_REF] Glass | Prediction of limit cycles in mathematical models of biological oscillations[END_REF]. The red stars represent the steady states listed in Table 1.

Proof:

At steady state, equation (2) implies:

(γ 1 + L jj )x 1,j = q j (x) + =j L 1, x1,
If x1,j = 0, then xk,j = 0 for all species k which implies q j (x) = 0. Therefore also x1, = 0 for all such that L 1, > 0. If L is irreducible, then x1,j = 0 for all j = 1, . . . , N and consequently all compartments are synchronized at the origin.

This result implies that, in any of the new steady states originated from the coupling, no compartment may have all species at zero. Our next results focus on a specific interconnection topology.

All to all interconnections

In this section we consider the special case of an all-to-all interconnection structure, where every pair of compartments are connected, and the diffusive parameter is identical to k for all connections. In this special case

L + Γ = k(N I -11 ) + γ 1 I. ( 9
)
By computing the inverse we get

(L + Γ) -1 = 1 N k + γ 1 I + 1 γ 1 (N k + γ 1 ) k11 ,
and by substituting in (5) we obtain

x1 = V 1 1 N k + γ 1 I + 1 γ 1 (N k + γ 1 ) k11 q(x n ).
We conclude that the first coordinates of the equilibria are

x1,i = V 1 1 N k + γ 1 q i + V 1 1 γ 1 (N k + γ 1 )
N k q [START_REF] Poignard | A stability result for periodic solutions of some nonmonotonic smooth negative feedback systems[END_REF] where q = ( N i=1 q i )/N , and for notational conveniency we have dropped the dependency on xn . In this case, if k is sufficiently large, then the only possible steady states of the coupled system are given by [START_REF] Gardner | Construction of a genetic toggle switch in escherichia coli[END_REF]. Define

k min = max γ 2 1 θ 1 |V 1 -N γ 1 θ 1 | , γ 1 (V 1 -γ 1 θ 1 ) |γ 1 θ 1 N -V 1 | . ( 11
)
Assume that

V 1 γ 1 θ 1 = N N -b , b = 1, . . . , N -1. ( 12
)
Proposition 3 Suppose L has an all-to-all interconnection structure and that the parameters satisfy condition [START_REF] Saini | The mammalian circadian timing system: synchronization of peripheral clocks[END_REF]. If k > max{γ 1 /N, k min }, then the coupled system of positive loops has only two steady states, given by ( 6).

Proof:

Consider L satisfying [START_REF] Grinstead | Introduction to Probability[END_REF]. Then the steady states of the coupled system satisfy [START_REF] Poignard | A stability result for periodic solutions of some nonmonotonic smooth negative feedback systems[END_REF]. We want to show that there is a sufficiently large k such that any equilibrium different from the synchronized ON or OFF states does not exist. Let x1 be the first components of such an equilibrium. The corresponding vector q(x n ) must have at least one component equal to zero and one component equal to one. Let q i (x n ) = 1 and q j (x n ) = 0 be any pair of those components. Then, for the equilibrium to not exist, it is required that either

V 1 1 N k + γ 1 + V 1 1 γ 1 (N k + γ 1 ) N k q < θ (13)
or

V 1 1 γ 1 (N k + γ 1 ) N k q > θ. (14)
Also note that, by definition, for the non synchronized equilibria q takes values of the form (N -b)/N with b = 1, . . . , N -1. Assumption [START_REF] Saini | The mammalian circadian timing system: synchronization of peripheral clocks[END_REF] implies that q = θ 1

γ 1 V 1 . Suppose first that q < θ 1 γ 1 V 1 . Then inequality (13) provides the bound k > γ 1 (V 1 -γ 1 θ 1 ) (γ 1 θ 1 -V 1 q) N
while inequality ( 14) is never satisfied. Now, suppose that q > θ 1

γ 1 V 1 .
Then inequality ( 13) is never satisfied, while inequality [START_REF] Yu-Chen | Robust, tunable biological oscillations from interlinked positive and negative feedback loops[END_REF] provides the bound

k > γ 2 1 θ 1 (V 1 q -γ 1 θ 1 ) N
Since the above must hold for all the possible q corresponding to non synchronized ON or OFF equilibria, the choice

k > k min = max γ 2 1 θ 1 |V 1 -N γ 1 θ 1 | , γ 1 (V 1 -γ 1 θ 1 ) |γ 1 θ 1 N -V 1 | .
guarantees the non existence of the non synchronized equilibrium. An immediate consequence of Propositions 2 and 3 is the synchronization of the positive loops for large diffusive parameter, either to the origin "all OFF" or to the "all ON" state.

Corollary 1 cr:synch-pos Suppose L has an all-to-all interconnection structure, the parameters satisfy condition [START_REF] Saini | The mammalian circadian timing system: synchronization of peripheral clocks[END_REF], and k > max{γ 1 /N, k min }. Then the individual systems synchronize to one of the states [START_REF] Gardner | Construction of a genetic toggle switch in escherichia coli[END_REF].

The corollary holds because: (a) for large k, Proposition 3 guarantees that there are only the OFF and ON steady states; (b) the trajectories remain in a bounded region of the state space; (c) Proposition 2 guarantees that the solutions in all compartments may converge to the origin or, if species x 1j in compartment j converges to a positive solution, then none of the other compartments can converge to the origin, hence all must converge to the ON equilibrium.

The lower bound on k

In this section, numerical simulations illustrate the behavior of the lower bound on k. For the simulations, a large number of positive loops with n = 3 species were generated: the parameters were randomly chosen from uniform distributions on the intervals

V r i ∈ (0, 3), θ r i ∈ 0, V r i γ i ,
for i = 1, 2, 3, r = 1, . . . , 1000 and fixed degradation rates, for simplicity,

γ 1 = γ 2 = γ 3 = 0.4.
For each of these 1000 sets of parameters, we considered N compartments with equal diffusive parameter k, in all possible combinations of N = 2, . . . , 15, k ∈ {0.2, 0.5, 1.0, 3.0, 7.0, 10.0, 15.0, 20.0}. [START_REF] Dassow | The segment polarity network is a robust developmental module[END_REF] For the coupling topology, we considered both an all-to-all and a ring-like interconnection. Below only the results of the all-to-all are shown, but the simulations corresponding to the ring-like interconnection give quite similar results. The first observation is that the number of distinct steady states increases with N and decreases with k, as expected (Fig. 2).

The question arises of analysing the lower bound for the diffusive coupling parameter. In the case of an all-to-all interconnection, the expression k min [START_REF] Thomas | Dynamical behaviour of biological regulatory networks -I. biological role of feedback loops and practical use of the concept of feedback loop-characteristic states[END_REF] provides a lower bound on k such that the coupled system admits only the ON and OFF states [START_REF] Gardner | Construction of a genetic toggle switch in escherichia coli[END_REF]. However, the proof of Prop. 3 provides a sharper analytical bound for each q = (N -b)/N , b = 1, . . . , N -1:

k bound > 1 |V 1 q -γ 1 θ 1 |N max{γ 2 1 θ 1 , γ 1 (V 1 -γ 1 θ 1 )} (16)
These lower bounds are shown in Fig. 3 (black curves), and illustrated by numerical results (red dots). Furthermore, for each fixed N , there are N -1 critical values of the ratio V 1 θ 1 γ 1 , given in [START_REF] Saini | The mammalian circadian timing system: synchronization of peripheral clocks[END_REF], for which the lower bound k bound is not defined. These values are shown in Fig. 4, ilustrated by numerical simulations (red dots).

Even though the set of critical parameters is a set of measure zero, it nevertheless introduces points of rupture in the regions of parameters, as seen in Fig. 3. This may indicate a loss of robustness of the compartmental system on the neighborhood of the critical parameters.

SYNCHRONIZATION OF NEGATIVE LOOPS

This last section briefly illustrates the case of a compartmental system of negative loops (case s * = s -). Proposition 1 states that no locally stable steady states may be generated for system [START_REF] Michael | A synthetic oscillatory network of transcriptional regulators[END_REF]. However, it is not possible to exclude the generation of new periodic orbits. Indeed, a random exploration of the parameters shows that it is also possible to find many examples where new periodic orbits appear, as in the example below. In this periodic orbit the compartments are not synchronized. Analysis of the effect of the diffusive parameter on the dynamics of the coupled system is the object of future work.

Consider N = 4 negative loops, of dimension n = 3, connected in a "line" topology:

L =     -2 1 0 0 1 -2 1 0 0 1 -2 1 0 0 1 -2     ,
with k = 0.5. The parameters are as follows, with randomly chosen V and θ:

γ = [0.4, 0.4, 0.4] , V = [1.3260, 1.0857, 1.8768] , (17)
Simulations show that each of the four systems may converge to a distinct periodic orbit (Fig. 5).

CONCLUSIONS

This paper investigates a compartmental model with N identical compartments, linked by a diffusive process. Each compartment is composed of n species whose dynamics is governed by a basic regulatory motif. The analysis focuses on a frequently occurring motif, the positive feedback loop. The model is represented in the piecewise linear formalism, which leads to explicit calculations in terms of the parameters of the model. Our main results show that a number of locally stable states are created for the compartmental model, that were not present in each single compartment. As expected, the number of these new steady states increases with the number of coupled compartments, but decreases with the diffusive strength. For the particular case of an all-to-all interconnection structure, we give an explicit lower bound on the diffusive parameter, to guarantee that the coupled system retains only the steady states common to all single compartments.

An unexpected result is the existence of a special relationship between the number of components N and the maximal concentration-to-activity threshold ratio V 1 /(γ 1 θ 1 ), which introduces sets of parameters for which the lower bound is not finite. This relationship implies that, for very specific parameter sets, the compartments cannot be guaranteed to synchronize.

Future work includes an extension of the lower bound on the diffusive strength for general interconnection topologies, as well as further analysis of the compartmental system with negative feedback loops.

Our analysis indicates that compartmental models composed of coupled positive loops are compatible with pattern formation mechanisms.

Figure 2 :

 2 Figure 2: Average number of distinct steady states for (2), as a function of N and k, under an all-to-all interconnection. Each point in a curve is an average over the 1000 systems with N compartments and diffusion k.

Figure 3 :

 3 Figure 3: The lower bound on the diffusive parameter, so that the coupled system admits only the ON and OFF steady states. Clockwise from upper left, the number of compartments are: N = 2, N = 3, N = 4, N = 5. The curves (11) for k min , b = 1, . . . N -1, are shown in black. The red points represent numerical data, showing the minimal k for each parameter V 1 /(θ 1 γ 1 ). Note that, in the simulations, the k values belong to a finite grid (15), hence the red dots do not always coincide with the black curves.

Figure 4 :

 4 Figure 4: The critical values of the ratio of parameters, V 1 θ 1 γ 1 . Each straight line corresponds to a fixed value of b. The numerical points correspond to those ratios where no k value within the grid (15) could reduce the number of steady states to 2.

θ

  = [0.9388, 0.2869, 0.2446] .

Figure 5 :

 5 Figure 5: Four coupled negative feedback loops. The trajectories of the four coupled negative feedback loops, each converging to a distinct periodic orbit.

Table 1 :

 1 The locally stable steady states for the coupled positive feedback loops with parameters[START_REF] Glass | Prediction of limit cycles in mathematical models of biological oscillations[END_REF] and k = 0.5, and their respective domain.

	x11	x21	x31	x41	q(x)
	0	0	0	0	0000
	0.5036 0.3597 0.5036 1.0505 0001
	0.3597 0.5036 1.0505 0.5036 0010
	0.8634 0.8634 1.5541 1.5541 0011
	0.5036 1.0505 0.5036 0.3597 0100
	0.8634 1.5541 1.5541 0.8634 0110
	1.0505 0.5036 0.3597 0.5036 1000
	1.5541 0.8634 0.8634 1.5541 1001
	1.5541 1.5541 0.8634 0.8634 1100
	2.4175 2.4175 2.4175 2.4175 1111

ACKNOWLEDGMENTS

M.C. and E.F. acknowledge the support of the French National Agency for Research through project ICycle ANR-16-CE33-0016-01. L. S. acknowledges the support of the National Sciences and Engineering Research Council of Canada and the Canada First Research Excellence Fund.