<figure><img src="image_1.png" /></figure>

HAL
open science

# Benchmark for quantitative characterization of circadian clock cycles

Odile Burckard, Michèle Teboul, Franck Delaunay, Madalena Chaves

► To cite this version:

Odile Burckard, Michèle Teboul, Franck Delaunay, Madalena Chaves. Benchmark for quantitative characterization of circadian clock cycles. BioSystems, 2025, 247, pp.105363. 10.1016/j.biosystems.2024.105363. hal-04834066

HAL Id: hal-04834066

https://inria.hal.science/hal-04834066v1

Submitted on 12 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Benchmark for quantitative characterization of circadian clock cycles

Odile Burckard<sup>a</sup>, Michèle Teboul<sup>b</sup>, Franck Delaunay<sup>b</sup> and Madalena Chaves<sup>a</sup>

## ARTICLE INFO

**Keywords:**

circadian clock
computational method
model comparison
cycle segmentation into stages
algorithm for circadian cycle characterization

## ABSTRACT

Understanding circadian clock mechanisms is fundamental in order to counteract the harmful effects of clock malfunctioning and associated diseases. Biochemical, genetic and systems biology approaches have provided invaluable information on the mechanisms of the circadian clock, from which many mathematical models have been developed to understand the dynamics and quantitative properties of the circadian oscillator. To better analyze and compare quantitatively all these circadian cycles, we propose a method based on a previously proposed circadian cycle segmentation into stages. We notably identify a sequence of eight stages that characterize the progress of the circadian cycle. Next, we apply our approach to an experimental dataset and to five different models, all built with ordinary differential equations. Our method permits to assess the agreement of mathematical model cycles with biological properties or to detect some inconsistencies. As another application of our method, we provide insights on how this segmentation into stages can help to analyze the effect of a clock gene loss of function on the dynamic of a genetic oscillator. The strength of our method is to provide a benchmark for characterization, comparison and improvement of new mathematical models of circadian oscillators in a wide variety of model systems.

## 1. Introduction

Since life first appeared on Earth, organisms have developed an internal clock mechanism called 'circadian clock' that they use to synchronize their physiology and behavior with the light/dark cycle and help them anticipate and prepare for daily changes in their environment. Circadian clocks form a robust network of periodic interactions between genes and proteins that provide a rhythmic regulation of a large number of biological processes. Circadian clocks drive oscillations at all scales, from molecular, subcellular, cellular, tissues and physiological processes, including for instance chromatin modification, mitochondrial dynamics, cell division cycle, sleep-wakes cycles, change in body temperature, hormonal secretion, etc. In mammals, the circadian clock mechanism (see Fig 1 for the mouse liver clock) is mainly based on two interlocked negative feedback loops (one between the complexes of proteins CLOCK:BMAL1 and PER:CRY and one between CLOCK:BMAL1 and REV-ERB), allowing Gouzé (1998); Novák and Tyson (2008); Richard and Comet (2011) circadian clocks to show sustainable oscillations over time in the absence of external synchronizers, with an approximated period of 24h. To accurately describe time within a free-running circadian system, that is in the absence of external synchronizers (e.g., constant darkness), the unit of time used is Circadian Time, such that 1CT is equal to the cycle period divided by 24 Jud, Schmutz, Hampp, Oster and Albrecht (2005). For mice previously entrained in a 12h/12h light/dark cycle before being kept in free-running condition, subjective day (sleep/rest/fasting phase in mice) extend from CT0 to CT12 while CT12 to CT24 correspond to subjective night (wake/activity/feeding phase in mice) Vitaterna, Takahashi and Turek (2001).

In the early subjective day, at CT0, the circadian cycle is in a poised state Takahashi (2017), while the protein CRY1 blocks CLOCK:BMAL1 activity (see top right of Fig 1). As CRY1 is gradually degraded, CLOCK:BMAL1, bound to the DNA response element Ebox (Enhancer box), promotes transcription of different genes (including the clock genes Period *Per1* and *Per2*, Cryptochrome *Cry1* and *Cry2* or *Rev-erbα*, *Rev-erbβ*, but also the clock controlled genes (Ccgs), needed for the regulation of various cellular, molecular, physiological and behavioral processes), with a peak observed around CT8 Takahashi (2017). *Rev-erb* mRNA expression peaks around CT6-10 Bugge et al. (2012). The resulting proteins REV-ERBα/β bind to the RRE (Rev Response Element) and block transcription of *Clock* and *Bmal1*, with a

<figure><img src="image_3.png" /></figure>

<footnote><sup>a</sup>Centre Inria d'Université Côte d'Azur, INRAE, CNRS, Macbes team, Sophia Antipolis, France</footnote>

<footnote><sup>b</sup>Université Côte d'Azur, CNRS, Inserm, iBV, Nice, France</footnote>

<footnote>*Corresponding author</footnote>

<footnote>odile.burckard@inria.fr (O. Burckard)</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1:** Representation of the core circadian clock mechanism of the mouse liver in the nucleus and temporal expression levels of the core clock components at the mRNA level (dashed line) and at the DNA binding level (solid line). White and light grey coloring background represent respectively subjective day (CT0-CT12) and night (CT12-CT24). REV-ERB, CRY and PER refer respectively to REV-ERBα/β, CRY1/2 and PER1/2. This figure combines a collection of data from various articles Bugge et al. (2012); Cao et al. (2021); Chiou et al. (2016); Gabriel et al. (2021); Jeong et al. (2022); Koike et al. (2012); Minami et al. (2013); Ripperger and Schibler (2006); Takahashi (2017); Ueda et al. (2005); Ukai-Tadenuma et al. (2011); Ye et al. (2014); Abe et al. (2022); Woller et al. (2016); Ko and Takahashi (2006); Ripperger et al. (2011); Atwood et al. (2011); Cho et al. (2012); Feng et al. (2011); Gatfield et al. (2009); Triqueneaux et al. (2004); Gréchez-Cassiau et al. (2008); Torra et al. (2000); Relógio et al. (2011); Balsalobre et al. (1998).</figcaption></figure>

binding peak observed around CT10 Solt, Wang, Banerjee, Hughes, Kojetin, Lundasen, Shin, Liu, Cameron, Noel et al. (2012); Banerjee, Wang, Solt, Griffett, Kazantzis, Amador, El-Gendy, Huitron-Resendiz, Roberts, Shin et al. (2014); Bugge et al. (2012).

The transition between subjective day and night, around CT12 is marked by PER and CRY nuclear translocation (see bottom of Fig 1), leading to CLOCK:BMAL1 repression during the subjective night. First, PER enhances local concentration of the casein kinase CK1δ at CLOCK:BMAL1-bound promoters in a CRY-dependent manner (interactions between CLOCK:BMAL1 and PER are facilitated by CRY), which leads to the phosphorylation of CLOCK. Hyperphosphorylation of CLOCK results in the displacement and sequestration of the complex CRY-CLOCK:BMAL1 from the Ebox and in the repression of transcriptional activity of CLOCK:BMAL1 Cao et al. (2021); Jeong et al. (2022). PER1, PER2 and CRY2 binding activity peak between CT15 and CT18 Takahashi (2017).

