# Analytic solutions for the circadian oscillator characterize cycle dynamics and its robustness

Odile Burckard, Madalena Chaves

► To cite this version:

Odile Burckard, Madalena Chaves. Analytic solutions for the circadian oscillator characterize cycle dynamics and its robustness. Journal of Mathematical Biology, 2025, 90 (1), pp.5. 10.1007/s00285-024-02164-y. hal-04391710v2

HAL Id: hal-04391710

https://hal.science/hal-04391710v2

Submitted on 18 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Analytic solutions for the circadian oscillator
characterize cycle dynamics and its robustness

Odile Burckard¹* and Madalena Chaves¹

¹*Centre Inria d'Université Côte d'Azur, INRAE, CNRS, Macbes team,
Sophia Antipolis, France.

*Corresponding author(s). E-mail(s): odile.burckard@inria.fr;

Contributing authors: madalena.chaves@inria.fr;

**Abstract**

Circadian clocks form a fundamental mechanism that promotes the correct behav-
ior of many cellular and molecular processes by synchronizing them on a 24 hour
period. However, the circadian cycles remain difficult to describe mathematically.
To overcome this problem, we first propose a segmentation of the circadian cycle
into eight stages based on the levels of expression of the core clock components
CLOCK:BMAL1, REV-ERB and PER:CRY. This cycle segmentation is next
characterized through a piecewise affine model, whose analytical study allows us
to propose an Algorithm to generate biologically-consistent circadian oscillators.
Our study provides a characterization of the cycle dynamics in terms of four fun-
damental threshold parameters and one scaling parameter, shows robustness of
the circadian system and its period, and identifies critical points for correct cycle
progression.

**Keywords:** Circadian clock cycle dynamics, piecewise affine model, analytic solutions, parameter regions, period robustness

# 1 Introduction

In mammals, every cell contains its own circadian clock network, a complex fam-
ily of interactions between genes and proteins that helps to control and synchronize
many cellular and molecular processes (such as heart beat, blood pressure, body tem-
perature,...). The main synchronizer of circadian clock is known to be the Earth's
light/dark cycle, forcing circadian clock period at around 24h. The core mechanism of
circadian clock is based on two negative feedback loops, between the complexes of pro-
teins PER:CRY and CLOCK:BMAL1 and between CLOCK:BMAL1 and REV-ERB,
allowing circadian clocks to show sustainable and rhythmic oscillations [1, 2].

The regulatory interactions of circadian clocks are frequently described with ordinary differential equations (ODEs). Throughout the years, several mathematical models have been developed using this framework, especially for the clock located in the Suprachiasmatic Nucleus (SCN) [3, 4] but also for peripheral organs such as the liver [5–7], or pancreatic beta-cells [8] or for more ‘generic’-type cells [2, 9], focusing on different circadian clock features (such as transcription, translation, import/export, degradation, phosphorylation...) and recovering important biological properties of circadian clocks. These models enable, among other points, to provide insights to show circadian period and oscillation robustness whether by studying it with a reduced number of clock molecules (mRNA or proteins) in the cell [10] or demonstrating the role of the two loops [4].

Modeling complex dynamic of genetic regulatory networks with ODEs often neces-
sitates the use of nonlinear terms, such as Hill functions or mass action terms, that
complicate the study of these models since analytical solutions cannot be deter-
mined. An alternative framework considers the switching behaviour of the gene-protein
interactions, by which Glass [11–14] first proposed a simplified representation based
on constant rates for species production and linear degradation rates. The net-
work dynamic is then governed by threshold-dependant regulations between variables:
depending on the position of a regulatory variable relative to a given threshold, the
constant production rate of a regulated variable can be switched either high or low.

This representation was generalized to a framework also known as piecewise affine
(PWA) systems that has been widely used to analyse genetic networks in computa-
tional and theoretical [15–18] ways. It has the advantage that the set of thresholds
allows to partition the state space into distinct domains or “boxes” in which solutions
can be easily calculated piecewise and then the global solution is a concatenation of
the different pieces [17]. More specifically, the global behavior of PWA models has been
widely studied in order to prove the existence and stability of periodic solutions, often
requiring restriction to specific systems and under some assumptions such as identical
degradation rates for each variable and a single threshold per variable [16, 19]. Subse-
quent research extended to systems with several thresholds per variable (i.e. a given
variable regulating several variables does so at different thresholds) [20, 21], or with
different degradation rates [18].

A new wave of generalization of these networks was more recently developed, with among others, the introduction of a more general concept of switching networks [22] associated with a method to characterize all the possible qualitative behaviors “contained” in a given network, by combinatorially identifying regions of parameters for which the qualitative dynamics don’t change. This method is implemented through the software for Dynamical Signatures Generated by Regulatory Networks (DSGRN) [22, 23]. Applications of this new DSGRN framework include model validation or rejection [24], with the particular example of the yeast cell cycle [25].
Continuous models that use only sigmoidal-type production terms and linear degradations can usually be associated to a corresponding PWA system by approximating the sigmoidal functions by step functions. However, in the case of models more generally described by functions of mass-action form or polynomial terms, these frameworks cannot in general be directly applied. In this work, we propose the construction of a piecewise affine model that approximates a reduced continuous nonlinear ODE model of circadian clock [26] including both a Hill function and mass action terms. Due to the polynomial terms, the previously cited methods cannot be applied, and the simplification leads to a system which is not of the classical PWA form, and which directly relies on the knowledge of the qualitative dynamics of the ODE model and of the general dynamics of circadian clocks.

More specifically, the qualitative behaviour refers to two features of circadian clock cycles experimentally observed, that are the phase opposition between the complexes CLOCK:BMAL1 and PER:CRY and the chronological order of protein peaks CLOCK:BMAL1, REV-ERB and PER:CRY. The idea of using the sequence of biological events, including minima and maxima, to qualitatively describe a system behavior is much used in Boolean systems, and presents the advantage of avoiding the requirement of quantitative experimental data that can be tricky to obtain. In [24], the authors propose a method to match models with datasets based on the sequence of biological events, that they next develop in [27] and apply in [25] to qualitatively match yeast cell-cycle networks with various datasets. Another class of models, called hybrid models, that are based on a formalism introduced in [28], also considers sequences of biological events to capture the qualitative dynamic of biological networks but they additionally include temporal features [29]. In this framework, celerities are assigned to each qualitative state, to take into account the time that the system spends in each state.

The construction of the piecewise affine model presented in this paper, although including peak order and temporal constraints, strongly differs from works previously mentioned. Indeed, the system is constructed so that its solutions follow a periodic orbit that approximates the cycle of an ODE model, and exhibit certain qualitative properties, such as the peak order. These requirements are guaranteed with a set of constraints derived from the analytical solutions of our PWA model, and have the form of limited dwell-time within some regions and of threshold inequality constraints. Based on these constraints, a dedicated Algorithm allows to estimate and explore the parameter space of threshold values, leading to a realistic region that ensures existence of the circadian periodic orbit.

Therefore, through the approximation of an ODE model, our approach proposes a characterization of the circadian clock cycle and shows that the circadian system is robust in the sense that a periodic orbit with period near 24 hours is observed for a wide range of threshold values. This robustness is determined by an adjustment between the minimal concentration of PER:CRY and the maximal concentration of CLOCK:BMAL1: the product of these two quantities may range over a large interval and sets a scaling for the other variables in the system.
# 2 Numerical and mathematical characterization of circadian time

## 2.1 A model reproducing circadian clock properties

Throughout this paper, we base our analysis on the continuous circadian clock model recently developed by Almeida et al. [26]. This model, calibrated with data from mouse fibroblast cells [30], focuses on the transcriptional details and describes the dynamics of the core clock protein complexes CLOCK:BMAL1, DBP, REV-ERB (representing both REV-ERBα and REV-ERBβ) and PER:CRY (complex formed by PER2 and CRY1). The proteins are respectively denoted B, D, R and P throughout this paper. The ordinary differential equations of this model are:

$$ 
\begin{aligned}
\dot{B} &= V_R h^{-}(R) - \gamma_B B P, \\
\dot{D} &= V_B B - \gamma_D D, \\
\dot{R} &= V_D D - \gamma_R R, \\
\dot{P} &= V_D D - \gamma_B B P,
\end{aligned} 
\qquad (1) $$

where $h^{-}(R) = \frac{k_{Rr}^{4}}{k_{Rr}^{4} + R^{4}}$, $V_{R} = 44.4\%.h^{-1}$, $k_{Rr} = 80.1\%$, $V_{B} = 0.142\%.h^{-1}$, $V_{D} = 19\%.h^{-1}$, $\gamma_{R} = 0.241h^{-1}$, $\gamma_{D} = 0.156h^{-1}$, $\gamma_{B} = 2.58h^{-1}$. $V_{B}$, $V_{D}$ and $V_{R}$ stand for synthesis rates and $\gamma_{B}$, $\gamma_{D}$ and $\gamma_{R}$ for degradation rates. This set of parameters comes from the original model, (see Table 2 in [26]) and is used throughout this paper, unless stated otherwise.

## 2.2 A segmentation of the circadian time into eight stages

We suggest a segmentation of the circadian time into stages based on the repression of CLOCK:BMAL1 by REV-ERB (represented by the term $V_{R}h^{-}(R)$ in model 1) and on the dynamics between CLOCK:BMAL1 and PER:CRY (represented by the mutual repression term $\gamma_{B}BP$ in model 1). We define the different regions by setting thresholds on the variables B, P and R in order to segment the circadian cycle and characterize the clock protein dynamics. Throughout this paper, we will use “stage i” to designate the time interval spent by a circadian cycle in “region i” of the state space. By abuse of notation, we will use “stage i” to designate both the time interval and (as a label of) “region i”.

