E Firippi

M Chaves

Topology-induced dynamics in coupled synthetic oscillators Topology-induced dynamics in a network of synthetic oscillators with piecewise affine approximation

Topology-induced dynamics in a network of synthetic oscillators with piecewise affine approximation

Elena Firippi, Madalena

I. INTRODUCTION

A deeper understanding of cellular complex processes is due to recent research in synthetic biology, which involves the design and implementation of small genetic circuits to mimic the function of more complex biochemical mechanisms. Analysis of these smaller circuits thus contributes to elucidate the origins of complex biological behaviour such as the bacteriophage lambda switch [START_REF] Gardner | Construction of a genetic toggle switch in escherichia coli[END_REF] or biological rhythms and oscillations (see the repressilator by [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF] and a review of synthetic oscillators by [START_REF] Purcell | A comparative analysis of synthetic genetic oscillators[END_REF]).

New research in synthetic biology is heading towards the use of basic modules (such as a toggle switch or a repressilator) as functional units to (re-)create more complex systems, followed by the control and regulation of these systems [START_REF] Del Vecchio | Control theory meets synthetic biology[END_REF]. In this context, the mathematical analysis of the coupled synthetic modules is essential to predict the behaviour of the new system, investigate the paa) Electronic mail: eleni.firippi@inria.fr b) Electronic mail: madalena.chaves@inria.fr rameter conditions leading to desired behaviour, and generally help experimental design. Examples of such studies include the coupling of different circuits to increase the dynamical range of the system [START_REF] Perez-Carrasco | Combining a toggle switch and a repressilator within the ac-dc circuit generates distinct dynamical behaviors[END_REF], form new spatial patterns [START_REF] Barbier | Controlling spatiotemporal pattern formation in a concentration gradient with a synthetic toggle switch[END_REF] or the regulation of a bacterial consortium [START_REF] Fiore | In-silico analysis and implementation of a multicellular feedback control strategy in a synthetic bacterial consortium[END_REF].

In this paper, we will consider the coupling of N identical oscillatory circuits, to study some questions relevant in biological rhythms [START_REF] Hafner | Effect of network architecture on synchronization and entrainment properties of the circadian oscillations in the suprachiasmatic nucleus[END_REF], including the synchronization of the individual circuits or the effect of the topology of interactions and their strength in the robustness of periodic oscillations.

A minimal two gene oscillator that was first proposed by [START_REF] Smolen | Frequency selectivity, multistability, and oscillations emerge from models of genetic regulatory systems[END_REF] and later implemented by [START_REF] Stricker | A fast, robust and tunable synthetic gene oscillator[END_REF], consists of two transcriptional factors competing with each other in a negative feedback loop, together with two self-regulatory loops, one positive and one negative. In previous work, we suggested an improved design for this two gene oscillator by showing that eliminating the self-regulatory loop of the repressor substantially increases the parameter region that admits oscillatory regimes [START_REF] Firippi | Period -control in a coupled system of two genetic oscillators for synthetic biology[END_REF]. To further analyze this improved circuit, we now write the continuous system in the piecewise affine model framework first introduced by [START_REF] Glass | The logical analysis of continuous, non-linear biochemical control networks[END_REF].

This formalism is based on the approximation of sigmoidaltype functions by step functions, and the phase space of the system is typically partitioned into rectangular regular domains, in which the system is affine and solutions can be explicitly computed [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF]. However, in our approximation of the two-gene oscillator, the regular domains are not rectangular but instead defined in terms of a quadratic expression relating the two coordinates. The solutions in each domain can still be explicitly computed, but the concatenation of the different pieces at the boundary between domains, has to be analyzed in a new way. In this new framework, we will compute the Poincaré map of the twodimensional system to prove the existence and uniqueness of a periodic orbit, and will also characterize the corresponding parameter region.

As an application of our two-gene oscillator analysis, we then study a network of N identical piecewise affine twogene oscillators, linked by diffusion of one of the variables. There are still few studies of coupled identical oscillators in a piecewise affine formalism: for instance, [START_REF] Edwards | On synchronization and cross-talk in parallel networks[END_REF] prove synchronization under a weak coupling, and [START_REF] Nicks | Clusters in nonsmooth oscillator networks[END_REF] study clusters of synchronized elements, with applications to neuron networks.

Our goal is to investigate the effect of the topology of interactions in the network's dynamics, so we consider three different topologies: one-to-all, all-to-all, and random connections. The main result is that all these topologies generate new "coupling" stable steady states, in addition to the periodic orbit, thus possibly contributing to a loss of synchronization and/or a loss of oscillatory behaviour in the coupled network. The generation of new steady states in the coupled system was also observed in [START_REF] Chaves | Coupling and synchronization of piecewise linear genetic regulatory systems[END_REF], for a network of coupled positive circuits. The number of new steady states depends on the number N of oscillators in the network and can be exactly computed for both the one-to-all and all-to-all topologies.

For all topologies, we use both numerical simulations and theoretical analysis to study the effect of the coupling parameter, the basins of attraction of the original periodic orbit, and other properties, to establish conditions that guarantee cell synchronization.

II. A TWO-GENE SYNTHETIC OSCILLATOR

The original model introduced by Smolen, Baxter, and Byrne (1998) describes a negative feedback circuit of two components, the transcription factors TF-A and TF-R: TF-A is a transcriptional activator that can bind to responsive elements DNA sequences; while TF-R is a protein that represses transcription by competing with TF-A for binding to the responsive elements. TF-A is an activator both for TF-R and itself and, likewise, TF-R is a repressor for both TF-A and itself, leading to a negative feedback loop with a self-activating activator and a self-repressing repressor. For simplicity, TF-A and TF-R are denoted by A and R respectively.

In a previous study [START_REF] Firippi | Period -control in a coupled system of two genetic oscillators for synthetic biology[END_REF], through bifurcation analysis, we have found that removing the selfrepressing loop significantly increases the parameter range where oscillations are observed. This leads to the circuit shown in Fig. 1 and described by the following equations:

dA dt = V A A 2 A 2 + θ 0 (1 + R/θ 2 ) -γ A A + r bas dR dt = V R A 2 A 2 + θ 1 -γ R R (1)
The parameters V A , V R > 0 express the synthesis rate of the model components and are measured in concentration.min -1 . The parameters γ A , γ R > 0 represent the degradation rates and they are measured in min -1 . The concentration and activity thresholds are described by parameters θ 0 , θ 1 , θ 2 > 0 that also have concentration units.

A R Figure 1. The improved design of the Smolen oscillator, proposed by [START_REF] Firippi | Period -control in a coupled system of two genetic oscillators for synthetic biology[END_REF].

A. Piecewise affine formalism

The piecewise affine (PWA) systems studied first by [START_REF] Glass | The logical analysis of continuous, non-linear biochemical control networks[END_REF] facilitate the system further study and analysis. This formalism is based on the approximation of the Hill functions that express the activation and inhibition (or repression), by increasing and decreasing step functions.

A synopsis of the PWA formalism based on the study by [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF] is recalled in this section: if x i is the protein concentration of the gene i and θ j i the concentration threshold that characterises the effect of x i on x j , the increasing step function s + : R + × R + → {0, 1} is defined:

s + (x i , θ j ) = 1, if x i > θ j i 0, if x i < θ j i (2)
Accordingly, the step function to express inhibition (or repression) will be:

s -: R + × R + → {0, 1}, s -(x i , θ j i ) = 1 - s + (x i , θ j i ).
The superscript j indicates the interaction of gene i on gene j of the network.

The production rate of variables x i ∈ [0, +∞), i = 1, . . . , n is given by the function f i : R n + → R + :

f i (x) = ∑ 1≤l≤m κ il D il (x) (3)
where κ il > 0 are the synthesis rate parameters and D il : R n + → {0, 1} are Boolean-valued functions that express the regulation of variable x i by other variables. The functions D il (x) may be products of increasing or decreasing step functions s +/-(x i , θ j i ). A general form of the system will be:

ẋi = κ i0 + f i (x) -γ i x i (4)
where κ i0 is a basal term, x = (x 1 , . . . , x n ) t > 0 contains the protein concentrations and γ i x i > 0 expresses the degradation rate of protein x i . The thresholds θ j i define hyperplanes and divide the phase space into regular domains. The regions such that x i = θ j i for some i and j are called switching domains and form the boundaries of the regular domains. In each regular domain B, the vector field is uniquely defined, the synthesis rate is constant (κ i0 + κ iB ):

ẋi = κ i0 + κ iB X -γ i x i
and the computation of the system solution, while x ∈ B, is straightforward, with trajectories moving towards the corresponding focal point φ iB = (κ i0 + κ iB )/γ i . As the trajectory reaches the boundary of B, the vector field changes and the solution must be computed in the new domain.

In PWA systems, a general property is that a focal point belonging to its own domain (i.e., φ iB ∈ B) is a locally stable steady state. The connection between continuous and piecewise affine models has been closely studied by [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: Equilibria and their stability[END_REF]. In the general approach the continuous activation and inhibition functions given by Hill functions are approximated by step functions and lead to a partition of the state space into rectangular domains. In contrast, for the twogene oscillator (1), this approximation leads to a partition defined by a quadratic curve.

B. The PWA system for the two-gene oscillator

In system (1), the function expressing synthesis of the activator h R : R + → R + is simply:

h R (A; θ 1 ) = A 2 A 2 + θ 1 , (5)
while the function expressing synthesis of the repressor is given by a composition of two Hill functions h A : R 2 + → R + :

h A ((A, R); (θ 0 , θ 2 )) = A 2 A 2 + θ 0 (1 + R/θ 2 ) . (6)
Synthesis of the activator (5) can be straightforwardly approximated by s + R : R + → {0, 1}:

s + R (A; θ 1 ) =      1, if A 2 > θ 1 0, if A 2 < θ 1 (7)
The case of ( 6) is more complex to deal with. A multiplicative approximation of the form s + (A; θ 1 )s -(R; θ 2 ) does not faithfully represent the dynamics of the continuous system. So, we chose to approximate the activation part by an increasing step function with a switching value that depends on R, as follows s + A : R 2 + → {0, 1}:

s + A ((A, R); (θ 0 , θ 2 )) =      1, if A 2 > θ 0 (1 + R/θ 2 ) 0, if A 2 < θ 0 (1 + R/θ 2 ). (8)
This is a step function with a coordinate-dependent switching value. This induces a partition of the state space with nonrectangular regular domains, as analyzed in the next section. Applying these approximations, the PWA associated with (1) becomes: The phase space of the PWA system (9) will be divided in domains by the two threshold hyperplanes A 2 = θ 1 and

dA dt = V A s + A ((A, R); (θ 0 , θ 2 )) -γ A A + r bas dR dt = V R s + R (A; θ 1 ) -γ R R. (9)
A 2 = θ 0 (1 + R/θ 2 ), R ∈ [0,V R /γ R ].
However, these domains are no longer rectangles, but instead are defined by a quadratic curve:

µ(A) = θ 2 (A 2 -θ 0 )/θ 0 . (10)
According to the original parameter set in Smolen, Baxter, and Byrne (1998) it holds that √ θ 1 = √ θ 0 , and with this condition the phase space is partitioned into three domains:

• D00 = {(A, R) ∈ R 2 + : A ∈ [0, √ θ 1 ], R ∈ [0, +∞)} • D10 = {(A, R) ∈ R 2 + : A ∈ [ √ θ 1 , θ 0 (1 + (V R /γ R )/θ 2 )], R > µ(A)} • D11 = {(A, R) ∈ R 2 + : A ∈ [ √ θ 1 , +∞), R ∈ [0, µ(A))}
It is not difficult to compute the focal points in each of these three domains and check that, for a large range of parameters, φ 00 , φ 01 ∈ D00 and φ 11 ∈ D11 (see also Fig. 2). Therefore, trajectories in the domain D00 converge to the stable solution φ 00 = (φ A 00 , φ R 00 ) = (r bas /γ A , 0). Similarly, trajectories in D01 are attracted by the focal point

φ 01 = (φ A 01 , φ R 01 ) = (r bas /γ A ,V R /γ R ) and trajectories in D11 converge to the stable solution φ 11 = (φ A 11 , φ R 11 ) = ((V A + r bas )/γ A ,V R /γ R ), see Fig. 2.
It is not difficult to see that this state space conformation prevents the existence of periodic oscillations.

A closer look shows that requiring √ θ 1 > √ θ 0 will shift the line A = √ θ 1 to the right of the plane and introduce a fourth region in the state space. Four domains D i j , i, j ∈ {0, 1} will be formed under the condition √ θ 1 > √ θ 0 :

• D 10 = {(A, R) ∈ R 2 + : A ∈ [0, √ θ 1 ], R ∈ [0, µ(A))} • D 11 = {(A, R) ∈ R 2 + : A ∈ [ √ θ 1 , +∞), R ∈ [0, µ(A))} Table I.
Focal points, and their location under conditions (a)-(c).

φ 10 = (φ A 10 , φ R 10 ) = ((V A + r bas )/γ A , 0) D 11 φ 11 = (φ A 11 , φ R 11 ) = ((V A + r bas )/γ A ,V R /γ R ) D 01 φ 01 = (φ A 01 , φ R 01 ) = (r bas /γ A ,V R /γ R ) D 00 φ 00 = (φ A 00 , φ R 00 ) = (r bas /γ A , 0) D 10 • D 01 = {(A, R) ∈ R 2 + : A ∈ [ √ θ 1 , +∞), R > µ(A)} • D 00 = {(A, R) ∈ R 2 + : A ∈ [0, √ θ 1 ], R > µ(A)}
The focal points in each of the four domains are given in Table I where, for an appropriate but still very large range of parameters, it is possible to place each focal point outside its domain. Indeed, if the focal point φ 00 in domain D 00 is shifted to the domain D 10 and similarly φ i j to the following domain in anti-clockwise sense, then the trajectories will be attracted to the corresponding focal point in the neighbor domain. In this way, trajectories are expected to cross the threshold hyperplanes and, since the vector field is uniquely defined in each domain, an oscillatory behavior will eventually be formed by the four trajectory parts.

III. ANALYSIS OF THE TWO-GENE PWA SYSTEM

We now choose suitable parameter conditions so that the focal points will be shifted as: φ 00 is in the domain D 10 , φ 10 ∈ D 01 , φ 11 ∈ D 01 and φ 01 ∈ D 00 , see Fig. 3 andTable I. For this conformation of the focal points, we will next show that system (9) does admit sustained periodic oscillations.

To prove the existence of a limit cycle we will compute the Poincaré map of system (9) and show that it has a unique fixed point. In each of the four domains D i j there is a map G i j that establishes a correspondence between an entry point at one boundary of D i j and a corresponding exit point after the solution crosses D i j . The Poincaré map is given by the composition of these four maps.

A. Parameter conditions

The configuration listed in Table I is obtained from the conditions on the parameters (a) to (f), stated below and assumed to hold throughout the remainder of this paper. For simplicity, let γ denote γ R /γ A and let δ > 0 be small. Condition (a) guarantees the existence of four regions (compare Figs. 2 and3), while (b) and (c) guarantee, respectively, that φ 11 ∈ D 01 and φ 00 ∈ D 10 :

(a) √ θ 0 < √ θ 1 and γ R < γ A , (b) √ θ 0 + δ < r bas γ A < √ θ 1 -δ (c) √ θ 1 + δ < V A + r bas γ A < θ 0 1 + V R /γ R θ 2 -δ 1 0 1 1.4
Figure 3. A trajectory of the PWA system following a periodic orbit with parameter values:

V A = 0.43 min -1 , V R = 0.3 min -1 , γ A = 0.2 min -1 , γ R = 0.18 min -1 , r bas = 0.3 min -1 , θ 0 = 0.7, θ 1 = 7.84, θ 2 = 0.06
The inequalities (d)-(f) will be applied below in Lemma 3 to determine the signs of the first and second derivatives of the crossing maps G i j , which will help to establish the existence of periodic orbits. Condition (d) guarantees that the map G 11 has a negative first derivative:

(d) -γ µ( √ θ 1 ) -φ R 11 √ θ 1 -φ A 11 > 2φ A 11 θ 2 θ 0 .
Conditions (e) establish the size of the first derivatives of the maps G i j at the point

x = √ θ 1 : (e) i. -γ µ( √ θ 1 ) -φ R i j √ θ 1 -φ A i j + 2 √ θ 1 θ 2 θ 0 > 1, i j ∈ {10, 01} ii. -γ µ( √ θ 1 ) -φ R i j √ θ 1 -φ A i j + 2 √ θ 1 θ 2 θ 0 < 1, i j ∈ {00, 11}.
Conditions (f) establish the sign of the second derivatives of G i j at the point

x = √ θ 1 : (f) i. (γ 2 + γ) µ( √ θ 1 ) φ A 11 - √ θ 1 > 2 θ 2 θ 0 (2γ √ θ 1 + φ A 11 - √ θ 1 ) ii. (γ 2 + γ) µ( √ θ 1 ) -φ R 01 ( √ θ 1 -φ A 01 ) 2 -4γ θ 2 θ 0 √ θ 1 √ θ 1 -φ A 01 > 2 θ 2 θ 0 iii. (γ 2 + γ) µ( √ θ 1 ) √ θ 1 -φ A 00 + 2 θ 2 θ 0 ( √ θ 1 -φ A 00 ) > 4γ √ θ 1 θ 2 θ 0 .
It can be easily checked that a set of parameters such as indicated in the caption of Fig. 3 verifies these conditions.

B. The first return map of the system

Lemma 1. The following functions characterize the correspondence between entry and exit points of D i j :

G i j (x) = √ θ 1 -φ A i j x -φ A i j γ R γ A (µ(x) -φ R i j ) + φ R i j (11)
for i, j ∈ {0, 1} referring to the corresponding domain D i j , and defined in the following sets:

• {(x, G 10 (x)) ∈ R 2 + : x ∈ (φ A 00 , √ θ 1 ), G 10 (x) < µ(x)} • {(x, G 11 (x)) ∈ R 2 + : x ∈ ( √ θ 1 , φ A 10 ), G 11 (x) < µ(x)} • {(x, G 01 (x)) ∈ R 2 + : x ∈ ( √ θ 1 , φ A 11 ), G 01 (x) > µ(x)} • {(x, G 00 (x)) ∈ R 2 + : x ∈ (φ A 01 , √ θ 1 ), G 00 (x) > µ(x)}
Proof. To obtain these functions, note that in each domain D i j , equations ( 9) can be written as

dA dt = γ A (φ A i j -A), dR dt = γ R (φ R i j -R).
Then, the trajectory crossing D i j from a point (A 0 , R 0 ) to a point (A, R) is given by:

A(t) = (A(t 0 ) -φ A i j )e -γ A (t-t 0 ) + φ A i j R(t) = (R(t 0 ) -φ R i j )e -γ R (t-t 0 ) + φ R i j (12)
and solving with respect to time yields the equalities

e -(t-t 0 ) = A -φ A i j A 0 -φ A i j 1 γ A = R -φ R i j R 0 -φ R i j 1 γ R . (13)
In addition, for i j ∈ {11, 00}, we have

A 0 = √ θ 1 and R 0 = µ(A 0 ). It follows that R 0 = G i j (A 0 ). Conversely, for i j ∈ {01, 10}, let R = µ(A) and A = √ θ 1 . Then it follows that R = G i j (A 0 ).
For a trajectory starting at the point P2 = ( Ã2 , R2 ) = ( Ã2 , µ( Ã2 )) and returning to another point P 2 , both at the boundary between domains D 11 and D 01 , the first return map of the system can be computed as follows.

Lemma 2. The Poincaré map of the system (9) is given by the function:

G(x) = G -1 11 • G 10 • G -1 00 • G 01 (x) (14) defined on the interval G : [ √ θ 1 , φ A 10 ) → [ √ θ 1 , φ A 10 ).
Proof. Applying the system's solution for each transition in the phase space from P2 to P 2 as it is described in the Proof of Lemma 1 and after some computation follows the Poincaré map of the system.

The first and second derivatives of the Poincaré map G, to be used in the next section, can be computed from those of the four crossing maps. The first derivative of each map G i j is given by:

dG i j dx = √ θ 1 -φ A i j x -φ A i j γ -γ µ(x) -φ R i j x -φ A i j + 2x θ 2 θ 0 (15)
and the second derivative is given by:

d 2 G i j dx 2 = √ θ 1 -φ A i j x -φ A i j γ (γ +γ 2 ) µ(x) -φ R i j (x -φ A i j ) 2 -4xγ θ 2 θ 0 +2 θ 2 θ 0 . (16)

C. Existence and uniqueness of limit cycle

By definition, the fixed points G(x * ) = x * of the Poincaré map indicate the number of periodic orbits of the PWA system (9), and their stability indicates the stability of the orbit. Our main result is thus:

Theorem 1. Assume conditions (a)-(f) hold. Then the first return map (14) of system (9) has a unique fixed point in the open interval ( √ θ 1 , φ A 0 ).
To analyse the fixed points of G, we will first compute the signs of the first and second derivatives of its components, and show that G is a strictly increasing and concave function. The following properties are easy to check: Lemma 3. Assume conditions (a)-(f) hold. For each function G i j , i, j ∈ {0, 1}, the following inequalities hold:

(i) dG i j dx > 0, i j ∈ {10, 01} (ii) dG i j dx < 0, i j ∈ {00, 11} (iii) d 2 G i j d 2 x > 0, i j ∈ {10, 00} (iv) d 2 G i j d 2 x < 0, i j ∈ {11, 01} (v) dG i j dx ( √ θ 1 ) > 1, i j ∈ {10, 01} (vi) dG i j dx ( √ θ 1 ) < 1, i j ∈ {00, 11}.
We can now state the monotonicity properties of G (see also Fig. 4): Lemma 4. Assume conditions (a)-(f) hold. Then map G defined as in ( 14) is an increasing and concave function. ), computed for the same parameters as in Fig. 3. Its fixed points are given by the intersection with the straight line y = x. The point √ θ 1 corresponds to an unstable fixed point of system (9) while x * corresponds to a locally stable periodic orbit.