During the end of the subjective night, RORc, in competition with REV-ERB for the RRE binding, promotes transcription of *Bmal1*. Around CT22-24, the translation of *Clock* and *Bmal1* mRNA generates the CLOCK and BMAL1 partner proteins (see top left of Fig 1) which upon dimerization will bind to Ebox containing promoters. However, while PER is degraded, CRY1 expression is delayed and CRY1 binds to the CLOCK:BMAL1-E-box (around CT0) blocking CLOCK:BMAL1 transcriptional activity until the next cycle begins. CLOCK:BMAL1 is thus repressed during all subjective night, by displacement, sequestration and blocking. Cao et al. (2021); Chiou et al. (2016); Ye et al. (2014); Jeong et al. (2022).

Systems of ordinary differential equations (ODEs) are relevant for modeling biological oscillators, and over the past decade several circadian clock models have been developed using this mathematical framework Almeida, Chaves
and Delaunay (2020a,b); Relógio et al. (2011); Woller et al. (2016); Brown and Doyle III (2020); Hesse, Martinelli, Aboumanify, Ballesta and Relógio (2021). To analyze the periodic behavior of circadian clocks, we previously proposed Burckard, Teboul, Delaunay and Chaves (2022); Burckard and Chaves (2024) to partition the circadian time into stages defined according to the levels of abundance of the core clock components CLOCK:BMAL1, REV-ERB and PER:CRY, both in the cytoplasm and in the nucleus. Based on this segmentation, our group developed a piecewise linear model Burckard and Chaves (2024) to represent the dynamics of the circadian cycle within each stage, generating numerous oscillators with different parameters, and giving insights on temporal estimation of some circadian events or critical points for correct cycle progression. In this paper, we expand this cycle segmentation into stages and propose an Algorithm and a quantitative method to partition, characterize and compare any circadian cycle from both experimental data and diverse ODEs models. Finally, considering experimental results showing that partial or complete loss function of Cry2 lengthens the circadian period *in vivo* and *in vitro* Horst, Muijtjens, Kobayashi, Takano, Kanno, Takao, Wit, Verkerk, Eker, Leenen et al. (1999); Vitaterna, Selby, Todo, Niwa, Thompson, Fruechte, Hitomi, Thresher, Ishikawa, Miyazaki et al. (1999), we modeled this by decreasing the transcription rate of Cry in each model, and showed how the cycle segmentation can be used to characterize the new dynamic of an oscillator after a genetic perturbation.

Our quantitative method facilitates the study of a cycle and its time dynamics (sequence and duration of visited stages, time spent above some thresholds), highlighting biological properties recovered by the models, or detecting inconsistencies in the cycle, thus creating a benchmark which allows for validation, improvement and comparison of different models.

## 2. Characterization of circadian cycles

In this Section, we develop an automated method to partition and then characterize, study and compare the dynamics of any circadian clock cycle including the three core components CLOCK:BMAL1, REV-ERB and PER:CRY.

### 2.1. Circadian time partition into stages

Our method uses a cycle segmentation (introduced in Burckard and Chaves (2024) and based on the model from Almeida et al. (2020a)), whose stages are defined according to the levels of abundance of the three fundamental clock components involved in the two main feedback loops: CLOCK:BMAL1, REV-ERB and PER:CRY (denoted respectively *B*, *R* and *P* throughout this paper).

By setting one threshold on each of the variables *B* and *R* (respectively *B*<sub>high</sub> and *R*<sub>int</sub>) and two thresholds on the variable *P* (*P*<sub>low</sub> and *P*<sub>high</sub>), the state space is thus subdivided into 12 regions. For biological reasons (such as the phase opposition between *B* and *P*, or the order of protein peaks - *B*, *R*, *P* -, both induced by the two negative feedback loops), only eight stages are considered as valid. Any circadian cycle is expected to successively evolve through these eight stages, in a well-defined order (in the order from stage 1 to stage 8), and to go through each stage only once during a same cycle (see Table 1 and Burckard and Chaves (2024) for more details).

To represent the transitional states when *P* is at an intermediate level, two thresholds are set on the variable *P*. However, in the case where a cycle fails to progress through the eight stages in the expected order, the partition can be adapted and enlarged by removing the threshold *P*<sub>low</sub>, which corresponds to merging stages 5 and 6 together -renamed as stage 5-, 7 and 8 together -renamed as stage 7-, 4 and 9 together -renamed as stage 4- and 1 and 10 together -renamed as stage 1-. The state space is then divided into eight stages: stages 1, 2, 3, 4, 5, 7, 11, 12, of which only stages 1 to 7 are expected from a biological point of view.

A more coarse-grained partition can still be considered by removing *R*<sub>int</sub>, leading to a partition of the cycle into only four stages and permitting at least to check the phase opposition between *B* and *P*. In this case the following stages are thus merged together: stages 1 and 10 -renamed as stage 1-, stages 2 and 3 -renamed as stage 2-, stages 4 and 9 -renamed as stage 4-, stages 5, 6, 7, 8 -renamed as stage 5- and stages 11 and 12 -renamed as stage 11-. The corresponding expected stages order is 1,2,4,5. Notice that the new stages 1 and 4 are identically defined, but for convenience, to keep the "ascending" order through the stages, the property such that each stage is crossed only once in a same cycle, and to make a difference between the time before the peak of *B* and the one of *P*, we distinguish them.

### 2.2. An adaptive algorithm for thresholds estimation

In order to detect whether a given model exhibits the desired sequence of stages described above, we propose Algorithm 1 to partition automatically any circadian clock cycle (a detailed description of each step of Algorithm 1 is given in Appendix). The goal of our algorithm is to find suitable threshold values while satisfying two conditions:
<table><caption>**Table 1**
A partition of the circadian time into stages.</caption><thead><tr><th>Stage</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><td rowspan="3">Stages conditions</td><td colspan="4">PER:CRY peak</td></tr><tr><td>R > R<sub>int</sub><br>B < B<sub>high</sub><br>P<sub>low</sub> &lt; P &lt; P<sub>high</sub></td><td>R > R<sub>int</sub><br>B < B<sub>high</sub><br>P > P<sub>high</sub></td><td>R < R<sub>int</sub><br>B < B<sub>high</sub><br>P > P<sub>high</sub></td><td>R < R<sub>int</sub><br>B < B<sub>high</sub><br>P<sub>low</sub> &lt; P &lt; P<sub>high</sub></td></tr><tr><td colspan="4">CLOCK:BMAL1 peak</td></tr><tr><td rowspan="3">Stages conditions</td><td>R < R<sub>int</sub><br>B > B<sub>high</sub><br>P<sub>low</sub> &lt; P &lt; P<sub>high</sub></td><td>R < R<sub>int</sub><br>B > B<sub>high</sub><br>P < P<sub>low</sub></td><td>R > R<sub>int</sub><br>B > B<sub>high</sub><br>P < P<sub>low</sub></td><td>R > R<sub>int</sub><br>B > B<sub>high</sub><br>P<sub>low</sub> &lt; P &lt; P<sub>high</sub></td></tr><tr><td colspan="4">Unexpected stages</td></tr><tr><td>R < R<sub>int</sub><br>B < B<sub>high</sub><br>P < P<sub>low</sub></td><td>R > R<sub>int</sub><br>B < B<sub>high</sub><br>P < P<sub>low</sub></td><td>R < R<sub>int</sub><br>B > B<sub>high</sub><br>P > P<sub>high</sub></td><td>R > R<sub>int</sub><br>B > B<sub>high</sub><br>P > P<sub>high</sub></td></tr></tbody></table>