Let's define $P_{low}$, a low level of reference for the variable P, $B_{high}$ and $P_{high}$, high levels of reference for the variables B and P and $R_{int}$ an intermediate level of reference for R. The state space can be partitioned using all the different combinations of the thresholds: as there are two thresholds for P, one for R and one for B, the state space can be divided into twelve stages (see Table 1).
<table><caption>**Table 1** Segmentation of the circadian time into twelve stages and equations describing the dynamics of the circadian clock in each of the biologically-consistent stages, with function trends expected.</caption><thead><tr><th>Stage</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><td>Stage conditions</td><td>R &gt; R<sub>int</sub><br>B &lt; B<sub>high</sub><br>P<sub>low</sub> &lt; P &lt; P<sub>high</sub></td><td>R &gt; R<sub>int</sub><br>B &lt; B<sub>high</sub><br>P &gt; P<sub>high</sub></td><td>R &lt; R<sub>int</sub><br>B &lt; B<sub>high</sub><br>P &gt; P<sub>high</sub></td><td>R &lt; R<sub>int</sub><br>B &lt; B<sub>high</sub><br>P<sub>low</sub> &lt; P &lt; P<sub>high</sub></td></tr><tr><td>CLOCK:BMAL1</td><td>→</td><td>→</td><td>→</td><td>✓</td></tr><tr><td>Ḃ =</td><td>0</td><td>0</td><td>0</td><td>V<sub>R</sub> - γ<sub>B</sub>B P<sub>low</sub> &gt; 0</td></tr><tr><td>DBP</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Ḋ =</td><td>V<sub>B</sub>B - γ<sub>D</sub>D &lt; 0</td><td>V<sub>B</sub>B - γ<sub>D</sub>D &lt; 0</td><td>V<sub>B</sub>B - γ<sub>D</sub>D &lt; 0</td><td>V<sub>B</sub>B - γ<sub>D</sub>D &lt; 0</td></tr><tr><td>REV</td><td>✓ or ↗</td><td>✓ or ↗</td><td>✓</td><td>✓</td></tr><tr><td>Ṙ =</td><td>V<sub>D</sub>D - γ<sub>R</sub>R</td><td>V<sub>D</sub>D - γ<sub>R</sub>R</td><td>V<sub>D</sub>D - γ<sub>R</sub>R &lt; 0</td><td>V<sub>D</sub>D - γ<sub>R</sub>R &lt; 0</td></tr><tr><td>PER:CRY</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Ṗ =</td><td>V<sub>D</sub>D &gt; 0</td><td>V<sub>D</sub>D &gt; 0</td><td>V<sub>D</sub>D - V<sub>R</sub> &lt; 0</td><td>V<sub>D</sub>D - γ<sub>B</sub>B P<sub>low</sub> &lt; 0</td></tr><tr><td>Stage</td><td>5</td><td>6</td><td>7</td><td>8</td></tr><tr><td>Stage conditions</td><td>R &lt; R<sub>int</sub><br>B &gt; B<sub>high</sub><br>P<sub>low</sub> &lt; P &lt; P<sub>high</sub></td><td>R &lt; R<sub>int</sub><br>B &gt; B<sub>high</sub><br>P &lt; P<sub>low</sub></td><td>R &gt; R<sub>int</sub><br>B &gt; B<sub>high</sub><br>P &lt; P<sub>low</sub></td><td>R &gt; R<sub>int</sub><br>B &gt; B<sub>high</sub><br>P<sub>low</sub> &lt; P &lt; P<sub>high</sub></td></tr><tr><td>CLOCK:BMAL1</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Ḃ =</td><td>V<sub>R</sub> - γ<sub>B</sub>B P<sub>low</sub> &gt; 0</td><td>V<sub>R</sub> - γ<sub>B</sub>δ &gt; 0</td><td>-γ<sub>B</sub>B<sub>eq</sub>P<sub>low</sub><br>= -V<sub>R</sub> &lt; 0</td><td>-γ<sub>B</sub>B<sub>eq</sub>P<sub>low</sub><br>= -V<sub>R</sub> &lt; 0</td></tr><tr><td>DBP</td><td>✓</td><td>✓</td><td>✓ or ↗</td><td>✓ or ↗</td></tr><tr><td>Ḋ =</td><td>V<sub>B</sub>B - γ<sub>D</sub>D &lt; 0</td><td>V<sub>B</sub>B - γ<sub>D</sub>D &gt; 0</td><td>V<sub>B</sub>B - γ<sub>D</sub>D</td><td>V<sub>B</sub>B - γ<sub>D</sub>D</td></tr><tr><td>REV</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Ṙ =</td><td>V<sub>D</sub>D - γ<sub>R</sub>R &lt; 0</td><td>V<sub>D</sub>D - γ<sub>R</sub>R &gt; 0</td><td>V<sub>D</sub>D - γ<sub>R</sub>R &gt; 0</td><td>V<sub>D</sub>D - γ<sub>R</sub>R &gt; 0</td></tr><tr><td>PER:CRY</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Ṗ =</td><td>V<sub>D</sub>D - γ<sub>B</sub>B P<sub>low</sub> &lt; 0</td><td>V<sub>D</sub>D - γ<sub>B</sub>δ</td><td>V<sub>D</sub>D - γ<sub>B</sub>B<sub>eq</sub>P<sub>low</sub><br>= V<sub>D</sub>D - V<sub>R</sub> &gt; 0</td><td>V<sub>D</sub>D - γ<sub>B</sub>B<sub>eq</sub>P<sub>low</sub><br>= V<sub>D</sub>D - V<sub>R</sub> &gt; 0</td></tr><tr><td>Stage</td><td>9 (stage refuted)</td><td>10 (stage refuted)</td><td>11 (stage refuted)</td><td>12 (stage refuted)</td></tr><tr><td>Stage conditions</td><td>R &lt; R<sub>int</sub><br>B &lt; B<sub>high</sub><br>P &lt; P<sub>low</sub></td><td>R &gt; R<sub>int</sub><br>B &lt; B<sub>high</sub><br>P &lt; P<sub>low</sub></td><td>R &lt; R<sub>int</sub><br>B &gt; B<sub>high</sub><br>P &gt; P<sub>high</sub></td><td>R &gt; R<sub>int</sub><br>B &gt; B<sub>high</sub><br>P &gt; P<sub>high</sub></td></tr></tbody></table>
However, as the circadian clock performs one cycle, the phase opposition property between CLOCK:BMAL1 and PER:CRY has to be respected [31], refuting some stages. This property implies that the oscillator cannot have $B < B_{high}$ and $P < P_{low}$ simultaneously (eliminating stages 9 and 10) nor $B > B_{high}$ and $P > P_{high}$ simultaneously (eliminating stages 11 and 12). The circadian cycle is then described by the first eight stages of Table 1.

## 2.3 Temporal sequence of the stages

The circadian clock shows a specific order of protein peaks: BMAL1 then REV-ERB, then PER:CRY. This property suggests a specific order among stages followed by any circadian oscillator. Notice that the probability that two variables cross their thresholds simultaneously is very small. So, from one stage to another, we assume that only one variable crosses a threshold. It's therefore not possible to pass from stage 1 to 3 for example. We define the peak of $B$ such that $B$ is above its threshold $B_{high}$, corresponding to stages 5, 6, 7 and 8. The peak of $R$ follows the peak of $B$ and is defined such that $R$ is above $R_{int}$ which is represented by stages 7, 8, 1 and 2. Finally, the peak of $P$ follows the peak of $R$ and corresponds to stages 2 and 3, as $P > P_{high}$. A cycle that respects physiological properties (such as the order of protein peaks) should not enter twice a same stage in a same cycle and should follow an ascending and periodic order of stages: 1-2-3-4-5-6-7-8-1. This sequence of transitions is represented by green arrows in Figure 1. Each stage is represented by a vertex in the space $(R, B, P)$ and is defined via the thresholds $R_{int}, B_{high}, P_{low}$ and $P_{high}$. The unlabelled vertices correspond to the discarded stages.

<figure><img src="image_2.png" /><figcaption>Fig. 1 Ideal physiological order of stages followed by an oscillator as it performs one cycle</figcaption></figure>

## 3 A piecewise affine clock model

Inspired by the piecewise affine model suggested before by our group [32] and based on the repression of CLOCK:BMAL1 by PER:CRY, on the dynamics between REV-ERB and CLOCK:BMAL1 and on the circadian time decomposition into eight stages, we simplified the nonlinear continuous model of Almeida et al. (1) into a piecewise affine one, in order to characterize the dynamic of the circadian clock cycle. A scheme of the regulatory interaction network of the continuous model, including the negative feedback loops and the nonlinear mechanisms, is shown in Figure 2.
The goal is to provide an analytic treatment of the circadian system, based on the three fundamental variables that are also represented in practically every model of the circadian clock, to obtain:

* analytic solutions in each region of the piecewise affine system, leading to sufficient conditions for a periodic trajectory of the circadian system;

* a characterization of the cycle dynamics in terms of four fundamental threshold parameters;

* robustness of the circadian cycle, in terms of threshold parameters.

<figure><img src="image_3.png" /><figcaption>**Fig. 2** Network representation of the circadian clock ODE model. The two negative feedback loops are highlighted by blue and red backgrounds. Activating interactions are shown by green arrows, and are described by linear terms in the model. CLOCK:BMAL1 inhibition by REV:ERB is represented by a red flat-head arrow and is described, in the ODE model, by the decreasing Hill function $h^{-}(R) = \frac{k_{Rr}^{4}}{k_{Rr}^{4}+R^{4}}$. Complex formation of PER:CRY and CLOCK:BMAL1 is shown by blue arrows and is described, in the model, by a mass action term, $g(B, P) = \gamma_{B}BP$. The piecewise affine model is based on this ODE model and has the same basic structure. The nonlinear functions $f$ and $g$ are approximated by affine terms whose definitions depend on the position of variables relative to their thresholds. See Table 1 for more details.</figcaption></figure>

## 3.1 Model design

The first nonlinear term of system (1) is the decreasing Hill function $V_{R}h^{-}(R)$ which can be simplified by a step function [17]. A straightforward choice for the threshold is $R_{int}$, which gives:

$$V_{R}h^{-}(R) \approx \begin{cases} 0, & R > R_{int} \\ V_{R}, & R < R_{int} \end{cases} \qquad (2)$$

The second nonlinear term is the function from the complex formation of PER:CRY and CLOCK:BMAL1, $g(B,P)=\gamma_{B}BP$. To linearize this polynomial term, we rely on the dynamic of the ODE model (see [32]) and propose a solution in which this term is either defined null, linear or constant according to the stages. The same thresholds as for the cycle segmentation are used: $B_{high}$, $P_{low}$ and $P_{high}$.