Proof. First, from Lemma 3 (i) we have that G 11 , G 00 are strictly decreasing functions in their domains, and since they are invertible, G -1 11 and G -1 00 are also strictly decreasing. Conversely, by Lemma 3 (ii), the functions G 10 , G 01 and their inverse functions are strictly increasing in their domains. Therefore, both G -1 11 • G 10 and G -1 00 • G 01 are strictly decreasing, meaning that G is itself strictly increasing as a composition of two decreasing functions.

For the concavity, recall that the inverse f -1 of a strictly monotone convex (resp., concave) function f in an open domain D is concave (resp., convex) in its corresponding domain.

Since G 10 , G 00 are convex functions and G 11 , G 01 are concave (by Lemma 3 (iii) and (iv)), set f 1 = G -1 00 • G 01 and note that

f 1 = (G -1 00 ) • G 01 (G 01 ) 2 + (G -1 00 ) • G 01 G 01 > 0
and so f 1 is convex and decreasing as a composition of an increasing and a decreasing function. Similarly

f 2 = [G 10 • f 1 ] > 0
so it follows that f 2 is a convex and decreasing function.

Lastly, we have that

G = [G -1 11 • f 2 ] < 0 implying that G is concave. 1. Extremal values G( √ θ 1 ) and G(φ A 10 )
Before proving our main result, we will evaluate map G at the two endpoints of its interval of definition, that is √ θ 1 and φ A 10 . Let us first study the map at

x = √ θ 1 , which is equivalent to the point A = √ θ 1 , R = µ( √ θ 1
) at the boundary of the four regions. At this point, the system (9) must be defined as a differential inclusion and should be written as a convex combination of the vector fields f i j in each of the four domains D i j :

    dA dt dR dt     = 1 ∑ i, j=0 α i j f i j ( θ 1 , µ( θ 1 )) (17)