1. the cycle progresses through the eight (or six or four) stages in the expected order;

2. the duration of protein peaks, defined as the time spent by each variable above its corresponding threshold, are biologically-consistent (see Fig 1 and Table 1).

First, notice that looking for thresholds such that the cycle progresses through the stages in the expected order amounts to saying that crossings of the different thresholds are expected in a well-defined order. Thus, the main idea of the Algorithm is to detect the errors occurring in the order of crossing thresholds and adjust the values of the thresholds in consequence.

Secondly, the peaks of CLOCK:BMAL1, REV-ERB and PER:CRY are defined such that *B*, *R* and *P* are respectively above *B*<sub>high</sub>, *R*<sub>int</sub> and *P*<sub>high</sub>. Let's call τ<sub>X<sub>φ</sub></sub> the time in hours that each variable (*X*) remains above its corresponding threshold (*X*<sub>φ</sub>) over one period (*T*). In order to obtain duration of protein peaks that can be compared with biological duration, expressed in Circadian Time (see Fig 1), we scale the results to a period of 24h (τ<sub>X<sub>φ</sub></sub> is thus called τ<sub>sX<sub>φ</sub></sub>):

$$ \tau_{sX_{\phi}} = \frac{t_{end} - t_{int}}{T} \cdot 24, \text{ such that } X(t) > X_{\phi}, \forall t \in [t_{int}, t_{end}] \quad (1) $$

Let's set τ<sub>sB<sub>high</sub></sub> = 10h, τ<sub>sR<sub>int</sub></sub> = 8h, τ<sub>sP<sub>low</sub></sub> = 20h and τ<sub>sP<sub>high</sub></sub> = 12h as initial peaks duration. By abuse of notation, we use through this paper "h" to refer to a "circadian time unit".

### 3. Cycle segmentation as a benchmark for circadian clock models

To better appreciate how circadian clock models are close to experimental data, in terms of cycle segmentation timing and dynamics, we both applied our segmentation on a biological dataset and on five ODEs models.

Since the variables of our segmentation represent (complexes of) core clock proteins in a cell (mixing nuclear and cytoplasmic values), we chose a dataset measuring the total amount of protein within a cell. Narumi *et al.* Narumi, Shimizu, Ukai-Tadenuma, Ode, Kanda, Shinohara, Sato, Matsumoto and Ueda (2016) recently proposed a method to quantify absolute amounts of proteins, using mass spectrometry techniques. They measured the number of copies per cell (see Fig 3A from Narumi *et al.* (2016)) of several proteins from mouse liver, including those of interest to us: BMAL1, CLOCK, REV-ERBα, REV-ERBβ, PER1, PER2, PER3, CRY1 and CRY2. However, these data don't distinguish unbound proteins from complexes and measurements are carried out at low temporal resolution with time points only every 4 hours. Thus, to be as close as possible to the conditions of application of our segmentation, we associated each variable of the segmentation to the data available as shown in Table S1 and we interpolated the data with a spline, generating semi-synthetic time series for *B*, *R* and *P*.
**Algorithm 1** Numerical estimation of the threshold parameters

1: **Input:** time series data of *B*, *R* and *P* over one period, a partition, $τ_{sX_φ}$ for each $X_φ$ and a maximal number of iterations (nb_iter_max)

2: Initialize thresholds: compute each threshold value ($X_φ$) such that the time spent above each threshold is $τ_{sX_φ}$

3: Initialize threshold boundaries: compute lower (resp. upper) boundaries for each threshold such that the time spent above each boundary is $τ_{sX_φ} + 1$ (resp. $τ_{sX_φ} - 1$)

4: Check stages order: verify if the sequence of stages corresponds to the expected order

5: **while** expected stages order is not reached and nb_iter < nb_iter_max **do**

6: Identify $X_{φ_{crossed}}$ and $X_{φ_{expected}}$, the threshold crossed instead of the threshold expected (determine the first stage transition error)

7: **if** $X_{φ_{crossed}}$ or $X_{φ_{expected}}$ reached its boundary **then**

8: Enlarge the interval of acceptable threshold values: increase or decrease the boundary(ies) reached as needed

9: **end if**

10: Tune the values of $X_{φ_{crossed}}$ and $X_{φ_{expected}}$ by increasing or decreasing them

11: Check stages order: verify if the sequence of stages corresponds to the expected order

12: **end while**

13: **Output:** stages order and threshold values ($X_φ$)

Then, we focused on five models well established in the literature: two models from Almeida *et al.* (a model with eight variables proposed in 2019 Almeida *et al.* (2020b) and its four variables reduced version Almeida *et al.* (2020a)), the Relógio *et al.* model Relógio *et al.* (2011) proposed in 2011, the Hesse *et al.* model Hesse *et al.* (2021) proposed in 2021 (which refines the Relógio *et al.* model) and the Brown *et al.* model Brown and Doyle III (2020) proposed in 2020. These models, covering earlier to quite recent examples, have different dimensions, but are built in a similar way: they focus on the core components of the circadian clock and on the two transcription-translation feedback loops, between PER:CRY and CLOCK:BMAL1 and between REV-ERB and CLOCK:BMAL1. Some models make the distinction between nuclear and cytoplasmic proteins or between phosphorylated and unphosphorylated proteins and others do not. Thus, similarly to the experimental dataset, we associated each variable of the segmentation with the variables of each model, as given by Table S1.

We next used Algorithm 1 to partition all these cycles, ideally into the eight previously defined stages but, otherwise, into six stages, or, as a last resort, into four stages. As shown in Fig 2(a-g), for the experimental data based cycle and for all models except the Brown *et al.* model, the cycles could be successfully partitioned into the eight specified stages. Only a six-stages segmentation was effective for the Brown *et al.* model.

Next, we quantitatively characterize and compare all these circadian cycles. We specifically focus on the duration and timing of protein peaks and on the duration of stages.

**3.1. Duration and timing of protein peaks**

Algorithm 1 is expected to set thresholds such that the time spent by each variable above its corresponding threshold (see Eq 1) is biologically-coherent (see Fig 1). Recall that the initial (and expected) peaks duration established were 10h for *B*, 8h for *R* and 12h for *P*. For all the cycles (see Table 2) the peak durations are relatively close to the ones expected, except for the peak of *B* of the Brown *et al.* model, which is 4h shorter than expected.

A temporal representation of the segmentation obtained for each cycle is illustrated in Fig 2(b-g). Fig 2(i) also illustrates the timing of *B*, *R* and *P* peaks across Circadian Time, for the eight-stages segmentations. Timing of *B* and *P* peaks are similar throughout all cycles: respectively 6 hours (from CT2 to CT8.1) and 8 hours (from CT14.9 to CT22.9) of the Circadian Time are shared by all the cycles partitioned into eight stages. This representation highlights the robustness of the relative phase of *B* and *P*. Only 3.5 hours of the Circadian Time (from CT10.1 to CT13.6) is shared by all the cycles for *R*, showing the variability in the timing of *R* peak which agrees with experimental results Li, Mohammad-Djafari, Dumitru, Dulong, Filipski, Siffroi-Fernandez, Mteyrek, Scaglione, Guettier, Delaunay *et al.* (2013).
<figure><img src="image_5.png" /><figcaption>**Figure 2:** Illustration of circadian cycles segmentation. (a) Cycle segmentation with data (see points) from Narumi *et al.* (times series are the result of data interpolation with a spline). Temporal representation (scaled to a period of 24h) of maximums of *B*, *R* and *P* (respectively in red, green and blue) and segmentation into stages (delimited between two orange lines) for (b) Narumi *et al.* dataset, (c) Almeida *et al.* model (reduced), (d) Almeida *et al.* model (original), (e) Relógio *et al.* model, (f) Hesse *et al.* model, (g) Brown *et al.* model and (h) Brown *et al.* model with $P$:=PER+CRY+PER:CRY. The maximum of *B* is set at CT6.4 for models segmentation (since observed at CT8 in Takahashi (2017) and at CT4.7 in our reference experimental dataset Narumi *et al.* (2016)) and maximums of *R* and *P* and stages are represented relatively to the maximum of *B*. (i) Duration and timing of *B*, *R* and *P* peaks and (j) stages duration associated to each segmentation.</figcaption></figure>
<table><caption>**Table 2**