First, $B$ remains close to its level $B_{high}$ (whilst $R$ and $P$ are large) implying that the two terms $V_{R}h^{-}(R)$ and $\gamma_{B}BP$ have similar values and add up to 0 in the $B$
equation, leading to the following simplification:

$$ \gamma_B B P \approx \begin{cases} 0, & R > R_{int}, B < B_{high}, P_{low} < P < P_{high} \\ 0, & R > R_{int}, B < B_{high}, P > P_{high} \\ V_R, & R < R_{int}, B < B_{high}, P > P_{high} \end{cases} \quad (3) $$

These 3 equations correspond to the first three stages.

Next, we consider the dynamics between B and P. P decreases and tends to reach its low level $P_{low}$. B increases and so the degradation term is mostly dominated by B. The degradation term of B and P during stages 4 and 5 is then defined by:

$$ \gamma_B B P \approx \begin{cases} \gamma_B B P_{low}, & R < R_{int}, B < B_{high}, P_{low} < P < P_{high} \\ \gamma_B B P_{low}, & R < R_{int}, B > B_{high}, P_{low} < P < P_{high} \end{cases} \quad (4) $$

At the end of stage 5, P reaches its low level $P_{low}$ and B becomes close to a reference state $B_{eq}$, corresponding to the equilibrium between its production and degradation terms, $\dot{B} = V_R - \gamma_B B_{eq} P_{low} = 0$ which gives:

$$ B_{eq} = \frac{V_R}{\gamma_B P_{low}}, \quad (5) $$

a relationship between $B_{eq}$ and $P_{low}$.

Then, during stage 6, B reaches its maximum and P its minimum. We define the degradation term by a constant $\delta$ such that during stage 6 there is a time $t$ where $\dot{P}(t) = 0$:

$$ \gamma_B B P \approx \begin{cases} \gamma_B \delta, & R < R_{int}, B > B_{high}, P < P_{low} \end{cases} \quad (6) $$

The value $\delta$ stands for a product between a maximal B and a minimal P.

During stages 7 and 8, B decreases and tends to reach $B_{eq}$ while P increases and crosses the threshold $P_{low}$. The degradation term is here described by:

$$ \gamma_B B P \approx \begin{cases} \gamma_B B_{eq} P_{low}, & R > R_{int}, B > B_{high}, P < P_{low} \\ \gamma_B B_{eq} P_{low}, & R > R_{int}, B > B_{high}, P_{low} < P < P_{high} \end{cases} \quad (7) $$

Finally, we define the time at which the cycle exits stage i by $t_i$ and the duration of each stage by the difference of exiting times:

$$ d_i = t_i - t_{i-1} \text{ with } i = 1..8 \text{ and } t_0 = t_8. \quad (8) $$

We set the numerical starting point of the cycle as a point reached by the oscillator at the end of stage 5 (so $t_5 = 0$), when P crosses its threshold, so $P_0 = P_{low}$. The other initial conditions are based on the balance between production and degradation rates at the end of stage 6:

$$ P_0 = P_{low}, \quad B_{high} < B_0 \le B_{eq}, \quad D_0 = \frac{V_B B_0}{\gamma_D}, \quad R_0 = \frac{V_D V_B B_0}{\gamma_D \gamma_R} \quad (9) $$
In what follows, the initial condition for the PWA system is a point at the boundary between stages 5 and 6, of the form $(B_0, D_0, R_0, P_0)$.

With these definitions, an analytical solution can be given for the dynamics of the model in each stage. Notice that these conditions and equations ensure positivity of B, D and R, provided that the threshold parameters are well-defined (see Section 3.2). As P reaches its minimum during stage 6, δ should also be defined carefully to ensure the positivity of P (in more detail below in Section 3.3). Table 1 summarizes the definitions of each stage, as well as the equations for each variable and the function trends expected to obtain oscillations and a consistent model, that is to say an oscillator that goes through all the stages in the specified physiological order.

## 3.2 Circadian cycle as a sequence of transitions between stages

The following assumptions and propositions summarize sufficient conditions on the values of the threshold parameters, of the scaling parameter δ and on the stage durations, to ensure that the analytical solutions represent an oscillator following the cycle transitions in the specified physiological order for one period (as defined in Section 2.3). Sketches of proofs are given in Appendix A.

First, some assumptions among the parameters are as follows:

$$A.1 \quad \frac{V_R}{3\gamma_B} < \delta < \frac{V_R}{\gamma_B},$$

$$A.2 \quad B_{high} < B_{eq} = \frac{V_R}{\gamma_B P_{low}},$$

$$A.3 \quad P_{low} > \frac{2V_R}{3\delta\gamma_B - V_R} \frac{V_D V_B}{\gamma_D \gamma_B},$$

$$A.4 \quad \frac{V_D V_B V_R}{\gamma_D \gamma_R \gamma_B P_{low}} < R_{int} < \frac{V_R}{\gamma_R}$$

Considering the requirements for the system to be in stage 6, Assumption A.1 ensures that B increases and remains above $B_{high}$ and Assumption A.4 ensures the initial starting point $R_0$ lower than $R_{int}$. In addition to A.1 - A.4, the following assumptions on variables $D(t)$ and $R(t)$ are also needed. They are all verified by numerical results (see in Appendix Figure A1):

$$A.5 \quad \frac{V_B B_{eq}}{\gamma_D} < D_{min},$$

$$A.6 \quad D(t_i) < \frac{V_R}{V_D}, \quad i \in \{2, 3, 4, 5\},$$

$$A.7 \quad \frac{V_R}{V_D} < D(t_6) < D(t_7),$$

$$A.8 \quad R_{int} < R(t_7) < R(t_8)$$

The next propositions state the conditions needed at each step for transition from one stage to the next. At the starting point of the cycle, only the values of $P_{low}$ and $\delta$ are set. Then, at each step, the idea is to define the value of a specific threshold by choosing the duration $d_i = t_i - t_{i-1}$ of stage $i$, such that the desired threshold is crossed first. The duration $d_i$ of stage $i$ also depends on the duration $d_{i-1}$ of stage $i-1$. Since the starting point of the system occurs at the end of stage 5 (at $t_5$),
the duration of stage 6 is the first one to be chosen. However, this duration directly impacts the good progress of the cycle, since a too short or too long duration blocks the beginning of the oscillations (see conditions of Proposition 1). Duration of stage 8 also requires precise conditions on the time intervals, to guarantee that the expected variable crosses its given threshold. Thus, for stages 6 and 8 (Propositions 1 and 2) an interval for the stage durations is provided, say $d_i \in (d_{i,min}, d_{i,max})$, to allow for flexibility in the threshold parameter choices.

Note that the transitions between two consecutive stages from stage 2 to 6 (see Proposition 4) follow naturally from the Assumptions A.1 - A.8 and from Propositions 1 to 3 and require no further conditions.

After one period, that is to say when the system reaches again the boundary between stage 5 and stage 6, $B(t_5) \approx B_0$ and $P(t_5) = P_0$. Indeed, during stage 5, $B(> B_{high})$ increases and converges quickly to $B_{eq}$, with $B_{eq} \ge B_0 > B_{high}$ and the end of stage 5 corresponds to the moment where $P$ equals $P_{low}$ (= $P_0$)). Then, although $D(t_5)$ and $R(t_5)$ are observed bigger than $D_0$ and $R_0$, the good progress of the cycles of the next periods through the stages is not impacted.

**Proposition 1** (Stage 6 to 7 to 8). Consider the PWA system as defined in Table 1. Assume that the following inequalities hold:

1(i) $d_{6min} < d_6$ such that $D(t_5 + d_{6min}) > \frac{V_R}{V_D}$ and $\exists d_7 < \frac{(V_R - \delta\gamma_B)d_{6min}}{V_R}$ such that $P(t_5 + d_{6min} + d_7) = P_{low}$;

1(ii) $d_6 < d_{6max}$ such that $P(t_5 + d_{6max}) < P_{low}$ and $R(t_5 + d_{6max}) < \frac{V_R}{\gamma_R}$.

Then, starting at the boundary between stages 5 and 6 with initial condition $(B_0, D_0, R_0, P_0)$ defined in (9) the system will successively cross stage 6 and enter stage 7 and cross stage 7 and enter stage 8.

Condition 1(i) ensures that $P$ increases at the end of stage 6 and that $d_6$ is long enough to ensure a solution of the equation $P(t_6 + d_7) = P_{low}$, with $d_7$ such that $P$ crosses $P_{low}$, before $B$ crosses $B_{high}$. Condition 1(ii) ensures that $R$ crosses $R_{int}$ before $P$ crosses $P_{low}$, and such that $R_{int} < \frac{V_R}{\gamma_R}$. These conditions guarantee transitions from stage 6 to 7 and from stage 7 to 8.

**Proposition 2** (Stage 8 to 1). Consider the PWA system as defined in Table 1. In addition to the inequalities of Proposition 1, assume that the following inequalities hold:

2(i) $d_{8min} < d_8 < d_{8max}$ with $d_{8min} = \frac{B(t_7) - \frac{V_R}{\gamma_B P_{low}}}{V_R}$ and $d_{8max} = \frac{B(t_7)}{V_R}$;

2(ii) $P(t_7 + d_{8max}) \ll P_{high}$.

Then, starting from stage 8, with initial conditions $B_0 = B(t_7)$, $D_0 = D(t_7)$, $R_0 = R(t_7)$, $P_0 = P(t_7) = P_{low}$, the system will cross stage 8 and enter stage 1.

Conditions 2(i) and 2(ii) ensure that $B$ crosses $B_{high}$ (such that $0 < B_{high} < B_{eq}$) before $P$ crosses $P_{high}$, which guarantees transition from stage 8 to 1.
**Proposition 3 (Stage 1 to 2).** Consider the PWA system as defined in Table 1. In addition to the inequalities of Propositions 1 and 2, assume that the following inequality holds.

$3(i) P_{high} < P(t_8 + d_{1max}),$ such that $R(t_8 + d_{1max}) = R_{int}$.