where ∑ 1 i, j=0 α i j = 1. It is easy to see that the zero vector is included in this convex combination, hence the point

( √ θ 1 , µ( √ θ 1 )
) is a Filippov-type equilibrium point of (9) [START_REF] Filippov | Differential equations with discontinuous right-hand side[END_REF]. Moreover, by definition of the partial crossing maps:

G i j (x) -φ R i j µ(x) -φ R i j = √ θ 1 -φ A i j x -φ A i j γ R γ A it follows that G i j ( √ θ 1 ) = µ( √ θ 1 ) for all i j, hence G( √ θ 1 ) = √ θ 1 and
√ θ 1 is a fixed point of the first return map. However, the combination of Lemma 3 (v) and (vi) implpies that G ( √ θ 1 ) > 0, which in turn implies that √ θ 1 corresponds to an unstable fixed point of (9).

At the other endpoint

x = φ A 10 , by definition of focal point, it follows that G(φ A 10 ) < φ A 10 , since φ A 11 = φ A 10 and for A 0 ∈ D 11 dA/dt = -γ A (φ A 11 -A) implies A(t; A 0 ) < φ A
10 for all t and A < A 0 .

Proof of Theorem 1

We want to show that there exists a unique

√ θ 1 < x * < φ A 10 such that x * = G(x * ). From the above discussion, G is strictly increasing and concave in [ √ θ 1 , φ A 0 ) and satisfies G( √ θ 1 ) = √ θ 1 and G(φ A 10 ) < φ A 10 .
Therefore G intersects the line y = x exactly once in the open interval.

To show that this fixed point x * is stable, i.e. dG/dx(x * ) < 1, consider the function

L(x) = G(x) -x on [ √ θ 1 , φ A 10 )
. This function has exactly two zeros on this interval, √ θ 1 and x * . By the Intermediate Value Theorem, L (x) has one zero at r ∈ ( √ θ 1 , x * ) and L is decreasing between r and x * + δ for some δ > 0 (since G(x) is continous and

G(x) < x in (x * , φ A 10 )). Therefore, L (x) = G (x) -1 < 0 and 0 < G (x) < 1 for all x ∈ (r, x * + δ ).

D. Focal points shift due to basal synthesis rate for the repressor

This Section generalizes system (9) by adding a basal production rate to the repressor equation. This result will be useful later on in Section IV C, where we need to characterize the shift in the focal points due to the network coupling. Define the following PWA system associated with (9):

dA dt = V A s + A ((A, R); (θ 0 , θ 2 )) -γ A A + r bas dR dt = V R s + R (A; θ 1 ) -(γ R + β )R + βC, (18)
with s + A , s + R , as defined in Section II B; β > 0 is a constant that controls the shift of φ R i j coordinates in the phase space and C > 0 expresses a basal term that remains constant. The other parameters are those of system (9).

The system domains are defined as in Section II C and the coordinates φ R i j can be computed as a shift of the corresponding focal points R-coordinates of system (9):

φ R 00 = φ R 10 = βC γ R + β φ R 01 = φ R 11 = V R + βC γ R + β (19)
Note that the φ A i j coordinates remain the same as in Table I. Additionally to parameter conditions (a)-(f), we next introduce conditions on β , for which the location of the focal points (the domain in which φ i j is located) is the same as in Table I); these conditions will imply that the system (18) also admits a unique periodic solution as shown in Section III C. Introduce first the notation:

a bas = r bas γ A , a max = V A + r bas γ A . (20)
For system (18), the following holds:

(i) φ 00 ∈ D 10 if and only if φ R 00 < µ(a bas ), (ii) φ 11 ∈ D 01 if and only if φ R 11 > µ(a max ). From ( 19), these inequalities become:

(i) β < β 1 (C) := µ(a bas )γ R /(C -µ(a bas )), whenever C > µ(a bas ), (ii) β < β 2 (C) := (µ(a max )γ R -V R )/(C -µ(a max )), when- ever C < µ(a max ).
Combining these two inequalities, yields the following condition on β :

β <      β 1 (C), if C > µ(a max ) β 2 (C), if C < µ(a bas ) min{β 1 (C), β 2 (C)}, if µ(a bas ) < C < µ(a max ) (21)
If β does not satisfy these conditions, it follows that the shift in Rcoordinates of the focal points induces a change in the dynamic of system (18) from sustained oscillations to stability: for C > µ(a max ), if the condition on β is not satisfied, then φ 00 shifts to its corresponding domain D 00 and becomes locally stable solution for the system (18). Analogously, for C < µ(a bas ) and the corresponding condition on β not satisfied, φ 11 changes location to D 11 and becomes steady state.

Lastly, for µ(a bas ) < C < µ(a max ), φ 00 shifts to D 00 for β > β 1 and φ 11 to D 11 for β > β 2 . In either of these three cases, the periodic solution disappears. This observation will be used in Section IV C to show stability of a network under some conditions.

IV. NETWORK DYNAMICS OF N COUPLED TWO-GENE OSCILLATORS

Motivated by the synchronization of mammalian circadian clocks in tissues and organs, we will consider a network of N two-gene oscillators interconnected according to three different topologies. Experiments have shown that the mammalian suprachiasmatic nucleus (SCN), that consists of neurons and is located in the brain, is the central synchronizer of the circadian rhythms, [START_REF] Ralph | Transplanted suprachiasmatic nucleus determines circadian period[END_REF]. In parallel, in cultured cells and tissue explants, there have been observed self-sustained and autonomous circadian clocks [START_REF] Balsalobre | A serum shock induces circadian gene expression in mammalian tissue culture cells[END_REF], [START_REF] Yamazaki | Resetting central and peripheral circadian oscillators in transgenic rats[END_REF]. The suprachiasmatic nucleus (SCN) is the central circadian synchroniser, sending signals to the rest of the brain and to the pheripheral clocks in the body. For the form of this interaction between the SCN and the peripheral clocks, we can consider diffusive exchanges between cells [START_REF] Silver | A diffusible coupling signal from the transplanted suprachiasmatic nucleus controlling circadian locomotor rhythms[END_REF].

In the context of the hierarchical organisation of the circadian system, we explore a network N two-gene oscillators coupled through diffusion terms. In a first case, we consider one PWA system (9) to represent the SCN and the other N -1 identical systems as peripheral clocks, leading to an "one-to-all" or "star" topology where oscillator i = 1 is connected to all others, as in Fig. 5. In a second case, we consider all oscillators to play identical roles, as in a homogeneous tissue, and study an "all-to-all" topology where each oscillator is connected to all others. As an intermediate case, we will consider randomly chosen connections, within certain parameters.

A. PWA oscillator network and location of focal points

Consider N identical oscillators of the form (9), with coordinates (A i , R i ) i ∈ {1, 2, . . . , N}, coupled through R i according to some topology L and a diffusion parameter κ:

dA i dt = V A s + A,R (A i , R i ; θ 0 , θ 2 ) -γ A A i + r bas dR i dt = V R s + A (A i ; θ 1 ) -γ R R i + κ ∑ j∈L (R j -R i ), i = 1, . . . , N, (22)
where the set L contains all systems j that are connected to i. It will be useful to define the vectors A = (A 1 , . . . , A N ) and R = (R 1 , . . . , R N ) , as well as the state of the coupled system as

X = (A 1 , R 1 , . . . , A i , R i , . . . , A N , R N ) ∈ R 2N ≥0 . We further define sh A,R = [s + A ((A 1 , R 1 ); (θ 0 , θ 2 )), . . . , s + A ((A n , R n ); (θ 0 , θ 2 ))]
and

sh A = [s + A (A 1 ; θ 1 ), . . . , s + A (A n ; θ 1 )]
to be N-dimensional Boolean vector functions, evaluated according to the state X. The coupled system ( 22) can be rewritten as:

dA dt = V A sh A,R -Γ A A + r bas dR dt = V R sh A -(Γ R + L)R, (23)
where Γ A , Γ R are diagonal N × N matrices with elements γ A , γ R > 0 respectively, and (by abuse of notation) r bas is an Ndimensional vector with elements r bas . L is the symmetric N × N Laplacian matrix of the connection topology, to be defined for each case.

