Odile Burckard

Michèle Teboul

Franck Delaunay

Madalena Chaves
email: madalena.chaves@inria.fr

Cycle dynamics and synchronization in a coupled network of peripheral circadian clocks

The intercellular interactions between peripheral circadian clocks, located in tissues and organs other than the suprachiasmatic nuclei of the hypothalamus, are still very poorly understood. We propose a theoretical and computational study of the coupling between two or more clocks, using a calibrated, reduced model of the circadian clock to describe some synchronization properties between peripheral cellular clocks. Based on a piecewise linearization of the dynamics of the mutual CLOCK:BMAL1 / PER:CRY inactivation term, we suggest a segmentation of the circadian cycle into six stages, to help analyse different types of synchronization between two clocks, including single stage duration, total period, and maximal amplitudes. Finally, our model reproduces some recent experimental results on the effects of different regimes of time-restricted feeding in liver circadian clocks of mice.

Introduction

Most organisms have evolved a circadian timing system to synchronize their physiology and behavior with the Earth's light/dark cycle. The core mechanism underlying circadian rhythms is a molecular clock governed by an oscillatory gene network present in virtually all cells and which translates the external time information into optimally phased molecular and cellular rhythms [START_REF] Takahashi | Transcriptional architecture of the mammalian circadian clock[END_REF]. Environmental or genetic disruption of the circadian clock leads to multiple metabolic alterations in animal models. In humans, misalignment between body time and social time caused for instance by rotating shiftwork or work at night is also a risk factor for cardio-metabolic, inflammatory and malignant diseases [START_REF] Mukherji | The circadian clock and liver function in health and disease[END_REF]. Reciprocally, diseases perturb the circadian coordination thus potentially leading to a vicious cycle. A large body of experimental and clinical evidence suggests that high amplitude and correctly phased oscillations of core clock and clock-controlled genes is a hallmark of health while disrupted or desynchronized rhythms are associated with pathology [START_REF] Kohsaka | High-fat diet disrupts behavioral and molecular circadian rhythms in mice[END_REF][START_REF] Weger | Systematic analysis of differential rhythmic liver gene expression mediated by the circadian clock and feeding rhythms[END_REF]. Whether and how cells temporally coordinate together at the scale of the tissue and organ is therefore a question of primary importance to later implement preventive or therapeutic strategies aiming at reinforcing circadian rhythms.

It is now well established that neurons of the central clock located in the suprachiasmatic nuclei (SCN) of the hypothalamus oscillate in synchrony so that a single phase is observed at the tissue level. The coupling mechanism within the SCN requires intercellular communication through gap junctions and neurotransmitters [START_REF] Dewoskin | Distinct roles for gaba across multiple timescales in mammalian circadian timekeeping[END_REF]. Intercellular coupling is an important determinant of oscillator amplitude, robustness and entrainment by synchronizers. The SCN encompasses a pacemaker clock that hierarchically controls other peripheral circadian clocks [START_REF] Finger | Mammalian circadian systems: Organization and modern life challenges[END_REF]. It is a largely accepted idea that peripheral clocks function mostly as "slave" oscillators that damp in the absence of internal systemic synchronizers such as body temperature and glucocorticoid hormones.

However, this view has been recently challenged by in vivo experiments performed in mice performed by Schibler's laboratory [START_REF] Sinturel | Circadian hepatocyte clocks keep synchrony in the absence of a master pacemaker in the suprachiasmatic nucleus or other extrahepatic clocks[END_REF]. By using a sophisticated technology allowing longitudinal long-term bioluminescence recordings in free living animals, the authors could demonstrate that oscillations of the core clock gene Rev-erbα persist in the livers of animals lacking a functional SCN or circadian oscillators in all extra-hepatic organs. Therefore, it appears that phase coherence between hepatocyte clocks can be achieved independently of external Zeitgebers, suggesting that intercellular coupling may occur in the liver. Another recent work [START_REF] Finger | Intercellular coupling between peripheral circadian oscillators by tgfβ signaling[END_REF] by Kramer's laboratory using the human osteosarcoma U2OS cell line reinforces this hypothesis by providing evidence that peripheral cellular oscillators couple in vitro via paracrine factors. This study further identified transforming growth factor β (TGFβ) as a candidate signaling molecule contributing to the intercellular coupling between circadian oscillators although this remains to be established in vivo.

In this paper, we propose a theoretical and computational study of the coupling between two or more clocks, to investigate the question of intercellular coupling between peripheral circadian oscillators. Our objective is to analyze the dynamical evolution and coordination between the main processes involved in intercellular coupling, to better understand the mechanisms leading to oscillator synchronization. For this study, we will use a reduced mathematical model of the circadian clock previously developed in our group by Almeida et al. [START_REF] Almeida | Transcription-based circadian mechanism controls the duration of molecular clock states in response to signaling inputs[END_REF][START_REF] Almeida | Control of synchronization ratios in clock/cell cycle coupling by growth factors and glucocorticoids[END_REF], which recovers important properties such as the interactions with the cell cycle. This model is calibrated with data from mouse fibroblast cells cultured in regular medium [START_REF] Feillet | Phase locking and multiple oscillating attractors for the coupled mammalian clock and cell cycle[END_REF], so we can assume that it is a faithful model of an autonomous peripheral clock.