Then, starting from stage 1, with initial conditions $B_0 = B(t_8) = B_{high}$, $D_0 = D(t_8)$, $R_0 = R(t_8)$, $P_0 = P(t_8)$, the system will cross stage 1 and enter stage 2.

Condition 3(i) ensures that $P$ crosses $P_{high}$ before $R$ crosses $R_{int}$, which guarantees transition from stage 1 to 2.

**Proposition 4 (Stage 2 to 6).** Consider the PWA system as defined in Table 1. Assume that inequalities of Propositions 1-3 hold. Then, starting from stage 2, with initial conditions $B_0 = B(t_1)$, $D_0 = D(t_1)$, $R_0 = R(t_1)$, $P_0 = P(t_1) = P_{high}$, the system will cross stage 2 and successively evolve through stages 3, 4, 5 and finally 6, in this order.

### 3.3 Algorithm for numerical estimation of circadian cycle threshold parameters

We use the conditions defined before in the propositions and propose Algorithm 1 to estimate the five-dimensional parameter space, composed of the four thresholds $B_{high}$, $R_{int}$, $P_{low}$ and $P_{high}$ and the scaling factor $\delta$, such that any oscillator of the PWA model follows the specified physiological stage order for several periods. Both definitions and conditions of the system ensure positivity of $B$, $D$, and $R$ but positivity of $P$ has to be verified, so we add a condition in this sense (see steps 3 and 4 of Algorithm 1).

As input of Algorithm 1, we take the same parameters as the ones of the continuous model 1, except for $\gamma_R$. This parameter has a big impact on the period and allows us to adjust the period around 24h for a better coherence with biological clocks (the initial value of $\gamma_R$ gave periods around 33h).
<caption>**Algorithm 1.** Numerical estimation of the threshold parameters.</caption>

0. **Input:** $V_R = 44.4$, $k_{Rr} = 80.1$, $V_B = 0.142$, $V_D = 19$, $\gamma_R = 0.7$, $\gamma_D = 0.156$, $\gamma_B = 2.58$.
1. Pick $\delta$ such that $\frac{V_R}{3\gamma_B} < \delta < \frac{V_R}{\gamma_B}$.
2. Pick $P_{low}$ such that $P_{low} > \frac{2V_R}{3\delta\gamma_B - V_R}\frac{V_D V_B}{\gamma_D \gamma_B}$, and define $B_{eq} = \frac{V_R}{\gamma_B P_{low}}$, $B_0 = B_{eq}$, $D_0 = \frac{V_B B_{eq}}{\gamma_D}$, $R_0 = \frac{V_D V_B B_{eq}}{\gamma_D \gamma_R}$, $P_0 = P_{low}$.
3. Compute $t_{P_{min}}$ solving $D(t_{P_{min}}) = \frac{\delta \gamma_B}{V_D}$, with initial conditions of 2.
4. Check:
    * if $P(t_{P_{min}}) > 0$ go to 5.
    * if $P(t_{P_{min}}) < 0$ go back to 1 and pick another $\delta$.
5. Pick $d_6$ such that $d_{6min} < d_6 < d_{6max}$ with $d_{6min}$ such that $D(t_5 + d_{6min}) > \frac{V_R}{V_D}$ and $d_{6max}$ such that $P(t_5 + d_{6max}) < P_{low}$ and $R(t_5 + d_{6max}) < \frac{V_R}{\gamma_R}$
6. Compute $d_7$ such that $P(t_6 + d_7) = P_{low}$.
7. Check:
    * if $d_7 < \frac{(V_R - \delta \gamma_B) d_6}{V_R}$ go to 8.
    * if there is no solution or if $d_7 > \frac{(V_R - \delta \gamma_B) d_6}{V_R}$ go back to 5 and pick another $d_6$.
8. Compute $R_{int}$ such that $R_{int} = R(t_5 + d_6)$.
9. Pick $d_8$ such that $\frac{B(t_7) - \frac{V_R}{\gamma_B P_{low}}}{V_R} < d_8 < \frac{B(t_7)}{V_R}$.
10. Compute $B_{high}$ such that $B(t_7 + d_8) = B_{high}$.
11. Compute $d_{1max}$ such that $R(t_8 + d_{1max}) = R_{int}$.
12. Compute $P_{high}$ such that $P(t_8) + p_{add} = P_{high}$, with $p_{add} = 200$ initially.
13. Check:
    * if $P(t_8 + d_{1max}) > P_{high}$ go to 14.
    * if $P(t_8 + d_{1max}) < P_{high}$ go back to 12 with $p_{add} = p_{add} - 10$.
14. Check if the oscillator goes through all the eight stages in the specified physiological order for several periods and save the threshold parameter values and $\delta$.
15. **Output:** set of parameters describing a circadian cycle: $\delta$, $P_{low}$, $R_{int}$, $B_{high}$, $P_{high}$.

## 4 Continuous and PWA models share qualitative properties

Since our goal is to faithfully approximate the dynamics of the continuous model (1) by piecewise affine equations, we theoretically show that the PWA model reproduces the same qualitative properties as the continuous model, that we next illustrate with numerical simulations. The five-dimensional parameter space of the PWA model is estimated by computing with Algorithm 1 around 2000 combinations of values for the thresholds and $\delta$ (see Section 5.1 for more details). The results show that the PWA
model is robust in the sense that the oscillators built from the set of threshold parame-
ters generated by Algorithm 1 all have the same qualitative behaviour, independently
of the exact values of the thresholds.

**4.1 Theoretical robustness of qualitative properties**

The construction of the PWA system provides a robust order of protein peaks and of
maximal values and ensures the phase opposition between *B* and *P* for any oscillator
whose threshold parameters were generated by Algorithm 1. Indeed Algorithm 1 satisfies a set of constraints on the threshold parameter values which ensures the progress
of the oscillator through the eight stages in the expected order, which itself guarantees the peak-ordering and phase opposition properties, as shown in Proposition 5.

**Proposition 5.** *Under Assumptions A.1-A.8 and Propositions 1 to 4, the following hold:*

(a) The maximum of *B* is reached at the boundary between stages 6 and 7 and *B* is at its minimum during stages 1, 2 and 3;

(b) The maximum of *P* is reached at the boundary between stages 2 and 3 and the minimum of *P* is reached during stage 6;

(c) The maximum of *R* is reached in between the maxima of *B* and *P*.

*Proof:* (a) The variable *B* peaks when *B* is above its threshold $B_{high}$, which is during stages 5 to 8. Especially, since during stages 5 and 6 *B* strictly increases (in stage 5, *B* increases and tends to $B_{eq}$ and in stage 6 $\dot{B} = V_R - \gamma_B \delta > 0$ with Assumption A.1) and during stages 7 and 8 *B* strictly decreases ($\dot{B} = -V_R < 0$), *B* reaches its maximal value at the boundary between stages 6 and 7 for any set of parameters returned by Algorithm 1. During stages 1, 2 and 3, *B* is constant and equal to its minimal value, which is just below $B_{high}$.

(b) *P* peaks when *P* is above $P_{high}$, which is during stages 2 and 3. The boundary between these two stages marks the point where *P* reaches its maximal value, with *P* strictly increasing during stage 2 ($\dot{P} = V_D D > 0$) and *P* strictly decreasing during stage 3 ($\dot{P} = V_D D - V_R < 0$ with Assumption A.6). *P* reaches its lowest values during stages 6 and 7, when *P* is below $P_{low}$. The system definition and Proposition 1 induce *P* minimal value during stage 6.

(c) *R* peaks when *R* is above its threshold, which is during stages 7, 8, 1 and 2. Since the maximal value of *B* (transition from stages 6 and 7) corresponds to the time where *R* increases above $R_{int}$, and since the maximal value of *P* (transition from stages 2 and 3) corresponds to the time where *R* decreases below $R_{int}$, the maximal value of *R* necessarily occurs between the maximal values of *B* and *P* (and is in fact always observed during stage 1). □

For any set of parameters returned by Algorithm 1, (a) and (b) demonstrate the
phase opposition between *B* and *P* with *B* (resp. *P*) reaching its maximal value in the
same stage as the minimal value of *P* (resp. *B*), and (c) shows the peak order *B*, *R*, *P*.
## 4.2 Cycle oscillations and phase portraits

From a numerical point of view, Figure 3(a) shows an example of an oscillator whose threshold values have been computed by Algorithm 1 while Figure 3(b) represents the oscillations of the continuous model. As might be expected, the maximal and minimal amplitudes of each protein vary between the two model formalisms, with the maximal value of *P* increasing by a factor of 1.5 and the maximal value *R* decreasing by about the same factor, while maximal value of *B* remains similar (for the PWA oscillator compared to the continuous one). In any case, the order among amplitudes and the phase opposition between *B* and *P* are maintained, as are the forms of the curves and their dynamics.

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Oscillations and cycle segmentation for (a) one oscillator generated by Algorithm 1 with $\delta=12.55$, $B_{high}=0.10$, $R_{int}=45.21$, $P_{low}=71.28$ and $P_{high}=290.55$ and (b) for the continuous model (1).</figcaption></figure>

These observations are confirmed by Figure 4, which shows the phase portraits of the continuous model with those of several oscillators from the PWA model, with different values of threshold parameters and of the scaling factor $\delta$. Figure 4(a) highlights the phase opposition property, recovered by both models: when *B* reaches its highest values, *P* is low and inversely. Similarly, Figure 4(b) shows the similar dynamic between *P* and *R*, shared by all the oscillators: successively, *R* increases and reaches its highest values while *P* is at its intermediate ones, and then, as *R* decreases, *P* increases, peaks and decreases. From these figures, it is clear that the simplified PWA model closely represents the dynamical properties of the continuous model over a large range of its threshold parameters.

# 5 Cycle robustness: quantitative analysis and critical elements

We first introduce the following elements which will next serve as comparative analysis tools between all the oscillators generated by Algorithm 1. Letting $X \in \{B, R, P\}$, define:
<figure><img src="image_5.png" /><figcaption>Fig. 4 Phase portrait projection on the planes: (a) B and P and (b) R and P for the continuous model (dashed black line) and for the PWA model (with different threshold parameters and one color stands for one δ value). Curves are normalized between 0 and 1 for comparison.</figcaption></figure>