Analogously to Section II C, we define the domains of system (23) as

D C = D c 1 × • • • × D c N ,
where c i belongs to {00, 01, 10, 11}, so that each D c i corresponds to a regular domain of the 2-dimensional system (9). The focal points for each domain D C are denoted by φ C = ( ÃC , RC ) and they are given by:

ÃC = V A Γ -1 A sh A,R [D C ] + Γ -1 A r bas RC = V R (Γ R + L) -1 sh A [D C ] (24)
where sh A,R [D C ] and sh A [D C ] are the synthesis rate vectors evaluated at domain D C . As before, the location of the focal points determines the dynamics of the system. In Section III C we have shown that the (single) two-gene oscillator (9) admits a unique periodic solution for certain parameter conditions. This periodic solution, which we call (ϕ A (t), ϕ R (t)), extends to naturally as a solution to the coupled system with ϕ i A (t) = ϕ A (t) and ϕ i R (t)) = ϕ R (t)) for all i. This can be seen by taking identical initial conditions for all systems i. Intriguingly, our analysis shows that diffusive coupling generates new locally stable steady states for the coupled system. Notice that only the coordinates RC depend on the coupling parameter κ, as well as on the number N of oscillators, while ÃC are independent of the coupling topology and can be computed for i = 1, . . . , N as follows

ÃC i = a bas , D C i ∈ {D 00 , D 10 } a max , D C i ∈ {D 11 , D 01 } (25)
where we used notation (20). This implies that regular domains D C whose product includes {D 01 , D 10 } for some i cannot contain their own focal point φ C . Conversely, for regular domains D C which are products of D 00 and D 11 only, the location of the focal points depends on the value of each RC i relative to the curve µ(A) (see Fig. 3), as stated in the next Lemma.

Lemma 5. For each oscillator i, the following hold: The inequalities (a) and (b) follow directly from the expressions (25) for ÃC and the geometry of the phase space. The result on φ C is an immediate consequence of (a) and (b).

(a) Let D c i = D 00 , RC i ∈ D 00 if and only if RC i > µ(a bas ), (b) Let D c i = D 11 , RC i ∈ D 11 if and only if RC i < µ(a max ). R 1 R 2 R 3 R 4 R 5 R 6 R N

B. Star Topology

This simplified interconnection scheme naturally assumes that SCN (represented by R 1 ) is connected to all other peripheral clocks, but the latter do not communicate among themselves, as in the network depicted in Fig. 5. The corresponding Laplacian matrix L is defined as

L = κ       N -1 -1 -1 . . . -1 -1 1 0 . . . 0 -1 0 1 . . . 0 . . . . . . -1 0 0 . . . 1       (26)
In each of the 4 N regular domains of the coupled system the coordinates of the focal points can be computed as in ( 24), where the matrix (Γ R + L) -1 is strictly positive and has the form

      Q S S • • • S S P T • • • T S T P • • • T . . . . . . . . . . . . S T T • • • P       (27)
where the entries are given by:

S = κ γ R (γ R + Nκ) (28) Q = γ R γ R (γ R + Nκ) + S T = κ γ R + κ S P = (γ R + κ) 2 + (N -2)γ R κ κ 2 T.
Based on Lemma 5, we can now show that new steady states may be generated due to the coupling scheme. In addition, using (28), we can estimate the number of these new steady states in terms of the coupling parameter and the number of oscillators in the network. Introduce the notation: Theorem 2. Consider system ( 23) in a star-coupled topology, with L defined as in ( 26). Assume that conditions (a)-(f) hold. Case (i) Define the positive integers:

D 00 × (D m 11 × D N-m-1 00 ) prm , D 11 × (D m-1 11 × D N-m 00 ) prm
m min = 1 TV R µ(a bas ) (29)
m max = 1 TV R µ(a max ) - (γ R + κ) 2 + (N -2)γ R κ κ 2 + 1 .
If m min < m max , then system (23) has exactly ∑ m max m=m min (N -1)/(m!(Nm -1)!) locally stable steady states corresponding to focal points

φ C ∈ D C = D 00 × (D m 11 × D N-m-1

00

) prm . Case (ii) Define the positive integers:

m min = 1 TV R µ(a bas ) - γ R κ (30) m max = 1 TV R µ(a max ) - γ R κ - 1 γ R + κ 1 T max γ R V R µ(a max ), 1 .
If m min < m max , then system (23) has exactly ∑ m max m=m min (N -1)!/((m -1)!(Nm)!) locally stable steady states corresponding to focal points

φ C ∈ D C = D 11 ×(D m-1 11 ×D N-m 00 ) prm . Case (iii) Let κ =              max γ R (µ(a max )γ R -V R ) V R -Nγ R µ(a max ) , γ 2 R µ(a bas ) V R -Nγ R µ(a bas ) , if N ≤ N * V R µ(a max ) -µ(a bas ) -γ R , if N > N * with N * = V R γ R µ(a bas ) -γ R (µ(a max )-µ(a bas )) V R -γ R (µ(a max )-µ(a bas ))
. Then no new steady states exist for any κ ≤ κ.

Proof. First, recall that a focal point φ C is a locally stable state if and only if it belongs to its own domain D C : since the eigenvalues of the matrices -Γ A and -(Γ R + L) have strictly negative real part, the Hurwitz condition holds, hence the focal points that lie in their own domains are locally asymptotically stable. Second, Lemma 5 provides the conditions required to guarantee φ C ∈ D C . Consider a general domain D C consisting of the product of m regions of the form D 11 and Nm of the form D 00 , in any order. Since oscillator 1 has a special role, we need to consider two cases, depending on whether c 1 = 00 or c 1 = 11:

Case (i) D C = D 00 × (D m 11 × D N-m-1

00

) prm . In this case, the vector sh A [D C ] has m entries equal to 1, and ( 27) and ( 24) lead to

RC 1 = mSV R RC i = (P + (m -1)T )V R , if sh A,i = 1, i = 2, . . . , N mTV R , if sh A,i = 0, i = 2, . . . , N.
By Lemma 5, the conditions to guarantee

φ C ∈ D C become mSV R > µ(a bas ) (P + (m -1)T )V R < µ(a max ) mTV R > µ(a bas )
From ( 28), it follows that S > T so, combining these inequalities leads to and (28) shows that P/T = ((γ R + κ) 2 + (N -2)γ R κ)/κ 2 which yields m min < m < m max . If there exists an integer m in this interval, this means that any choice of m oscillators in D 11 implies that the corresponding focal point belongs to D C . Now, there are several ways of choosing m oscillators among the N. In this first case, D c 1 = D 00 is fixed, so the formula is "N -1 choose m", that is (N -1)!/(m!(N -1m)!) ways for each m. Since there is one focal point for each of these choices and this is valid for any integer m in the interval [m min , m max ], obtain the given sum for the total number of new steady states.

1 TV R µ(a bas ) < m < 1 + 1 TV R µ(a max ) - P T ,
Case (ii) D C = D 11 × (D m-1 11 × D N-m 00 ) prm .
The coordinates RC of the focal points are now given by

RC 1 = (Q + (m -1)S)V R RC i = (S + P + (m -2)T )V R , if sh A,i = 1, i = 2, . . . , N (S + (m -1)T )V R , if sh A,i = 0, i = 2, . . . , N.
Applying Lemma 5 yields,

(Q + (m -1)S)V R < µ(a max ) (S + P + (m -2)T )V R < µ(a max ) (S + (m -1)T )V R > µ(a bas ).
Combining the first two inequalities

m < min 1 + 1 SV R µ(a max ) - Q S , 2 + 1 TV R µ(a max ) - P + S T , using (28) to obtain Q/S = 1 + γ R /κ, (S + P)/T = 2 + γ R /κ + 1/(κS) m < min 1 SV R µ(a max ) - γ R κ , 1 TV R µ(a max ) - γ R κ - 1 κS ,
and noticing that

T /S = κ/(γ R + κ) = 1 -γ R /(γ R + κ) yields m min < m < m max .
In this case, the first oscillator is fixed at D c 1 = D 11 so the formula becomes "N -1 choose m -1", that is (N -1)!/((m -1)!(Nm)!) for each m. We note that if the inequalities