Mathematical models have been developed specifically for the SCN clock [START_REF] Forger | A detailed predictive model of the mammalian circadian clock[END_REF], for other organs such as the liver [START_REF] Woller | A mathematical model of the liver circadian clock linking feeding and fasting cycles to clock function[END_REF] or pancreatic beta-cells [START_REF] Woller | Modeling clock-related metabolic syndrome due to conflicting light and food cues[END_REF], but more often for "generic" type cells, blending literature data from various sources [START_REF] Relógio | Tuning the mammalian circadian clock: robust synergy of two loops[END_REF][START_REF] Korenčič | The interplay of cis-regulatory elements rules circadian rhythms in mouse liver[END_REF][START_REF] Almeida | Transcription-based circadian mechanism controls the duration of molecular clock states in response to signaling inputs[END_REF]. Synchronization of neurons in the SCN is studied through a general model in [START_REF] Kim | Molecular mechanisms that regulate the coupled period of the mammalian circadian clock[END_REF] and further detailed in [START_REF] Dewoskin | Distinct roles for gaba across multiple timescales in mammalian circadian timekeeping[END_REF], showing the coupling mechanism requires intercellular communication through gap junctions and neurotransmitters.

Synchronization of clocks in other organs or tissues has been studied essentially as a result of SCN control or by entrainment to external cues such as light or food intake [START_REF] Woller | A mathematical model of the liver circadian clock linking feeding and fasting cycles to clock function[END_REF][START_REF] Woller | Modeling clock-related metabolic syndrome due to conflicting light and food cues[END_REF]. In this study, we suggest a segmentation of the circadian cycle into six stages based on a piecewise linearization of our model [START_REF] Almeida | Control of synchronization ratios in clock/cell cycle coupling by growth factors and glucocorticoids[END_REF], to help characterize the effect of coupling on the synchronization mechanism. We find that clock synchronization happens at distinct levels: for instance, clocks with slightly differing periods can synchronize both their amplitudes and periods but maintain distinct duration for each of the six cycle stages. For a cluster of neighbouring cells, our coupling model reproduces some of the recent results on mouse hepatocytes [START_REF] Sinturel | Circadian hepatocyte clocks keep synchrony in the absence of a master pacemaker in the suprachiasmatic nucleus or other extrahepatic clocks[END_REF], suggesting a relevant role for peripheral clock coupling.

Coupling and synchronization of two clock models

The circadian clock model used in this paper was developed in our group [START_REF] Almeida | Transcription-based circadian mechanism controls the duration of molecular clock states in response to signaling inputs[END_REF][START_REF] Almeida | Control of synchronization ratios in clock/cell cycle coupling by growth factors and glucocorticoids[END_REF] and focuses on the transcriptional details of the two main feedback loops present in the circadian clock mechanism [START_REF] Takahashi | Transcriptional architecture of the mammalian circadian clock[END_REF][START_REF] Relógio | Tuning the mammalian circadian clock: robust synergy of two loops[END_REF][START_REF] Korenčič | The interplay of cis-regulatory elements rules circadian rhythms in mouse liver[END_REF], both centered on the protein complex CLOCK:BMAL1, next briefly described. In this reduced model, paralogous proteins are not distinguished so that for instance the REV variable accounts for both REV-ERBα and REV-ERBβ. The complex CLOCK:BMAL1 promotes transcription of the genes Per and Cry, whose corresponding proteins PER and CRY dimerize to form another complex, PER:CRY. Subsequently, upon nuclear translocation, PER:CRY interacts with CLOCK:BMAL1 and the resulting complex is exported from the nucleus, thus blocking CLOCK:BMAL1 transcriptional activity and closing the first negative feedback loop. In parallel, the complex CLOCK:BMAL1 also promotes transcription of the genes Dbp and Rev-erbα. The protein DBP binds to the D-box clock controlled element located in the Per and Rev-erbα promoters to activate their transcription. Finally, protein REV-ERBα binds to the R-box in the Bmal1 promoter, thereby blocking synthesis of the protein BMAL1 and closing the second negative feedback loop.

The reduced form of the model studied in [START_REF] Almeida | Control of synchronization ratios in clock/cell cycle coupling by growth factors and glucocorticoids[END_REF] describes the dynamics of the protein complexes CLOCK:BMAL1 (denoted B), DBP (denoted D), REV (denoted R), and PER:CRY (denoted P ). The simplified interaction diagram is depicted in Fig. 1 and the ordinary differential equations model from [START_REF] Almeida | Control of synchronization ratios in clock/cell cycle coupling by growth factors and glucocorticoids[END_REF] can be written as: 4 is a decreasing function. The model was calibrated with REV-ERBα experimental data from [START_REF] Feillet | Phase locking and multiple oscillating attractors for the coupled mammalian clock and cell cycle[END_REF] and its concentration units are in percentages of the REV-ERBα data mean value, which was used to normalize the data (see [START_REF] Almeida | Transcription-based circadian mechanism controls the duration of molecular clock states in response to signaling inputs[END_REF]). The set of estimated parameters can be found in Table 2 of [START_REF] Almeida | Control of synchronization ratios in clock/cell cycle coupling by growth factors and glucocorticoids[END_REF] and is used throughout this paper (unless stated otherwise):

Ḃ = V R h -(R) -γ B BP, Ḋ = V B B -γ D D, (1)
Ṙ = V D D -γ R R, Ṗ = V D D -γ B BP, where h -(R) = k 4 R k 4 R +R
V R = 44.4%.h -1 , k R = 80.1%, V B = 0.142%.h -1 , V D = 19%.h -1 γ R = 0.241h -1 , γ D = 0.156h -1 , γ B = 2.58h -1 , (2)
with a period of τ 1 = 24.5 hours.

Coupling between two peripheral clocks

To study the coupling between two or more oscillators, we will include two new observations into our model (1): the variation of the oscillator's period depending on expression level of gene Cry2 [START_REF] Zhang | A genome-wide rnai screen for modifiers of the circadian clock in human cells[END_REF] and the role of gene Per2 in intercellular clock communication [START_REF] Finger | Intercellular coupling between peripheral circadian oscillators by tgfβ signaling[END_REF]. Although the proteins CRY and PER are not specifically represented in system (1), their complex PER:CRY is described by variable P . Therefore, both new observations will be modeled by modifying equation Ṗ in two different ways, as follows.