1. $A_{X_{\phi}}$: threshold value (called $X_{\phi}$) relative to the amplitude ($X_{max} - X_{min}$) of the corresponding variable (X):

$$A_{X_{\phi}} = \frac{X_{\phi} - X_{min}}{X_{max} - X_{min}} \cdot 100 \qquad (10)$$

2. $\tau_{X_{\phi}}$: time in hours that each variable (X) remains above its corresponding threshold ($X_{\phi}$) over one period (T), corresponding also to protein peak duration:

$$\tau_{X_{\phi}} = \frac{t_{end} - t_{in}}{T}, \text{ such that } X(t) > X_{\phi}, \forall t \in [t_{in}, t_{end}] \qquad (11)$$

## 5.1 Five-dimensional parameter space of the PWA model

The parameter space of the PWA model is composed of the four threshold parameters: $B_{high}$, $R_{int}$, $P_{low}$ and $P_{high}$, and of one scaling factor: $\delta$. Figure 5 represents this five-dimensional parameter space as estimated by Algorithm 1 but, for a clearer comparison, the figure represents the relative amplitudes for each threshold parameter, as defined in equation (10). Thus, relatively to the amplitude of their corresponding variables, the values of $R_{int}$, $P_{low}$ and $P_{high}$ allowing the expected oscillations can be taken in a large range of concentrations. Notice nevertheless that $A_{R_{int}}$ and $A_{P_{low}}$ increase or decrease directly with $\delta$ while the relation between $\delta$ and $A_{P_{high}}$ is less clear. The way the system is constructed implies that the minimal value of B will always be close to $B_{high}$ and thus $A_{B_{high}}$ is small.

Although the region of threshold parameter values is wide, the periods (called T) of the oscillations of the PWA model remain in a small interval and are consistent with a circadian clock period, from 22 to 26h.

## 5.2 Duration of protein peaks

Figure 6(a) represents the peak durations, as defined in equation (11), for each set of thresholds computed by Algorithm 1. Notice first that the duration of REV-ERB
<figure><img src="image_6.png" /><figcaption>**Fig. 5** (a), (b) Sets of threshold values computed by Algorithm 1 ensuring a correct representation of circadian clocks oscillations relatively to the amplitude of their corresponding variable (in %), and their associated period (T). One hexagon corresponds to one oscillator and its associated set of threshold and scaling parameters. One color for each $δ$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6** Duration of protein peaks (a) and stage durations (b) in hours. One hexagon corresponds to one set of parameters. One color for each $δ$.</figcaption></figure>

peak is inversely proportional to CLOCK:BMAL1 peak and to the scaling factor $δ$, highlighting one negative feedback loop characterizing circadian clocks: the more REV-ERB is expressed, the more CLOCK:BMAL1 is repressed, and inversely. The relation between $τ_{Bhigh}$ and $τ_{Phigh}$ is less clear since the value of $P_{high}$ is assessed last by Algorithm 1 and does not have a major influence on the correct progress of the cycle through the stages.
Then, the peak of *B* lasts between 7h and 10h and the peak of *P* lasts between 10h and 14h. The peak of *R* shows a larger range of durations, from 7 to 13h. Although no constraint about protein peak durations were given, the durations of *B* and *P* peaks reasonably correspond to biological results [31]. Concerning *R*, the oscillators with the biggest $\delta$ ($\delta \gtrsim 11$) have a shorter peak of *R* (between 7 and 10h) and seem more consistent with biological results [33].

### 5.3 Duration of stages and relation to daylight/night intervals

As the circadian oscillator performs one cycle, it goes through well-marked steps highlighted by our stages. Figure 7 represents minimal and maximal values of stage durations after one period according to conditions from Propositions 1 to 4, while Figure 6(b) represent stage durations for the set of oscillators generated by Algorithm 1, after several periods.

As explained in Section 3.2, the duration of stage 6 is the first one to be defined, and directly depends on the values of $\delta$ and $P_{low}$. It varies from 1 to 3h and increases with $\delta$ and $P_{low}$. Then, stages 7 and 8 are very short and last less than one hour. Duration of stages 1 and 2 evolve inversely with $\delta$. During stage 1, which lasts from 2.5 to 5h, *P* progressively increases, repressing *B* while *R* is still expressed. Stage 2 lasts from 3.5 to 8h and a short stage 2 represents an earlier degradation of *R* compared to a long stage 2. Concerning stage 3, its duration seems to evolve similarly to $\delta$, and goes from 3 to 8h. Stage 4 is short and lasts less than half an hour. Finally, duration of stage 5 is around 6 hours and stays consistent from an oscillator to another.

Next, we draw a closer comparison between our analytic stages and the experimental Circadian Time phases as described by Takahashi *et al.* [31] in the mouse liver and establish a correspondence with daylight and night hours. At Circadian Time CT0, the circadian cycle is in a poised state, until transcriptional activity of CLOCK:BMAL1 is derepressed. Then, until CT12, (subjective daylight hours), both CLOCK:BMAL1 transcriptional activity and REV-ERB expression are observed, with a respective peak at approximately CT8 [31] and CT8-10 [31, 33–35]. Subjective night hours extend from CT12 to CT24 and are characterized by high levels of PER:CRY.

In our segmentation, stage 4 corresponds to a weak expression of all the variables and the cycle enters stage 4 when *P* decreases under its threshold $P_{high}$, while *R* and *B* are low and leaves stage 4 when *B* increases above $B_{high}$. 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 and can be associated at the end of the poised state [31].

Then, since the peak of *B* is represented by stages 5 to 8 (when $B > B_{high}$) and since we observed the maximal value of *R* during stage 1, we assume that daylight hours are represented by stages 5 to 1, with stage 1 marking the transition between daylight and night hours, around CT12. The peak of *P* occurs during stages 2 and 3 (when $P > P_{high}$), and these two stages can be associated to night hours.
<figure><img src="image_8.png" /><figcaption>**Fig. 7** (a) Minimal and maximal possible duration of stage 6 over $\delta$ according to conditions of Proposition 1 (lines in blue tones for $d_{6min}$ and in orange tones for $d_{6max}$). In each color, a darker line corresponds to a bigger $P_{low}$. (b) Durations of stage 7 over $\delta$ according to the minimal (lines in blue tones) and maximal (lines in orange tones) duration of stage 6. (c) Minimal and maximal possible duration of stage 8 over $\delta$ according to conditions of Proposition 2 (dashed lines for $d_{8min}$ and solid lines for $d_{8max}$) and according to the minimal and maximal duration of stage 6. (d-h) Durations of stages 1 to 5 are computed based on the minimal and maximal durations of stages 6 and 8. The lines use the same color and style code as for (a-c).</figcaption></figure>

We can thus associate night hours to stages 2 and 3 and daylight hours to stages 4, 5, 6, 7, 8 and 1. In our model, night hours last around 10 to 14h and daylight hours are in a narrower interval, and are around 11-12h.
## 5.4 Two critical dynamical elements

Our analysis of stage durations and transitions identifies two critical moments during the circadian cycle: the first is the estimation of the duration of stage 6, where the system is defined by the scaling factor $\delta$ and which allows for $P_{low}$ and $R_{int}$ assessment, and the second is the evaluation of the duration of stage 8, which allows for $B_{high}$ and $P_{high}$ assessment.

In our piecewise affine model, the scaling factor $\delta$ represents the product of the concentrations of PER:CRY and CLOCK:BMAL1, when they are in phase opposition with PER:CRY at its minimum value. In fact, during this interval, it is difficult to approximate the dynamics of the degradation term $\gamma_B P$ by a linear term in the continuous model. The best (and simpler) option was to consider that the product is constant and equal to some $\delta$ for a suitable time interval. And, indeed, this value defines a kind of scaling of the cycle that directly impacts the first threshold parameters to be defined ($P_{low}$ and $R_{int}$), varying them in direct proportion to $\delta$ (see Figure 5(a)).

It is interesting to note that, as the factor $\delta$ ranges over a large interval, it sets the scale for the other threshold parameters while keeping the total period close to 24 hours (in an interval between 22 and 26 hours).

Another critical point in the cycle appears at stage 8, when extra conditions on the crossing times of the variables $B$ and $P$ are needed to guarantee the correct order of transition (see Proposition 2). In fact, stage 8 marks the last moments before CLOCK:BMAL1 repression and is of short duration (less than 1h). $P$ is at an intermediate value and increasing, and should effectively repress $B$ before $P$ itself crosses $P_{high}$ (since $P$ and $B$ cannot be simultaneously above their highest thresholds). The variable $R$ should wait for $B$ repression before it also starts decreasing, which is verified by Assumption A.8. In other words, the dynamics at stage 8 potentially allows for different outcomes, so time control needs to be more strict in order to rule out the undesirable outcomes.

Globally, a well-defined duration of the stages 6-8 allows for a good assessment of all the threshold parameters, ensuring the expected progress of the oscillator through the stages and sustained oscillations.

## 6 Discussion

Our study is one of the first to propose a more precise, analytical and quantitative characterization of circadian clock cycle dynamics based on the approximation of a continuous model [26] including mass-action terms by piecewise affine (PWA) equations. In this perspective, we proposed a segmentation of the circadian cycle into stages and to focus on a periodic sequence of transitions, consistent with qualitative properties of circadian cycles. Our main contribution is to propose a method to identify the parameter space of real threshold values leading to oscillators with the same periodic dynamics, defined by the specific sequence of transitions. This sequence is ensured by a set of inequality constraints on the threshold parameters and on some stage durations established from analytical solutions. An Algorithm satisfying all these conditions allows to explore, according to the scaling value $\delta$, the region of threshold values. Our analysis shows that the PWA model allows for a large range of threshold
parameters that exhibit the same dynamical behaviour, with a periodic orbit of about 24h.