Peak durations scaled at a 24h-period for each dataset/model of Fig 2(b-h)</caption><thead><tr><th>Dataset / model</th><th>Fig 2</th><th>τ<sub>s</sub>B<sub>high</sub></th><th>τ<sub>s</sub>R<sub>int</sub></th><th>τ<sub>s</sub>P<sub>high</sub></th></tr></thead><tbody><tr><td>Narumi et al. (2016)</td><td>b)</td><td>9.22h</td><td>9.74h</td><td>12.32h</td></tr><tr><td>Almeida et al. (2020a) (reduced)</td><td>c)</td><td>9.23h</td><td>8.16h</td><td>10.18h</td></tr><tr><td>Almeida et al. (2020b) (original)</td><td>d)</td><td>9.52h</td><td>10.63h</td><td>10.93h</td></tr><tr><td>Relógio et al. (2011)</td><td>e)</td><td>9.06h</td><td>8.01h</td><td>11.37h</td></tr><tr><td>Hesse et al. (2021)</td><td>f)</td><td>8.76h</td><td>8.29h</td><td>12.24h</td></tr><tr><td>Brown and Doyle III (2020)</td><td>g)</td><td>6.05h</td><td>9.18h</td><td>9.55h</td></tr></tbody></table>

## 3.2. Duration of stages

Fig 2(j) summarizes the stages duration obtained from Fig 2(b-g).

### 3.2.1. Cycles with an eight-stages segmentation.

First note that, regardless of the cycle, stage 4 is very short and lasts just a few seconds. This stage corresponds to a weak expression level of all the variables and the cycle enters stage 4 when *P* decreases under its threshold *P*<sub>high</sub>, while *R* and *B* are low and leaves stage 4 when *B* increases above *B*<sub>high</sub>. From a biological point of view, as soon as PER and CRY are degraded, CLOCK:BMAL1 is derepressed and starts its transcriptional activity. A short stage 4 could represent a fast derepression of *B* as soon as the levels of *P* are low enough. Then, the intrinsic duration of stages 5, 6 and 7 is highly variable according to the cycles (from 0 to 7h) and the duration of stage 8 remains short (around 1h), even though the sum of these four stages, corresponding to the time where *B* is above *B*<sub>high</sub>, stays consistent across the cycles, with a duration of about 9h. During these stages, levels of *P* are low while levels of *B* are high. Finally, the sum of stages 2 and 3 (equivalent to τ<sub>s</sub>P<sub>high</sub>) lasts around 11h for all the cycles, but the intrinsic duration of stages 2 and 3 depends on the cycles. Considering experimental observations (see Fig 2(b)), the period during which *P* and *R* peak simultaneously, corresponding to stage 2, is expected to be short (close to 0), while the period during which *B* and *R* peak simultaneously, corresponding to stage 7, is expected to be long (around 6h).

### 3.2.2. Cycle segmentation of the Brown et al. model.

For this model, the duration of stage 1 is twice to three times longer than the one of the other models. The cycle enters stage 1 when *B* decreases below *B*<sub>high</sub> and exits stage 1 when *P* increases above *P*<sub>high</sub>. Thus, relatively to *B* repression, *P* increases above its threshold *P*<sub>high</sub> tardily. *B* and *P* are breaking phase opposition, due to a delay and a low activity of *P*. In the Brown *et al.* model, only the nuclear PER:CRY is taken into account while in the other models, both cytoplasmic and nuclear are considered. Remarkably, if we model the complex PER:CRY in the cytoplasm by summing up the variables CRY1, CRY2 and PER (proteins in the cytoplasm, see Table S1) and apply Algorithm 1, we obtain an eight-stages segmentation with consistent durations of peaks of proteins and durations of stages similar to those of other models (see Fig 2(h,j)). Our cycle segmentation analysis thus suggests that the addition of a new variable representing cytoplasmic PER:CRY to the model would lead to a partition into eight stages and more realistic peak durations.

## 3.3. Stages related to Circadian Time phases in the literature

Using extensive DNA binding profiling in the mouse liver, Takahashi *et al.* (2017) proposed six distinctive phases characterizing the circadian cycle, as summarized in Fig 3.

According to their description, the circadian cycle begins with a poised state, at CT0, during which "CLOCK:BMAL1 and CRY1 bind to E-box". A derepression phase follows, at CT4, "as the levels of CRY1 decline". Thus, between CT0 and CT4, *B* is gradually derepressed as levels of *P* are decreasing, corresponding in our segmentation to stage 4 and the beginning of stage 5.

The third phase consists in "a major activation phase at CT8 during which CLOCK:BMAL1 occupancy peak", and is followed by a "transcriptional activation phase at CT12-15 during which pre-mRNA transcript expression peak". CLOCK and BMAL1 binding activities are observed throughout subjective daytime hours, until CT12, while REV-ERB peak is observed in the late daytime. In our segmentation, stages 5 to 8 characterize high levels of *B*. Thus, stages 5 to 8 spread out subjective daytime hours, from CT2 to CT12. Then, since (i) we observed the maximum of *R* either in
<figure><img src="image_6.png" /><figcaption>**Figure 3: Temporal analogy between the phases described by Takahashi in Takahashi (2017) and our cycle segmentation stages.** The dashed line approximates the temporal expression of each stage.</figcaption></figure>

stage 7, 1 or 2, (ii) the transition between stage 8 and stage 1 represents *B* repression, and (iii) levels of *P* are increasing during stage 1, we assumed that stage 1 marks the transition between daytime and night hours, around CT12.

Then, they identified "a repression phase at CT15–CT18 during which PER1, PER2 and CRY2 occupancy peaks and transcription declines". In other words, levels of *P* are high and *P* strongly represses *B*, as characterized by stages 2 and 3 of our segmentation.

Finally, the cycle is in "a transition between the end of repression and the beginning of the next cycle", marking the end of subjective night hours, that we associated with the end of our stage 3.

Roughly, confirming the results obtained in Burckard and Chaves (2024), we thus could associate subjective night hours to stages 2 and 3 and daytime hours to stages 4, 5, 6, 7, 8, while stage 1 is at the day-to-night transition.

## 4. Knockdown impacts on oscillators dynamic