To represent circadian period variability among a population of cells, we will follow experimental results showing that Cry2 knockout or knockdown cells have a longer period than wild type cells [START_REF] Börding | Generation of human cry1 and cry2 knockout cells using duplex crispr/cas9 technology[END_REF][START_REF] Zhang | A genome-wide rnai screen for modifiers of the circadian clock in human cells[END_REF]. To represent this in our model, we add a factor ε in the equation of PER:CRY, Ṗ = εV D D -γ B BP , where ε < 1 describes a lower availability of protein CRY2 Figure 1: Diffusive coupling between two circadian clocks through the complex PER:CRY. and thus a lower rate (εV D ) of complex PER:CRY formation (Table 2 in Appendix shows the correspondence between ε values and periods.)

Several earlier studies have shown that the Per2 gene is responsive to multiple signaling pathways both in the SCN and in peripheral organs ( [START_REF] Travnickova-Bendova | Bimodal regulation of mperiod promoters by creb-dependent signaling and clock/bmal1 activity[END_REF][START_REF] Balsalobre | Multiple signaling pathways elicit circadian gene expression in cultured rat-1 fibroblasts[END_REF]). Interestingly, recent experimental data [START_REF] Finger | Intercellular coupling between peripheral circadian oscillators by tgfβ signaling[END_REF] shows that paracrine signaling downregulates Dbp and upregulates Per2 expression, suggesting that gene Per2 is an early responder to intercellular clock communication. In this paper, we will assume that this communication is represented by a direct diffusive action on the complex PER:CRY variable P . Although a highly simplified representation, it avoids the introduction of new uncalibrated variables into the model, and allows us to continue our analytical study. The resulting coupled model is:

Ḃi = V R h -(R i ) -γ B B i P i , Ḋi = V B B i -γ D D i , (3) Ṙi = V D D i -γ R R i , Ṗi = ε i V D D i -γ B B i P i + k(P j -P i ),
In the next sections, we will analyze the effect of ε and k in the coupled network, assuming a weak coupling rate [START_REF] Guenthner | Circadian rhythms of PER2::LUC in individual primary mouse hepatocytes and cultures[END_REF] and a wider interval for the circadian periods ε min < ε < ε max .

Synchronization of two systems

In a first step we will study the synchronization of two systems (using ε 1 = 1 and ε 2 = ε), as a function of the factor ε and the coupling strength k. To do this, we compute the total difference between the variables,

|r| = r 2 B + r 2 D + r 2 R + r 2 P , with r B = B 1 -B 2 , r D = D 1 -D 2 , r R = R 1 -R 2 , r P = P 1 -P 2 .
The total difference |r| will increase with the perturbation (1 -ε) (mutation "strength") and decrease with the coupling strength k. Fig. 2 illustrates the opposite effects of period differences and coupling strength on the error. In the case of equal periods (i.e., 1 -ε = 0), the error converges to zero and the two clocks synchronize exactly, while the coupling strength cannot guarantee exact synchronization.

To further study the role of coupling in maintaining cell synchronization and the capacity of the cells to synchronize, we performed two types of simulations to calculate the amount of time Desynchronization (black curve): the two oscillators start from a point were their phase difference is smallest and the coupling is turned off. Synchronization (blue curves): the two oscillators start from a point where their phase difference is either 0.2τ 1 = 4.9h, 0.25τ 1 = 6.1h, 0.3τ 1 = 7.3h, or 0.4τ 1 = 9.8h, with τ 1 = 24.5h and the coupling is turned to k = 0.04.

In both cases, the oscillators are considered synchronized while their phase difference is lower than 0.15τ 1 = 3.6h.

needed to desynchronize two cells once coupling is disabled or, conversely, to calculate the time needed to synchronize cells that start with some phase difference. (Throughout this Section, we compute the oscillators phases relative to the crossing of threshold value k R by variable R, as it decreases.) Fig. 3 shows the number of circadian cycles needed to desynchronize or synchronize two cells, as one of the periods varies in terms of ε. Unsurprisingly, the time needed to desynchronize cells increases in an exponential form as the two oscillator periods become closer (black curve). More interestingly, the form of the synchronization curve highly depends on the initial phase difference between the two oscillators. For phase differences up to ≈ 8 hours, the number of cycles to synchronization decreases with the difference between periods. For low initial phase differences (about 5 to 6 hours), the time to synchronization is less than six circadian cycles for two cells with periods τ 1 = 24.5h and 24.8 < τ 2 < 25.5 (that is, 0.85 < ε < 0.95). For a higher initial phase difference of about 7.3h, the two cells need from 20 to 5 cycles to synchronize. Finally, for very high initial phase differences (around 10 hours, see lighter blue curve) the trend is inverted, since synchronization times are around five circadian cycles for 25.2 < τ 2 < 25.5 (that is, 0.85 < ε < 0.9) and then increase to about seven cycles for τ 2 = 24.8 (or ε = 0.95).

Our simulations evoke an observation on SCN clocks, for which it has been reported that the role of neuropeptide VIP (a factor implicated in cell communication) in synchronization is dose dependent, with a low VIP leading to synchronization but a high VIP leading to desynchronization [START_REF] An | A neuropeptide speeds circadian entrainment by reducing intercellular synchrony[END_REF].

Characterization of circadian cycle timing