An interesting feature of this model is the multiple roles of the thresholds. First used to divide the state space into well-defined regions, their real values turn out to be important to guarantee a periodic trajectory for the system. They finally allow to assess and show robustness of our system. Another interesting feature is the factor $\delta$ representing the product of the minimal concentration of PER:CRY and maximal of CLOCK:BMAL1. This parameter provides a scaling for the PWA model such that, for each value of $\delta$, a corresponding set of threshold parameters exists, that generates an ordered circadian cycle with period around 24 hours. In this way, a large range of allowed concentrations is established for the PWA, supporting high robustness of the system based on an inverse proportionality between PER:CRY and CLOCK:BMAL1, $B = \frac{\delta}{\gamma_B P}$, as the maximal concentration of CLOCK:BMAL1 becomes high, the minimum of PER:CRY should become smaller. Thus, among other quantitative points, our analysis indicates that the two concentrations minimal PER:CRY and maximal CLOCK:BMAL1 play a large role in setting the oscillatory regime of the cycle, and guarantee its robustness over a wide range of threshold values.

Our characterization of the circadian cycle involved the approximation of a mass-action term into piecewise affine equations, which was made possible by the knowledge of qualitative properties of the ODE model we approximated. An open direction remains the development of a general method for approximating mass-action terms through piecewise affine expressions in a suitable partition of the state space. Nevertheless, our analysis of the term $\gamma_B P$ may be useful for application in other clock models, where a similar term represents protein sequestering through binding. Some examples include models for *Neurospora crassa* [36], where other proteins play roles similar to BMAL1, PER and CRY.

**Acknowledgments.** We thank Franck Delaunay for many enlightening discussions. We also thank one of the referees for suggesting several interesting points of view and generalizations for our method.

## Declarations

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

* Competing interests. The authors declare no competing interests.

* Availability of data and materials. This paper generates a new piecewise affine model which is based on a previously published ODE model (cited in the bibliography, and described in the text). An algorithm is then proposed to explore and characterize the threshold parameter space of the new model. Numerical simulations of this algorithm generate the data analysed in the paper. All parameters and conditions are provided to allow for computational reproduction of the data. The codes underlying the results of this study are available from gitlab.
Appendix A Illustration of Assumptions A.6-A.8
and proofs of Propositions 1 to 4

<figure><img src="image_9.png" /><figcaption>**Fig. A1** Numerical simulations justifying Assumptions A.6-A.8 (see Section 3.2). (a-e) Values of $D(t_i)$, $i \in \{2, 3, 5, 6, 7\}$ to compare with $\frac{V_R}{V_D}$: illustration of Assumptions A.6 and A.7. (f) $R(t_8)-R(t_7)$ to compare with 0: illustration of Assumption A.8. Each ending stage value depends on $\delta$, $P_{low}$ (a darker line corresponds to a bigger $P_{low}$), the duration of stage 6 (lines in blue tones for $d_{6min}$ and in orange tones for $d_{6max}$) and the duration of stage 8 (dashed lines for $d_{8min}$ and solid lines for $d_{8max}$). See caption of Figure 7 for more details.</figcaption></figure>

*Proof of Proposition 1:* Consider the system during stage 6, that is to say for $t_5 < t < t_6$. Since $\dot{B} = V_R - \gamma_B \delta$, Assumption A.1 implies that $B(t)$ strictly increases, so $B(t)$ remains higher than $B_{high}$.
Then, we can show that in $[t_5, t_6]$, $D(t)$ increases. Indeed, $\dot{D}(t) = V_B B(t) - \gamma_D D(t)$. Replacing $B(t)$ and $D(t)$ by their analytical solutions during stage 6, with initial conditions $B_{high} < B_0 \le B_{eq}$ and $D_0 = \frac{V_B B_0}{\gamma_D}$,

$$ \dot{D}(t) = \frac{V_B(V_R - \delta\gamma_B)}{\gamma_D} (1 - e^{-\gamma_D(t-t_5)}). $$

Assumption A.1 gives $\frac{V_B(V_R - \delta\gamma_B)}{\gamma_D} > 0$. Since $t - t_5 > 0$, $\dot{D}(t) > 0$ and $D(t)$ increases. Conditions A.1 and A.3 imply a decrease of $P(t)$ at the beginning of stage 6:

$$ \dot{P}(t_5) = V_D D(t_5) - \delta\gamma_B \le \frac{V_D V_B V_R}{\gamma_D \gamma_B P_{low}} - \delta\gamma_B < 0. $$

Next, set $\sigma = \frac{V_D V_B (V_R - \gamma_B \delta)}{\gamma_D}$ and notice that,

$$ \dot{R}(t) = V_D D(t) - \gamma_R R(t) = \frac{\sigma}{\gamma_D - \gamma_R} (e^{-\gamma_D(t-t_5)} - e^{-\gamma_R(t-t_5)}) + \frac{\sigma}{\gamma_R} (1 - e^{-\gamma_R(t-t_5)}) $$

which satisfies $\dot{R}(t_5) > 0$, so $\dot{R}(t) > 0$ for some time interval after entering stage 6. To guarantee that the system evolves to stage 7, it is necessary that $P(t)$ increases at the end of stage 6 and that $R(t)$ reaches $R_{int}$ before $P(t)$ reaches $P_{low}$. This is assured by conditions 1(i) and 1(ii), which provide appropriate conditions on the duration $d_6 = t_6 - t_5$. The system enters stage 7.

Consider the system during stage 7, that is $t_6 < t < t_7$. Since $\dot{B}(t) = -V_R$ in this interval, $B(t)$ decreases according to the expression $B(t) = B(t_6) - V_R(t - t_6)$. By Assumption A.7, $D(t_6) > \frac{V_R}{V_D}$. At the beginning of stage 7 we have $R(t_6) = R_{int} < \frac{V_R}{\gamma_R}$ (by Assumption A.4). Then, for an interval $t \in (t_6, t_6 + o_7)$, with $o_7 > 0$ both $R(t)$ and $P(t)$ increase:

$$ \dot{R}(t) = V_D D(t) - \gamma_R R(t) > V_D \frac{V_R}{V_D} - \gamma_R R_{int} > 0 $$

$$ \dot{P}(t) = V_D D(t) - V_R \approx V_D D(t_6) - V_R > 0. $$

Now, we want to find a duration $d_7$ such that $P(t)$ crosses its threshold $P_{low}$ before $B(t)$ decreases to $B_{high}$. To achieve this, set $P(t_6 + d_7) = P_{low}$ and $B(t_6 + d_7) = B(t_6) - V_R d_7 > B_{high}$. Substituting for $B(t_6)$ obtains:

$$ B(t_5) + (V_R - \gamma_B \delta) d_6 - V_R d_7 > B_{high} \Leftrightarrow V_R d_7 < B(t_5) - B_{high} + (V_R - \gamma_B \delta) d_6 $$

This last inequality holds by condition 1(i) and since $B(t_5) = B_0 > B_{high}$, ensuring that the system enters stage 8. $\square$

*Proof of Proposition 2:* Consider the system during stage 8, that is $t_7 < t < t_8$. Since $\dot{B}(t) = -V_R$ in this interval, $B(t)$ decreases according to the expression $B(t) = B(t_7) - V_R(t - t_7)$. By Assumption A.8, $R(t_8) > R_{int}$. By Assumption A.7,
$D(t_7) > \frac{V_R}{V_D}$, which implies $\dot{P}(t_7) = V_D D(t_7) - V_R > 0$ and $P(t)$ increases in an interval $[t_7, t_7 + o_8]$, with $o_7 > 0$.

By condition 2(ii), $P(t_7 + d_{8max}) < P_{high}$. To check that $B$ is the first variable to cross its threshold, compute the values of $B(t)$ at the instants $t_7+d_{8min}$ and $t_7+d_{8max}$, as given by equalities 2(i):

$$B(t_7 + d_{8min}) = B(t_7) - V_R(t_7 + d_{8min} - t_7) = B(t_7) - (B(t_7) - B_{eq}) = B_{eq},$$ $$B(t_7 + d_{8max}) = B(t_7) - V_R(t_7 + d_{8max} - t_7) = B(t_7) - B(t_7) = 0,$$

and recall that $B_{eq} > B_{high}$. So, since $B(t)$ is continuous, there exists some $d_8$ with $d_{8min} < d_8 < d_{8max}$, such that $B(t_7 + d_8) = B_{high}$. The system enters stage 1. □

*Proof of Proposition 3:* Consider the system during stage 1, that is $t_8 < t < t_1$. During this interval, $\dot{B}(t) = 0$, implying that $B(t)$ remains constant. Since $B(t)$ enters stage 1 as $B(t)$ crosses below $B_{high}$, we can write $B(t) \approx B_{high}$.

The equation for $P$ is $\dot{P}(t) = V_D D(t) > 0$ implying that $P(t)$ increases. Then, condition 3(i) implies that $P(t)$ will cross its threshold $P_{high}$ before $R(t)$ crosses $R_{int}$. The system enters stage 2. □

*Proof of Proposition 4:* In $[t_8, t_3]$, $\dot{B}(t) = 0$ implies that $B(t)$ is constant, with $B(t) \approx B_{high}$. While $B(t)$ is constant, $D(t)$ satisfies $\dot{D}(t) = V_B B_{high} - \gamma_D D(t)$ and has the form:

$$D(t) = \left( D(t_8) - \frac{V_B B_{high}}{\gamma_D} \right) e^{-\gamma_D(t-t_8)} + \frac{V_B B_{high}}{\gamma_D}. \quad (A1)$$

Since $D_{min} > \frac{V_B B_{eq}}{\gamma_D} > \frac{V_B B_{high}}{\gamma_D}$ (by Assumption A.5), $D(t)$ decreases and approaches $D_{min}$ as $t$ increases. Eventually, $R(t)$ also decreases for $t < t_3$ and approaches $R(t_3) = \frac{V_D V_B B_{high}}{\gamma_R \gamma_D}$.

In $[t_1, t_2]$, $\dot{P}(t) = V_D D(t) > 0$ implies that $P(t)$ increases and will remain above $P_{high}$. Therefore, the next variable to cross a threshold is $R(t)$, it decreases and will reach $R(t_2) = R_{int}$, so the system enters stage 3.

In $[t_2, t_3]$, $D(t_2) < \frac{V_R}{V_D}$ by Assumption A.6, and $D(t)$ is still decreasing. Now this implies $\dot{P} = V_D D - V_R < 0$, so $P$ decreases until it crosses $P_{high} = P(t_3)$, and the system enters stage 4.

