# Circadian Clock Model with Sequestration Repression Motif: Existence of Periodic Orbits and Entrainment Properties

Benjamin Böbel, Madalena Chaves, Jean-Luc Gouzé

► To cite this version:

Benjamin Böbel, Madalena Chaves, Jean-Luc Gouzé. Circadian Clock Model with Sequestration Repression Motif: Existence of Periodic Orbits and Entrainment Properties. Interface Focus, 2025, 15 (3), 10.1098/rsfs.2025.0008. hal-05060049

HAL Id: hal-05060049

https://hal.science/hal-05060049v1

Submitted on 27 May 2025

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
# Circadian Clock Model with Sequestration Repression Motif:
Existence of Periodic Orbits and Entrainment Properties

Böbel, Benjamin ¹, Chaves, Madalena ¹, and Gouzé, Jean-Luc ¹

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

## Abstract

Protein sequestration motifs appear in many biological regulatory networks and introduce special properties into the network dynamics. Sequestration can be described as a mode of inactivation of a given protein by its binding to a second protein to form a new complex. In this complexed form, the original protein is prevented from performing its specific functions and is thus rendered inactive. We study a mathematical model of the mammalian circadian clock with a protein sequestration motif, which generates one of the negative feedback loops that guarantees periodic behavior. First, the motif permits a time-scale separation, which can be used to simplify the model. We show that the simplified model admits a periodic orbit over a fairly large region of parameters. Second, we are able to show that the sequestration motif induces a phase response curve with a very specific form, implying that external perturbations can affect the system only in a narrow window of the periodic orbit. Finally, we compare our model with a classic Goodwin oscillator, to illustrate the dynamical and robustness properties induced by the protein sequestration motif.

## 1 Introduction

Virtually all living systems display internal rhythms that are coordinated with the geophysical day-night schedule [1]. These circadian rhythms are generated on the cellular level by a complex and circular interplay of synthesis, degradation and sequestration of transcription factors, a group of regulatory proteins [2]. Cellular circadian oscillators, or clocks, have a freerunning period of about 24h and possess the ability to vary this period with respect to external input.

Each type of mammalian cell exhibits the same core genetic oscillator, while its external inputs and also downstream outputs vary with its tissue type [3]. Common inputs for circadian clocks include light information from the eye in the central clock of the brain, body temperature, insulin levels or exercise for circadian clocks in peripheral tissues.

Circadian systems in mammals are characterized

by a hierarchical structure where the suprachiasmatic nucleus (SCN) in the hypothalamus plays the central role and cells of peripheral organs are downstream [3]. In the healthy state, each organ exhibits internal coherence between the rhythms of its cells. This has been most prominently shown for the cells in the SCN. In vitro slices of SCN tissue show sustained coherent oscillations which do not dampen or dephase with time [4]. In vitro samples of other tissues do not share this behaviour as their oscillations dampen or dephase. In mammals, time-of-day information from the SCN is passed down to oscillators in peripheral tissues through a variety of pathways such as body temperature or endocrine factors [3].

In the last 30 years a variety of models have been proposed to describe the transcriptional-translational feedback loops that constitute the core of the circadian system in every cell [5, 6,
7, 8, 9]. Among them are approaches to include all known processes in order to make quantitative predictions for the medical community, as well as reduced models which take into account only the minimum amount of constituents needed to produce self-sustained oscillations.

The transcriptional-translational feedback loops in the circadian clock are based on transcription factors, i.e., proteins in the nucleus which recognize specific DNA sequences, bind to these, and influence the transcription rate of the corresponding gene [10, p. 402]. The DNA sequences which are bound are regulatory elements, or *boxes*, that lie in the promoter region of a gene. There are three boxes which are bound by transcription factors of the core circadian system, the clock-controlled elements (CCEs): E-box, D-box and R-box (RORE). In turn, all of the genes coding for the core circadian proteins have at least one of these CCEs in their promoter region [11].

For example, the complex CLOCK:BMAL1 is a transcriptional activator that binds to E-boxes and contributes to the expression of the gene *Dbp* (among others), increasing the production of the protein DBP. DBP, in turn, is an activator of the D-box, increasing the expression of proteins REV-ERB and PER. Finally, REV-ERB is a transcriptional repressor binding to the R-box, repressing the expression of CLOCK and BMAL1.

CLOCK:BMAL1 is the core protein complex of the circadian system. Its concentration peaks at dawn and it is the major activator of clock-controlled genes, both downstream genes which are rhythmically expressed and other genes of the clock system [2, 12]. In mice it binds over 2000 sites and, in particular, upregulates genes which are important for metabolism [12].

The mammalian circadian system is made up of two interlocked negative feedback loops that affect the complex CLOCK:BMAL1. The first negative feedback loop via REV represses transcription of gene *Bmal1* and eventually contributes to a decrease in the CLOCK:BMAL1 complex. This effect is often represented by a decreasing Hill-type function of REV on the synthesis of CLOCK:BMAL1.

The second feedback is one of the major processes in the circadian clock involving the tight binding of the complex PER:CRY to the transcription factor CLOCK:BMAL1. We follow [13, 14, 15, 16] in calling this process protein sequestration by complex-

ation, or simply *sequestration*. The sequestration process inactivates the transcriptional activity of CLOCK:BMAL1.

Most models have focused on the transcriptional mechanisms of the circadian clock, in which repression can be modeled with a Michaelis-Menten or Hill type function. Sequestration of CLOCK:BMAL1 is a post-translational process which can be more suitably described by mass action kinetics. The process of protein sequestration as a repression mechanism in biological systems is well-known [13, 17]. It is known to generate ultrasensitive responses in genetic networks, comparable to strongly cooperative processes (high Hill-coefficients) [18, 13]. In [16] the sequestration motif is used as the key process in the design of a biological feedback controller.

Sequestration type repression as a central process of the circadian system has been studied in the Kim-Forger model [19, 15]. A derivation of their model can be found in [20, pp. 46-49]. In their model the activator in the sequestration has a constant concentration which simplifies the mass-action expression into a static repression function. This static repression function has a big impact on the behaviour of a network of coupled oscillators, notably, the coupled period stays close to the freerunning period [21].

In [22], the author models the *Neurospora crassa* circadian clock, where a dynamic sequestration repression is the only negative feedback loop and finds that this repression mechanism suffices to produce oscillations in ODE models. Stable oscillations in an ODE model with only sequestration repression are also found in [16]. However, in the latter example the authors show that the mean of stochastic simulations of the same model does not produce oscillations.

In this paper, we study sequestration as a dynamical process, where both the dynamics of the activator concentration and the repressor concentration are modeled dynamically. We use a model that takes into account both major feedback loops and combine Hill type for transcriptional repression and mass action kinetics for sequestration repression. This model has been previously developed in our team and used to describe the interaction of circadian clock and cell cycle [23] and to study the effect of time-restricted feeding on clusters of coupled cells [24].

We now mathematically show the existence of oscillations using the symmetry of the sequestration
process. The sequestration process is symmetrical in the sense, that it has the same effect on the concentrations of the two involved chemical species. In our case, one PER:CRY complex inactivates exactly one CLOCK:BMAL1. Mathematically, sequestration is described by the same mass-action term in the dynamics of two variables and introduces therefore a symmetrical term in the model equations.

Furthermore, we employ the theory of phase responses to explain some entrainment properties of the oscillations. The Phase Response Curve (PRC) is a useful measure in theoretical models as well as real biological systems [20, 25, 26]. It tabulates resulting phase shifts ($\Delta\phi$) from a localized external perturbation, a pulse, as a function of perturbation timing ($\phi$) and is thus an input-output measure of the periodic steady state of a dynamical system [27]. An external perturbation changes amplitude and phase of an oscillator. For stable periodic orbits, amplitude changes are transient, however, while phase changes persist [28]. The PRC $\Delta\phi(\phi)$ therefore captures the persistent response of an oscillator, relying on the assumption that trajectories rapidly fall back to the periodic orbit.

Section 2 introduces a reduced piecewise affine model for circadian oscillations exploiting the sequestration process. Section 3 contains a proof for the existence of a periodic behaviour by a cyclic transition of subsystems. Section 4 shows the phase response-curve of the model. Finally, Section 5 studies entrainment properties that result from the phase responses.

## 2 A Reduced Circadian Clock Model

<figure><img src="image_2.png" /><figcaption>Figure 1: Scheme of the circadian clock model in (1). Complexes CLOCK:BMAL1 and PER:CRY and proteins DBP and REV-ERB interact forming two negative feedback loops.</figcaption></figure>