As another application, we used our segmentation into stages to investigate the changes in the dynamics and the reorganization of the circadian cycle upon the knockdown of a core clock component. The experimental deletion or silencing of *Cry2* reproducibly leads to a long period phenotype under free-running conditions Horst et al. (1999); Zhang, Liu, Hirota, Miraglia, Welch, Pongsawakul, Liu, Atwood, Huss, Janes et al. (2009); Ramanathan, Xu, Khan, Shen, Gitis, Welsh, Hogenesch and Liu (2014). Building on these experimental results, we mimicked the *Cry2* knockdown in the previously studied models by multiplying the transcription rate of *Cry* (or, if not applicable, the transcription rate of PER:CRY Burckard et al. (2022)) by a factor $\epsilon$ lower than 1. This leads to a decrease of the transcription rate of *Cry* and to an increase of the period as $\epsilon$ decreases below 1.

Using the segmentation into stages, we aimed at studying and comparing the cycle dynamic between a reference oscillator (i.e. $\epsilon = 1$) and long-period oscillators (i.e. $\epsilon < 1$). Considering the set of constraints required and the quantity of information returned by the eight-stages segmentation, as well as the variability in the timing expression of REV:ERB, it appeared that the eight-stages segmentation was not the most suitable cycle partition. For an easier comparison and for better highlighting some observations, we chose the segmentation into four stages. In this case, only two thresholds are set: $B_{high}$ on *B* and $P_{high}$ on *P*. The resulting stages order is, as explained in the Section 2.1, 1, 2, 4, 5. Stage 2 (resp. 5) represents *P* (resp. *B*) above $P_{high}$ (resp. $B_{high}$) and *B* (resp. *P*) below $B_{high}$ (resp. $P_{high}$). Stages 1 and 4 are identical, with *B* and *P* lower than their respective thresholds. We nevertheless distinguish them to differentiate the time before the peak of *B* and the peak of *P*. Finally, stage 11 may also be part of the four-stages segmentation and corresponds to *B* and *P* above their respective thresholds simultaneously.

To compare and characterize each cycle (varying according to $\epsilon$) within and across models, we adopted the following strategy:

1. Normalize the amplitudes of all the cycles between 0 and 1

2. Use Algorithm 1 to compute the values of $B_{high}$ and $P_{high}$ for a four-stages segmentation applied to the reference oscillator (i.e. $\epsilon = 1$)

3. Fix the thresholds for each oscillator (with the values determined during the previous step)

4. Partition the cycles with the fixed threshold values
## 4.1. Stages order and duration

In the reduced model of Almeida et al., the initial order of stages (1, 2, 4, 5, 1) persists regardless the value of $\epsilon$. For the Relógio et al., Hesse et al. and the original Almeida et al. models, stage 4 disappears and is replaced by stage 11 giving rise to the order 1, 2, 11, 5, 1. Fig 4 shows the stages duration, as a function of $\epsilon$.

<figure><img src="image_7.png" /><figcaption>**Figure 4:** Stages duration as a function of $\epsilon$. (a) Almeida et al. (reduced), (b) Almeida et al. (original), (c) Relógio et al. and (d) Hesse et al. models. The range of $\epsilon$ in each model corresponds to a change of 1.3 hours in the total period (except in case (b), where only a small period variation is observed with $\epsilon$).</figcaption></figure>

A common observation is the increase of stage 1 as $\epsilon$ decreases due to a delay in the time at which P crosses $P_{high}$, relatively to the time where B crosses $B_{high}$. Then, two different behaviors arise according to the models: either stages 2 and 4 increase as $\epsilon$ decreases (Fig 4(a)), or stage 2 decreases and stage 11 increases (Fig 4(b-d)). In the first case, an increase of stage 4 represents a delay of the increase of B above $B_{high}$ relatively to the decrease of P below $P_{high}$. This delay results in a shorter peak of B as represented by the decrease of stage 5. In the second case, stage 2 decreases with the decrease of $\epsilon$ and is followed by stage 11. Stage 11 is defined by both B and P above their thresholds. Thus, stage 2 followed by stage 11 highlights that B increases above $B_{high}$ while P is still above $P_{high}$: the phase opposition between B and P is misaligned.

## 4.2. Phase Difference Curves