In order to better understand how the coupling mechanism works, the effect of the period differences (characterized by ε) and the coupling strength k on the synchronization, we will divide the circadian cycle into six different stages, based on the dynamics of the CLOCK:BMAL1 / PER:CRY mutual inactivation. To do this, the nonlinear terms of model ( 1) are simplified to create a piecewise linear model [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: equilibria and their stability[END_REF] that partitions the state space into six different regions. The time spent by the orbit in each of these regions can then be explicitly calculated in terms of the system's parameters and switching conditions, providing estimates for the duration of each stage. A computational study shows the duration of each stage and their re-organization in the coupled network.

A piecewise linear clock model

There are two nonlinear terms in system [START_REF] Almeida | Control of synchronization ratios in clock/cell cycle coupling by growth factors and glucocorticoids[END_REF]. First, let's consider the decreasing Hill function 4 . This class of functions is often simplified by a step function [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: equilibria and their stability[END_REF]:

f (R) = V R k 4 R k 4 R +R
f (R) ≈ 0, R > k R V R , R < k R (4)
The second nonlinear term is the common degradation or inactivation of the complexes CLOCK:BMAL1 and PER:CRY, that is, -γ B BP . As the oscillator performs one full cycle, this term goes through several well marked conditions (see Fig. 4). Define B * and P * to be low levels of reference for the variables B and P , respectively. Consider also B high and P high to be high levels of reference for those variables. Let B max and P max denote their maximal values. These reference values define a family of intervals in B and P ([0, B * ],[B * , B max ] and [0, P * ], [P * , P high ], [P high , P max ]) where the CLOCK:BMAL1 and PER:CRY common inactivation term can be simplified by setting either B or P to be at their low or high reference values. Notice that while B is approximately constant at its low level B * , or when P is very large (P > P high ), the two terms γ B BP and f (R) have similar values, so they add up approximately to 0 in the B equation. This suggests the following simplification: To approximate γ B BP in the other cases, we adopt the following criterion: for intermediate levels of P , we assume that the inactivation term is dominated by the B dynamics with P constant, while for very low P , the opposite scenario holds. Taking into account the phase opposition between B and P , we expect B to be close to its maximum value when P < P * . We thus propose the following piecewise linear approximation for γ B BP :

γ B BP ≈ 0, R > k R and B < B * V R , R < k R and P > P high (5)
γ B BP ≈                0, R > k R , B < B * Stage 1 V R , R < k R , P > P high Stage 2 γ B BP * , R < k R , P * < P < P high Stage 3 γ B B max P, R < k R , P < P * Stage 4 γ B B max P, R > k R , B > B * , P < P * Stage 5 γ B BP * , R > k R , B > B * , P > P * Stage 6 (6)
These conditions cover all the state space and ensure positivity of the piecewise linear system composed of ( 1)-( 4)- [START_REF] Börding | Generation of human cry1 and cry2 knockout cells using duplex crispr/cas9 technology[END_REF]. To simulate this system we need to provide suitable values for the reference levels.

To estimate B max , consider the steady state obtained when transcription is activated and the degradation term is fixed at γ B BP * :

Ḃ ≈ V R -γ B BP * ⇒ V R -γ B B max P * = 0 ⇔ B max = V R γ B P * .
To estimate the low level P * , consider the cycle's phase opposition property and assume that both B and D reached their maximum values, and P is produced at a rate εV D D max :

Ḋ ≈ V B B max -γ D D ⇒ V B B max -γ D D max = 0 ⇔ D max = V B B max γ D Ṗ ≈ εV D D max -γ D B max P ⇒ εV D D max -γ D B max P = 0 ⇔ P = εV D D max γ B B max = εV D V B γ B γ D .
For a safer estimate we consider P * = P /4, so that

P * = εV D V B 4γ B γ D and B max = 4V R γ D εV D V B
In contrast, there is no obvious way to estimate a suitable non-zero value for B * , so we choose a fixed low value, B * = 0.01 for the simulations in this paper. Set B high = B max since a distinction is not required between them. Finally, P max is not used explicitly but an intermediate value for P must be fixed. In the current simulations, we use an appropriately high value within the system's parameters and set P high = V R . Fig. 5 shows oscillations of the 4 proteins of the piecewise linear model ( 1)-( 4)-(6, as well as the duration of each stage (black straight lines), with parameter set [START_REF] Almeida | Transcription-based circadian mechanism controls the duration of molecular clock states in response to signaling inputs[END_REF]. We observe that, depending on the parameter set, the system may not go through all six stages.

A partition of the circadian cycle

Through the simplification of the continuous model (1) into a piecewise affine model we analyse the circadian cycle as a sequence of six main stages, based on the repression of CLOCK:BMAL1 by REV and the dynamics of the mutual degradation/sequestration between CLOCK:BMAL1 and PER:CRY. Our stages can be compared to the common circadian time (CT) description [START_REF] Takahashi | Transcriptional architecture of the mammalian circadian clock[END_REF], which is based on the circadian transcriptional landscape in the mouse liver. This means that night hours (from CT12 to CT0/24) correspond to the mouse activity and feeding while the daylight hours (from CT0 to CT12) correspond to the resting and fasting phase. In this context, we find that the circadian night is divided into two stages, both characterized by a decreasing or low concentration of CLOCK:BMAL1 and the common degradation term either modulated by CLOCK:BMAL1 dynamics (stage 6) or at a very low level (stage 1).

The circadian daylight interval is divided into approximately three stages. An early day (stage 2) corresponding to a very high PER:CRY level, which imposes a high but saturated common degradation term, and CLOCK:BMAL1 de-repressed. The beginning of this stage would approximately correspond to the beginning of the circadian time count CT0. The second daylight stage (3) corresponds to CLOCK:BMAL1 augmentation as PER:CRY decreases with their degradation modulated by the dynamics of CLOCK:BMAL1. In the third daylight stage (4) PER:CRY is now at a low level while CLOCK:BMAL1 is still increasing, but their common degradation term is now modulated by PER:CRY dynamics. The end of this stage would correspond to middle circadian time (CT12). Finally, there is a day-to-night transition (stage 5), where REV has reached its switching value and repression of CLOCK:BMAL1 is turned on.

R > K R R < K R R < K R R < K R R > K R R > K R

Duration of each cycle stage in terms of the system's parameters

The advantage of the piecewise linear approximation is to allow the explicit computation of a solution for four of the six cycle stages and hence also an expression for these four stages duration in terms of the parameters. We can express the duration of stages 6,1,2, and 3 in terms of the parameters, but expressions for stages 4 and 5 are more complex to derive. However, in all our examples, stages 4 and 5 have very short duration (lower than about 2 hours).

For stage 6, we have the following equation for B: Ḃ = -γ B BP * whose analytical solution is: B(t) = e -γ B P * (t-t 0 ) B(t 0 ). Letting (t 3 , t 6 ) denote the times at which the cycle end stage 3 and 6, respectively, we assume that B(t 3 ) = B max = V R /(γ B P * ) and B(t 6 ) = B * . So, recalling

P * = εV D V B /(4γ B γ D ), obtain B(t 6 ) = e -γ B P * (t 6 -t 3 ) B(t 3 ) ⇒ t 6 -t 3 = - 4γ D εV D V B ln εB * V D V B 4V R γ D (7)
Since ε represents Cry2 knockout, as ε decreases the period of the system increases. This expression shows that duration of stage 6 also increases with the effectiveness of Cry2 knockout. In a similar way, we can compute the duration of stages 1, 2 and 3:

t 1 -t 6 = - 1 γ D ln D(t 1 ) -V B B * γ D D(t 6 ) -V B B * γ D (8)
t 2 -t 1 = - 1 γ D ln D(t 2 ) -V B B * γ D D(t 1 ) -V B B * γ D (9)
t 3 -t 2 = - 1 γ B P * ln B(t 3 ) -V R γ B P * B * -V R γ B P * (10)
Duration of stage 6 depends only on the parameters whereas duration of the other stages (k = 1, . . . , 5) depends also on the values reached by D and B at the beginning (D(t k-1 ), B(t k-1 )) or at the end (D(t k ), B(t k )) of stage k, that we don't necessarily know how to estimate.

Notice that stages 6, 1 and 2 follow each other in this order. Since the equation of D is the same for all stages, during stages 1 and 2, its solution is

D(t) = e -γ D (t-t 0 ) D(t 0 ) - V B B * γ D + V B B * γ D (11)
The lower level D * of D is reached when B = B * so we can estimate it :

Ḋ ≈ V B B * -γ D D ⇒ V B B * -γ D D * = 0 ⇔ D * = V B B * γ D So, as D(t 0 ) > V B B * γ D , D(t)
is decreasing during stages 1 and 2, so that D(t 6 ) > D(t 1 ) > D(t 2 ). Fig. 6 shows the duration of each stage as a function of ε. For parameter set (2), the piecewise linear system doesn't go through the stages 4 and 5. With respect to ε, stage 6 decreases strongly then stabilizes, stage 3 increases then stabilizes, stage 1 increases for very low ε then decreases. Stage 2 seems to increase close to linear. Let's focus now on the duration of the period: because of the strong decrease of stage 6, it decreases, then because of the increase of stage 2 and the stabilization of stages 3 and 5, it increases. More specifically, it is interesting to compare the model behavior for ε < 1 with Cry2 knockout cells: lower ε values represent more effective knockouts, with less available protein CRY2. Cry2 knockouts have higher periods amplitudes [START_REF] Guenthner | Circadian rhythms of PER2::LUC in individual primary mouse hepatocytes and cultures[END_REF] and our model shows that the period increases with the effectiveness of the knockout. Moreover, our model shows that the longer period is due to a longer duration of stage 6 corresponding to a longer end of the day, in agreement with the property that Cry2 deficient mice represent models of "eveningness" [START_REF] Destici | Altered phase-relationship between peripheral oscillators and environmental time in cry1 or cry2 deficient mouse models for early and late chronotypes[END_REF] (showing activity patterns later in the day).

Circadian stages and period duration in a mixed cell population

To study the coupling dynamics in a cell population with mixed periods, we now consider the continuous system (3) with ε 1 = 1 and ε 2 = ε, to reproduce two coupled oscillators with different periods. To further analyse the dynamics, we measure the duration of the six stages for each oscillator, according to the conditions summarized in Table 1. Fig. 7 shows the stage and period durations of the coupled system, as a function of ε, for two values of the coupling strength k.

As expected, a higher coupling constant allows synchronization of the periods for a wider interval around ε = 1. For ε < 1, this simulation evokes the mixture between Cry2 knockout and wild type cell populations performed for instance in [START_REF] Guenthner | Circadian rhythms of PER2::LUC in individual primary mouse hepatocytes and cultures[END_REF] or [START_REF] Finger | Intercellular coupling between peripheral circadian oscillators by tgfβ signaling[END_REF]. The period becomes longer as ε < 1 decreases (check thin dashed lines, uncoupled oscillator 2), as is the case in Cry2 knockout cells. The mixture of wild type and Cry2 knockout cells is represented by the (a) coupled oscillators solid and dashed thick lines, respectively (ε < 1). In this case, we see that the longer periods are shortened by up to 2 hours (dashed line) while the wild type cells remain at approximately the same period (solid line) in a small neighborhood of ε = 1, in good agreement with the recent observations by [START_REF] Finger | Intercellular coupling between peripheral circadian oscillators by tgfβ signaling[END_REF] (Figure 1 H,J in this reference). Although amplitudes and periods of both oscillators remain similar, notice that the duration of each stage may differ substantially from one oscillator to the other (compare same-color solid and dashed lines in Fig. 7): the two oscillators spend different amounts of time in different stages, while remaining synchronized on the full period. Fig. 8 confirms that, in a close neighborhood of ε = 1, the coupled oscillators have essentially the same periods but with a phase difference, in particular the variables' peaks are not reached simultaneously.

To further understand the mechanisms leading to the distinct duration of stages 1 and 2 in the two coupled oscillators, we first identified the entry and exit points from each stage and then identified and compared the dominant processes involved in each of the stages. In agreement with (6), stage 1 is completely defined by the value of BMAL1, since both entry and exit points from this stage correspond to the crossing of threshold B * by BMAL1, as this protein is either decreasing (entry to stage 1) or increasing (exit to stage 1) (see Fig. 9 (a)). The entry and exit to stage 2 are defined by REV and PER:CRY, respectively. The oscillators enter into stage 2 when REV crosses its activation level k R , while decreasing (see Fig. 9 (b)) and exit stage 2 when PER:CRY crosses its high level P high while decreasing (see Fig. 9 (c)). (c) Times series of variables B (red), R (green), and P (blue). The solid and dashed black horizontal lines denoted the length of the six stages. The vertical lines denote the beginning and ending of stages 1 and 2. (d)-(e) Times series of the main processes at play in the equation of the systems: synthesis of B and P (V R h -1 (R j ) and V D D j ), their mutual degradation term (γ B B j P j ), and the effect of coupling (plotted as V D D -k(P i -P j )).

(d) V R h -(R 1 ) B B 1 P 1 V R h -(R
V D D 1 V D D 1 -k(P 1 -P 2 ) B B 1 P 1 V D D 2 V D D 2 -k(P 2 -P 1 ) B B 2 P 2
We then compared the strength of the different processes during the two stages (see Fig. 9 (d)-(e)) It follows that during stage 1, repression of BMAL1 by REV is high and the mutual degradation of BMAL1 and PER:CRY is at a low level. The dominant process during stage 1 is clearly the production term of PER:CRY, V D D. Comparison of the dark and light blue curves of Fig. 9 (e), shows that the coupling term is also negligible during this stage. This last observation indicates that the two oscillators are operating in a close to uncoupled mode and therefore the longer period oscillator naturally shows a longer duration of stage 1.

Conversely, during stage 2, the opposite happens, with low repression of BMAL1 by REV and the mutual degradation of BMAL1 and PER:CRY at high level is now the dominant process (see Fig. 9 (d)). Also, during stage 2, the coupling term is at its highest strength, as can be seen by its effect on P:C dynamics: the coupling term increases (resp., decreases) the production rate of the shorter (resp., longer) period oscillator. This leads to the longer period oscillator decreasing PER:CRY levels faster and hence reaching the threshold P high faster than expected relative to the uncoupled oscillator. Inversely, the shorter period oscillator reaches P high slower than expected. Thus, during stage 2, the coupling term is largely responsible for decreasing the difference between the two periods by imposing a shorter duration of this stage for the longer period oscillator and a longer duration of this stage for the shorter period oscillator.

These conclusions are confirmed by Fig. 10, which show that the duration of stage 1 is roughly the same in the case of coupled or uncoupled oscillators in a small neighborhood of ε = 1. In contrast, the duration of stage 2 for oscillator 2 decreases significantly from the uncoupled to coupled situation (for ε < 1), while it increases slightly for oscillator 1. The synchronization dynamics between two oscillators along the circadian cycle is summarized in Fig. 11, which highlights the dominant processes at each stage.

Clock synchronization in response to time-restricted feeding

So far, we have studied the coupling between two oscillators. However, the study of cell periods in a mixture of wild type and Cry2 knockout cells suggests that hepatocytes are coupled in clusters of 5 to 10 cells (see [START_REF] Guenthner | Circadian rhythms of PER2::LUC in individual primary mouse hepatocytes and cultures[END_REF]). In liver lobules, cells are organised in parallel rows. We tested 5 distinct topologies (star, all-to-all, ring, rows and chain) to represent different patterns of exchanges between the cells. We report our results on star topology (see Fig. 12).

Figure 12: Star topology

In this section, we use our clock model to study the relevance of peripheral clock coupling in the response to different regimes of time-restricted feeding in liver circadian clocks of mice. Recent in vivo experiments by Sinturel et al. [START_REF] Sinturel | Circadian hepatocyte clocks keep synchrony in the absence of a master pacemaker in the suprachiasmatic nucleus or other extrahepatic clocks[END_REF] show that coherent oscillations are maintained in the liver of mice devoid of a functional SCN or having a clock exclusively in hepatocytes. We designed a computational study to reproduce the experimental protocol in these experiments, consisting of feeding mice (with and without SCN) ad libitum for about 20 days, then alternating periods of fasting and feeding (of 12 hours each) for about 20 days before feeding them again ad libitum. Circadian oscillations were longitudinally free-living animals thanks to PER2::luciferase reporter. The show that during ad libitum intervals, oscillations mice without SCN are damped. However, a 12 hour alternance between eating and fasting intervals permits rhythmic and increasing in peripheral organs.

food intake, increases and causes the pancreatic β-cells to secrete insulin [5]. Glucose and insulin levels are strongly linked and that is why we model the food intake with insulin (I(t)), added in the CLOCK:BMAL1 equation:

Ḃ = V R h -(R) -γ B BP + I(t) (12)
We consider a cluster of 8 cells with initial periods taken randomly between 23.01h and 26.32h (0.7 ≤ ε ≤ 1.2), coupled in a diffusive way through PER:CRY (k = 0.04) in a star topology. To model ad libitum food intake in SCN lesioned mice (arrhythmic feeding situation), we let insulin be at a constant but low level I(t) = 2%/h, throughout the first 20 days. To model time-restricted feeding, we add 2-hour pulses of insulin high amplitude I(t) = 13%/h, once every circadian cycle. As the peak of insulin lasts about 2 hours and using results of [START_REF] Almeida | Transcription-based circadian mechanism controls the duration of molecular clock states in response to signaling inputs[END_REF], the 2 hour pulse is applied from CT5 to CT7 during ten days and then from CT17 to CT19 during the next ten days. This schedule corresponds to the restricted access to the food during either the resting or activity phase respectively as in [START_REF] Sinturel | Circadian hepatocyte clocks keep synchrony in the absence of a master pacemaker in the suprachiasmatic nucleus or other extrahepatic clocks[END_REF]. Fig. 13 reports the average concentration of PER:CRY among the 8 cells, P (t) = The FFT analysis highlights several aspects. First, in the intervals of time-restricted feeding (see Fig 14 (b)), the cluster of uncoupled cells show two small amplitude peaks around 22h and 25h. Conversely, the cluster of coupled cells shows only one well defined peak around 23h, with at least a two-fold higher magnitude. This observation indicates that the coupled cluster has the capacity to produce a significant coherent signal under the time-restricted feeding regime, corresponding well to the experimental results of [START_REF] Sinturel | Circadian hepatocyte clocks keep synchrony in the absence of a master pacemaker in the suprachiasmatic nucleus or other extrahepatic clocks[END_REF]. A clear contrast is presented by the uncoupled cluster, unable to generate coherent oscillations and showing two low amplitude periods, conform to those present in the 8 individual cells.

A second observation is that the main peak representing ad libitum feeding (see Fig 14 (a)) is centered around 25h for the uncoupled cluster, while it has a lower period for the coupled cluster.

To complete these observations, we analyzed the amplitudes of the 3 main peaks of the FFT in terms of different values of insulin, for ad libitum feeding periods (values of insulin between 0.5 and 5) and time-restricted feeding periods values of insulin between 11 and 15.5) separately (see Fig 15), with and without coupling.

In the intervals of time-restricted feeding, the coupled cluster shows a peak with a significantly higher amplitude (see Fig 15(b), dashed red line) than the other peaks, suggesting that the cells synchronize over the same period. At the opposite, for uncoupled clusters all peaks have comparable amplitudes indicating that several different periods are similarly represented, thus uncoupled cells are weakly or not synchronized. In the intervals of ad libitum feeding, peaks have similar amplitudes, no matter whether there is coupling or not, suggesting weak synchronization.

Comparison between the amplitude of the peaks of the FFT between ad libitum feeding and time-restricted feeding for a coupled system (see Fig 15) shows that the main peak of timerestricted feeding is significantly higher than the main peak of ad libitum feeding, in agreement with the experimental results of [START_REF] Sinturel | Circadian hepatocyte clocks keep synchrony in the absence of a master pacemaker in the suprachiasmatic nucleus or other extrahepatic clocks[END_REF] (Fig. 2,F). There is not such difference between the peaks of ad libitum feeding and time restricted feeding for an uncoupled system, again suggesting a lack of synchronization in uncoupled cells.

Remarkably, our simulations are in agreement with experimental results for mice without SCN [START_REF] Sinturel | Circadian hepatocyte clocks keep synchrony in the absence of a master pacemaker in the suprachiasmatic nucleus or other extrahepatic clocks[END_REF], suggesting that small clusters of weakly coupled cells are important elements for achieving coherence and synchronization of peripheral clocks.

Discussion

This theoretical and computational study was motivated by recent experimental observations suggesting that intercellular coupling between circadian clocks occurs in peripheral organs such as the liver and contributes to maintain coherent oscillations in the absence of the SCN master clock. We constructed a network of coupled peripheral clocks to provide an early analysis of some aspects of the coupling dynamics.

Using a validated transcriptional model of the clock. Our model of the mammalian circadian clock was developed in previous work, calibrated with data from mouse fibroblasts, and validated on a number of applications (for instance, on the interactions between the mammalian clock and cell cycle). This model, with its validated set of physiological parameters, constituted a solid basis for our coupled network and its analysis, a clear advantage relative to generalist phase models used in many studies. Namely, this model allowed us to represent populations of cells with either longer or shorter periods, by tuning only one parameter (ε).

A partition of circadian cycle into six stages to characterize synchronization dynamics. To more closely analyze the coupling dynamics, we devised a partition of the circadian cycle based on the repression of BMAL1 and the common inactivation between the complexes CLOCK:BMAL1 and PER:CRY. In this way we obtain a circadian cycle with six stages, each representing one part of the circadian day and its length characterized by the system's parameters. This partition was used to study the coupling process in a mixture of wild type and longer period Cry2 -/-cells. Comparison of the clock periods with or without coupling shows that the longer period Cry2 -/--type cells shorten their periods in a coupled network, while wild type cells do not greatly change their periods, a result in very good agreement with the recent work of [START_REF] Finger | Intercellular coupling between peripheral circadian oscillators by tgfβ signaling[END_REF].

Our results highlight the fact that two coupled oscillators may be synchronized in period but have quite different duration for each stage (Fig. 9). Analysis of the main processes (Fig. 11) shows that, during stage 1 (circadian night) the coupling term is negligible, and the oscillators operate near to uncoupled mode, with the longer period oscillator naturally exhibiting a longer stage 1. Conversely, during stage 2 (circadian early day), the coupling is at its strongest level and succeeds in balancing the two periods by shortening the longer period oscillator.

Synchronization in a cluster of cells in response to time-restricted feeding. In the liver, communication between hepatocytes spreads to clusters of less than 10 neighbor cells, under weak coupling [START_REF] Guenthner | Circadian rhythms of PER2::LUC in individual primary mouse hepatocytes and cultures[END_REF]. We used our model to represent these cell clusters as observed in liver lobules, and simulated both uncoupled or coupled clusters with high period variability (between 23 and 26.5 hours), to study the average PER:CRY signal in response to either ad libitum or time-restricted feeding. Our simulations indicate that both timed food intake and intercellular coupling are determinant factors for the PER:CRY amplitude and period of oscillations (Fig. 13). Moreover, both factors are needed to reproduce the experimental observations, as shown by comparison of results with and without network coupling (Fig. 14 (a) and (b) and Fig. 15). The model with network coupling reproduces more faithfully than the one without coupling the recent experiments showing that hepatocyte clocks still synchronize in mice devoid of a functional SCN. This comparison of model and data thus suggests that coherence and synchronization of peripheral clocks is possible thanks to coupling between clocks.

Based on a physiological model of peripheral circadian clocks, our theoretical and computational results are among the first to analyze peripheral clock coupling. Our modeling approach provides a dynamical view of the coupling mechanism and illustrates how the coordination between the different stages leads to synchronization of two peripheral oscillators. Our results also show that intercellular coupling of peripheral clocks is a major factor in obtaining a coherent signal in response to feeding schedules.

Figure 3 :

 3 Figure3: Number of circadian cycles needed for desynchronization and synchronization of two oscillators. Desynchronization (black curve): the two oscillators start from a point were their phase difference is smallest and the coupling is turned off. Synchronization (blue curves): the two oscillators start from a point where their phase difference is either 0.2τ 1 = 4.9h, 0.25τ 1 = 6.1h, 0.3τ 1 = 7.3h, or 0.4τ 1 = 9.8h, with τ 1 = 24.5h and the coupling is turned to k = 0.04. In both cases, the oscillators are considered synchronized while their phase difference is lower than 0.15τ 1 = 3.6h.

Figure 4 :

 4 Figure 4: Protein concentrations over time and γ B BP and f (R) curves.

Figure 5 :

 5 Figure 5: Protein concentrations over time and stages duration, for the piecewise linear model. The stages are represented by the black straight lines, corresponding to the ordinates 50, 100, 150, 175, 200 and 225, respectively.

Figure 6 :

 6 Figure 6: Duration of each stage as a function of ε

Figure 7 :

 7 Figure 7: Duration of stages and period for system (3) as a function of ε. (a) k=0.04; (b) k=0.09. In both panels, the solid lines represent oscillator 1 (ε = 1) and the dashed lines represent oscillator 2 (ε = 1). The thin solid and dashed yellow lines represent the periods of the oscillators before coupling.

Figure 8 :

 8 Figure 8: Phase advance of oscillator 2 (ε = 1) over oscillator 1 (ε = 1), as an average over about 50 periods. (a) k = 0.04; (b) k = 0.09. Thin lines represent average plus standard deviation.

Figure 9 :

 9 Figure 9: Dynamics of two coupled oscillators with k = 0.04 and ε = 0.82. Oscillator 1 is shown in solid curves, oscillator 2 is shown in dashed curves. The uncoupled system periods are τ u 1 = 24.5h and τ u 2 = 25.5h, while the coupled periods are τ c 1 = τ c 2 = 24.4h. (a)-(c) Times series of variables B (red), R (green), and P (blue). The solid and dashed black horizontal lines denoted the length of the six stages. The vertical lines denote the beginning and ending of stages 1 and 2. (d)-(e) Times series of the main processes at play in the equation of the systems: synthesis of B and P (V R h -1 (R j ) and V D D j ), their mutual degradation term (γ B B j P j ), and the effect of coupling (plotted as V D D -k(P i -P j )).

Figure 10 :

 10 Figure 10: Comparison of duration of stages 1 and 2 for the coupled and uncoupled oscillators.

Figure 11 :

 11 Figure 11: Dominant processes throughout one circadian cycle. During stage 1 (circadian night) the two oscillators evolve in a weakly coupled mode, each oscillator following its individual natural period. Concentrations of PER:CRY and REV increase overall. During stage 2 (circadian early day), coupling is at its strongest level, leading to the adjustment of the two oscillators periods, with a shorter stage 2 duration for the longer period oscillator.

8 i=1 P i (t) 8 ,Figure 13 :

 8813 Figure 13: Average PER:CRY concentrations over time for ad libitum feeding (black bands) or time-restricted feeding (white bands), over a cluster of 8 cells with different periods (τ 1 = 24.11h, τ 2 = 23.48h, τ 3 = 23.55h, τ 4 = 26.29h, τ 5 = 26.17h, τ 6 = 26.15h, τ 7 = 24.95h and τ 8 = 24.58h). (a) Without coupling and (b) Star topology and k = 0.04.

Figure 14 :

 14 Figure 14: FFT analysis of Fig. 13. Blue curve: Without coupling. Orange curve: Star topology and k = 0.04. (a) Periods during ad libitum feeding; (b) Periods during time-restricted feeding.

Figure 15 :

 15 Figure 15: Amplitudes of the 3 main peaks of FFT for k=0 (solid lines) and k=0.04 (dashed lines) in terms of insulin values. (a) Ad libitum feeding; (b) Time-restricted feeding.

Table 1 :

 1 Characterization of the six circadian time stages based on the mutual CLOCK:BMAL1 and PER:CRY inactivation dynamics.

	Stage	1	2	3	4	5	6
	Circadian	night	early day	day	day	day/night night
	time		CT0→		→ CT12	transition
	CLOCK:BMAL1 low	low	increasing	high	high	decreasing
	(B)	B < B *	B > 0	B > 0	B > 0	B > B *	B > B *
	PER:CRY	increasing	high	decreasing	low	low	increasing
	(P)	P > 0	P > P high P * < P < P high P < P *	P < P *	P > P *
	Common	negligeable satured	modulated	modulated modulated modulated
	degradation		by P	by B	by P	by P	by B
	Repression of	on	off	off	off	on	on
	CLOCK:BMAL1

Acknowledgments

We are grateful to the referees for their constructive remarks and for a suggestion that led to the analysis of the dominant processes in Section 4. This work was supported in part by project UCA JEDI ANR-15-IDEX-01 and the LABEX SIGNALIFE (ANR-11-LABX-0028-01) from the program Investments for the Future of the French National Agency for Research and by project ICycle ANR-16-CE33-0016-01.