We study a four variable ODE model that describes the transcriptional regulation and major post-translational effects of circadian rhythms in mammalian peripheral cells. It has been proposed in [23] as a skeleton model which describes only the main transcriptional-translational feedback loops (TTFL). Our team built a larger model of the core circadian clock based on the occupation of the CCEs E-box, D-box and R-box and their role in the expression of the genes *Bmal1*, *Ror*, *Rev*, *Dbp*, *E4bp4*, *Per* and *Cry* [9]. The larger 8-variable model uses Michaelis-Menten type activation of the three CCEs which are incorporated in the ODE description of protein production. The complex CLOCK:BMAL1 was chosen to be represented by the rate-limiting BMAL1 dynamics, since CLOCK and BMAL1 share the same CCE (R-box). The complex PER:CRY was modeled through the complex formation of PER and CRY. By quasi-steady state assumptions on several protein concentrations, a 4-variable skeleton model was obtained. Because of the steady-state calculations, the CCE activation of PER percolates to the complex PER:CRY. As a consequence, the production term for the protein REV-ERB and the complex PER:CRY become identical in the skeleton model. We choose to work with the skeleton model, as it emphasizes the main dynamical elements of the system and is amenable to detailed analytical study, while retaining the biological properties. See Fig. 1 for the scheme of our model.

The model variables are concentrations of the protein complexes CLOCK:BMAL1 and PER:CRY and proteins DBP and REV. We group proteins of homologous genes into one variable such that variable CRY takes into account both CRY1 and CRY2. The model equation read as follows:

$$ 
\begin{array}{l}
\text{[CLOCK:BMAL1]} \quad \dot{B} = V_R h(R) - \gamma_B BP \\
\text{[DBP]} \quad \dot{D} = V_B B - \gamma_D D \\
\text{[REV-ERB]} \quad \dot{R} = V_D D - \gamma_R R \\
\text{[PER:CRY]} \quad \dot{P} = V_D D - \gamma_B BP
\end{array}
\qquad (1) \\
\text{with } h(R) = \frac{k_R^2}{k_R^2 + R^2}
$$

where every variable describes a concentration of a protein or protein complex and has a synthesis and a degradation process. B stands for [BMAL1:CLOCK], D for [DBP], R for [REV:ERB] and P for [PER:CRY]. The degradation parameters
are named $\gamma$ and the synthesis parameters $V$.

There are two negative feedback loops which center around CLOCK:BMAL1 and imply a delayed repression of the complex after each increase in concentration. The first negative feedback loop via REV acts on CLOCK:BMAL1 using a Hill-type repression. The second feedback loop includes a protein sequestration type repression where the complex PER:CRY binds to CLOCK:BMAL1 and renders it inactive. By law of mass action, the rate of the sequestration process is proportional to both concentrations. Note, that the synthesis terms stem from a reduction of the larger 8-variable model, as described above.

Equations for $B$, $D$, $R$ in (1) resemble a Goodwin oscillator, since they invoke the nonlinear degradation of $B$ and then a downstream cascade of linear terms. The difference in our system is the degradation part of $B$, which is not a natural degradation but the mutual sequestration of $B$ and $P$ described above.

The standard parameter values used in this article are direct descendants of the fitted 8-variable model used in [9]. All parameters of their large model were fitted to experimental timeseries of fluorescence intensity of a REV-ERB reporter. Afterwards the 8-variable model was reduced to (1) as described above.

## 2.1 Parameter Reduction

To analyse the dynamical behaviour of model (1) and, in particular, to investigate the contribution of the sequestration phenomenon, we first reduce the dimensionality of the parameter space. To this end, we introduce the change of variables: $B = V_R b$, $D = V_B V_R d$, $R = k_R r$, $P = V_R p$. This brings about the parameter-reduced system:

$$ 
\begin{aligned}
\dot{b} &= h_d(r) - \alpha bp \\
\dot{d} &= b - \beta d \\
\dot{r} &= \gamma d - \delta r \\
\dot{p} &= \epsilon d - \alpha bp
\end{aligned}
\qquad (2) \\
\text{with } h_d(r) = \frac{1}{1+r^2}
$$

where the new parameters are given by: $\alpha = V_R \gamma_B$, $\beta = \gamma_D$, $\gamma = V_B V_D V_R / k_R$, $\delta = \gamma_R$, $\epsilon = V_B V_D$.

Secondly, the sequestration process introduces a symmetrical degradation term in $b$ and $p$. This leads

to phase-opposition in the dynamics of variables $b$ and $p$, where $b \approx 0$ during the peak of $p$ and vice versa. The following change of variables

$$ x = b - p \qquad y = b + p \qquad (3) $$

with its inverse

$$ b = \frac{1}{2}(x + y) \qquad p = \frac{1}{2}(y - x) \qquad (4) $$

makes use of this symmetrical process and separates slow and fast elements. The dynamics of the transformed system read as follows:

$$ \dot{x} = h_d(r) - \epsilon d \qquad (5a) $$

$$ \dot{d} = \frac{1}{2}(x + y) - \beta d \qquad (5b) $$

$$ \dot{r} = \gamma d - \delta r \qquad (5c) $$

$$ \dot{y} = h_d(r) + \epsilon d + \frac{\alpha}{2}(x^2 - y^2) \qquad (5d) $$

Variables $x$ and $y$ no longer represent individual concentrations. As a difference of positive quantities, $x$ can take on positive and negative values.

By change of variables, the two mass-action expressions describing the sequestration are integrated into the new variable $y$. This is now the only variable containing the dynamical nonlinearity of the sequestration process.

Simulation of this system shows oscillations which are remarkably insensitive in period and amplitude to increasing the sequestration parameter $\alpha$. This stems from the fact, that the quadratic term in $\dot{y}$ dominates the other terms when the condition $y^2 = x^2$ is violated.

Note, that without synthesis processes, the difference between the two concentrations involved in the sequestration (called $x$) remains constant, as previously used in [16]. In the article, the authors also remark, that this reaction is on a faster timescale than the rest, given a suitably high sequestration parameter. This separation of timescales is central to the following model reduction.

## 2.2 The large sequestration rate limit

As described in Section 2.1, the sequestration term dominates other terms in (5d) when the timescales are appropriately separated. This is a consequence of the expression involving the difference of squares and the fact, that $x$ and $y$ take values on the same
<figure><img src="image_3.png" /><figcaption>Figure 2: Influence of sequestration parameter $\alpha$ on system dynamics and period. Top: Solutions of variables $b$ and $p$ for two different values of parameter $\alpha$, after transients. Bottom: Period of periodic solution of model (2) as a function of parameter $\alpha$. The standard value of $\alpha = 114.6$, determined by fitting to experimental data, is shown in grey. Other parameters: $\beta = 0.156$, $\gamma = 0.15$, $\delta = 0.241$, $\epsilon = 2.698$.</figcaption></figure>

order of magnitude as the other variables. For the small time-window when $x \approx 0 \Leftrightarrow b \approx p$, the time-scale separation does not hold for $\alpha$ below some threshold, $\alpha_{\text{thres}}$. As long as the sequestration parameter $\alpha$ is sufficiently large, timescales remain separated and (5d) will be controlled by the sequestration term. Conversely, once timescales are separated, increasing $\alpha$ above $\alpha_{\text{thres}}$ no longer changes the dynamics. To visualize this, we show two time-series and the period of the limit cycle solution depending on $\alpha$ (see Fig. 2).

For convenience, in the following we will consider the limit $\alpha \to \infty$ as a proxy for sufficiently separated timescales. However, the separation is inherent in the model dynamics and does not stem from the choice of a limit sequestration rate.

In the limit $\alpha \to \infty$ the sequestration process is infinitely fast, so $h(r)/\alpha \ll 1$ and $(\epsilon d)/\alpha \ll 1$. We can use a quasi steady-state argument to pose that the dynamics of $y$ are dictated by: $|y| = |x|$ (see [29]). Since $y$ is the sum of positive quantities, we can write: $y = |x|$.

This limit effectively means that there can never be a coexistence of both species $b$ and $p$, since they sequester with an infinite rate. One of the concentra-

tions has to be zero and any production of the zero-concentration species is immediately subtracted from the other one. Therefore, in the large sequestration rate limit, the composed variable $x = b - p$ describes the interlocked dynamics of $b$ and $p$ completely, with the two cases ($p = 0, b > 0 \Rightarrow x > 0$ and $b = 0, p > 0 \Rightarrow x < 0$).

Because of (4) and $y = |x|$, it is easy to see that:

$$b = \frac{1}{2}(x + y) = \frac{1}{2}(x + |x|) = \begin{cases} x, & x > 0 \\ 0, & \text{else} \end{cases} \quad := f_b(x). \tag{6}$$

Function $f_b(x)$ is used as the synthesis term for variable $d$ from now on.

## 2.3 Piecewise affine model

In order to further simplify the system, we replace the nonlinearly decreasing Hill-function $h_d(r)$ with a switch function $h_s(r)$ that is constant 1 below its threshold and constant 0 above it. The switch function $h_s(r)$ is the limit $n \to \infty$ of decreasing Hill-functions. With these two approximations, our model reduces to the following piecewise affine (PWA) system with linear time-invariant subsystems:

$$\begin{align} \dot{x} &= h_s(r) - \epsilon d \nonumber \\ \dot{d} &= f_b(x) - \beta d \nonumber \\ \dot{r} &= \gamma d - \delta r \tag{7} \end{align}$$

$$\text{with } h_s(r) = \begin{cases} 1, & r < 1 \\ 0, & \text{else} \end{cases},$$

where $f_b(x)$ as defined in (6).