To go deeper into this analysis, we introduce a tool to compare the dynamics of two oscillators throughout their cycle. The phase response curve (PRC) is commonly used to capture the effect of a perturbation on the phase of a periodic oscillator. Considering a reference position/phase in a cycle (or an observable event in experiments such as the onset of locomotor activity of mice), the phase response curve shows the phase shift of the perturbed oscillator according to the phase at which the stimulus (light pulse for example) is applied (Sacre and Sepulchre (2014). Here, we introduce the phase difference curve to compare the dynamic between two oscillators, of which one is affected by a knockdown. Thus, instead of analyzing the phase shift of one reference position following a stimulus at different times, we propose to compute the phase shift between a reference oscillator and an oscillator with a knockdown at each amplitude point of the cycle. For this, let oscillator 1 be the reference oscillator, and oscillator 2 be the one impacted by a knockdown. The amplitudes of each oscillator are normalized between 0 and 1. Let the time of the minimum of R be the initial starting point. Then, starting from this initial time, the phase difference curve corresponds to the phase shift between oscillator 2 and oscillator 1 for each amplitude point (each phase) of each variable throughout a cycle. In practice, this amounts to computing the time that variable X of oscillator 2 needs to reach the same amplitude as oscillator 1 (see Fig S2), as given by Eq 2:
$$ \Delta t = t_2 - t_1, \text{ such that } X_1(t_1) = X_2(t_2) \text{ and} \qquad (2) $$

$$ \operatorname{sign}(X_1(t_1 - dt_1) - X_1(t_1)) = \operatorname{sign}(X_2(t_2 - dt_2) - X_2(t_2)) \forall t_1 \in [0, T_1] \qquad (3) $$

where $X_i$ refers to variable $X$ of oscillator $i$, and $T_1$ to the period of oscillator 1.

We thus computed the phase difference curves of a long-period oscillator compared to the reference oscillator (see Fig 5). For each model, $\epsilon$ is tuned such that the period of oscillator 2 ($\epsilon < 1$) increases by approximately 0.7h (0.32h for the original model of Almeida et al.) compared to oscillator 1 ($\epsilon = 1$). This period augmentation is in the observed range of period differences between WT mice and Cry2 deficient mice (e.g., 0.86h in Horst et al. (1999)).

The PDCs are represented over one cycle, whose duration corresponds to the period of oscillator 1, and the difference between $\Delta t$ computed at the initial time and $\Delta t$ computed at the end of the cycle is equal to the period difference between the two oscillators.

<figure><img src="image_8.png" /><figcaption>Figure 5: Phase difference curves between a long-period oscillator and the reference oscillator for B (dotted curves), R (dashed curves) and P (tri-star curves). (a) Almeida et al. (reduced) (b) Almeida et al. (original) (c) Relógio et al. (d) Hesse et al. models. The marker color indicates the stage of the long-period oscillator at the time where a given variable of this oscillator reaches the same amplitude point as the same variable of the reference oscillator. Dashed vertical lines correspond to the time at which oscillator 2 transitions from one stage to another according to $P_2$.</figcaption></figure>

If the period augmentation is accompanied by a proportional augmentation across all stages of the cycle, one may expect a linear PDC. However, this is not the case for any of the oscillators. In fact, even within each oscillator the impact of $\epsilon$ on the dynamic of each variable differs, with distinct accelerations (phase advance) or slowdowns (delays) occurring at different cycle times according to the variables. For instance, a decrease in the PDC of variable $X$ corresponds to an acceleration of oscillator 2 relative to the reference oscillator, that is, $X_2(t_2)$ reaches $X_1(t_1)$ “faster”, leading to a smaller $\Delta t = t_2 - t_1$. Similarly, an increase in the PDC corresponds to a delay in oscillator 2, with larger $\Delta t$. In general, the PDC for variable $R$ tends to evolve monotonically and is closer to a linear form, for all models. Conversely, the PDC curves for $B$ and $P$ go through regions of acceleration and slowing-down throughout the cycle.

Fig 5 also shows the stage of oscillator 2 at instant $t_2$ of the given variable (recall $X_2(t_2) = X_1(t_1)$). The dashed vertical lines mark the transitions of oscillator 2 from one stage to another according to the amplitude points reached by variable $P_2$. Note that, along the horizontal axis (time), the transitions between two stages do not happen simultaneously among the three variables, and especially for $P$ compared to $B$ (see also Fig S2). Misaligned transitions highlight new
behaviors in the dynamic of the three variables for oscillator 2 compared to oscillator 1: consider the transition between stage 5 and 1 (transition from pink to green) of each model. This transition occurs earlier for the PDC of *P* (tri-star marker) than for *B* (dot marker). It corresponds to a delay in the phase of *P*₂ relatively to *B*₂ as compared to *P*₁ with *B*₁: at the time at which *B*₂ decreases below *B*<sub>high</sub> (transition from stage 5 to 1), *P*₂ is phase delayed compared to the phase at which *P*₁ was when *B*₁ decreased below *B*<sub>high</sub>. Moreover, during stage 1, the PDCs of *P* keep increasing: oscillator 2 slows down, which increases the delay of *P*₂ compared to *P*₁. Thus, due to the delay of *P*₂ at the beginning of stage 1 and its subsequent slowing down, the duration of stage 1 increases for oscillator 2, corroborating Fig 4. This illustrates how a longer period is obtained by stretching stage 1 (which corresponds to the subjective day-to-night transition) through a slow down of *P* and a misalignment between the temporal evolution of *B* and *P*.

**5. Discussion**

Given the extensive biological knowledge on the circadian clock and the many mathematical models derived from it, our study aims to provide an explicit and general method for characterizing the dynamics of the circadian clock cycle.

Based on an analytic circadian clock cycle segmentation into stages, we first proposed an Algorithm to automati-
cally partition the cycle of any circadian clock model built with ODEs and including the three core clock components:
the protein REV-ERB and complexes CLOCK:BMAL1 and PER:CRY. The goal of this Algorithm is to find suitable
values of thresholds such that a given cycle progresses through the stages of the partition in an expected order. Starting
with initial threshold values, the Algorithm detects the errors occurring in the order of crossing thresholds (and thus
in the order of stages), and automatically corrects the errors by adjusting the values of the thresholds.

We applied this Algorithm to partition the cycles of an experimental dataset and of five different models. We analyzed the resulting segmentations, characterized the circadian cycles and checked their progression and dynamics along the cycle. One advantage of our method is to detect some problems that may appear and indicate ways for improvement. For example, a delay in the peak of PER:CRY in the Brown *et al*. model Brown and Doyle III (2020) seems due to a late and low activity of this complex relatively to the other proteins. To solve this problem, application of our Algorithm suggested the addition of a variable representing PER:CRY in the cytoplasm. Moreover, our method facilitated the analysis of circadian cycles and highlighted some features shared by all models, showing for example robustness of CLOCK:BMAL1 and PER:CRY phases. Indeed, it appeared that intrinsic duration of stages 5, 6 and 7 (and of stages 2 and 3) may differ according to the models, but the timing and duration of their sum (corresponding respectively to *B* and *P* peaks) remained similar through the different cycles studied. Moreover, all models exhibited a short duration of stage 4, which could be associated to a fast derepression of CLOCK:BMAL1 as soon as levels of PER and CRY are low enough. The properties shared by all models in term of timing of protein peaks or stages duration allowed us to establish a closer correspondence between our analytic stages and experimental circadian time phases.

Our segmentation analysis therefore provides a benchmark for quantitative comparison of different models of the mammalian clock: independently of their absolute amplitudes, units and total period, the time dynamics and progress of the cycle can be evaluated by the relative durations of protein peaks and stages. By observing both the progress of the circadian cycle through the stages and the relative duration of stages and protein peaks associated, this benchmark model will thus permit (in)validation of a clock model and provide insights for potential enhancements.

Finally, we showed how the segmentation into stages can be applied to investigate the new dynamics of an oscillator (throughout its entire cycle) following a perturbation such as a gene knockdown. As a mechanism for the long period phenotype caused by *Cry2* silencing, our analysis suggested a misalignment between the temporal evolution of BMAL1 and PER:CRY, through a delay in the augmentation of PER:CRY above its high threshold.

While our approach used the well documented mouse circadian clock as a benchmark, it could be easily imple-
mented with the clock of many other model systems by proposing a segmentation into stages adapted to the circadian
cycle of the organism of interest. Indeed, clock mechanisms and components are highly conserved across metazoans
from Drosophila Leloup and Goldbeter (1998) to humans. Further, even when components differ between Phyla for
instance between fungi, plants and animal cells, very similar design principle including transcriptional/translational
negative feedback loops are shared despite the evolutionary divergence Clock (1996); Baker, Loros and Dunlap (2012);
François (2005); Paranjpe and Kumar Sharma (2005). Our method provides a tool for quantifying the characteristics
of the circadian clock cycle, allowing for instance an estimation of the time duration of the observed biological phases
Takahashi (2017). Moreover, it reveals those phases that are time-restricted and need to be finely adjusted (cf fast
derepression, stage 4) or, conversely, those phases that can sustain higher variability. Finally, with its potential for
application to other circadian clocks, this tool provides a means for studying and quantitatively comparing the cycle dynamics from various clocks, and to investigate whether and how quantitative temporal features are shared by different organisms.

## Funding

This work was supported in part by project InSync (ANR-22-CE45-0012) and by LABEX SIGNALIFE (ANR-11-LABX-0028) and IDEX UCA Jedi (ANR-15-IDEX-01) from the program Investments for the Future of the French National Agency for Research.

## Data availability

The codes underlying the results presented in the study are available from https://gitlab.inria.fr/oburckar/cycle_partition.git.

## A. Appendix

### Description of each step of Algorithm 1.

**Steps 1 to 3.** Algorithm 1 takes as input time series data of B, R and P over one period -from the output of a model or from experimental measures-, a target partition (see above) with an associated specified order through the stages and initial peaks duration ($\tau_{sX_{\phi}}$) and generates an initial estimation of the thresholds values (as well as a tolerance range for threshold values), such that the value spent by each variable X above its threshold $X_{\phi}$ corresponds to $\tau_{sX_{\phi}}$ (to $\tau_{s_{max}X_{\phi}} = \tau_{sX_{\phi}} + 1$ and to $\tau_{s_{min}X_{\phi}} = \tau_{sX_{\phi}} - 1$ for the respective lower and upper threshold boundaries).

**Step 4.** Algorithm 1 verifies if the sequence of stages corresponds to the expected order.

**Steps 5 to 12.** While both the expected stages order and a maximal number of iteration are not reached, Algorithm 1 detects and corrects the errors occurring in the transitions between stages, by notably adjusting the thresholds $X_{\phi}$.

**Steps 6.** To begin, Algorithm 1 determines the first stage transition error and identifies which threshold is crossed instead of the desired one.

In what follows, the idea is to adjust the values of each of the two thresholds identified as responsible for the error by increasing or decreasing them such that the oscillator crosses the expected threshold before the threshold crossed in error.

**Steps 7 to 9.** If a threshold has already reached its upper (resp. lower) boundary and that it must be increased (resp. decreased), the value of the upper (resp. lower) boundary is widened such that $\tau_{s_{min}X_{\phi}} = \tau_{s_{min}X_{\phi}} - 1$ (resp. $\tau_{s_{max}X_{\phi}} = \tau_{s_{max}X_{\phi}} + 1$).

**Step 10.** Algorithm 1 tunes the values of the threshold currently crossed (called $X_{\phi_{crossed}}$) and of the threshold expected ($X_{\phi_{expected}}$). The oscillator crosses $X_{\phi_{crossed}}$ at time $t_{X_{\phi_{crossed}}}$ and $X_{\phi_{expected}}$ at time $t_{X_{\phi_{expected}}}$, with $t_{X_{\phi_{crossed}}} < t_{X_{\phi_{expected}}}$. To ensure that $X_{\phi_{expected}}$ is crossed before $X_{\phi_{crossed}}$, the idea is to reverse the order in which the thresholds are crossed (see also Fig S1):

$$X_{\phi_{crossed}} := X\left(\frac{t_{X_{\phi_{crossed}}} + t_{X_{\phi_{expected}}}}{2} + 2dt\right) \quad (4)$$
<table><caption>**Table S1**
Association between the variables of the segmentation and the ones of each model or dataset</caption><thead><tr><th>Variable</th><th>B :=</th><th>R :=</th><th>P :=</th></tr></thead><tbody><tr><td>Represents</td><td>Complex CLOCK:BMAL1</td><td>REV-ERB and paralogs</td><td>Proteins and complexes related to PER and CRY and blocking the transcriptional activity of CLOCK:BMAL1</td></tr><tr><td>Narumi et al. (2016)</td><td>min(CLOCK, BMAL1)</td><td>REV-ERB<sub>α</sub> + REV-ERB<sub>β</sub></td><td>min(PER1+PER2+PER3, CRY2) + CRY1</td></tr><tr><td>Almeida et al. (2020a) (reduced)</td><td>CLOCK:BMAL1</td><td>REV-ERB</td><td>PER:CRY</td></tr><tr><td>Almeida et al. (2020b) (original)</td><td>CLOCK:BMAL1</td><td>REV-ERB</td><td>PER:CRY</td></tr><tr><td>Relógio et al. (2011)</td><td>CLOCK:BMAL</td><td>REV-ERB<sub>C</sub> + REV-ERB<sub>N</sub></td><td>PER:CRY<sub>C</sub> + PER*:CRY<sub>C</sub> + PER:CRY<sub>N</sub> + PER*:CRY<sub>N</sub></td></tr><tr><td>Hesse et al. (2021)</td><td>CLOCK:BMAL<sub>C</sub> + CLOCK:BMAL<sub>N</sub></td><td>REV-ERB<sub>C</sub> + REV-ERB<sub>N</sub></td><td>PER:CRY<sub>C</sub> + PER:CRY<sub>N</sub></td></tr><tr><td>Brown and Doyle III (2020)</td><td>BMAL1</td><td>REV-ERB<sub>α</sub></td><td>PER:CRY<sub>1</sub><sub>N</sub> + PER:CRY<sub>2</sub><sub>N</sub></td></tr><tr><td>Brown and Doyle III (2020)</td><td>BMAL1</td><td>REV-ERB<sub>α</sub></td><td>CRY1 + CRY2 + PER + PER:CRY<sub>1</sub><sub>N</sub> + PER:CRY<sub>2</sub><sub>N</sub> (Version 2)</td></tr></tbody></table>

$$X_{\phi_{\text{expected}}} := X \left( \frac{t X_{\phi_{\text{crossed}}} + t X_{\phi_{\text{expected}}}}{2} - 2dt \right), \qquad (5)$$

with $dt = t_{i+1} - t_i$.

If the updated values of the thresholds $X_{\phi_{crossed}}$ or $X_{\phi_{expected}}$ exceed the range of their boundary, the threshold(s) value(s) concerned is (are) brought back to their respective limits.

**Step 11.** Algorithm 1 verifies if the sequence of stages corresponds to the expected order.

**Step 13.** Algorithm 1 returns as output the first set of threshold values encountered allowing for the expected partition or the last stages order obtained if the Algorithm could not converge to a suitable set of thresholds within the maximal number of iterations.

## References

Abe, Y.O., Yoshitane, H., Kim, D.W., Kawakami, S., Koebis, M., Nakao, K., Aiba, A., Kim, J.K., Fukada, Y., 2022. Rhythmic transcription of bmal1 stabilizes the circadian timekeeping system in mammals. Nature Communications 13, 4652.

Almeida, S., Chaves, M., Delaunay, F., 2020a. Control of synchronization ratios in clock/cell cycle coupling by growth factors and glucocorticoids. Royal Society Open Science 7, 192054.

Almeida, S., Chaves, M., Delaunay, F., 2020b. Transcription-based circadian mechanism controls the duration of molecular clock states in response to signaling inputs. Journal of theoretical biology 484, 110015.

Atwood, A., DeConde, R., Wang, S.S., Mockler, T.C., Sabir, J.S., Ideker, T., Kay, S.A., 2011. Cell-autonomous circadian clock of hepatocytes drives rhythms in transcription and polyamine synthesis. Proceedings of the National Academy of Sciences 108, 18560–18565.

Baker, C.L., Loros, J.J., Dunlap, J.C., 2012. The circadian clock of neurospora crassa. FEMS microbiology reviews 36, 95–110.

Balsalobre, A., Damiola, F., Schibler, U., 1998. A serum shock induces circadian gene expression in mammalian tissue culture cells. Cell 93, 929–937.

Banerjee, S., Wang, Y., Solt, L.A., Griffett, K., Kazantzis, M., Amador, A., El-Gendy, B.M., Huitron-Resendiz, S., Roberts, A.J., Shin, Y., et al., 2014. Pharmacological targeting of the mammalian clock regulates sleep architecture and emotional behaviour. Nature communications 5, 5759.

Brown, L.S., Doyle III, F.J., 2020. A dual-feedback loop model of the mammalian circadian clock for multi-input control of circadian phase. PLoS Computational Biology 16, e1008459.

Bugge, A., Feng, D., Everett, L.J., Briggs, E.R., Mullican, S.E., Wang, F., Jager, J., Lazar, M.A., 2012. Rev-erbα and rev-erbβ coordinately protect the circadian clock and normal metabolic function. Genes & development 26, 657–667.

Burckard, O., Chaves, M., 2024. Analytic solutions for the circadian oscillator characterize cycle dynamics and its robustness. URL: https://hal.science/hal-04391710.
<figure><img src="image_9.png" /><figcaption>**Figure S1:** Illustration of Step 10 of Algorithm 1. Each row corresponds to time series values of respectively *B*, *R* and *P* (in bold when the variable is above its threshold), with *dt* = 0.5. Dashed vertical lines represent the times of thresholds crossing and thus mark the transition between two stages. Down (resp. up) arrow means that a given threshold is crossed while the corresponding variable is decreasing (resp. increasing). Figures (d), (e) and (f) illustrate respectively tables (a), (b) and (c). (a) Stage 6 is followed by stage 9 instead of stage 7, because *B*<sub>high</sub> is crossed instead of *R*<sub>int</sub>: the values of these two thresholds are updated according to Eq 2 and 3, and shown in (b). (b) Stage 7 is followed by stage 10 instead of stage 8, because *B*<sub>high</sub> is crossed instead of *P*<sub>low</sub>: the values of these two thresholds are updated according to Eq 2 and 3, and shown in (c). (c) The expected order of stages on this time period is reached.</figcaption></figure>

Burckard, O., Teboul, M., Delaunay, F., Chaves, M., 2022. Cycle dynamics and synchronization in a coupled network of peripheral circadian clocks. Interface Focus 12, 20210087.

Cao, X., Yang, Y., Selby, C.P., Liu, Z., Sancar, A., 2021. Molecular mechanism of the repressive phase of the mammalian circadian clock. Proceedings of the National Academy of Sciences 118, e2021174118.

Chiou, Y.Y., Yang, Y., Rashid, N., Ye, R., Selby, C.P., Sancar, A., 2016. Mammalian period represses and de-represses transcription by displacing clock–bmal1 from promoters in a cryptochrome-dependent manner. Proceedings of the National Academy of Sciences 113, E6072–E6079.

Cho, H., Zhao, X., Hatori, M., Yu, R.T., Barish, G.D., Lam, M.T., Chong, L.W., DiTacchio, L., Atkins, A.R., Glass, C.K., et al., 2012. Regulation of circadian behaviour and metabolism by rev-erb-α and rev-erb-β. Nature 485, 123–127.

Clock, T.P.C., 1996. Ordered phosphorylation governs oscillation of. Proc. Natl. Acad. Sci. USA 93, 6947.

Feng, D., Liu, T., Sun, Z., Bugge, A., Mullican, S.E., Alenghat, T., Liu, X.S., Lazar, M.A., 2011. A circadian rhythm orchestrated by histone deacetylase 3 controls hepatic lipid metabolism. Science 331, 1315–1319.

François, P., 2005. A model for the neurospora circadian clock. Biophysical journal 88, 2369–2383.

Gabriel, C.H., Del Olmo, M., Zehtabian, A., Jäger, M., Reischl, S., van Dijk, H., Ulbricht, C., Rakhymzhan, A., Korte, T., Koller, B., et al., 2021. Live-cell imaging of circadian clock protein dynamics in crispr-generated knock-in cells. Nature communications 12, 3796.
<figure><img src="image_10.png" /><figcaption>Figure S2: Illustration of a phase difference computation at one time point between oscillator 1 (solid line) and oscillator 2 (dashed line) for the Hesse *et al.* model.</figcaption></figure>

Gatfield, D., Le Martelot, G., Vejnar, C.E., Gerlach, D., Schaad, O., Fleury-Olela, F., Ruskeepää, A.L., Oresic, M., Esau, C.C., Zdobnov, E.M., et al., 2009. Integration of microrna mir-122 in hepatic circadian gene expression. *Genes & development* 23, 1313–1326.

Gouzé, J.L., 1998. Positive and negative circuits in dynamical systems. *Journal of Biological Systems* 6, 11–15.

Gréchez-Cassiau, A., Rayet, B., Guillaumond, F., Teboul, M., Delaunay, F., 2008. The circadian clock component bmal1 is a critical regulator of p21waf1/cip1 expression and hepatocyte proliferation. *Journal of Biological Chemistry* 283, 4535–4542.

Hesse, J., Martinelli, J., Aboumanify, O., Ballesta, A., Relógio, A., 2021. A mathematical model of the circadian clock and drug pharmacology to optimize irinotecan administration timing in colorectal cancer. *Computational and Structural Biotechnology Journal* 19, 5170–5183.

Horst, G.T.v.d., Muijtjens, M., Kobayashi, K., Takano, R., Kanno, S.i., Takao, M., Wit, J.d., Verkerk, A., Eker, A.P., Leenen, D.v., et al., 1999. Mammalian *cry1* and *cry2* are essential for maintenance of circadian rhythms. *Nature* 398, 627–630.

Jeong, E.M., Song, Y.M., Kim, J.K., 2022. Combined multiple transcriptional repression mechanisms generate ultrasensitivity and oscillations. *Interface Focus* 12, 20210084.

Jud, C., Schmutz, I., Hampp, G., Oster, H., Albrecht, U., 2005. A guideline for analyzing circadian wheel-running behavior in rodents under different lighting conditions. *Biological procedures online* 7, 101–116.

Ko, C.H., Takahashi, J.S., 2006. Molecular components of the mammalian circadian clock. *Human molecular genetics* 15, R271–R277.

Koike, N., Yoo, S.H., Huang, H.C., Kumar, V., Lee, C., Kim, T.K., Takahashi, J.S., 2012. Transcriptional architecture and chromatin landscape of the core circadian clock in mammals. *science* 338, 349–354.

Leloup, J.C., Goldbeter, A., 1998. A model for circadian rhythms in drosophila incorporating the formation of a complex between the per and tim proteins. *Journal of biological rhythms* 13, 70–87.

Li, X.M., Mohammad-Djafari, A., Dumitru, M., Dulong, S., Filipski, E., Sifroi-Fernandez, S., Mteyrek, A., Scaglione, F., Guettier, C., Delaunay, F., et al., 2013. A circadian clock transcription model for the personalization of cancer chronotherapy. *Cancer research* 73, 7176–7188.

Minami, Y., Ode, K.L., Ueda, H.R., 2013. Mammalian circadian clock: the roles of transcriptional repression and delay. *Circadian clocks*, 359–377.

Narumi, R., Shimizu, Y., Ukai-Tadenuma, M., Ode, K.L., Kanda, G.N., Shinohara, Y., Sato, A., Matsumoto, K., Ueda, H.R., 2016. Mass spectrometry-based absolute quantification reveals rhythmic variation of mouse circadian clock proteins. *Proceedings of the National Academy of Sciences* 113, E3461–E3467.

Novák, B., Tyson, J.J., 2008. Design principles of biochemical oscillators. *Nature reviews Molecular cell biology* 9, 981–991.

Paranjpe, D.A., Kumar Sharma, V., 2005. Evolution of temporal order in living organisms. *Journal of Circadian Rhythms* 3, 1–13.

Ramanathan, C., Xu, H., Khan, S.K., Shen, Y., Gitis, P.J., Welsh, D.K., Hogenesch, J.B., Liu, A.C., 2014. Cell type-specific functions of period genes revealed by novel adipocyte and hepatocyte circadian clock models. *PLoS genetics* 10, e1004244.

Relógio, A., Westermark, P.O., Wallach, T., Schellenberg, K., Kramer, A., Herzel, H., 2011. Tuning the mammalian circadian clock: robust synergy of two loops. *PLoS computational biology* 7, e1002309.

Richard, A., Comet, J.P., 2011. Stable periodicity and negative circuits in differential systems. *Journal of mathematical biology* 63, 593–600.

Ripperger, J.A., Jud, C., Albrecht, U., 2011. The daily rhythm of mice. *FEBS letters* 585, 1384–1392.

Ripperger, J.A., Schibler, U., 2006. Rhythmic clock-bmal1 binding to multiple e-box motifs drives circadian dbp transcription and chromatin transitions. *Nature genetics* 38, 369–374.