Consider next the system during stage 4, that is to say $t_3 < t < t_4$. As the oscillator enters stage 4, $B(t_3) \lesssim B_{high}$, $R(t_3) \ll R_{int}$ and $P(t_3) = P_{high}$.

Since $B_{high} < B_{eq} = \frac{V_R}{\gamma_B P_{low}}$ (see Assumption A.2), in an interval $[t_3, t_3 + o_4]$, with $o_4 > 0$, $B(t)$ will increase since $\dot{B}(t) = V_R - \gamma_B B(t) P_{low} > 0$. But, notice that $B(t_3) \lesssim B_{high}$ at the end of stage 3, which implies that $B(t_3 + \varepsilon) > B_{high}$ for any $\varepsilon > 0$. Therefore, $B(t)$ will cross its threshold with $B(t_4) = B_{high}$, the system enters stage 5, and $d_4 = t_4 - t_3$ is necessarily a very short interval.
Finally, consider the system during stage 5, that is to say $t_4 < t < t_5$. For $t_8 < t < t_4$ we had $B(t) \approx B_{high}$ so, at the beginning of stage 5, B satisfies $\dot{B}(t) = V_R - \gamma_B B(t) P_{low} > 0$. $B(t)$ increases and converges quickly to $B_{eq}$. Thus, we can write $B(t) \approx B_{eq}$ during stage 5.

On some interval $[t_4, t_4+o_5)$, with $o_5 > 0$ both $D(t)$ and $R(t)$ continue to decrease, as indicated by equation (A1). Namely $R(t) < R_{int}$.

In stage 5, $P(t)$ is governed by the equation $\dot{P}(t) = V_D D(t) - \gamma_B B(t) P_{low}$. By Assumption A.6, during stage 4 and 5, $D(t) < \frac{V_R}{V_D}$. $\dot{P}(t) \approx V_D D(t) - \gamma_B B_{eq} P_{low} \approx V_D D(t) - V_R < 0$ $P$ decreases until it crosses $P_{low}$ at $P(t_5)$, defining the entry of the system into stage 6. $\square$

# Appendix B Transitions from refuted stages

In this Section, to complete our analysis, we provide a piecewise affine model for the biologically refuted stages 9-12 and draw the possible transitions. These stages are characterized by combinations of values for the complexes CLOCK:BMAL1 and PER:CRY that are contradictory to the phase opposition typically observed between those two complexes, since the variables $B$ and $P$ are both below their minimal thresholds in stages 9 and 10, and both above their maximal thresholds in stages 11 and 12. In this context, the equations governing the dynamics in stages 9-12 do not play a role in the algorithm which prevents transitions into these stages. However, it may be conceived that the system might be temporarily shifted to one of these stages, as in the case of a perturbation. To treat these cases, we next show that all trajectories starting in either of these four stages eventually leave them to enter one of the healthy cycle stages.

To define the PWA model in stages 9-12, we need the corresponding approximation of the nonlinear inactivation term $\gamma_B BP$. For low $B$ and low $P$ (stages 9 and 10), the inactivation term is quite small and we use the same approximations as in the adjacent stages with low $B$ and intermediate $P$ (stages 4 and 1, respectively). For high $B$ and $P$ (stages 11 and 12), the inactivation term should be significant and we assume that it is proportional to $P$ while $B$ is fixed at large value:

$$ \gamma_B BP \approx \begin{cases} \gamma_B BP_{low}, & R < R_{int}, B < B_{high}, P < P_{low} \\ 0, & R > R_{int}, B < B_{high}, P < P_{low} \\ \gamma_B B_{eq} P, & R < R_{int}, B > B_{high}, P > P_{high} \\ \gamma_B B_{eq} P, & R > R_{int}, B > B_{high}, P > P_{high}. \end{cases} \quad (B2) $$

We can now check that any trajectory starting in one of these four stages eventually enters a stage in the healthy circadian cycle.

**Proposition 6.** Consider the PWA system as defined in Tables 1 and B1, and satisfying Assumptions A.1 to A.4. Then, a trajectory starting in one of stages 9 - 12 eventually crosses to another stage, as follows:

* Stage 9 to 4, 6 or 10;

* Stage 10 to 1 or 9;

* Stage 11 to 3, 5, or 12;
• Stage 12 to 2, 8, or 11.

*Proof:* In stage 9, $\dot{B} = V_R - \gamma_B B P_{low}$ and $\dot{P} = V_D D - \gamma_B B P_{low}$. From Assumption A.2, $B < B_{high} < \frac{V_R}{\gamma_B P_{low}}$, so $\dot{B} > 0$ and B will increase until it crosses its threshold. Since no specific information is available on D or R, we assume that, depending on their initial conditions, trajectories starting in 9 may also cross either the $P_{low}$ or $R_{int}$ thresholds first. That is from 9 to 4, 6, or 10.

In stage 10, $\dot{B} = 0$ and $\dot{P} = V_D D$ so B remains constant while P increases and trajectories will eventually cross the $P_{low}$ threshold. Depending on the B value, some trajectories may cross $R_{int}$ first. From 10 to 1 or 9.

We analyse stages 11 and 12 as a single state: the equations for P, D and R are the same, while that for B is $\dot{B} = V_R - \gamma_B B_{eq} P$ in stage 11 and $\dot{B} = -\gamma_B B_{eq} P$ in stage 12. Since $P_{high} > \frac{V_R}{\gamma_B B_{eq}} = P_{low}$ (by definition of $B_{eq}$), it follows that $\dot{B} < V_R - \gamma_B B_{eq} P_{high}$ strictly decreases on stages 11 and 12. However, we do not have enough knowledge of D and R to establish a specific order in crossing for all trajectories. Since trajectories may depend on initial condition, we consider that either of the three thresholds can be crossed from 11 and 12, enabling three transitions from each of these stages (from 11 to 3, 5 or 12 and from 12 to 2, 8 or 11). □

Note that this set of transitions theoretically allows the existence of (artificial) cycles of the form $9 \rightleftharpoons 10$ or $11 \rightleftharpoons 12$. But we next show that such cycles are unstable, in the sense that the trajectories always leave them after finite time.

**Proposition 7.** Consider the PWA system as defined in Tables 1 and B1, and satisfying Assumptions A.1 to A.4. Any cycles of the form $9 \rightleftharpoons 10$ or $11 \rightleftharpoons 12$ are unstable, that is, any trajectory leaves in finite time.

*Proof:* An infinite cycle of the form $9 \rightleftharpoons 10$ may happen, if $R_{int}$ is always the first threshold to be crossed for all times. On crossing between these two stages, the equations of R and D remain unchanged. Hence, by continuity of solutions, the derivative of R doesn't immediately change sign upon switching, and we may assume that $R(t)$ continues increasing (for transition $9 \to 10$) or decreasing (for transition $10 \to 9$) for some time $\epsilon$ after the switch. Suppose that there exists a sequence of $R_{int}$ crossings between 9 and 10, starting by the transition $9 \to 10$: $t_{9,0} < t_{10,1} < t_{9,2} < t_{10,3} < \dots < t_{i,j} < \dots$, tending to infinity, and where $t_{i,j}$ denotes the exit time of stage i at the $j^{th}$ transition between stages 9 and 10 ($i = 9, 10$ and $j$ tends to infinity). To get a contradiction and show that trajectories must leave such a cycle after a finite time interval, note that B stays constant in 10 but increases and tends to $\frac{V_R}{\gamma_B P_{low}}$ in 9 according to:

$$B(t) = e^{-\gamma_B P_{low}(t-t_0)}(B(t_0) - \frac{V_R}{\gamma_B P_{low}}) + \frac{V_R}{\gamma_B P_{low}}$$

with $B_{high} < \frac{V_R}{\gamma_B P_{low}}$ (from Assumption A.2).
We have that: $B(t_{9,0}) = B(t_{10,1}) < B(t_{9,2}) = B(t_{10,3}) < \cdots < B(t_{9,j}) = B(t_{10,j+1}) < B_{high} < B(t_{9,j+2}) < \frac{V_R}{\gamma_B P_{low}}$. Therefore, there exists a finite instant $t_{9,*}$ such that: $B(t_{9,*}) > B_{high}$.

At $t_{9,*}$, $B_{high}$ is the first threshold to be crossed, and hence the trajectory crosses from 9 to 6, breaking the artificial cycle and entering the healthy cycle. Notice that since in stage 10 $P$ increases ($\dot{P} = V_D D > 0$), it may be possible that $P$ increases above its threshold $P_{low}$ before $B$ crosses $B_{high}$ but, since we don't know the tendency of $P$ in stage 9, it is easier to guarantee the exit from stages 9-10 with the crossing of $B_{high}$ by $B$.

Similarly, suppose that a trajectory follows a cycle of the form 11 $\rightleftharpoons$ 12, with a sequence of switches at times $t_{11,0} < t_{12,1} < t_{11,2} < t_{12,3} < \cdots < t_{i,j} < \cdots$, with $i = 11, 12$ and $j \to \infty$. Recall that $B > B_{high}$ and $P > P_{high}$ in stages 11 and 12, recall also $B_{eq} = \frac{V_R}{\gamma_B P_{low}}$ and note that $B$ equation satisfies $\dot{B} = V_R - \gamma_B B_{eq} P < V_R - \gamma_B B_{eq} P_{high} < 0$ in stage 11 and $\dot{B} = -\gamma_B B_{eq} P < -\gamma_B B_{eq} P_{high} < V_R - \gamma_B B_{eq} P_{high} < 0$ in stage 12.

Therefore, $B$ decreases along both stages 11 and 12, with $B(t) < B(t_0) + V_R - \gamma_B B_{eq} P_{high}t$. Thus, starting from the transition 11 $\to$ 12 we have that: $B(t_{11,0}) > B(t_{12,1}) > B(t_{11,3}) > \cdots > B_{high} > B(t_{i,j})$. There is a finite instant $t_{*,*} = \frac{B_{high} - B(t_0)}{V_R - \gamma_B B_{eq} P_{high}}$ at which the trajectory will cross to either stage 2 or 3, breaking the 11 $\rightleftharpoons$ 12 cycle and entering the healthy cycle. Similarly to the first part of the proof, notice that determining the tendency of $P$ in stages 11 and 12 is difficult, and $P$ may cross $P_{high}$ first but to guarantee the exit from stages 11-12, it is easier to consider the crossing of $B_{high}$ by $B$. $\square$