The advantages of this model are that 4 regions can be clearly identified and an explicit solution for the dynamics in each region can be written. From these solutions, we can infer possible transitions between the regions and eventually detect global solutions such as periodic orbits.

A sample trajectory of the PWA model is shown in Fig. 3. Note how the dynamics of $x$ change abruptly at its trough due to the switch function $h_s(r)$.

# 3 Existence of a Periodic Orbit

It is in general hard to explicitly show the existence of a periodic orbit in systems of nonlinear ODEs
<figure><img src="image_4.png" /><figcaption>Figure 3: Timeseries of the PWA model. Parameters: $\beta = 0.156, \gamma = 0.15, \delta = 0.241, \epsilon = 2.698$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Scheme of the transitions for the piecewise-linear model. Possible transitions are shown as arrows between subsystems. Unwanted transitions are dashed. Transitions with a question mark will be proved in order to follow sequence (8). For clarity, not all possible transitions are shown.</figcaption></figure>

of dimension larger than 2 [30, p. 43]. To study the oscillatory behavior of circadian systems with two feedback loops, we resort to the PWA system (7). The PWA system captures both feedback loops while still being analytically tractable. In the following we show that, under certain conditions, a specific cyclic sequence of transitions between subsystems is the only possible solution of the system dynamics. The existence of a periodic orbit follows via the Brouwer fixed point theorem.

First, we define three thresholds $r^*, x^*, d^*$ to naturally divide the state space into 8 regions. In each region, the dynamics is given by one of four affine subsystems. The regions follow a notation based on whether the variables are above or below their respective thresholds, e.g., the region rXd is the subsystem where $r < r^*, x > x^*$ and $d < d^*$.

The threshold in $x$ is $x^* = 0$, the non-differentiable point in $f_b(x)$ (see Eq. (6)). The

threshold in $r$ is $r^* = 1$, which is the discontinuity of the switch function $h_s(r)$. Whether $r$ can transition its threshold is defined by $\dot{r}(r = r^*) = \gamma d - \delta$. The value of variable $d$ determines whether $\dot{r}$ on the plane $r = r^*$ is bigger or smaller than 0. Therefore the value of $d$ completely determines whether the variable $r$ can transition. To define a threshold in variable $d$, we introduce $d^* = \delta/\gamma$ such that $d < d^* \Rightarrow \dot{r}(r = 1) < 0$ and likewise $d > d^* \Rightarrow \dot{r}(r = 1) > 0$.

In the following, we give sufficient conditions for a specific sequence of subsystems. Since this sequence is cyclic, solutions must exhibit oscillations.

**Assumption 3.1.** Assume $\gamma > \epsilon\delta$, or equivalently $1 > \epsilon d^*$.

This assumption constrains the dynamics of $x$ with respect to $d$ (whether $d$ is above or below its threshold). Under this assumption $d < d^* \Rightarrow \dot{x} > 0$.

Under Assumption 3.1 it is easy to verify that the following transitions are unambiguous such that any trajectory in the first region will necessarily end up in the second region: $RxD \to Rxd, Rxd \to rxd, rxd \to rXd, rXd \to rXD$. Note that $1/\epsilon > d^* = \delta/\gamma$.

There are several possible cyclic sequences under Assumption 3.1 (see Fig. 4). We decided to investigate the sequence

$$RxD \to Rxd \to rxd \to rXd \to rXD \to RXD \to RxD \quad (8)$$

because of a numerical analysis, showing that this sequence is the most common one for random parameters. In other words, we find a large volume in parameter space that produces periodic orbits which follow this sequence. In order to prove that trajectories follow sequence (8), we will provide conditions and proofs to avoid the dashed transitions in Fig. 4 and keep instead transitions $rXD \to RXD$ and $RXD \to RxD$ (labeled by "?" in Fig. 4).

## 3.1 Transition rXD to RXD

In the region rXD, the dynamics are $\dot{x} = 1 - \epsilon d$, $\dot{d} = x - \beta d$, $\dot{r} = \gamma d - \delta r$, which allows for three transitions: $rXD \to rxD, rXD \to RXD, rXD \to rXd$. We want to constrain the parameters in such a way, that only the transition $rXD \to RXD$ remains possible. Preliminaries, such as the upper and lower bounds on the initial conditions of $x$ are given in Appendix A.
**Lemma 3.2.** Suppose stage rXD was reached from stage rXd. The time it takes for variable r to reach its threshold 1 is bounded from above by

$$T_r = \sqrt{\frac{2}{\gamma(\underline{x} - \beta/\epsilon)}},$$

where $\underline{x}$ is a lower bound on $x$ at the transition rXD $\rightarrow$ rXD (see Lemma A.2).

*Proof.* The trajectory has crossed the threshold $d = d^*$ from below. Therefore, condition $\dot{d} > 0$ holds at the transition, which implies $x > \beta d^*$. $\dot{x} > 1 - \epsilon d^* > 0$ holds by Assumption 3.1. Furthermore, $\dot{r} > \gamma d^* - \delta r = \gamma \delta / \gamma - \delta r > 0$ since $r < 1$. Consider

$$ 
\begin{aligned}
\dot{x} \ge 0 &\Rightarrow x(t) \ge \underline{x} \\
\dot{d} \ge \underline{x} - \beta \sup d &= \underline{x} - \beta/\epsilon \\
&\Rightarrow d(t) \ge (\underline{x} - \beta/\epsilon)t + d^* \\
\dot{r} \ge \gamma d - \delta \sup r &\ge \gamma d - \delta \\
&\Rightarrow r(t) \ge \frac{\gamma}{2}(\underline{x} - \beta/\epsilon)t^2 + \underbrace{(\gamma d^* - \delta)}_{0} t + r_0,
\end{aligned}
$$

where $\underline{x}$ is a lower bound on the IC of $x$ (see Lemma A.2), $r = 1$ is $r$'s threshold, $r_0 = 0$ is the lower bound on $r$'s initial condition, and $d^*$ is the IC of $d$. Denote $\tau_r$ the time at which $r = 1$. Time $\tau_r$ is then bounded from above by $T_r$:

$$ 
\begin{aligned}
r(\tau_r) &= 1 \ge \frac{\gamma}{2}(\underline{x} - \beta/\epsilon)\tau_r^2 \\
\tau_r &\le \sqrt{\frac{2}{\gamma(\underline{x} - \beta/\epsilon)}} =: T_r.
\end{aligned}
$$

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

**Lemma 3.3** (threshold conditions rXD). Suppose stage rXD was reached from stage rXd. The condition $d = 1/\epsilon$ must be fulfilled before either $d$ or $x$ reach their respective threshold.

*Proof.* The trajectory has crossed the threshold $d = d^*$ from below. Since stages rXd and rXD are governed by the same equations, condition $\dot{d} > 0$ holds on both sides of the transition, which implies $x > \beta\delta/\gamma$. We also know, that $\dot{x} > 1 - \epsilon d^* > 0$. In summary, $\dot{d} > 0$, $\dot{x} > 0$ at the time of threshold crossing. This holds true at least until $d = 1/\epsilon$ at which point $\dot{x} = 0$ and a subsequent decrease in $x$ is followed by decreasing $d$. □

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

**Lemma 3.4.** Suppose stage rXD was reached from stage rXd. Suppose we stay in stage rXD. The time it takes for variable $d$ to reach $1/\epsilon$ is bounded from below by

$$T_d = \frac{\sqrt{G^2 + 2H^2/\epsilon} - G}{H},$$

where $G = \bar{x} - \beta d^*$, $H = 1 - \epsilon d^*$ and $\bar{x}$ is an upper bound on $x$ at the transition rXd rightarrow rXD (see Lemma A.1).

*Proof.*

$$ 
\begin{aligned}
\dot{x} &\le 1 - \epsilon \inf d = 1 - \epsilon d^* \Rightarrow x(t) \le (1 - \epsilon d^*)t + \bar{x} \\
\dot{d} &\le (1 - \epsilon d^*)t + \bar{x} - \beta \inf d = (1 - \epsilon d^*)t + \bar{x} - \beta d^* \\
&\Rightarrow d(t) \le \frac{1}{2}(1 - \epsilon d^*)t^2 + (\bar{x} - \beta d^*)t + d^*,
\end{aligned}
$$

where $\bar{x}$ of Lemma A.1 was used as upper bound. We get an expression for the time $\tau_d$ at which $d$ reaches $1/\epsilon$:

$$d(\tau_d) = 1/\epsilon \le \frac{1}{2}(1 - \epsilon d^*)\tau_d^2 + (\bar{x} - \beta d^*)\tau_D + d^*$$

which is a quadratic inequality in $\tau_d$.

$$0 \le \frac{1}{2} \underbrace{(1 - \epsilon d^*)}_{H} \tau_d^2 + \underbrace{(\bar{x} - \beta d^*)}_{G} \tau_d - \underbrace{(1/\epsilon - d^*)}_{H/\epsilon}$$

with $H > 0$ because of Assumption 3.1 and $G > 0$ from the transition condition on $\bar{x}$ above. Since $\tau_d \ge 0$, there is only one solution:

$$\tau_d \ge \frac{\sqrt{G^2 + 2H^2/\epsilon} - G}{H} =: T_d.$$

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

**Theorem 3.5.** Suppose the stage rXD was reached by stage rXd. If the following inequality is fulfilled, the system can only transition to stage RXD:

$$\sqrt{\frac{2}{\gamma(\underline{x} - \beta/\epsilon)}} < \frac{\sqrt{G^2 + 2H^2/\epsilon} - G}{H}.$$

*Proof.* By Lemma 3.2, the time it takes to reach $r$'s threshold is bounded from above by

$$T_r = \sqrt{\frac{2}{\gamma(\underline{x} - \beta/\epsilon)}}.$$
<figure><img src="image_9.png" /><figcaption>Figure 5: Sample trajectory and Jordan curve in region RXD. This is a projection of phase space onto the (x, d)-plane. A sample trajectory following the dynamics of subsytem RXD is shown in solid. In dashed, we show the Jordan curve described in Lemma 3.6. The diagonal line separates the two subregions A and B of Lemma 3.6.</figcaption></figure>

Lemma 3.3 shows that the threshold $d = 1/\epsilon$ must be crossed before the system reaches either $d$'s or $x$'s threshold. By Lemma 3.4, the time it takes to reach $d = 1/\epsilon$ is bounded from below by

$$T_d = \frac{\sqrt{G^2 + 2H^2/\epsilon} - G}{H}.$$

If $T_r < T_d$, the only possible transition is rXD → RxD. □

## 3.2 Transition RXD to RxD

In the region RXD, the dynamics are $\dot{x} = -\epsilon d < 0$, $\dot{d} = x - \beta d$, $\dot{r} = \gamma d - \delta r$, which allows for two transitions: RXD → RXD, RXD → RxD. We want to constrain the parameters in such a way, that only the transition RXD → RxD remains possible. The variable $r$ is above its threshold and cannot transition below during this stage, since $\dot{r}(r = 1, d > d^*) > 0$. Furthermore, variable $r$ does not influence the dynamics of the other variables. Therefore, we can safely consider the dynamics of the 2d-system of $x$, and $d$. In 2d, the proof is based on the construction of a so-called Jordan curve [31, p. 376], a closed curve which separates the plane into two components. In our case one component contains the trajectory and the other component contains the unwanted transition $d < d^*$ (see Fig. 5). Since we are dealing with the region $d^* = d^*, x \ge 0$, there need only be a simple curve $\mathcal{J}$ connecting the two axes $x = 0, d = d^*$ to create a separation.

**Lemma 3.6 (Jordan curve).** In region RXD, consider the 2d-system:

$$\begin{aligned} \dot{x} &= -\epsilon d \\ \dot{d} &= x - \beta d. \end{aligned}$$

There exists a curve $\mathcal{J}$ which separates any trajectory starting at $(x > \underline{x}, d = d^*)$ and the axis $(0 \le x < \underline{x}, d = d^*)$.

*Proof.* Consider two subregions: subregion $A$ where $x > \beta d$ and subregion $B$, where $x < \beta d$. In subregion $A$, we create a bounding trajectory by using a new set of dynamics which has, for every coordinate $(x, d)$, a flatter negative slope, than the original system. It can therefore never be crossed by the original trajectory. Let

$$\begin{aligned} \dot{x}_m &= -\frac{\epsilon}{\beta}x_m \\ \dot{d}_m &= x_m - \beta d_m \end{aligned}$$

which has the following slope

$$\frac{d d_m}{d x_m} = \frac{x_m - \beta d_m}{-\frac{\epsilon}{\beta} x_m}.$$

The slope of the original system is

$$\frac{d d}{d x} = \frac{x - \beta d}{-\epsilon d}.$$

Since $-d > -\frac{\epsilon}{\beta}x$ in subregion $A$, it holds:

$$\frac{d d_m}{d x_m} > \frac{d d}{d x}.$$

In subregion $B$, we create a bounding trajectory by using another set of dynamics which has, for every coordinate $(x, d)$, a less steep positive slope than the original system. It can therefore never be crossed by the original trajectory. Let

$$\begin{aligned} \dot{x}_p &= -\epsilon d^* \\ \dot{d}_p &= x_p - \beta d_p \end{aligned} \qquad (9)$$

which has the following slope:

$$\frac{d d_p}{d x_p} = \frac{x_p - \beta d_p}{-\epsilon d^*}.$$

Since $-d < -d^*$, it holds:

$$\frac{d d_p}{d x_p} < \frac{d d}{d x}.$$
Let $\mathcal{J}$ be the curve that is described by the trajectory starting at IC $(\underline{x}, d^*)$, and following dynamics of subregion $\mathcal{A}$ until $x = \beta d$ at time $t_1$ and following dynamics of subregion $\mathcal{B}$ until $x = 0$ at time $t_2$. By the arguments above $\mathcal{J}$ separates the phase space into two disconnected components, one containing all trajectories of the real system and the other one containing $(0 \le \underline{x} < \underline{x}, d = d^*)$. $\square$

**Theorem 3.7 (final condition).** Under the following conditions, the bounding Jordan curve $\mathcal{J}$ will terminate at a point $O = (x = 0, d > d^*)$:

$$x_m(t_1) \left( \frac{1}{\beta} - \frac{1}{2\epsilon d^*} x_m(t_1) \right) > d^* \quad (10)$$

where

$$x_m(t_1) = \underline{x} \exp \left( -\frac{\epsilon}{\beta^2 - \epsilon} \ln \left( \beta \frac{d^* + F}{\underline{x} + \beta F} \right) \right) \\ F = \frac{\beta \underline{x}}{\epsilon - \beta^2}.$$

*Proof.* The curve $\mathcal{J}$ is described in Theorem 3.6. It is produced by the trajectory of the system:

$$\left\{ \begin{array}{ll} (x_m, d_m) & \text{if } (x, d) \in \mathcal{A} \\ (x_p, d_p) & \text{if } (x, d) \in \mathcal{B} \end{array} \right.$$

with

$$\begin{aligned} \dot{x}_m &= -\frac{\epsilon}{\beta}x_m \\ \dot{d}_m &= x_m - \beta d_m \\ \dot{x}_p &= -\epsilon d^* \\ \dot{d}_p &= x_p - \beta d_p . \end{aligned}$$

First, We will integrate the trajectory starting from $(\underline{x}, d^*)$ and reaching $x = \beta d$.

$$\begin{aligned} x_m(t) &= \underline{x} \exp \left( -\frac{\epsilon}{\beta} t \right) \\ d_m(t) &= -F \exp \left( -\frac{\epsilon}{\beta} t \right) + (d^* + F) \exp(-\beta t) \end{aligned}$$

$$\text{where } F = \frac{\beta \underline{x}}{\epsilon - \beta^2}.$$

Using the ending condition $x(t_1) = \beta d(t_1)$:

$$\begin{aligned} \left(F + \frac{\underline{x}}{\beta}\right) \exp\left(-\frac{\epsilon}{\beta}t_1\right) &= (d^* + F) \exp(-\beta t_1) \\ t_1 &= \frac{1}{\beta - \epsilon/\beta} \ln\left(\frac{d^* + F}{\underline{x}/\beta + F}\right). \end{aligned}$$

Second, we will integrate the trajectory starting from $(x_p(0) = x_m(t_1), d_p(0) = d_m(t_1))$ and reaching $x_p(t_2) = 0$:

$$\begin{aligned} x_p(t) &= x_m(t_1) - \epsilon d^* t \\ d_p(t) &= b_0 + b_1 t + b_2 \exp(-\beta t) \end{aligned} \quad (11)$$

$$\begin{aligned} \text{where } b_0 &= \frac{x_m(t_1)}{\beta} + \frac{\epsilon d^*}{\beta^2} \\ b_1 &= \frac{\epsilon d^*}{\beta} \\ b_2 &= -\frac{\epsilon d^*}{\beta^2}. \end{aligned}$$

The final condition yields:

$$\begin{aligned} x_p(t_2) &= x_m(t_1) - \epsilon d^* t_2 = 0 \\ \Rightarrow t_2 &= \frac{x_m(t_1)}{\epsilon d^*}. \end{aligned} \quad (12)$$

Last, we need to show that $d_p(t_2) > d^*$. Keeping in mind that $\exp(-x) < (1 - x + x^2/2) \quad \forall x > 0$ and $b_2 < 0$, we can calculate a lower bound:

$$\begin{aligned} d_p(t_2) &= b_0 + b_1 t_2 + b_2 \exp(-\beta t_2) \\ &> b_0 + b_1 t_2 + b_2 \left(1 - \beta t_2 + \frac{\beta^2 t_2^2}{2}\right) \\ &= \frac{x_m(t_1)}{\beta} + t_2(b_1 - \beta b_2) + t_2^2 \left(\frac{1}{2}\beta^2 b_2\right) \\ &> \frac{x_m(t_1)}{\beta} - \frac{1}{2}\epsilon d^* t_2^2 \\ &= x_m(t_1) \left(\frac{1}{\beta} - \frac{x_m(t_1)}{2\epsilon d^*}\right) \end{aligned}$$

where

$$x_m(t_1) = \underline{x} \exp \left( -\frac{\epsilon}{\beta^2 - \epsilon} \ln \left( \beta \frac{d^* + F}{\underline{x} + \beta F} \right) \right).$$

Now, if

$$x_m(t_1) \left( \frac{1}{\beta} - \frac{x_m(t_1)}{2\epsilon d^*} \right) > d^* \quad (13)$$

it follows that $d_p(t_2) > d^*$. Therefore if a set of parameters fulfills (13), the Jordan curve will terminate at $O = (x_p = 0, d_p(t_2) > d^*)$. $\square$

**Remark.** If the Jordan curve intersects $x = 0$ with $d > d^*$, a transition of variable $d$ is impossible and only the threshold at $x = 0$ can be attained.
<figure><img src="image_10.png" /><figcaption>Figure 6: Region in parameter space admitting oscillations. Shown are 2d projections of individual parameter sets fulfilling the sufficient conditions for the given cyclic sequence of subsystems. Histograms of the number of valid values per parameter are shown on the diagonal.</figcaption></figure>

Recall the definitions of $t_2$ in Lemma 3.6 and $d_p(t_2)$ in Eq. (9).

**Corollary 3.8.** *Under Assumption 3.1, Theorem 3.5 and Theorem 3.7, the only possible transition of subsystems is the cyclic one given by: $RxD \rightarrow Rxd \rightarrow rxd \rightarrow rXd \rightarrow rXD \rightarrow RXD \rightarrow RxD$. If all solutions to Eq. (7) remain bounded (as they do in our numerical simulations), using the $(r, d \ge d_p(t_2))$-plane at the transition $RXD \rightarrow RxD$ (which does not include the point $(r^*, x^*, d^*)$) by the Brouwer fixed point theorem [32, p. 51], there exists a periodic orbit.*

By sampling the parameter space, we find a dense region which fulfills all conditions: Assumption 3.1, Theorem 3.5 and Theorem 3.7 (see Fig. 6).

In this section we showed the existence of a periodic orbit in our piecewise affine model for a fairly large parameter region. Now, that we established the existence of a periodic orbit, we will study its responses to external inputs.

<figure><img src="image_11.png" /><figcaption>Figure 7: Phase Response-Curve of the PWA System (7). The phase responses were obtained by direct numerical integration. Negative stimuli were applied on variable *x* as proxy for positive stimuli on *p*. Stimulus amplitude *A* = -0.5, stimulus length *L* = 0.05 h. The upward and downward transition times of variable *x* are indicated. Parameters as in Fig. 3.</figcaption></figure>

# 4 Dynamics under external inputs

For a deeper understanding of the system, we generated the Phase Response-Curve (PRC) of our piecewise affine model. The PRC is an input-output measure that tabulates the measured phase responses against the timing of a localized input. Phase responses can be phase advances or phase delays of the oscillator leading to a period shortening or lengthening, respectively.

We constructed the Phase Response-Curve using direct numerical integration (see Fig. 7). A pulse was applied to variable *x* on 100 different phases in the cycle. The pulse is a square pulse of length *L* = 0.05 h and amplitude *A* = -0.5. A negative pulse on *x* was used, since this is equivalent to a positive perturbation of variable *p* ($x = b - p$, see Section 2.2).

Well-known PRCs of canonical models close to bifurcations that produce periodic orbits follow one of two forms depending on their classification. Class I oscillators follow a $\Delta\phi(\phi) = 1 - \cos(\phi)$ form and Class II oscillators a $\Delta\phi(\phi) = -\sin(\phi)$ form [27]. Remarkably, the PRC of our model is close to zero almost everywhere. This is a so-called *dead-zone* [20, p. 171]. It indicates, that the oscillation is very insensitive to external perturbations during a time window which comprises the greater part of the cycle.

The stimulus phases, where the response is non-
zero, lie exactly in the region where $x > 0$ (transitions of $x$ shown in Fig. 7). To understand the form of this PRC, note that in the subsystems where $x < 0$, the variable $x$ has no influence on other variables or even itself, since $f_b(x) = 0$. Any perturbation will therefore only change the state of $x$, but will not propagate through the system. As a consequence, while $x < 0$, an external perturbation applied on variable $x$, does not produce differential responses, but a constant response and thus a flat PRC.

Variable $x$ in model (7) is the difference between the concentrations of the two complexes $b$ and $p$ that mutually inactivate each other. $x > 0$ corresponds to the interval where the complex CLOCK:BMAL1 is not repressed and can trigger the activation of all clock-controlled genes.

In terms of concentrations, the near zero responses are a result of the fact that external perturbations on [PER:CRY] do not have an influence on the concentrations [DBP] and [REV] while CLOCK:BMAL1 is repressed. This is a consequence of the fast repression by sequestration. Conversely, perturbations on [PER:CRY] have an immediate influence on [DBP] and [REV] when CLOCK:BMAL is not repressed. Every increase in [PER:CRY] is immediately followed by the mutual sequestration of PER:CRY and CLOCK:BMAL1 and thus transformed into a negative perturbation in [CLOCK:BMAL1]. It is therefore the sequestration process which interlocks the dynamics of [CLOCK:BMAL1] and [PER:CRY], that is at the heart of this lack of delays as phase response. Note that in the limiting case $\alpha \to \infty$ sequestration is instantaneous. As explained in Section 2.2, sequestration is in general faster than other processes.

## 4.1 Comparison with Goodwin oscillator

As a comparison, we also generated the PRC for a periodic system with no sequestration term, the Goodwin oscillator (see Fig. 8). We assume external perturbations to influence the protein concentration and thus variable $Z$. Model equations and parameters are taken from Fig. 2 in [33]. The equations

<figure><img src="image_12.png" /><figcaption>Figure 8: PRC of the Goodwin oscillator (Eq. (14)). A single stimulus was applied at 100 phases to the $z$ variable of the system. Stimulus amplitude $A = 0.5$, stimulus length $L = 0.05$ h. Parameters: $\alpha_1 = \alpha_2 = \alpha_3 = 5.0$, $\gamma_1 = \gamma_2 = \gamma_3 = 0.5$, $K = 1.0$, $n = 10.0$.</figcaption></figure>

are as follows:

$$ 
\begin{aligned}
\dot{X} &= \alpha_1 \frac{K^n}{K^n + Z^n} - \gamma_1 X \\
\dot{Y} &= \alpha_2 X - \gamma_2 Y \\
\dot{Z} &= \alpha_3 Y - \gamma_3 Z.
\end{aligned} 
\qquad (14) 
$$

The PRC of the Goodwin oscillator was constructed with a square pulse of length $L = 0.05$ h and amplitude $A = 0.5$. Perturbations on variable $Z$ trigger a delay in one part of the cycle and an advance in another part, like a class II oscillator. The Phase Response-Curve of our oscillator, on the other hand, shows negligible phase delays for a big part of the cycle and comparatively large phase advances in a small time window, restricted by $x$. For a more detailed interpretation, see the Discussion.

# 5 Entrainment Properties of Cellular Oscillators

Whether a periodic stimulus can entrain an oscillator depends on the difference of stimulus period $T_{st}$ and freerunning period of the oscillator $T_{fr}$ and also on the stimulus amplitude. In the range of weak perturbations, an external stimulus only triggers a phase shift in the perturbed oscillator. For this case, a stable phase relationship between an oscillator of period $T_{fr}$ and a periodic stimulus of period $T_{st}$ can only be maintained if the period difference can be compensated by the phase shift induced by the external stimulus.

Therefore, the Phase Response-Curve, which quantifies the possible phase shifts of a single stim-
<figure><img src="image_13.png" /><figcaption>Figure 9: Arnold tongue of the piecewise-affine model. A periodic square-wave stimulus with stimulus length = 0.05 h was added to the x-variable of the system. Amplitude A and stimulus period T<sub>st</sub> are varied. Color coded is the standard deviation of the Kuramoto order parameter k calculated from the phases of variable x and the stimulus. Parameters as in Fig. 3.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 10: Arnold tongue of Goodwin model. A periodic square-wave stimulus with stimulus length = 0.05h was added to the Z-variable of the system. Amplitude A and stimulus period T<sub>st</sub> are varied. Color coded is the standard deviation of the Kuramoto order parameter k calculated from the phases of variable x and the stimulus. Parameters as in Fig. 8</figcaption></figure>

ulus, is able to qualitatively and quantitatively predict the range of entrainment.

We use Arnold tongues to visualize the entrain-
ment region of an oscillator to a specific periodic
stimulus (see [34] for an introduction). With increas-
ing stimulus amplitude, a greater period difference
of freerunning and stimulus period can be compen-
sated. We quantify the entrainment of an oscillator
via the Kuramoto order paramter. The Kuramoto
order parameter k(t) is the amplitude of the com-
plex mean of all phases in a network of N coupled
oscillators:

$$
k(t) = \left| \frac{1}{N} \sum_{j}^{N} e^{i\phi_{j}(t)} \right| .
$$

k stays constant if the phase difference between the
concerned phases φⱼ does not change over time. We
use the standard deviation std(k) after a transient
as the numerical measure for entrainment. For the
PWA system, we simulated the oscillator for 25000 h
and analyzed the last 25% of the time series. We
transformed the timeseries of the stimulus and the
oscillatior into phase variables. If std(k) ≈ 0, k
stays constant and the periodic stimulus and the
oscillator have a stable phase relationship. The
stimulus we use, is the same used to generate the
Phase Response-Curves in Section 4, but in a peri-
odic fashion. It is a square pulse of fixed length of

L = 0.05 h with amplitude A. The stimulus period
T<sub>st</sub> determines the distance of two successive square
pulses.

In Fig. 9, the entrainment region of the PWA-
oscillator is shown. It exhibits a well-entrained
region for stimulus periods which are smaller than
the freerunning period of T<sub>fr</sub> ≈ 27.2 (T<sub>st</sub> < T<sub>fr</sub>), but
a much smaller entrainment region for T<sub>st</sub> > T<sub>fr</sub>.

This behaviour is reflected in the Phase Response-
Curve of the PWA oscillator. Note, that the stimu-
lus amplitude used for the PRC is A = −0.5 and the
the maximal phase advance (period shortening) was
0.4 h. This corresponds to the minimal stimulus
period that entrains the oscillator in Fig. 9, which
is ≈ 0.4 h shorter than the freerunning period.

A single stimulus can advance the phase of the
oscillator (and thus shorten its period) if it arrives
in the window $x > 0$. But the possible phase de-
lays are much smaller, leading only to a marginal
lengthening of the oscillator period.

As we have described in Section 4, this structure
of large phase advances but marginal phase delays
in the PRC is a direct result of the structure of
system (7).

The largest possible phase delay translates di-
rectly into the largest T<sub>st</sub> that is able to entrain the
oscillator. It is therefore the direct consequence of
the interlocking of variables b and p due to the se-
questration process, that gives rise to this one-sided
Arnold's tongue.

To compare, Fig. 10 shows the Arnold's tongue of the Goodwin oscillator. Here we simulated the oscillator for 5000 h and analyzed the last 25% of the time series. It exhibits entrainment for stimuli with $T_{st} < T_{fr}$ and $T_{st} > T_{fr}$. Again, this is due to the structure of the phase response-curve, which shows that a square wave stimulus can advance or delay the phase of the Goodwin oscillator, depending on the phase relationship between the two. Consequently, the Goodwin oscillator can be entrained by periodic stimuli with shorter but also with longer periods than its freerunning period.

Comparison of the entrainment properties for the two oscillators suggests that the role of the sequestration motif may be to restrict the range of responses of the system to external perturbations, and thus introduce a form of robustness of the oscillator by averting deviations from the original period.

# 6 Discussion

We study a model of the major transcriptional and post-translational processes in mammalian circadian clocks.

The model describes the interacting dynamics of core clock proteins DBP, REV-ERB and protein complexes CLOCK:BMAL1, PER:CRY. The sequestering of complexes CLOCK:BMAL1 and PER:CRY is a major post-translational process in the core circadian clock network. It has been modeled in different fashions in the past. In [7] the process is not modeled directly but represented as multiplicative competition in the synthesis terms in which BMAL1 plays an activating role. Kim and Forger [15] model the sequestration as a static process starting from the laws of mass action with the additional assumption of constant activator concentration. François and Hakim [14] use mass-action dynamics for their description of an oscillator with one feedback loop. In addition, they consider an irreversible sequestration, so no dissociation is taken into account. We use an oscillator with two feedback loops to describe the core circadian clock, where the sequestration is modelled dynamically, like in [14].

We show that the mass-action process describing the sequestration separates two timescales in the model, when appropriately transformed. The timescale separation is a result of the tightly inter-

locked dynamics of the sequestering species. Because of the interlocking, it is the difference of concentrations [CLOCK:BMAL1] and [PER:CRY] that moves on the slow timescale, together with the dynamics of [DBP] and [PER]. We use this timescale separation to show that our system can be easily reduced and transformed into a 3 dimensional piecewise affine system with linear time-invariant subsystems.

By systematically excluding specific transitions between the subsystems, we show that under certain conditions, only one sequence of periodic transitions remains possible. Using the Brouwer fixed point theorem, we then prove the existence of a periodic orbit. Furthermore, we show numerically that all posed conditions are fulfilled for a rather large region in parameter space.

Next, the Phase Response-Curve for the system is studied. It depends highly on the state of the aforementioned difference [CLOCK:BMAL1] - [PER:CRY]. It shows a characteristic shape, being almost zero almost everywhere and positive in the region where the difference variable is positive. In terms of concentrations, the oscillator is, paradoxically, almost insensitive to stimuli in [PER:CRY] when [PER:CRY] is nonzero. It is, however, very sensitive to stimuli, when [PER:CRY] is zero. This is due to the fast sequestering of any surplus [PER:CRY], which decreases [CLOCK:BMAL1] and has downstream effects in [DBP] and [REV].

In contrast, perturbation on the inhibiting factor in a model with strongly cooperative Hill-type repression (Goodwin oscillator), produces both phase advances and phase delays. We exemplify this notion further by showing the regions of entrainment for two different oscillators, our sequestering-type repression oscillator and the Goodwin oscillator with Hill-type repression. The Goodwin oscillator can be entrained by a periodic stimulus which has a longer period than the freerunning oscillator period. The periodic stimulus can force the Goodwin oscillator to share the stimulus period and have a stable phase relationship to the stimulus, given that its amplitude is high enough. The sequestration type oscillator, on the other hand, cannot be entrained by a stimulus with a period significantly higher than its own freerunning period. The sequestration motif hinders the entrainment of longer period stimuli.

There is to our knowledge no experimental evidence for a PRC with large constant zone and rather
tight window of nontrivial phase responses. Usually, experimental phase response curves show zones of phase advances and phase delays [35]. However, we show that if the repressor in a sequestration-based oscillator is subjected to external input, the tight binding of repressor to activator should lead to phase-response curves of this kind.

In the context of biological clocks, this article investigates thoroughly the sequestration repression motif. We have shown that the dynamical modelling of the sequestration process introduces a symmetry in the concentration dynamics. Because of the nonlinearity of the mass-action process, a timescale separation independent of the sequestration rate is possible after change of variables. Furthermore, we use this separation to create a piecewise affine system and prove the existence of a periodic orbit. Finally, we show that the tight interlocking of the activator and repressor species gives rise to a unique phase-response curve. Altogether, our results highlight the unique dynamics due to sequestration repression compared to highly cooperative repression types.

## Acknowledgments

This work was supported in part by the French National Agency for Research through project InSync ANR-22-CE45-0012-01.

## References

[1] Alina Patke, Michael W. Young, and Sofia Axelrod. "Molecular Mechanisms and Physiological Importance of Circadian Rhythms". In: *Nature Reviews Molecular Cell Biology* 21.2 (Feb. 2020), pp. 67–84. DOI: 10.1038/s41580-019-0179-2.

[2] Joseph S. Takahashi. "Transcriptional Architecture of the Mammalian Circadian Clock". In: *Nature Reviews Genetics* 18.3 (Mar. 2017), pp. 164–179. DOI: 10.1038/nrg.2016.150.

[3] Ueli Schibler et al. "Clock-Talk: Interactions between Central and Peripheral Circadian Oscillators in Mammals". In: *Cold Spring Harbor Symposia on Quantitative Biology* 80 (2015), pp. 223–232. DOI: 10.1101/sqb.2015.80.027490.

[4] Michael H. Hastings, Elizabeth S. Maywood, and Marco Brancaccio. "Generation of Circadian Rhythms in the Suprachiasmatic Nucleus". In: *Nature Reviews Neuroscience* 19.8 (Aug. 2018), pp. 453–469. DOI: 10.1038/s41583-018-0026-z.

[5] Albert Goldbeter. "A Model for Circadian Oscillations in the *Drosophila* Period Protein (PER)". In: *Proceedings of the Royal Society of London. Series B: Biological Sciences* 261.1362 (Sept. 22, 1995), pp. 319–324. DOI: 10.1098/rspb.1995.0153.

[6] Lindsey S. Brown and Francis J. Doyle. "A Dual-Feedback Loop Model of the Mammalian Circadian Clock for Multi-Input Control of Circadian Phase". In: *PLOS Computational Biology* 16.11 (Nov. 23, 2020). Ed. by Bard Ementrout, e1008459. DOI: 10.1371/journal.pcbi.1008459.

[7] Janina Hesse et al. "A Mathematical Model of the Circadian Clock and Drug Pharmacology to Optimize Irinotecan Administration Timing in Colorectal Cancer". In: *Computational and Structural Biotechnology Journal* 19 (2021), pp. 5170–5183. DOI: 10.1016/j.csbj.2021.08.051.

[8] Didier Gonze et al. "Spontaneous Synchronization of Coupled Circadian Oscillators". In: *Biophysical Journal* 89.1 (July 2005), pp. 120–129. DOI: 10.1529/biophysj.104.058388.

[9] Sofia Almeida, Madalena Chaves, and Franck Delaunay. "Transcription-Based Circadian Mechanism Controls the Duration of Molecular Clock States in Response to Signaling Inputs". In: *Journal of Theoretical Biology* 484 (Jan. 2020), p. 110015. DOI: 10.1016/j.jtbi.2019.110015.

[10] Bruce Alberts. *Molecular Biology of the Cell*. Seventh edition. New York: W. W. Norton & Company, 2022.

[11] Hiroki R Ueda et al. "System-Level Identification of Transcriptional Circuits Underlying Mammalian Circadian Clocks". In: *Nature Genetics* 37.2 (Feb. 2005), pp. 187–192. DOI: 10.1038/ng1504.
[12] Guillaume Rey et al. "Genome-Wide and Phase-Specific DNA-Binding Rhythms of BMAL1 Control Circadian Output Functions in Mouse Liver". In: *PLoS Biology* 9.2 (Feb. 22, 2011). Ed. by Achim Kramer, e1000595. DOI: 10.1371/journal.pbio.1000595.

[13] Nicolas E Buchler and Frederick R Cross. "Protein Sequestration Generates a Flexible Ultrasensitive Response in a Genetic Network". In: *Molecular Systems Biology* 5.1 (Jan. 2009), p. 272. DOI: 10.1038/msb.2009.30.

[14] Paul François and Vincent Hakim. "Core Genetic Module: The Mixed Feedback Loop". In: *Physical Review E* 72.3 (Sept. 16, 2005), p. 031908. DOI: 10.1103/PhysRevE.72.031908.

[15] Jae Kyoung Kim and Daniel B Forger. "A Mechanism for Robust Circadian Timekeeping via Stoichiometric Balance". In: *Molecular Systems Biology* 8.1 (Jan. 2012), p. 630. DOI: 10.1038/msb.2012.62.

[16] Corentin Briat, Ankit Gupta, and Mustafa Khammash. "Antithetic Integral Feedback Ensures Robust Perfect Adaptation in Noisy Biomolecular Networks". In: *Cell Systems* 2.1 (Jan. 2016), pp. 15–26. DOI: 10.1016/j.cels.2016.01.004.

[17] Luis Treviño-Quintanilla, Julio Freyre-González, and Irma Martínez-Flores. "Anti-Sigma Factors in E. Coli: Common Regulatory Mechanisms Controlling Sigma Factors Availability". In: *Current Genomics* 14.6 (Sept. 1, 2013), pp. 378–387. DOI: 10.2174/1389202911314060007.

[18] Nicolas E. Buchler and Matthieu Louis. "Molecular Titration and Ultrasensitivity in Regulatory Networks". In: *Journal of Molecular Biology* 384.5 (Dec. 2008), pp. 1106–1119. DOI: 10.1016/j.jmb.2008.09.079.

[19] Jae Kyoung Kim et al. "Molecular Mechanisms That Regulate the Coupled Period of the Mammalian Circadian Clock". In: *Biophysical Journal* 106.9 (May 2014), pp. 2071–2081. DOI: 10.1016/j.bpj.2014.02.039.

[20] Daniel B Forger. *Biological Clocks, Rhythms, and Oscillations*. Cambridge, Massachusetts, USA: MIT Press, 2017.

[21] Jae Kyoung Kim. "Protein Sequestration versus Hill-type Repression in Circadian Clock Models". In: *IET Systems Biology* 10.4 (Aug. 2016), pp. 125–135. DOI: 10.1049/iet-syb.2015.0090.

[22] Paul François. "A Model for the Neurospora Circadian Clock". In: *Biophysical Journal* 88.4 (Apr. 2005), pp. 2369–2383. DOI: 10.1529/biophysj.104.053975.

[23] S. Almeida, M. Chaves, and F. Delaunay. "Control of Synchronization Ratios in Clock/Cell Cycle Coupling by Growth Factors and Glucocorticoids". In: *Royal Society Open Science* 7.2 (Feb. 2020), p. 192054. DOI: 10.1098/rsos.192054.

[24] Odile Burckard et al. "Cycle Dynamics and Synchronization in a Coupled Network of Peripheral Circadian Clocks". In: *Interface Focus* 12.3 (June 6, 2022), p. 20210087. DOI: 10.1098/rsfs.2021.0087.

[25] Pierre Sacré and Rodolphe Sepulchre. "Sensitivity Analysis of Circadian Entrainment in the Space of Phase Response Curves". In: *A Systems Theoretic Approach to Systems and Synthetic Biology II: Analysis and Design of Cellular Systems*. Ed. by Vishwesh V. Kulkaneri, Guy-Bart Stan, and Karthik Raman. Dordrecht: Springer Netherlands, 2014, pp. 59–81. DOI: 10.1007/978-94-017-9047-5_3.

[26] Arthur T. Winfree. *The Geometry of Biological Time*. Red. by J. E. Marsden, S. Wiggins, and L. Sirovich. Vol. 12. Interdisciplinary Applied Mathematics. New York, NY: Springer New York, 2001. DOI: 10.1007/978-1-4757-3484-3.

[27] Pierre Sacre and Rodolphe Sepulchre. "Sensitivity Analysis of Oscillator Models in the Space of Phase-Response Curves: Oscillators As Open Systems". In: *IEEE Control Systems Magazine* 34.2 (Apr. 2014), pp. 50–74. DOI: 10.1109/MCS.2013.2295710.

[28] A. Granada et al. "Chapter 1 Phase Response Curves". In: *Methods in Enzymology*. Vol. 454. Elsevier, 2009, pp. 1–27. DOI: 10.1016/S0076-6879(08)03801-9.
[29] Lee A. Segel and Marshall Slemrod. “The Quasi-Steady-State Assumption: A Case Study in Perturbation”. In: *SIAM Review* 31.3 (Sept. 1989), pp. 446–477. doi: 10.1137/1031091.

[30] John Guckenheimer and Philip Holmes. *Non-linear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields*. Vol. 42. Applied Mathematical Sciences. New York, NY: Springer New York, 1983. doi: 10.1007/978-1-4612-1140-2.

[31] James R. Munkres. *Topology*. 2nd ed. Upper Saddle River, NJ: Prentice Hall, Inc, 2000. 537 pp.

[32] Eberhard Zeidler. *Nonlinear Functional Analysis and its Applications*. Springer-Verlag New York, Inc., 1985.

[33] Didier Gonze and Peter Ruoff. “The Goodwin Oscillator and Its Legacy”. In: *Acta Biotheretica* 69.4 (Dec. 2021), pp. 857–874. doi: 10.1007/s10441-020-09379-8.

[34] Stefano Boccaletti et al. *Synchronization: From Coupled Systems to Complex Networks*. 1st ed. Cambridge University Press, Mar. 29, 2018. doi: 10.1017/9781107297111.

[35] Vijay Kumar Sharmad. “Light-Induced Phase Response Curves of the Circadian Activity Rhythm in Individual Field Mice, Mus Booduga”. In: *Chronobiology International* 13.6 (Jan. 1996), pp. 401–409. doi: 10.3109/07420529609020911.

# A Proofs for cyclic transitions

## A.1 Preliminaries

**Lemma A.1** (upper bound on $x$ after stage rXd). Suppose that the previous stage was rxd. Suppose that the stage rXd finishes by transitioning to stage rXD. The maximal value variable $x$ can assume is

$$ \bar{x} := \frac{\sqrt{\beta^2 d^{*2} + 2Hd^* - \beta d^*}}{H}, $$

where $H = 1 - \epsilon d^*$.

*Proof.* $x$ crosses 0 during the transition rxd → rXd. We can therefore calculate an upper bound on $x$ at the end of stage rXd. In order to do this, we need

to find an upper bound on the time it takes for $d$ to reach its threshold since this marks the end of the stage.

$$ \begin{aligned} \dot{x} &\ge 1 - \epsilon \sup d = 1 - \epsilon d^* \\ &\Rightarrow x(t) \ge (1 - \epsilon d^*)t \end{aligned} $$

$$ \begin{aligned} \dot{d} &\ge x - \beta \sup d = (1 - \epsilon d^*)t - \beta d^* \\ &\Rightarrow d(t) \ge \frac{1}{2}(1 - \epsilon d^*)t^2 - \beta d^*t \end{aligned} $$

We used $\sup d = d^*$ because we are in stage rXd. Inserting the time $\tau_d$ when $d$ reaches its threshold $d^*$, we get an upper bound for it:

$$ \begin{aligned} d(\tau_d) &= d^* \ge \frac{1}{2}(1 - \epsilon d^*)\tau_d^2 - \beta d^* \tau_d \\ &\ge \frac{1}{2}\underbrace{(1 - \epsilon d^*)}_{H}\tau_d^2 - \beta d^* \tau_d - d^* \\ &\le \frac{\sqrt{\beta^2 d^{*2} + 2Hd^* - \beta d^*}}{H}. \end{aligned} $$

Note, that $\dot{x} \le 1 - \epsilon \inf d = 1$. This can be used to calculate a maximal value of $d$ at the end of the stage:

$$ x(\tau_d) \le \tau_d \le \frac{\sqrt{\beta^2 d^{*2} + 2Hd^* - \beta d^*}}{H} =: \bar{x}. $$

□

**Lemma A.2** (lower bound on $x$ after stage rXd). Suppose that the previous stage was rxd. Suppose that the stage rXd finishes by transitioning to stage rXD. The minimal value variable $x$ can assume is

$$ \underline{x} := H\sqrt{d^*}, $$

where $H = 1 - \epsilon d^*$.

*Proof.* $x$ crosses 0 during the transition rxd → rXd. We can therefore calculate a lower bound on $x$ at the end of stage rXd. In order to do this, we need to find a lower bound on the time it takes for $d$ to reach its threshold since this marks the end of the stage.

$$ \begin{aligned} \dot{x} &\le 1 - \epsilon \inf d = 1 \\ &\Rightarrow x(t) \le t \end{aligned} $$

$$ \begin{aligned} \dot{d} &\le x - \beta \inf d = t \\ &\Rightarrow d(t) \le \frac{1}{2}t^2 \end{aligned} $$
Using the time $\tau_d$ when $d$ reaches its threshold $d^*$:

$$d(\tau_d) = d^* \le \frac{1}{2}\tau_d^2$$
$$\tau_d \ge \sqrt{d^*}.$$

Note, that $\dot{x} \ge 1 - \epsilon \sup d = 1 - \epsilon d^*$. This can be used to calculate a lower bound of $x$ at the end of the stage:

$$x(\tau_d) \ge (1 - \epsilon d^*) \tau_d \ge (1 - \epsilon d^*) \sqrt{d^*} =: \underline{x}.$$

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

**Assumption A.3.** Assume $\bar{x} > \underline{x}$.

**Remark.** Assumption A.3 implies that the upper bound on $x$ is greater than the lower bound on $x$ in stage rXD, since $\bar{x}' > \bar{x}$. The bounds on $X$ are therefore $\underline{x} \le x \le \bar{x}$ after stage rXd and $\underline{x} \le x \le \bar{x}'$ after stage rXD.

## A.2 Transition rXD to RXD

**Lemma A.4 (threshold conditions rXD).** Suppose stage rXD was reached from stage rXd. The condition $d = 1/\epsilon$ must be fulfilled before either $d$ or $x$ reach their respective threshold.

*Proof.* The trajectory has crossed the threshold $d = d^*$ from below. Since stages rXd and rXD are governed by the same equations, condition $\dot{d} > 0$ holds on both sides of the transition, which implies $x > \beta\delta/\gamma$. We also know, that $\dot{x} > 1 - \epsilon d^* > 0$. In summary, $\dot{d} > 0, \dot{x} > 0$ at the time of threshold crossing. This holds true at least until $d = 1/\epsilon$ at which point $\dot{x} = 0$ and a subsequent decrease in $x$ is followed by decreasing $d$. □

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

**Theorem A.5.** Suppose the stage rXD was reached by stage rXd. If the following inequality is fulfilled, the system can only transition to stage RXD:

$$\sqrt{\frac{2}{\gamma(\underline{x} - \beta/\epsilon)}} < \frac{\sqrt{G^2 + 2H^2/\epsilon} - G}{H}.$$

*Proof.* By lemma 3.2, the time it takes to reach $r$'s threshold is bounded from above by

$$T_r = \sqrt{\frac{2}{\gamma(\underline{x} - \beta/\epsilon)}}.$$

Lemma 3.3 shows that the threshold $d = 1/\epsilon$ must be crossed before the system reaches either $d$'s or $x$'s threshold. By lemma 3.4, the time it takes to reach $d = 1/\epsilon$ is bounded from below by

$$T_d = \frac{\sqrt{G^2 + 2H^2/\epsilon} - G}{H}.$$

If $T_r < T_d$, the only possible transition is rXD $\rightarrow$ RXD. □

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

## A.3 Transition RXD to RxD

**Lemma A.6 (initial condition).** Suppose that a trajectory followed the stage progression rXd, RXD and is now in stage RXD. The trajectory traversing stage RXD is bounded from below by the trajectory starting with the initial condition $x = \underline{x}, d = d^*$.

*Proof.* In stage rXD, the following holds by Theorem 3.3: $\dot{x} > 0, \dot{d} > 0$. Any initial condition in RXD, will therefore have an initial condition $(x_{IC}, d_{IC}) \ge (\underline{x}, d^*)$, where equality is attained only if $T_r = 0$. Since the dynamics are linear and trajectories in 2-d cannot cross because of uniqueness, the real trajectory is bounded from below by the one starting from $(\underline{x}, d^*)$. □

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

**Theorem A.7 (final condition).** Under the following conditions, the bounding Jordan curve $\mathcal{J}$ will terminate at a point $O = (x=0, d > d^*)$:

$$x_m(t_1) \left( \frac{1}{\beta} - \frac{1}{2\epsilon d^*} x_m(t_1) \right) > d^* \quad (15)$$

where

$$x_m(t_1) = \underline{x} \exp \left( - \frac{\epsilon}{\beta^2 - \epsilon} \ln \left( \beta \frac{d^* + F}{\underline{x} + \beta F} \right) \right)$$
$$F = \frac{\beta \underline{x}}{\epsilon - \beta^2}.$$

*Proof.* The curve $\mathcal{J}$ is described in Theorem 3.6. It is produced by the trajectory of the system:

$$\begin{cases} (x_m, d_m) & \text{if } (x, d) \in \mathcal{A} \\ (x_p, d_p) & \text{if } (x, d) \in \mathcal{B} \end{cases}$$

with

$$\dot{x}_m = -\frac{\epsilon}{\beta}x_m$$

$$\dot{d}_m = x_m - \beta d_m$$

$$\dot{x}_p = -\epsilon d^*$$

$$\dot{d}_p = x_p - \beta d_p.$$
First, We will integrate the trajectory starting from $(\underline{x}, d^*)$ and reaching $x = \beta d$.

$$x_m(t) = \underline{x} \exp \left( -\frac{\epsilon}{\beta} t \right)$$

$$d_m(t) = -F \exp \left( -\frac{\epsilon}{\beta} t \right) + (d^* + F) \exp(-\beta t)$$

$$\text{where } F = \frac{\beta \underline{x}}{\epsilon - \beta^2}.$$

Using the ending condition $x(t_1) = \beta d(t_1)$:

$$(F + \frac{\underline{x}}{\beta}) \exp \left( -\frac{\epsilon}{\beta} t_1 \right) = (d^* + F) \exp(-\beta t_1)$$

$$t_1 = \frac{1}{\beta - \epsilon/\beta} \ln \left( \frac{d^* + F}{\underline{x}/\beta + F} \right).$$

Second, we will integrate the trajectory starting from $(x_p(0) = x_m(t_1), d_p(0) = d_m(t_1))$ and reaching $x_p(t_2) = 0$:

$$x_p(t) = x_m(t_1) - \epsilon d^* t$$

$$d_p(t) = b_0 + b_1 t + b_2 \exp(-\beta t) \quad (16)$$

$$\text{where } b_0 = \frac{x_m(t_1)}{\beta} + \frac{\epsilon d^*}{\beta^2}$$

$$b_1 = \frac{\epsilon d^*}{\beta}$$

$$b_2 = -\frac{\epsilon d^*}{\beta^2}.$$

The final condition yields:

$$\begin{align}
x_p(t_2) &= x_m(t_1) - \epsilon d^* t_2 = 0 \nonumber \\
\Rightarrow \quad t_2 &= \frac{x_m(t_1)}{\epsilon d^*}. \tag{17}
\end{align}$$

Last, we need to show that $d_p(t_2) > d^*$. Keeping in mind that $\exp(-x) < (1 - x + x^2/2) \quad \forall x > 0$ and $b_2 < 0$, we can calculate a lower bound:

$$d_p(t_2) = b_0 + b_1 t_2 + b_2 \exp(-\beta t_2)$$

$$\begin{align*}
&> b_0 + b_1 t_2 + b_2 \left( 1 - \beta t_2 + \frac{\beta^2 t_2^2}{2} \right) \\
&= \frac{x_m(t_1)}{\beta} + t_2(b_1 - \beta b_2) + t_2^2 \left( \frac{1}{2} \beta^2 b_2 \right)
\end{align*}$$

$$> \frac{x_m(t_1)}{\beta} - \frac{1}{2} \epsilon d^* t_2^2$$

$$= x_m(t_1) \left( \frac{1}{\beta} - \frac{x_m(t_1)}{2\epsilon d^*} \right)$$

where

$$x_m(t_1) = \underline{x} \exp \left( -\frac{\epsilon}{\beta^2 - \epsilon} \ln \left( \beta \frac{d^* + F}{\underline{x} + \beta F} \right) \right).$$

Now, if

$$x_m(t_1) \left( \frac{1}{\beta} - \frac{x_m(t_1)}{2\epsilon d^*} \right) > d^* \quad (18)$$

it follows that $d_p(t_2) > d^*$. Therefore if a set of parameters fulfills (18), the jordan curve will terminate at $O = (x_p = 0, d_p(t_2) > d^*)$. □