QV R < µ(a max ) (31) (S + P + T )V R > µ(a max ) SV R > µ(a bas )
are satisfied and m min ≥ m max for the values defined in (29), then there is exactly one locally stable steady state for the system (23), that is of the form D C = D 11 × (D 00 ) N-1 . This derives from the inequalities satisfied in Case (ii), for m = 1.

Case (iii)

To prove the last statement of Theorem 2, observe that κ is the value of κ that renders m min > m max for the values as defined by ( 29) and ( 31), that is the intersection of

1 TV R µ(a max ) - P T + 1 < 1 TV R µ(a bas ) (32)
∩ QV R > µ(a max ) ∩ SV R < µ(a bas )
gives an interval for κ such that no focal point belongs to its own domain. Solving µ(a max )/TV R -P/T + 1 < µ(a bas )/TV R with respect to κ and using (28) yields the value of κc = V R /(µ(a max )µ(a bas ))γ R , which is constant ∀N. Solving SV R < µ(a bas ) with respect to N for κc to obtain

N * = V R γ R µ(a bas ) - γ R (µ(a max ) -µ(a bas )) V R -γ R (µ(a max ) -µ(a bas ))
, that is the maximum value of N for which, the value of κ is by solving {QV R > µ(a max )} ∩ {SV R < µ(a bas )} with respect to κ.

Example 1: N = 2. System (23) has 16 regular domains and the corresponding focal points are given in Table (II).

In this case, for a coupling constant κ = 0.25, and the other parameters as in Fig. 3, we have S = 2.0425, T = 1.1875, Q = 3.5131, P = 3.5131, µ(a bas ) = 0.1336, µ(a max ) = 1.0865, κ = 0.21 and the bounds (29) and (30): 0.375 < m < 1.0914 QV R = 1.05 < µ(a max ) and SV R = 0.6127 > µ(a bas ).

According to Theorem 2, m = 1 thus implying at most 1!/(1!0!) = 1 new steady state in D 0011 and 1!/(0!1!) = 1 new steady state in D 1100 , in agreement with Table II. The value of κ also agrees with simulations, see Fig. 7.

Example 2: N = 10. Similarly, the computations of the m bounds in Theorem 2 case (i) yield 1.4783 < m < 4.3013, thus Table II. Focal points of 2 systems coupled and their corresponding location in the domains for κ = 0.25.

Φ 1010 D 1111 Φ 1011 D 1111 Φ 1001 D 1100 Φ 1000 D 1110 Φ 1110 D 1111 Φ 1111 D 0101 Φ 1101 D 0100 Φ 1100 D 1100 Φ 0110 D 0011 Φ 0111 D 0001 Φ 0101 D 0000 Φ 0100 D 0000 Φ 0010 D 1011 Φ 0011 D 0011 Φ 0001 D 0000 Φ 0000 D 1010
m min = 2, and m max = 4. We expect ∑ 4 j=2 (9)!/(( j)!(9-j)!) = 246 steady states in domains of the form D 00 × (D 11 ) 4 × (D 00 ) N-5 . Applying the inequalities of case (ii) to obtain 0.7583 < m < 3.5813, with m min = 1 and m max = 3, we expect ∑ 3 j=1 (9)!/(( j -1)!(10j)!) = 46 steady states in domains of the form D 11 × (D 11 ) 3 × (D 00 ) N-3 . Our estimations are exact as it is illustrated by comparison with numerical computations in Fig. 6.

C. Star-coupled network stabilization

In this section we show that the star topology does not admit "mixed" dynamics: if any of the two-dimensional systems converges to a steady state, then all the network stabilizes.

Theorem 3. Consider system (23) in a star-coupled topology, with L defined as in (26). Assume that conditions (a)-(f) hold. If any 2D system k converges to a steady state (away from the domain boundaries), then the full network converges to a steady state.

Proof. Suppose that 2D system k = 1, converges a to steady state in a regular domain, Rk = constant. Then, from (25), there are two possible values for Ãk : Ãk = a bas , which implies sh A,R = 0 and sh A = 0 a max , which implies sh A,R = 1 and sh A = 1.

For simplicity we denote δ j = (sh A ) j ∈ {0, 1}.

By assumption, Rk must remain constant (or within an interval of length ε as small as desired), δ k remains fixed, and Ṙk ≈ 0, thus (22) implies:

κL kk Rk + γ R Rk ≈ V R δ k (t) + κ ∑ j∈L R j , (33)
Since oscillator k is only connected to system j = 1, from ( 24) and ( 26), it follows that L kk = 1 and ∑ j∈L R j = R1 . Solving (33) with respect to R1 gives :

R1 ≈ κ + γ R κ Rk - V R κ δ k (t) (34)
Since both Rk and δ k (t) are fixed, the latter implies that R1 = constant, and so (A 1 (t), R 1 (t)) converges to a steady state.

Next, for an oscillator i of the coupled network, with i = 1, k, from (22) it holds that

dR i dt = V R δ i (t) -(γ R + κ)R i + κ R1 . (35)
Since R1 remains constant and the coupling parameter κ is fixed, the term κ R1 acts like a basal term for all components R i , i = 2, . . . , N. This implies that systems (A i , R i ) become decoupled, and are described by the following equations:

dA i dt = V A s + A,R (A i , R i ; θ 0 , θ 2 ) -γ A A i + r bas dR i dt = V R s + A (A i ; θ 1 ) -(γ R + κ)R i + κ R1 , (36)
This implies that each i-th system ( 36) is equivalent to system ( 18) with β = κ, C = R1 and other parameters as in Fig. 3.

Since there are at least two oscillators in the network (k and 1) that converge to steady states, it holds that κ > κ (Theorem 2). From Lemma 5, it follows that either µ(a bas ) < R1 or R1 < µ(a max ). Recalling (21), we have that

R1 > µ(a max ) ⇒ β 1 ( R1 ) < γ R µ(a bas ) µ(a max ) -µ(a bas ) < κ, R1 < µ(a bas )⇒ β 2 ( R1 ) < V R -γ R µ(a max ) µ(a max ) -µ(a bas ) < κ,
with the κ expression for N sufficiently large.