The state transition graph corresponding to the PWA system in Table 1 together with the equations for stages 9-12 in Table B1, is shown in Figure B2. The healthy cycle is represented by green arrows and each transition from the next stage is ensured by Assumptions A.1-A.8, by conditions from Propositions 1-4 and provided that the cycle started at the point $(B_0, D_0, R_0, P_0)$ at the boundary between stages 5 and 6 (as defined by Eq. 9). The arrows in black are the transitions eliminated by Propositions 1-4. The gray arrows represent transitions from the biologically refuted stages and assume no hypothesis on the initial conditions of trajectories. Therefore, a trajectory evolving from one of the refuted stages 9-12 into one of the healthy cycle may not immediately follow the healthy cycle sequence of transitions. However, simulations show that trajectories do follow this sequence after a finite time.

## References

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

[2] Almeida, S., Chaves, M., Delaunay, F.: Transcription-based circadian mechanism controls the duration of molecular clock states in response to signaling inputs.
<table><caption>**Table B1** Equations describing the dynamics of the circadian clock in stages 9-12.</caption><thead><tr><th>Stage</th><th>9 (stage refuted)</th><th>10 (stage refuted)</th><th>11 (stage refuted)</th><th>12 (stage refuted)</th></tr></thead><tbody><tr><td rowspan="3">Stage conditions</td><td>R &lt; R<sub>int</sub></td><td>R &gt; R<sub>int</sub></td><td>R &lt; R<sub>int</sub></td><td>R &gt; R<sub>int</sub></td></tr><tr><td>B &lt; B<sub>high</sub></td><td>B &lt; B<sub>high</sub></td><td>B &gt; B<sub>high</sub></td><td>B &gt; B<sub>high</sub></td></tr><tr><td>P &lt; P<sub>low</sub></td><td>P &lt; P<sub>low</sub></td><td>P &gt; P<sub>high</sub></td><td>P &gt; P<sub>high</sub></td></tr><tr><td>CLOCK:BMAL1</td><td>↑</td><td>→</td><td>↓</td><td>↓</td></tr><tr><td>˙B =</td><td>V<sub>R</sub> - γ<sub>B</sub>BP<sub>low</sub> &gt; 0</td><td>0</td><td>V<sub>R</sub> - γ<sub>B</sub>B<sub>eq</sub>P &lt; 0</td><td>-γ<sub>B</sub>B<sub>eq</sub>P &lt; 0</td></tr><tr><td>˙D =</td><td>V<sub>B</sub>B - γ<sub>D</sub>D</td><td>V<sub>B</sub>B - γ<sub>D</sub>D</td><td>V<sub>B</sub>B - γ<sub>D</sub>D</td><td>V<sub>B</sub>B - γ<sub>D</sub>D</td></tr><tr><td>˙R =</td><td>V<sub>D</sub>D - γ<sub>R</sub>R</td><td>V<sub>D</sub>D - γ<sub>R</sub>R</td><td>V<sub>D</sub>D - γ<sub>R</sub>R</td><td>V<sub>D</sub>D - γ<sub>R</sub>R</td></tr><tr><td>PER:CRY</td><td>↑</td><td>↑</td><td>↓</td><td>↓</td></tr><tr><td>˙P =</td><td>V<sub>D</sub>D - γ<sub>B</sub>BP<sub>low</sub></td><td>V<sub>D</sub>D &gt; 0</td><td>V<sub>D</sub>D - γ<sub>B</sub>B<sub>eq</sub>P</td><td>V<sub>D</sub>D - γ<sub>B</sub>B<sub>eq</sub>P</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Fig. B2** Considering that the starting point of the cycle is at the boundary between stages 5 and stage 6 (in red): green arrows show the healthy cycle, black arrows show the unexpected transitions from the healthy cycle but which are prevented by conditions from Propositions 1 to 4. Otherwise, other possible transitions from the healthy cycle are excluded by definition of the model equations and by Assumptions A.1-A.8. Considering that the starting point is in a biologically refuted stage (stages 9 to 12): grey arrows show transitions from these stages.</figcaption></figure>

Journal of theoretical biology **484**, 110015 (2020)

[3] Forger, D.B., Peskin, C.S.: A detailed predictive model of the mammalian circadian clock. Proceedings of the National Academy of Sciences **100**(25), 14806-14811 (2003)

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

[5] Woller, A., Duez, H., Staels, B., Lefranc, M.: A mathematical model of the liver circadian clock linking feeding and fasting cycles to clock function. Cell reports **17**(4), 1087-1097 (2016)

[6] Korenčič, A., Bordyugov, G., Košir, R., Rozman, D., Goličnik, M., Herzel, H.: The interplay of cis-regulatory elements rules circadian rhythms in mouse liver. PLoS one **7**(11), 46835 (2012)
[7] Hesse, J., Martinelli, J., Aboumanify, O., Ballesta, A., Relógio, A.: 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 (2021)

[8] Woller, A., Gonze, D.: Modeling clock-related metabolic syndrome due to conflicting light and food cues. Scientific reports **8**(1), 1–10 (2018)

[9] Brown, L.S., Doyle III, F.J.: A dual-feedback loop model of the mammalian circadian clock for multi-input control of circadian phase. PLoS Computational Biology **16**(11), 1008459 (2020)

[10] Gonze, D., Halloy, J., Goldbeter, A.: Robustness of circadian rhythms with respect to molecular noise. Proceedings of the National Academy of Sciences **99**(2), 673–678 (2002)

[11] Glass, L., Kauffman, S.A.: The logical analysis of continuous, non-linear biochemical control networks. Journal of theoretical Biology **39**(1), 103–129 (1973)

[12] Glass, L.: Combinatorial and topological methods in nonlinear chemical kinetics. The Journal of chemical physics **63**(4), 1325–1335 (1975)

[13] Glass, L., Pasternack, J.S.: Stable oscillations in mathematical models of biological control systems. Journal of Mathematical Biology **6**, 207–223 (1978)

[14] Glass, L., Pasternack, J.S.: Prediction of limit cycles in mathematical models of biological oscillations. Bulletin of mathematical biology **40**(1), 27–44 (1978)

[15] De Jong, H., Gouzé, J.-L., Hernandez, C., Page, M., Sari, T., Geiselmann, J.: Qualitative simulation of genetic regulatory networks using piecewise-linear models. Bulletin of mathematical biology **66**(2), 301–340 (2004)

[16] Edwards, R.: Analysis of continuous-time switching networks. Physica D: Nonlinear Phenomena **146**(1-4), 165–199 (2000)

[17] Casey, R., Jong, H.d., Gouzé, J.-L.: Piecewise-linear models of genetic regulatory networks: equilibria and their stability. Journal of mathematical biology **52**(1), 27–56 (2006)

[18] Farcot, E., Gouzé, J.-L.: Periodic solutions of piecewise affine gene network models with non uniform decay rates: the case of a negative feedback loop. Acta biotheoretica **57**(4), 429–455 (2009)

[19] Mestl, T., Plahte, E., Omholt, S.W.: Periodic solutions in systems of piecewise-linear differential equations. Dynamics and stability of systems **10**(2), 179–193 (1995)

[20] Farcot, E.: Geometric properties of a class of piecewise affine biological network
models. Journal of Mathematical Biology **52**, 373-418 (2006)

[21] Edwards, R., Farcot, E., Foxall, E.: Explicit construction of chaotic attractors in glass networks. Chaos, Solitons & Fractals **45**(5), 666-680 (2012)

[22] Cummins, B., Gedeon, T., Harker, S., Mischaikow, K., Mok, K.: Combinatorial representation of parameter space for switching networks. SIAM journal on applied dynamical systems **15**(4), 2176-2212 (2016)

[23] Gameiro, M.: Dynamic signatures generated by regulatory networks (2018). https://github.com/marciogameiro/DSGRN

[24] Cummins, B., Gedeon, T., Harker, S., Mischaikow, K.: Model rejection and parameter reduction via time series. SIAM journal on applied dynamical systems **17**(2), 1589-1616 (2018)

[25] Fox, J., Cummins, B., Moseley, R.C., Gameiro, M., Haase, S.B.: A yeast cell cycle pulse generator model shows consistency with multiple oscillatory and checkpoint mutant datasets. Mathematical Biosciences **367**, 109102 (2024)

[26] Almeida, S., Chaves, M., Delaunay, F.: Control of synchronization ratios in clock/cell cycle coupling by growth factors and glucocorticoids. Royal Society Open Science **7**(2), 192054 (2020)

[27] Berry, E., Cummins, B., Nerem, R.R., Smith, L.M., Haase, S.B., Gedeon, T.: Using extremal events to characterize noisy time series. Journal of mathematical biology **80**(5), 1523-1557 (2020)

[28] Thomas, R.: Boolean formalization of genetic control circuits. Journal of theoretical biology **42**(3), 563-585 (1973)

[29] Behaegel, J., Comet, J.-P., Bernot, G., Cornillon, E., Delaunay, F.: A hybrid model of cell cycle in mammals. Journal of bioinformatics and computational biology **14**(01), 1640001 (2016)

[30] Feillet, C., Krusche, P., Tamanini, F., Janssens, R.C., Downey, M.J., Martin, P., Teboul, M., Saito, S., Lévi, F.A., Bretschneider, T., et al.: Phase locking and multiple oscillating attractors for the coupled mammalian clock and cell cycle. Proceedings of the National Academy of Sciences **111**(27), 9828-9833 (2014)

[31] Takahashi, J.S.: Transcriptional architecture of the mammalian circadian clock. Nature Reviews Genetics **18**(3), 164-179 (2017)

[32] Burckard, O., Teboul, M., Delaunay, F., Chaves, M.: Cycle dynamics and synchronization in a coupled network of peripheral circadian clocks. Interface Focus **12**(3), 20210087 (2022)

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

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

[35] 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.*: Regulation of circadian behaviour and metabolism by rev-erb-α and rev-erb-β. *Nature* **485**(7396), 123–127 (2012)

[36] François, P.: A model for the neurospora circadian clock. *Biophysical journal* **88**(4), 2369–2383 (2005)