For the range R1 ∈ (µ(a bas ), µ(a max )), note that β 1 (C) is a decreasing function in the interval (µ(a bas ), +∞), with lim C→µ(a bas ) β 1 (C) = +∞ and, conversely, β 2 (C) is an increasing function in the interval (0, µ(a max ) with lim C→µ(a max ) β 2 (C) = +∞. Therefore, in the interval C ∈ (µ(a bas ), µ(a max )), the two functions intersect at a single point Ĉ and it is easy to check that:

min{β 1 ( R1 ), β 2 ( R1 )} < β 1 ( Ĉ) = β 2 ( Ĉ) with β 1 ( Ĉ) = β 2 ( Ĉ) = V R -γ R (µ(a bas ) -µ(a max )) mu(a bas ) -µ(a max ) ≡ κ.
Thus, from ( 21), the i-th system (36) does not admit a periodic solution, meaning that all systems in the star-coupled network converge to a (locally) stable steady state.

D. All-to-all topology

An alternative configuration mimics the network of interconnections that might be found between the cells in a living tissue, where each cell communicate with all the others. This is a fully symmetrical network, since no oscillator can be distinguished from another. The corresponding Γ R + L matrix and its inverse are:

Γ R + L = κ(NI -11') + γ R I (37) (Γ R + L) -1 = 1 Nκ + γ R I + κ γ R (Nκ + γ R ) 11', (38)
where I is the N × N identity matrix, 1 is the vector with all entries equal to 1 and 11' is the matrix with all entries 1. As for the star topology, we can count the number of new steady states due to the coupling. Since all oscillators have an equal role, let (D m 11 × D N-m 00 ) prm denote any permutation of m D 11 and Nm D 00 domains.

Theorem 4. Consider system (23) coupled in an all-to-all topology, with L defined as in (37). Assume that conditions (a)-(f) hold. Define:

m min = γ R κ Nκ + γ R V R µ(a bas ) (39)
m max = γ R κ Nκ + γ R V R µ(a max ) -1 ,
If m min < m max , then system (23) has exactly ∑ m max m=m min N!/(m!(Nm)!) locally stable steady states corresponding to focal points

φ C ∈ D C = (D m 11 × D N-m 00 ) prm . Furthermore, let κ = 1 N V R µ(a max ) -µ(a bas ) -γ R ( 40
)
then none of the focal points corresponds to a steady state for any κ ≤ κ.

Proof. As in the proof of Theorem 2, we will provide conditions such that regular domains of the form D C = (D m 11 × D N-m 00 ) prm contain their own focal points. In this case, all oscillators play a similar role. By ( 24) and (37), the RC coordinates of the focal point of D C are given by

RC i =        V R κ γ R (Nκ + γ) m, if sh A,i = 0 V R 1 Nκ + γ R +V R κ γ R (Nκ + γ R ) m, if sh A,i = 1.
By Lemma 5, the conditions to guarantee that φ

C ∈ D C V R κ γ R (Nκ + γ R ) m > µ(a bas ) V R 1 Nκ + γ R +V R κ γ R (Nκ + γ R ) m < µ(a max ),
immediately yield the values (39). Since the interconnection scheme is fully symmetric, the maximal number of new steady states equals all the possible combinations of the form D C , that is N!/(m!(Nm)!), for each m min ≤ m ≤ m max . The value κ is the highest value that implies m min > m max in (39). It follows by solving the inequality with respect to κ: A fundamental difference arises between the star and the all-to-all topology: the upper bound on the diffusion constant decreases to zero with N in the all-to-all topology, while it remains constant in the star topology (Fig. 7). For large N, we can expect to avoid new steady states in the star topology by assigning a small but finite and reasonable diffusion constant (on the order of the degradation constants); however, steady states will appear for sufficiently large N in all-to-all topology, possibly contributing to disrupt synchronization of the oscillators to the periodic orbit. In any case, the periodic orbit is preserved: in the last section we give the basin of attraction for the periodic solution. Therefore, considering clusters of cells connected in all-to-all topology, for a reasonable number of cells in each cluster the oscillations are preserved regardless the large number of the steady states generated by the coupling, see Fig. 8.

γ R (γ R + Nκ)(µ(a max ) -µ(a bas )) -V R γ R < 0.

E. Random topology

Star and all-to-all topologies are standard network architectures in graph theory with numerous applications. When it comes to biology, in particular to intercellular coupling, it becomes important to explore a larger class of interconnection structures.

In the work of Hafner, Koeppl, and Gonze (2012), a comprehensive mathematical model of the circadian clock is used to simulate the cells inside the SCN with different types of interconnections. The results of this study suggest that cells with problematic oscillatory behaviour synchronize more efficiently in a random topology than in other topologies tested. Motivated by this result, we introduce a random topology for the network (23), with degree of connectivity d. The cell i is randomly bidirectionally connected to cell j, with i, j = 1, . . . , N, i = j, and each cell (node) has exactly d connections (edges). We will denote this bidirectional connection as : i ↔ j. The elements of the Laplacian matrix L are defined as

l k j = κ(N -1 -d), j = k -κ, if j ↔ k, j = k,
The numerical results for the total number of new locally stable steady states in random connected networks of N oscillators, with degree of connectivity fixed for each network at 25% -50% of N, are illustrated in Fig. 6, along with the corresponding numbers for the two other topologies.

The comparison of the total number of steady states for the three topologies, shows that random connections generate more steady states than star topology but less than all-to-all. This is an intuitively reasonable result since, for the average degree of connectivity δ , for the three types of connections, we have that δ star < δ random < δ all-to-all .

F. Initial conditions to guarantee convergence to the periodic solution

The generation of new stable steady states due to the diffusive coupling may be a disrupting factor in a network where the goal is synchronization towards the biological rhythm present in each single cell. Thus, in addition to the lower bound on the coupling constant, which prevents the generation of new stable steady states, it is also crucial to define initial conditions that guarantee convergence to the periodic orbit. Due to the specific form of the domains of system (9), and consequently of system (22) (that are defined by the quadratic curve µ(A)), the explicit definition of the basins of attraction for the limit cycle and the locally stable steady states is not straightforward. Instead, numerical simulations were performed to identify the basins of attraction for the two types of attractors in three regions (zones) of the phase space, as defined in Table III. Table IV shows the percentage of initial conditions that converge to the steady states in each zone.

To facilitate the characterization of our results, first we define the line l θ (A) = cA + b, that passes through the point ( √ θ 0 , 0) and intersects the R-nullcline A = √ θ 1 and the R Table III. Phase space of system ( 22) partitioned in three regions to test convergence of initial conditions. Parameters as in Fig. 3,

l θ (A) = cA + b, c = 0.8482, b = -0.7082. Zones 1 {R i ∈ [0,V R /γ R ), A i ≤ l -1 θ (R i )} 2 {R i ∈ [0,V R /γ R ), A i > l -1 θ (R i ) ∩ ( √ θ 1 + ε, µ -1 (R i ))} 3 A i ∈ ( √ θ 1 + ε, µ -1 (V R /γ R )), R i ∈ [0, µ(A i ))
Table IV. Percentage of the initial conditions that lead the trajectories to the steady states, N = 10, κ = 0.25, and the initial states for each system to slightly vary as:

(A i , R i ) = α(A 1 , R 1 ), with (A 1 , R 1 ) ∈ zone j, α ∈ [0.8, 1.2], i = 2, . . . N and j = 1, 2, 3.
-line R = V R /γ R , with c = V R /γ R ( √ θ 1 - √ θ 0 ) and b = -V R √ θ 0 /γ R ( √ θ 1 - √ θ 0 )
. This line is used to define the zones of Table III.

The results of Table IV imply that trajectories starting in zone 1 have the highest probability to oscillate, with 100 % of the initial conditions tested leading to the periodic orbit. We note that the other two zones both intersect the basins of attraction of more than one attractor. Initial conditions in the region the three topologies tested. Indicatively for star topology, initial conditions in zone 1 are depicted in Fig. 9; trajectories starting from the marked superimposed points (A i , R i ) converge to the periodic orbit.

C = {(A, R) ∈ R 2N ≥0 }, R i ∈ [0,V R /γ R ), A i ≤ l -1 θ (R i ), i =

V. DISCUSSION

An improved synthetic oscillator with PWA analysis In this work we study an improved model of a synthetic biology oscillator using a piecewise affine approximation and prove existence and uniqueness of periodic solution for this PWA system, in a suitable region of parameters. This qualitative characterization of the system allows us to better understand the influence of each parameter on the dynamical behaviour, crucial factors to help guide system design and synthetic implementation. The two most important steps for the improvement of the system design are: first, separate the activation thresholds. The activator should act much faster on itself than on the repressor (θ 0 < θ 1 ). Second, the maximal rate for the activator synthesis (V A + r bas )/γ A should be smaller than the maximal repressor modulating factor 1 +V R /(γ R θ 2 ).

Interconnection of identical oscillators generates new dynamics A network of synthetic oscillators was next constructed by connecting N identical elements through diffusive coupling, organized according to different interconnection structures: one-to-all or star, all-to-all and random. All three forms of interconnection preserve the periodic orbit but acquire additional attractors: new locally stable steady states. We give the exact total number of the new stable states, which increases with the number N of the systems in the network. In addition, we provide a lower bound on the diffusion constant that prevents the generation of new steady states due to the coupling. We also identify a region of the state space which prevents convergence to any of the steady states, for the three forms of interconnections applied. Numerical simulations show that, in the case of convergence to the periodic orbit, the N systems synchronize both in period and in phase, even for a small value of the coupling constant. An interesting question to be addressed in future work is whether the coupled systems still synchronize if the individual systems are not identical but have different parameters and periods.

Network topologies and robustness of oscillations The three interconnection topologies lead to similar qualitative dynamical behaviours and raise the problem of the robustness of the periodic behavior, when the goal is, for instance, the synchronization of a cluster of cells to the same periodic rhythm. The analysis of the new steady states and their basins of attraction indicates that, for a large amount of initial conditions, the network elements will not synchronize to a similar asymptotic behaviour. Star topology has fewer new steady states and a finite lower bound on the diffusion constant that prevents locally stable steady states and guarantees the existence of periodic oscillations. Moreover, this lower bound has physiologically reasonable values and is independent of the number of cells in the network. In constrast, the all-to-all topology creates locally stable steady states for any value of diffusion constant, for sufficiently large N. This analysis indicates that a star interconnection, resembling a central synchronizer with peripheral clocks, has the capacity to more robustly maintain a synchronized oscillatory behaviour, while in an all-to-all network, representing identical cells in a homogeneous tissue, it will be difficult to synchronize to a periodic orbit as the number of cells increases. In random interconnections, reminiscent of a network of circadian oscillators inside the SCN, the total number of locally stable states lies in between the corresponding numbers of the two other topologies, making this type of connection more robust to oscillations than all-to-all topology.

Conclusion Our results show that the architecture of interconnections induces new dynamics in a network of identical oscillators, creating new steady states and interfering with synchronization properties, which can be undesirable in synthetic biology applications. Our analysis provides useful guiding lines for choosing network design and implementation showing that, by carefully adjusting the coupling strength or tuning initial conditions, it is still possible to synchronize the network of N biological oscillators, in all architectures. However, a star-shaped network has the potential to more robustly maintain oscillators synchronized.

Figure 2 .

 2 Figure2. An initial PWA system with activity thresholds θ 0 and θ 1 equal.

Figure 4 .

 4 Figure 4. The map G, defined in the interval [ √ θ 1 , φ A 10), computed for the same parameters as in Fig.3. Its fixed points are given by the intersection with the straight line y = x. The point √ θ 1 corresponds to an unstable fixed point of system (9) while x * corresponds to a locally stable periodic orbit.

Figure 5 .

 5 Figure 5. Star topology scheme: coupling N two-gene oscillators of the form (9), through variable R.

  to denote products of m domains of the form D 11 and N -m of the form D 00 , where only the first component D c1 is fixed, and the other components are any permutation of the remaining domains. For example, in the case N = 3, m = 2 we have: D 00 × D 11 × D 11 and D 11 × D 11 × D 00 , D 11 × D 00 × D 11 .

Figure 6 .

 6 Figure 6. The number of new steady states (log scale) generated by each topology.

Figure 7 .

 7 Figure7. The upper bound κ on the diffusion constant that guarantees no new steady states for κ < κ.

  Figure 8. Activity in time of N = 200 coupled two-gene oscillators of the form (9) in all-to-all topology, initial conditions in zone 1, κ = 0.25 and other parameters as in Fig. 3. Trajectories converge to the periodic solution.

Figure 9 .

 9 Figure 9. For all initial conditions tested in zone 1 (blue + symbols), trajectories converge to the periodic orbit. Here are shown the numerical tests for star topology with N = 10.

ACKNOWLEDGMENTS

We gratefully thank Luca Scardovi for many interesting and useful discussions on the coupling and synchronization of oscillators, during his extended visit to Inria. The authors were partly supported by the French national agency for research through project ICycle ANR-16-CE33-0016-01.

DATA AVAILABILITY

The data that supports the findings of this study are available within the article.


Rémi Vezy
email: remi.vezy@cirad.fr

Raphaël Perez

François Grand

Jean Dauzat

Light exchanges in discrete directions as an alternative to raytracing and radiosity

Keywords: ARCHIMED, light interception, scattering, photosynthesis, FSPM

Introduction

Light modelling at the scale of organs is essential to account accurately for the complex interactions between biophysical processes such as photosynthesis, stomatal conductance and energy balance. Yet, the calculation of radiative exchanges at fine scales is computationally-intensive and it remains a hindrance to a widespread use of FSPMs despite advances in light modelling using either radiosity (Chelle and Andrieu, 1998) or raytracing (Bailey, 2018). This study shows that simplifications based on the discretization of radiative fluxes allow processing radiative exchanges in a natural environment while maintaining good accuracy on the simulation of biophysical processes such as carbon assimilation.

Material and Methods

The present study is based on biophysical simulations performed using the ARCHIMED model. Incident radiation is depicted as a set of specular fluxes (i.e. parallel rays) in discrete directions using the sun direction for direct radiation and predefined "turtle" directions for the diffuse radiation. The "turtle" directions are obtained by splitting the sky hemisphere into sectors of equal solid angle (Dauzat et al, 2001). Optionally, direct radiation can be distributed in neighboring "turtle" sectors (turtle only). For each direction, the scene is projected on an image plane and the interception of incident light is deduced from rasterized pixel projections. Additionally, Z-Buffering gives the overlay of scene objects and, in this regard, pixels can be viewed as rays traced from outside down to the ground level. Light scattering can thus be processed similarly to raytracing. In the case of Lambertian objects, we further assume that all rays scattered by an object carry the same energy whatever the "turtle" direction. Net assimilation (An) is calculated with Farquhar's model (Farquhar et al. 1980), stomatal conductance with Medlyn's model (Medlyn et al. 2011) and the leaf temperature is found by solving the energy balance of the system. Simulations are run on a dense three-dimensional scene including two palms (Elaeis guineensis) with the following configuration: latitude= 15°, Day of year 71, time steps of 30mn, clearness index Kt= 0.5. A "toricity" option is used to generate a virtually infinite canopy. The number of "turtle" directions is set to 6, 16, 46 or 136. The sun position is either integrated into the turtle or separately computed. The pixel density ranges from 341 to 6821 pixels m -2 . The reference outputs are obtained with the highest number of directions and pixels.  Scene metrics: plot= 15.9m*9.2m,meshes= 24 863,triangles= 571 934,LAI= 3.2,leaflets= 24 493 Results and Discussions Figure 1 (left) illustrates the effect of the number of discrete light directions on the estimation of biophysical processes in comparison with the reference of 136 directions. Sampling the sun direction provides best results since direct radiation largely contributes to the PAR irradiance, the energy load of leaflets and, finally, their assimilation. Bias remain low when the sun direction is not sampled except when the number of "turtle" directions is decreased to six. The dispersion of residuals remains quite limited for 46 directions, meaning that reliable values can be obtained at leaflet scale for such configuration. Figure 1 (right) shows that a low pixel density (682 pixels m -2 , i.e. 50 000 pixels) is sufficient to get a relatively unbiased estimation of carbon assimilation at plot level, but a higher density is necessary to get reliable estimation at leaflet scale.

Figure 1: Evaluation of the error induced by a reduction in number of directions (Left), or a reduction of the number of pixels (Right) for the intercepted photosynthetically active radiation (PAR), absorbed energy (PAR + near infrared) and net carbon assimilation (An) at the leaflet scale for a palm plot. Values are presented relative to the reference simulation shown as the first value on left, i.e. 136 directions on the left plot (500 000 pixels), 1000 pixels on the right (46 directions, turtle only). Red color is used for a simulation with a precise computation of the sun position, and blue for an integration of the sun position in the turtle.

The reference configuration in the left pane of Fig. 1 generates 68.5M rays for each time step and, since several hits are recorded per ray (6 on average) this generates about 5 sub-rays that are used for the calculation of light scattering. Running the complete simulation with the reference configuration from the right pane of Fig. 1 lasts ~3.4 min for each time step1 (23M rays). This time can be decreased to only 2 seconds per step by storing partial scene illumination for each direction, but this preliminary step can be time-consuming, mainly during the multiple scattering for the PAR and NIR ranges. A considerable shortening is expected by treating light exchanges using directional form factors between pairs of objects instead of propagating scattered light by individual rays.

Conclusion

Using discrete ordinates allows performing accurate and unbiased simulations of light interception. Biases arise when decreasing the number of directions but with limited consequences on carbon assimilation. Larger biases occur when pixel density is too low to sample correctly individual leaflets. A configuration with 46 turtle directions for depicting both direct and diffuse radiation and a pixel density of 682 pixels m -2 allows fast computations while providing sufficient information to get precise light budget at fine scales.

Tests on a computer with 6 cores, Intel Xeon W2133 3.60 GHz, RAM 32Go.