# Feedback-delay dependence of the stability of cluster periodic orbits in
populations of degrade-and-fire oscillators with common activator

Bastien Fernandez  Laboratoire de Probabilités, Statistique et Modélisation
CNRS - Univ. Paris Cité - Sorbonne Univ.
Paris, France  Matteo Tanzi  Laboratoire de Probabilités, Statistique et
Modélisation
CNRS - Univ. Paris Cité - Sorbonne Univ.
Paris, France  Department of Mathematics, King’s College London
Strand Building, The Strand
London, UK

###### Abstract

Feedback delay has been identified as a key ingredient in the quorum sensing
synchronization of synthetic gene oscillators. While this influence has been
evidenced at the theoretical level in a simplified system of degrade-and-fire
oscillators coupled via a common activator protein, full mathematical
certifications remained to be provided. Here, we prove from a rigorous
mathematical viewpoint that, for the very same model, the synchronized
degrade-and-fire oscillations are 1/ unstable with respect to out-of-sync
perturbations in absence of delay, and 2/ are otherwise asymptotically stable
in presence of delay, no matter how small is its amplitude. To that goal, we
proceed to an extensive study of the population dynamics in this system, which
in particular identifies the mechanisms of, and related criteria for, the
delay-dependent stability of periodic orbits with respect to out-of-sync
perturbations. As an additional outcome, the analysis also reveals that,
depending on the parameters, multiple stable partially synchronized periodic
orbits can coexist with the fully synchronized one.

.

##  1 Introduction

Starting with the toggle-switch and repressilator [2, 3], elementary
regulatory circuits have been proposed, and implemented, as basic building
blocks of gene networks in Synthetic Biology. By understanding (and by
controlling) the functioning of simple representative examples, this now
popular field of research intends to yield advances in bio-engineering and
medical applications, beyond intrinsic interest to fundamental biology [16].

A more recent development of Synthetic Biology aims at investigating the
collective dimension of the regulatory dynamics in populations composed by
many individuals equipped with simple genetic circuits. In particular, in
suitably designed populations of quorum sensing oscillators, stunning evidence
of fully synchronized oscillations has been obtained, bringing a standard
notion in Physics into the realm of micro-biological colonies [1].

This experimental phenomenology has called for theoretical conceptualization
based on (simplified) mathematical models for the underlying systems dynamics.
In particular, the state of each individual in [1] is represented by the
concentration of an auto-repressor protein. Oscillatory behaviour then results
from a negative feedback loop for the protein concentration level. The system
has been designed so that the production occurs via short but large bursts.
Fast production is followed by a slow degradation over longer time intervals.
Under appropriate considerations, this degrade-and-fire (DF) mechanism can be
represented by a differential equation whose vector field is negative constant
(corresponding to degradation) and where the concentration is instantaneously
reset when it reaches 0 (firing) [17].

In simple models of coupled DF oscillators, rigorous proofs of a sharp
transition from a virtually uncoupled regime to massive clustering upon
increase of the interaction strength, have been established [4, 5] and confirm
appropriate modelling of the phenomenology.

In a more elaborated model [18], the quorum sensing mechanism that favours
synchrony relies on the presence of a common activator protein that increases
the amplitude of the firings. The activator concentration is encoded in an
additional variable that is coupled to the mean concentration of the
individual repressors. Moreover, as a simplification of the detailed
mechanistic model in [1], a delay has been introduced into the dynamics, which
affects the activator concentration involved in the firings. Numerics and
theoretical investigations have revealed that this delay plays a crucial role
in stabilizing the synchronized oscillations.

Delayed interactions are known to have a significant impact of the functioning
of biological systems, for instance in the regulation of the synchronization
of oscillations in gene networks during development [10, 13, 14]. There exists
a large literature on the Lyapunov stability of synchronized oscillations in
systems of delayed differential equations [9, 11, 19]. The literature also
includes instances of systems whose characteristics are close (but distinct)
to the degrade-and-fire oscillators and for which stability is shown to
sharply depend on the parameter(s) [7, 15]. The crucial role played by the
delay to favor synchronization has been highlighted in different situations.
For example in [12], it has been shown how a time-delayed coupling between two
oscillators can result into synchronization for arbitrarily small size of the
coupling strength (while perturbation arguments show that synchrony cannot be
achieved in absence of delay).

Back to [18], the authors have carried out a numerical investigation of the
typical orbits arising in their model under various conditions. In particular,
different values of the degradation rate of the concentration of the activator
protein, different responses of the reset value on the concentration of the
activator, different delays, and different amplitude of the noise (added at
each firing event) have been considered. The system exhibits a plethora of
behaviours, especially synchronization, stable and metastable clustering, and
absence of synchronization. One of the main outcomes of these investigations
is that the delay in the coupling has a crucial role in favoring the
synchronization of the concentrations.

The goal of the present paper is to mathematically certify the observations in
[18]. More precisely, we prove in particular that for the deterministic
dynamics of the model therein, the synchronized oscillations are unstable
(with respect to out-of-sync perturbations) in absence of delay, and that they
are otherwise asymptotically stable for every positive delay.

The paper is organized as follows. Firstly, we recall the definition of the
model and we provide the basic properties that are useful for the analysis of
its dynamics. Then, we proceed in Section 3 to the study of the existence and
the delay-dependent stability of the fully synchronized periodic orbit
(Proposition 3.1). This analysis identifies the main mechanisms involved in
that phenomenology. These elements are further developed in Section 4 which
presents a systematic approach to the stability of arbitrary periodic orbits
with partially synchronized repressor concentrations. In particular, a
criterion for instability in absence of delay and another criterion for
stability in presence of delay are established. The theory culminates with a
delay-dependent existence and stability statement (Theorem 4.4) which is based
on the corresponding properties inside the partially synchronized subspace, in
absence of delay. Theorem 4.4 is the extension of Proposition 3.1 to arbitrary
partially synchronized orbits. Finally, an example of application to periodic
orbits with several clusters of equi-distributed repressor concentrations is
given, which shows in particular that several forms of asymptotically stable
(partially) synchronized oscillations can coexist in this system.

##### Aknowledgments:

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Sklodowska-Curie grant
agreement No 843880.

##  2 The dynamical system and its basic properties

###  2.1 Definition of the dynamics

Following [18], we consider a population of $N$ DF oscillators
($N\in\mathbb{N}$) represented by the variable
$(\mathbf{x},A)\in(\mathbb{R}^{+})^{N+1}$ where
$\mathbf{x}=(x_{1},\cdots,x_{N})\in(\mathbb{R}^{+})^{N}$ collects the
concentrations $x_{i}$ of the repressor proteins and $A$ denotes the activator
concentration.

The dynamics can be depicted as follows. Each repressor concentration decays
independently at constant speed -1. When it reaches zero, it is
instantaneously reset to a value that depends on the activator concentration.
We shall refer to reset events as firings. In addition, the repressor proteins
contribute to the synthesis of activator, which itself also degrades at
constant rate. In formal terms, time variations of the variable
$(\mathbf{x},A)$ are governed by the following coupled equations

$\displaystyle\left\\{\begin{array}[]{lcc}\dot{x}_{i}(t)=-1&\text{if}&x_{i}(t)>0\\\
x_{i}(t^{+})=R+\nu
A(t-\tau)&\text{if}&x_{i}(t)=0\end{array}\right.\quad\quad\forall
i\in\\{1,\cdots,N\\}$ (3)
$\displaystyle\begin{array}[]{l}\dot{A}(t)=m(t)-\beta A(t)\quad\text{if}\quad\
m(t^{+})=m(t)\quad\text{where}\quad
m(t)=\frac{1}{N}\sum_{i=1}^{N}x_{i}(t)\end{array}$ (5)

with initial repressor concentration vector
$\mathbf{x}(0)=\mathbf{x}^{0}\in\mathbb{R}^{+}$ and initial activator
concentration profile $A|_{[-\tau,0]}\in(\mathbb{R}^{+})^{[-\tau,0]}$.

The dynamics depends on four parameters,111The definition in [18] has an
additional parameter for the repressors decay rate. However, this rate can be
set to 1 by an appropriate rescaling of the variables and other parameters.
namely $R,\beta,\nu\in\mathbb{R}^{+}_{\ast}$ and $\tau\in\mathbb{R}^{+}$ for
which, for the sake of the analysis, we impose the conditions222Notice that
the inequalities $\nu<\beta$ and $\tau<R$ hold in all numerical results in
[18].

$\nu<\beta\quad\text{and}\quad\tau<R,$

and the inequalities (11) and (14) below. Some of the formal statements below
explicitly express the dependence on certain parameters, especially $\beta$
and $\tau$. In these cases, all other parameters are implictly assumed to be
given beforehand.

As norms are concerned, both in $\mathbb{R}^{N}$ and for real functions, we
shall use the following notations

$\|\mathbf{x}\|_{N}=\max_{i\in\\{1,\cdots,N\\}}|x_{i}|\quad\text{and}\quad\|A|_{I}\|_{0}:=\sup_{t\in
I}|A(t)|,$

where $N\in\mathbb{N}$ and the interval $I$ are arbitrary.

###  2.2 Basic considerations and elementary properties

####  2.2.1 Existence of global trajectories and conditions for well-posedness
in $(\mathbb{R}^{+})^{N+1}$

##### Existence and uniqueness of global solutions.

Given an arbitrary function $A:[-\tau,+\infty)\to\mathbb{R}^{+}$, an index
$i\in\\{1,\cdots,N\\}$ and $x_{i}^{0}\in\mathbb{R}^{+}$, equation (3) for the
single real variable $x_{i}$, with initial condition $x_{i}(0)=x_{i}^{0}$,
trivially admits a unique piecewise linear left continuous solution
$x_{i}:[0,+\infty)\to\mathbb{R}^{+}$ with constant slope $-1$ and positive
jump discontinuities at firings.

Independently, let $m:[0,+\infty)\to\mathbb{R}^{+}$ be an arbitrary left
continuous piecewise affine function with finitely many discontinuities in
every bounded interval and let $A^{0}\in\mathbb{R}^{+}$ be arbitrary. Equation
(5) with initial condition $A(0)=A^{0}$ admits a unique continuous solution
$A(t)=\phi^{t}_{m}(A^{0})$ where
$\phi^{\cdot}_{m}:[0,+\infty)\to\mathbb{R}^{+}$ is defined by the following
variation of constant formula333Indeed, equation (5) determines
$\phi_{m}^{t}(A^{0})$ on a dense subset of $\mathbb{R}^{+}$, which can then be
uniquely extended to the entire $\mathbb{R}^{+}$ by continuity.

$\phi^{t}_{m}(A):=\left(A+\int_{0}^{t}e^{\beta s}m(s)ds\right)e^{-\beta t}.$
(6)

In addition, this expression implies that $A(t)$ is Lipschitz continuous on
every bounded interval. Its Lipschitz constant is controlled by the supremum
of $m$ on the interval under consideration together with the value of $A$.

Put together, and since the mean value $\frac{1}{N}\sum_{i=1}^{N}x_{i}(t)$
associated with the solution must be left continuous and has finitely many
jumps on every interval $[0,t]$, the arguments above imply that the coupled
equations (3)-(5) have, given any initial datum
$(\mathbf{x}^{0},A|_{[-\tau,0]})$, a unique global solution
$(\mathbf{x}(t),A(t))$ such that $\mathbf{x}(0)=\mathbf{x}^{0}$ and $t\mapsto
A(t)$ is continuous on $t\in\mathbb{R}^{+}$.

From equation (3) and expression (6), one easily deduces that each oscillator
must fire infinitely often in every trajectory. Moreover, the time duration
between two consecutive firings of any given oscillator is equal to the reset
concentration at the previous firing; in particular this time cannot be
smaller than $R$. In addition, the first firing time of oscillator $i$ is
$x^{0}_{i}$.

##### Condition for well-posedness in $(\mathbb{R}^{+})^{N+1}$.

To provide the profile $A|_{[-\tau,0)}$ in the initial datum only serves to
specify the reset concentration(s) at any firing that would occur in the time
interval $[0,\tau)$. Actually, since $\tau<R$, only the first firing of
some/all oscillators - those firings at the times $x^{0}_{i}<\tau$ \- can
occur in this interval, because any reset value is at least $R$. Therefore, it
suffices to provide ($\mathbf{x}^{0},A^{0}$ and) the values
$A(x^{0}_{i}-\tau)$ for $x^{0}_{i}<\tau$, in order to define the trajectory,
ie. the dynamics is indeed well-posed in finite dimension.

In particular, for $\tau=0$, the knowledge of $(\mathbf{x}^{0},A^{0})$
suffices to define the trajectory; hence the dynamics is well-posed in
$(\mathbb{R}^{+})^{N+1}$, which is particularly convenient for the stability
analysis of periodic orbits.

When $\tau>0$, if

$\text{min}_{\mathbf{x}^{0}}:=\min_{i}x_{i}^{0}\geq\tau,$

then no firing can occur in the time interval $[0,\tau)$; hence the subsequent
trajectory is again well-defined given only $(\mathbf{x}^{0},A^{0})$. This is
also the case when $\text{min}_{\mathbf{x}^{0}}<\tau$ if it is impossible that
a firing occurs in the past time interval
$[\text{min}_{\mathbf{x}^{0}}-\tau,0)$. Indeed, one can reverse the time
direction in equation (5) in order to compute the values $A(x_{i}^{0}-\tau)$
for $x_{i}^{0}<\tau$ using $(\mathbf{x}^{0},A^{0})$. More precisely, we have
$A(x_{i}^{0}-\tau)=\phi_{m}^{x_{i}^{0}-\tau}(A^{0})$ where for
$t\in\mathbb{R}^{+}$, $\phi^{-t}_{m}(A)$ is defined by the backward time
variation of constant formula

$\phi^{-t}_{m}(A):=\left(A-\int_{0}^{t}e^{-\beta s}m(-s)ds\right)e^{\beta t},$
(7)

with $m(-s)=m^{0}+s$ for $s\in[0,\tau-\text{min}_{\mathbf{x}^{0}}]$ and
$m^{0}:=\frac{1}{N}\sum_{i=1}^{N}x^{0}_{i}$. We must also make sure that the
values $A(x_{i}^{0}-\tau)$ are non-negative, viz. $A^{0}\geq
A_{\mathbf{x}^{0}}(\tau)$, where

$A_{\mathbf{x}^{0}}(\tau):=\int_{0}^{\tau}e^{-\beta s}(m^{0}+s)ds,$

(which tends to 0 as $\tau\to 0$).

As expression (6), the backward time formula (7) is Lipschitz continuous on
every bounded interval. Given $\mathbf{x},A,\beta$ and $\tau$, let
$K:=K_{\mathbf{x},A,\beta,\tau}$ be its Lipschitz constant on $[0,\tau]$ when
computed with $m(-s)=\frac{1}{N}\sum_{i=1}^{N}x_{i}+s$.

######  Claim 2.1.

Assume that $\text{min}_{\mathbf{x}^{0}}\geq\tau$ or
$\text{min}_{\mathbf{x}^{0}}<\tau$ and

$\text{max}_{\mathbf{x}^{0}}:=\max_{i}x^{0}_{i}<R+\nu A^{0}+(\nu
K+1)\text{min}_{\mathbf{x}^{0}}-(2\nu K+1)\tau\quad\text{and}\quad A^{0}\geq
A_{\mathbf{x}}(\tau).$

Then the trajectory is well-defined given
$(\mathbf{x}^{0},A^{0})\in(\mathbb{R}^{+})^{N+1}$ and lies in
$(\mathbb{R}^{+})^{N+1}$.

Proof: The proof is immediate. Together with the Lipschitz continuity
mentioned above, the condition in the statement implies

$x_{i}(t)=\text{max}_{\mathbf{x}^{0}}+t<R+\nu A(t-\tau),\ \forall
t\in[\text{min}_{\mathbf{x}^{0}}-\tau,0],$

ie. no reset value can be attained by the backward flow; hence no firing can
occur in the time interval $[\text{min}_{\mathbf{x}^{0}}-\tau,0)$. $\Box$

Notice finally that the last inequality indicates that the following stronger
condition

$\text{max}_{\mathbf{x}^{0}}+\text{min}_{\mathbf{x}^{0}}<R-\tau\quad\text{and}\quad
A^{0}\geq A_{\mathbf{x}}(\tau),$ (8)

suffices to obtain the same conclusion. We shall use this sronger condition in
the stability analysis of the synchronized periodic orbit in Section 3.3.
Moreover, notice also that the condition in Claim 2.1 will be necessary for
the stability analysis of periodic orbits with equidistributed repressor
concentrations (see Section 4.6).

####  2.2.2 Attracting invariant set

The dissipative term in equation (5) suggests that the trajectories should
asymptotically approach a bounded forward invariant set. This property is
formally expressed in the next statement. Let

$A_{\text{max}}:=\frac{R}{\beta-\nu}\quad\text{and}\quad Q:=[0,R+\nu
A_{\text{max}}]^{N}\times[0,A_{\text{max}}].$ (9)

######  Lemma 2.2.

Assume that $\mathbf{x}^{0}\in[0,R+\nu A_{\text{max}}]^{N}$ and
$A|_{[-\tau,0]}\in[0,A_{\text{max}}]^{[-\tau,0]}$. Then the solution satisfies

$(\mathbf{x}(t),A(t))\in Q,\quad\forall t\in\mathbb{R}^{+}.$

In addition, given any initial datum $(\mathbf{x}^{0},A|_{[-\tau,0]})$, the
subsequent trajectory satisfies

$\limsup_{t\rightarrow+\infty}A(t)\leq A_{\text{max}},\quad\text{and
then}\quad\limsup_{t\rightarrow+\infty}x_{i}(t)\leq R+\nu
A_{\text{max}}\quad\forall i\in\\{1,\cdots,N\\}.$

In the rest of the paper, we always assume that $(\mathbf{x}(t),A(t))\in Q$
for all $t\in\mathbb{R}^{+}$, even when this is not explicitly stated.

###### Proof.

Assume that $\mathbf{x}^{0}\in[0,R+\nu A_{\text{max}}]^{N}$ and
$A|_{[-\tau,0]}\in[0,A_{\text{max}}]^{[-\tau,0]}$. Then, expression (6)
implies that $A(t)\geq 0$ for all $t\in\mathbb{R}^{+}$. We prove that
$A(t)\leq A_{\text{max}}$ for all $t\in\mathbb{R}^{+}$ (and then
$\max_{i}x_{i}(t)\leq R+\nu A_{\text{max}}$ for all $t\in\mathbb{R}^{+}$) by
contradiction. Since $t\mapsto\phi_{m}^{t}(A(0))$ is continuous, assume
otherwise the existence of $\delta,t_{\delta}\in\mathbb{R}^{+}_{\ast}$ such
that

$A(t_{\delta})=A_{\text{max}}+\delta\quad\text{and}\quad
A(t)<A_{\text{max}}+\delta\ \text{for}\ t\in[-\tau,t_{\delta}).$

Then, we certainly have $m(t)<R+\nu(A_{\text{max}}+\delta)$ for
$t\in[0,t_{\delta})$. Using expression (6), the definition of $A_{\text{max}}$
and the condition $\nu<\beta$ successively imply

$A(t_{\delta})\leq\left(A_{\text{max}}+(R+\nu(A_{\text{max}}+\delta))\frac{e^{\beta
t_{\delta}}-1}{\beta}\right)e^{-\beta
t_{\delta}}=A_{\text{max}}+\frac{\nu}{\beta}\delta(1-e^{-\beta
t_{\delta}})<A_{\text{max}}+\delta$

which is impossible.

Considering now an arbitrary trajectory, we first prove that for every
$t_{\ast}>\text{max}_{\mathbf{x}^{0}}$ such that

$A(t_{\ast})=\max_{t\in[-\tau,t_{\ast}]}A(t),$

we must have $A(t_{\ast})\leq A_{\text{max}}$. By contradiction, the fact that
all oscillators must have been reset at least once before time $t_{\ast}$ and
the definition of $t_{\ast}$ imply that we must have $m(t_{\ast})\leq R+\nu
A(t_{\ast})$ and hence

$\dot{A}(t_{\ast})\leq R+\nu A(t_{\ast})-\beta A(t_{\ast})<0$

where the second inequality follows from $A(t_{\ast})>A_{\text{max}}$.
Moreover, the derivative $\dot{A}(t)$ is left continuous. Hence $A$ must be
decreasing in the left neighbourhood of $t_{\ast}$, which is impossible from
the definition of $t_{\ast}$. Therefore, we must have
$\sup_{t\in\mathbb{R}^{+}}A(t)<+\infty$ in every trajectory.

In order to prove that $\limsup_{t\to+\infty}A(t)\leq A_{\text{max}}$, we use
a bootstrap argument. Let $\rho\in(0,1)$ be sufficiently large so that
$\nu<\rho\beta$. We claim that for every $\delta>0$ such that

$A(t)\leq(1+\delta)A_{\text{max}},\ \forall t\in[-\tau,+\infty),$ (10)

there exists $t_{\delta}\in\mathbb{R}^{+}$ such that
$A(t)\leq(1+\rho\delta)A_{\text{max}}$ for all $t\in[t_{\delta},+\infty)$.
Indeed, assume firstly that $A(t)>(1+\rho\delta)A_{\text{max}}$ for all
$t\in\mathbb{R}$. Then we would have

$\dot{A}(t)<-\beta(1+\rho\delta)A_{\text{max}}+R+\nu(1+\delta)A_{\text{max}}=(\nu-\rho\beta)A_{\text{max}}\delta<0,\
\forall t\in\mathbb{R}^{+},$

ie. $A(t)$ would have to decrease at least linearly. Given the inequality
(10), it would be impossible that it remains above
$(1+\rho\delta)A_{\text{max}}$ forever. Moreover, a similar reasoning implies
that if $A(t_{\delta})=(1+\rho\delta)A_{\text{max}}$, then we must have
$A(t)\leq A(t_{\delta})$ for all $t>t_{\delta}$. ∎

####  2.2.3 Preservation of the order in which the oscillators fire

Since all degradation rates are equal, in every trajectory, the initial
ordering of the repressor concentrations is preserved until the first firing.
Under the assumption

$\nu<\frac{\beta}{1+\beta R},$ (11)

we are going to show that, when inside $Q$, between any two consecutive
firings of a given oscillator, every reset concentration must lie above the
current concentration of that oscillator. By induction, this implies that,
after the last of the first firing times of each oscillator, the order of the
repressor concentrations is cyclically permuted at each firing (and evidently
remains constant in time between firings), implying that the order in which
the oscillators fire is preserved forever. That property will make the
analysis of the dynamics simpler.

Given $i\in\\{1,\cdots,N\\}$ and $k\in\mathbb{N}$, let $t_{i}^{k}$ be the
instant of the $k^{\text{th}}$ firing of oscillator $i$.

######  Claim 2.3.

Assume that inequality (11) holds and consider a trajectory for which
$(\mathbf{x}(t),A(t))\in Q$ for all $t\in\mathbb{R}^{+}$. Let
$i\in\\{1,\cdots,N\\}$ and $k\in\mathbb{N}$ be arbitrary. If an oscillator
$j\in\\{1,\cdots,N\\}$ fires between the $k$th and $(k+1)$th firings of $i$
(ie. if there exists $\ell\in\mathbb{N}$ such that
$t^{\ell}_{j}\in(t^{k}_{i},t^{k+1}_{i})$), then we have

$x_{i}\left((t^{\ell}_{j})^{+}\right)<x_{j}\left((t^{\ell}_{j})^{+}\right).$

###### Proof.

If $(\mathbf{x}(t),A(t))\in Q$ then we have $|\dot{A}(t)|\leq\beta
A_{\text{max}}$. The assumption (11) then implies $|\dot{A}(t)|<\frac{1}{\nu}$
which in turn yields

$R+\nu A(t_{1}-\tau)-(t_{2}-t_{1})<R+\nu A(t_{2}-\tau),\ \forall t_{2}>t_{1},$

from where the conclusion is immediate. ∎

####  2.2.4 Return map

In every trajectory, firings must occur infinitely often and their consecutive
occurrences are separated by positive time intervals. Moreover, the order
preservation obtained in the previous section implies that (for
$t>\max_{i}t_{i}^{1}=\text{max}_{\mathbf{x}^{0}}$) between any two consecutive
firings of oscillator $i$, all other oscillators having repressor
concentration distinct from $x_{i}$ must fire exactly once. Accordingly, in
order to analyse the dynamics, it suffices to study the iterations of a return
map that acts on data immediately before the firing of a given oscillator, say
oscillator $N$.

An expression of the return map can be computed as follows. Given
$\mathbf{x}^{0}$ with $x^{0}_{N}=0$, for convenience in the sequel, we denote
by $t_{\text{R}}$, the time $t_{N}^{2}=x_{N}(0^{+})=R+\nu A(-\tau)$ of the
second firing of $x_{N}$. We also assume that
$\text{max}_{\mathbf{x}^{0}}<x_{N}(0^{+})$ so that the order in which the
oscillators fire is preserved from $t=0$. Accordingly, for any oscillator for
which $x^{0}_{i}>0$, the corresponding repressor concentration between
$t_{i}^{1}=x^{0}_{i}$ and $t_{\text{R}}$ is given by

$x_{i}(t)=R+\nu A(x^{0}_{i}-\tau)-(t-x^{0}_{i}),\ \forall
t\in(x^{0}_{i},t_{\text{R}}].$

Hence, we have

$x_{i}(t_{\text{R}})=x^{0}_{i}+\nu\left(A(x^{0}_{i}-\tau)-A(-\tau)\right),\
\forall i\in\\{1,\cdots,N-1\\}.$

In particular, if $(\mathbf{x},A)\in Q$ lies in the Poincaré section $x_{N}=0$
and satisfies the conditions of Claim 2.1, then the return map $F_{N}$ writes
$(\mathbf{x}^{\prime},A^{\prime})=F_{N}(\mathbf{x},A)$ where

$\left\\{\begin{array}[]{l}x^{\prime}_{i}=x_{i}+\nu(\phi_{m}^{x_{i}-\tau}(A)-\phi_{m}^{-\tau}(A))\
\text{for}\ i\in\\{1,\cdots,N-1\\}\\\
A^{\prime}=\phi^{R+\nu\phi^{-\tau}_{m}(A)}_{m}(A)\end{array}\right.$ (12)

In the stability analyis of periodic orbits (either fully synchronized or only
partially synchronized) in the various sections below, we shall ensure that
all iterates $(\mathbf{x}^{k},A^{k}):=F_{N}^{k}(\mathbf{x},A)$ of sufficiently
small initial perturbations satisfy the conditions of Claim 2.1 (or even
condition (8) in the case of the full synchronized orbit), so that the
stability analysis actually reduces to the study of their (linearized)
dynamics in $(\mathbb{R}^{+})^{N}$.

Of note, Appendix A states and proves a certain property of Lipschitz-
continuous dependence of the return map on its input datum, not only when the
map reduces to one of $(\mathbb{R}^{+})^{N}$ but also in the case of an
arbitrary datum $(\mathbf{x},A|_{[-\tau,0]})$. This property will be employed
in the proof of a stability criterion for partially synchronized periodic
orbits in Section 4.

##  3 Existence and stability analysis of the synchronized periodic orbit

As a system with mean-field interactions, the equations (3)-(5) commute with
every permutation of the repressor indices. This suggests to study the
synchronized dynamics inside the invariant subspace
$x_{1}=x_{2}=\cdots=x_{N}$. The synchronized dynamics of the population of $N$
oscillators reduces to that of the $N=1$ system. In this section, we first
investigate the corresponding return map $F_{1}$, a one-dimensional map that
acts on the variable $A$ (since we always have $x_{1}=0$ along every orbit).
We show that all orbits asymptotically converge to a unique fixed point
$A_{\text{FP}}$. Then, we study the Lyapunov stability in
$(\mathbb{R}^{+})^{N}$ for $N\geq 2$, of the fixed point
$(0,\cdots,0,A_{\text{FP}})$ of the return map $F_{N}$. We show that this
stability depends on whether the delay-parameter $\tau$ vanishes or it is
positive; however and remarkably, it is independent of $N$. The results are
collected in the following statement. Letting

$R_{\tau}=R+\nu\frac{1+\beta\tau-e^{\beta\tau}}{\beta^{2}}\quad\text{and}\quad\nu_{\tau}=\nu
e^{\beta\tau},$ (13)

we will assume the following conditions on the parameters

$R_{\tau}>0\quad\text{and}\quad\nu_{\tau}<\beta.$ (14)

######  Proposition 3.1.

Assume that (14) holds. Then the one-dimensional return map $F_{1}$ has a
unique globally attracting fixed point $A_{\text{FP}}$.

Assume that (11) also holds.

(i) If $\tau=0$, then for every $N\geq 2$, the corresponding fixed point
$(0,\cdots,0,A_{\text{FP}})$ of $F_{N}$ is unstable in $(\mathbb{R}^{+})^{N}$.

(ii) There exists $\tau_{0}\in\mathbb{R}^{+}_{\ast}$ such that for every
$\tau\in(0,\tau_{0})$ and every $N\geq 2$, the fixed point
$(0,\cdots,0,A_{\text{FP}})$ of $F_{N}$ is locally asymptotically stable in
$(\mathbb{R}^{+})^{N}$.

This statement confirms the numerical observations reported in [18] about the
delay-dependence of the stability of the synchronized periodic orbit
associated with $A_{\text{FP}}$.444In [18], the repressor resets $R+\nu
A(t-\tau)$ are perturbed by a (small) additive random noise. Our result shows
that the stability of the synchronized periodic orbit is not affected by these
random fluctuations. An illustration is given in Fig. 1.

Figure 1: Time series of two trajectories for $N=10,t\in[0,20]$ and $\tau=0$
(left)/$\tau=0.2$ (right). The other parameters are $R=2,\beta=1$, and
$\nu=0.2$. The sawtooth series correspond to the repressor concentrations
$x_{i}(t)$ ($i\in\\{1,\cdots,10\\}$) and the central series in red color
corresponds to the activator concentration $A(t)$. On the left picture, the
instability of the synchronized periodic orbit is marked although rather weak.
On the right picture, the asymptotic stability is more evident.

The rest of this section is devoted to the proof of Proposition 3.1.

###  3.1 Existence of a globally attracting fixed point of the map $F_{1}$

Figure 2: _Left._ Illustration of the vector field (light blue arrows) and
some segments of trajectories (solid black curves) of the synchronized
dynamics/system with one oscillator ($N=1$), under the assumptions (11) and
(14). Initial conditions: $(A_{0}(\tau),R_{\tau}+\nu_{\tau}A_{0}(\tau))$,
$(A_{\text{FP}},R_{\tau}+\nu_{\tau}A_{\text{FP}})$, $(A_{\tau},R_{\tau}+\nu
A_{\tau})$; where the orbit of
$(A_{\text{FP}},R_{\tau}+\nu_{\tau}A_{\text{FP}})$ is the unique attracting
periodic orbit. _Right._ Graph of the corresponding return map $F_{1}$.

For $N=1$, the condition (8) reduces to

$A\geq A_{0}(\tau)=\frac{1-(1+\beta\tau)e^{-\beta\tau}}{\beta^{2}}.$

Moreover, using (7) with $m(-s)=s$ for $s\in[0,\tau]$, we obtain that the
reset concentration at any firing is given by (NB: see (13) for the definition
of $R_{\tau}$ and $\nu_{\tau}$ and see Fig. 2 for an illustration)

$R+\nu\phi^{-\tau}_{m}(A)=R_{\tau}+\nu_{\tau}A,$

Therefore, using (6) with $m(s)=R_{\tau}+\nu_{\tau}A-s$ for all
$s\in(0,R_{\tau}+\nu_{\tau}A]$ yields the following explicit expression for
the return map $F_{1}(A)=\phi^{R+\nu\phi^{-\tau}_{m}(A)}_{m}(A)$

$F_{1}(A)=\frac{1}{\beta^{2}}+\left(A(1-\frac{\nu_{\tau}}{\beta})-\frac{R_{\tau}}{\beta}-\frac{1}{\beta^{2}}\right)e^{-\beta(R_{\tau}+\nu_{\tau}A)}.$

Letting $A_{\tau}:=\frac{R_{\tau}}{\beta-\nu_{\tau}}$, the assumption (14)
implies $A_{0}(\tau)<A_{\tau}$. Clearly, this assumption also implies that the
asymptotic dynamics of $F_{1}$ must lie in the interval
$[A_{0}(\tau),A_{\tau}]$, again see Fig. 2 for illustration. The following
features then yield that $F_{1}$ must have a unique globally attracting fixed
point inside this interval, proving the preliminary claim in Proposition 3.1:

  * •

$F_{1}(A_{0}(\tau))=\frac{1}{\beta^{2}}-\left(\frac{(1+\beta\tau)}{\beta^{2}}e^{-\beta\tau}+\frac{R}{\beta}\right)e^{-\beta
R}>A_{0}(\tau)$ for every $\tau\in[0,R)$.

  * •

$F_{1}(A_{\tau})=\frac{1}{\beta^{2}}\left(1-e^{-\beta^{2}A_{\tau}}\right)<A_{\tau}$,
so that altogether
$F_{1}([A_{0}(\tau),A_{\tau}])\subsetneq[A_{0}(\tau),A_{\tau}]$.

  * •

$F_{1}^{\prime}(A_{\tau})>0$ and $F_{1}^{\prime\prime}|_{[0,A_{\tau}]}<0$,
viz. $F_{1}|_{[0,A_{\tau}]}$ is increasing with decreasing derivative.

###  3.2 Proof of instability for $\tau=0$

In this section, we prove item (i) in Proposition 3.1. Consider the
restriction of $F_{N}$ to the subspace of 2-cluster states for which $N-1$
repressor concentrations are equal, more precisely, the restriction to those
$\mathbf{x}$ for which $x_{1}>0$ and $x_{2}=\cdots=x_{N}=0$. We are going to
show that $(0,\cdots,0,A_{\text{FP}})$ is linearly repelling along the
direction of $x_{1}$.

Letting $\epsilon>0$ sufficiently small, assume that

$x_{1}\in(0,\epsilon)\ \text{and}\ |A-A_{\text{FP}}|<\epsilon.$

In particular, $\epsilon$ must be small enough so that such $(\mathbf{x},A)$
satisfies the condition (8).

Given that $x_{2}=\cdots=x_{N}=0$, we have for the continuous time trajectory

$m(t)=\frac{(N-1)(R+\nu A)+x_{1}}{N}-t,\ \forall t\in(0,x_{1}],$

which yields using (6) and after simple algebra

$\phi^{x_{1}}_{m}(A)=A+x_{1}\left(\frac{N-1}{N}(R+\nu A)-\beta A\right)+{\cal
O}(x_{1}^{2}).$

From (12), we obtain the following expansion for the coordinate
$x^{\prime}_{1}$ of the first iterate
$(\mathbf{x}^{\prime},A^{\prime})=F_{N}(\mathbf{x},A)$

$x^{\prime}_{1}=x_{1}+\nu(\phi^{x_{1}}_{m}(A)-A)=x_{1}\left(1+\nu\left(\frac{N-1}{N}(R+\nu
A)-\beta A\right)\right)+{\cal O}(x_{1}^{2}).$

Therefore, in order to prove the desired instability, all we have to show is
$\frac{N-1}{N}(R+\nu A)-\beta A>0$ for all $|A-A_{\text{FP}}|<\epsilon$ with
$\epsilon$ small. By continuity, it suffices to show that $\frac{N-1}{N}(R+\nu
A_{\text{FP}})-\beta A_{\text{FP}}>0$. Since $\frac{N-1}{N}\geq\frac{1}{2}$
for all $N\geq 2$, it suffices to show that this inequality holds for $N=2$,
viz.

$A_{\text{FP}}<\frac{R}{2\beta-\nu}.$

Given the properties of $F_{1}$ described in the previous subsection, in order
to prove that inequality, it suffices to verify that
$F_{1}(\frac{R}{2\beta-\nu})<\frac{R}{2\beta-\nu}$. Explicit computations show
that the sign of $F_{1}(\frac{R}{2\beta-\nu})-\frac{R}{2\beta-\nu}$ is the
same as the one of
$\tanh\left(\frac{\beta^{2}R}{2\beta-\nu}\right)-\frac{\beta^{2}R}{2\beta-\nu}$,
which is negative for every $R,\beta\in\mathbb{R}^{+}$ and $\nu\in(0,\beta)$.
This concludes the proof of instability for $\tau=0$. $\Box$

###  3.3 Proof of stability for $\tau>0$

We first give some heuristic for the synchronization mechanism when $\tau>0$.
Given a small initial perturbation $(\mathbf{x},A)$ of the fixed point, the
vector field in (5), when computed at instants immediately before the first
firing (which occurs at time 0) must be negative because $m(t)$ must be close
to $0$ and $A(t)$ must be close to $A_{\text{FP}}$, which is positive.
Therefore for $\tau>0$ small enough and $x_{i}\in(0,\tau)$, since there cannot
be any firing in the interval $[-\tau,x_{i}-\tau]$, we must have

$\phi_{m}^{x_{i}-\tau}(A)=\phi_{m}^{-\tau}(A)-C^{\prime}x_{i}+\
\text{h.o.t.},$

for some $C^{\prime}\in\mathbb{R}^{+}_{\ast}$, which implies that the image
$x^{\prime}_{i}$ of $x_{i}$ under $F_{N}$ is given by

$x^{\prime}_{i}=x_{i}+\nu(\phi_{m}^{x_{i}-\tau}(A)-\phi_{m}^{-\tau}(A))=x_{i}(1-\nu
C^{\prime})+\ \text{h.o.t.}$

Thus, we must have $x^{\prime}_{i}<x_{i}$ when $(\mathbf{x},A)$ is
sufficiently close to $(0,\cdots,0,A_{\text{FP}})$, which together with the
fact that $x^{\prime}_{i}>0$ (from order preservation), implies the desired
synchronisation.

##### Proof of item (ii) in Proposition 3.1

We are going to show that for every $N\geq 2$, the map $F_{N}$ is a
contraction in a small neighbourhood of $(0,\cdots,0,A_{\text{FP}})$ in
$(\mathbb{R}^{+})^{N}$ when $\tau>0$ is small enough. Given $\tau\in(0,R)$,
let $(\mathbf{x},A)$ with $x_{N}=0$, and
$\|\mathbf{x}\|_{N-1},|A-A_{\text{FP}}|$ small. In particular, we assume that
$\|\mathbf{x}\|_{N-1}<\tau$ and that $(\mathbf{x},A)$ satisfies the condition
(8).

Let $i\in\\{1,\cdots,N-1\\}$. The condition (8) ensures that no firing occurs
in the time interval $[-\tau,0)$; hence

$m(t)=m(0)-t\quad\text{for}\quad t\in[-\tau,x_{i}-\tau],$

from where we obtain using (6), after simple algebra

$\phi^{x_{i}-\tau}_{m}(A)=\left(\phi^{-\tau}_{m}(A)+\int_{0}^{x_{i}}e^{\beta
s}m(s-\tau)ds\right)e^{-\beta
x_{i}}=\phi^{-\tau}_{m}(A)-x_{i}\left(\beta\phi^{-\tau}_{m}(A)-\tau\right)+{\cal
O}(\|\mathbf{x}\|_{N-1}^{2})$

and thus, the following expansion results for the coordinate $x^{\prime}_{i}$
of the first iterate $(\mathbf{x}^{\prime},A^{\prime})=F_{N}(\mathbf{x},A)$

$\displaystyle x^{\prime}_{i}$
$\displaystyle=x_{i}\left(1-\nu\left(\beta\phi^{-\tau}_{m}(A)-\tau\right)\right)+{\cal
O}(\|\mathbf{x}\|_{N-1}^{2})$ $\displaystyle=x_{i}(1-\nu\beta
A_{\text{FP}}+(\beta(A-\phi^{-\tau}_{m}(A))-\tau))+\nu\beta
x_{i}(A_{\text{FP}}-A)+{\cal O}(\|\mathbf{x}\|_{N-1}^{2})$

That no firing occurs in the interval $[-\tau,0)$ implies that
$\phi^{-\tau}_{m}$ is close to the identity when $\tau$ is small. Therefore,
there exists $K_{\tau}>0$ with $\lim_{\tau\to 0}K_{\tau}=0$ such that

$|\beta(A-\phi^{-\tau}_{m}(A))-\tau|\leq K_{\tau}$

for all $A$ uniformly bounded, and in particular when $|A-A_{\text{FP}}|$ is
small. In addition, order preservation, which must hold for the initial
conditions $(\mathbf{x},A)$ under consideration here, implies that
$x^{\prime}_{i}\geq 0$ for all $i$. Letting $\tau\to 0$ and $A\to
A_{\text{FP}}$, we must have $1-\nu\beta A_{\text{FP}}\geq 0$. Altogether,
this implies the existence of $\gamma_{1}\in(0,1)$ such that

$\|\mathbf{x}^{\prime}\|_{N-1}\leq\left(\gamma_{1}+\nu\beta|A-A_{\text{FP}}|\right)\|\mathbf{x}\|_{N-1},$

when $\tau,\|\mathbf{x}\|_{N-1}$ and $|A-A_{\text{FP}}|$ are sufficiently
small.

In order to control the dynamics of the activator variable, we are going to
show that $A^{\prime}$ is close to $F_{1}(A)$ when $\|\mathbf{x}\|_{N-1}$ is
small and then use that $F_{1}$ is a contraction in the neighbourhood of
$A_{\text{FP}}$ in $\mathbb{R}$. Recalling the notation
$t_{\text{R}}=R+\nu\phi_{m}^{-\tau}(A)$ for the return time and letting
$t_{\text{R}}^{\text{sync}}=R+\nu\phi_{m_{\text{sync}}}^{-\tau}(A)$ and
$m_{\text{sync}}$ for the quantities associated with the synchronized
trajectory issued from $(0,A)$ at $t=0$, we observe that together with (7),
the expressions

$m(t)=\frac{1}{N}\sum_{i=1}^{N-1}x_{i}-t\quad\text{and}\quad
m_{\text{sync}}(t)=-t\quad\text{for}\ t\in[-\tau,0],$

result in

$t_{\text{R}}-t_{\text{R}}^{\text{sync}}=-\nu\frac{e^{\beta\tau}-1}{\beta
N}\sum_{i=1}^{N-1}x_{i}<0.$

Accordingly, we obtain from (6)

$\displaystyle A^{\prime}=$ $\displaystyle
F_{1}(A)+\left(A+\int_{0}^{t_{\text{R}}}e^{\beta
s}m_{\text{sync}}(s)ds\right)(e^{-\beta t_{\text{R}}}-e^{-\beta
t_{\text{R}}^{\text{sync}}})+\int_{0}^{t_{\text{R}}}e^{\beta(s-t_{\text{R}})}(m(s)-m_{\text{sync}}(s))ds$
$\displaystyle+\int_{t_{\text{R}}}^{t_{\text{R}}^{\text{sync}}}e^{\beta(s-t_{\text{R}}^{\text{sync}})}m_{\text{sync}}(s)ds$

The expression of $t_{\text{R}}-t_{\text{R}}^{\text{sync}}$ above implies that
the second and last terms in the RHS are controlled by $\|\mathbf{x}\|_{N-1}$,
also because all quantities $A,t_{\text{R}},t_{\text{R}}^{\text{sync}}$ and
$m_{\text{sync}}$ are bounded when in (or close to) the attracting set.

In order to control the third term, we first observe that it suffices to
provide an estimate of the integral between time $t=\|x\|_{N-1}$ (namely the
time of the last first firing of all oscillators) and $t=t_{\text{R}}$ because
the remaining integral can be controlled by the same arguments as before.
Moreover, we have

$m(s)=R+\frac{\nu}{N}\sum_{i=1}^{N}\phi_{m}^{x_{i}-\tau}(A)+\frac{1}{N}\sum_{i=1}^{N-1}x_{i}-s\quad\text{and}\quad
m_{\text{sync}}(s)=R+\nu\phi_{m_{\text{sync}}}^{-\tau}(A)-s\quad\text{for}\
s\in[\|x\|_{N-1},t_{\text{R}}]$

and again from (6)

$|\phi_{m}^{x_{i}-\tau}(A)-\phi_{m}^{-\tau}(A)|\leq\phi_{m}^{-\tau}(A)(1-e^{-\beta
x_{i}})+\|m|_{[-\tau,x_{i}-\tau]}\|_{0}\frac{e^{\beta x_{i}}-1}{\beta}.$

In addition, from the expressions of $m$ and $m_{\text{sync}}$ above, we have

$\phi_{m}^{-\tau}(A)=\phi_{m_{\text{sync}}}^{-\tau}(A)-\frac{e^{\beta\tau}-1}{\beta
N}\sum_{i=1}^{N-1}x_{i}.$

Combining all the estimates above, we finally conclude about the existence of
$K\in\mathbb{R}^{+}$ such that

$|A^{\prime}-A_{\text{FP}}|\leq|F_{1}(A)-A_{\text{FP}}|+K\|\mathbf{x}\|_{N-1}+{\cal
O}(\|\mathbf{x}\|_{N-1}^{2})$

when $\|\mathbf{x}\|_{N-1}$ is sufficiently small. Besides, the analysis of
the map $F_{1}$ in a previous section implies the existence of
$\gamma_{2}\in(0,1)$ such that

$|F_{1}(A)-A_{\text{FP}}|\leq\gamma_{2}|A-A_{\text{FP}}|$

when $A$ is sufficiently close to $A_{\text{FP}}$.

Altogether, when $\tau,\|\mathbf{x}\|_{N-1}$ and $|A-A_{\text{FP}}|$ are
sufficiently small, the return dynamics of $\|\mathbf{x}\|_{N-1}$ and
$|A-A_{\text{FP}}|$ is dominated by a matrix whose eigenvalues are close to
$\gamma_{1}$ and $\gamma_{2}$. Hence, the map $F_{N}$ must be a contraction
(for an appropriate norm in $\mathbb{R}^{N}$) in the neighbourhood of
$(0,\cdots,0,A_{\text{FP}})$, as desired. In particular, we are sure that
$(\mathbf{x}^{\prime},A^{\prime})$ also satisfies condition (8) provided that
$\|\mathbf{x}\|_{N-1},|A-A_{\text{FP}}|$ are sufficiently small; hence the
argument can be repeated to conclude about asymptotic stability of
$(0,\cdots,0,A_{\text{FP}})$. $\Box$

##  4 Stability analysis for partially synchronized periodic orbits

The analysis in the previous section and its arguments are not limited to
synchronized trajectories. They extend to partially synchronized periodic
orbits. A partially synchronized trajectory (whether it is periodic or not) is
a solution of the equations (3)-(5) for which two or more repressor
concentrations are equal at all times (NB: this is the case iff the
concentrations are equal at $t=0$).

For an arbitrary fixed point of the return map in a given partially
synchronized subspace (see next section for an accurate definition), we
establish a stability criterion for $\tau>0$ (Lemma 4.2) and an instability
criterion for $\tau=0$ (Lemma 4.3). For simplicity, we only consider stability
with respect to initial perturbations that affect a single cluster of the
periodic orbit under investigation. The stability with respect to
perturbations that affect several clusters are direct extensions that are left
to the interested reader. In addition, we also consider initial conditions
that are given by the datum $(\mathbf{x},A|_{[-\tau,0]})$, even though we have
argued that only finite many values of $A$ suffice in order to define any
trajectory.

In a second step, we apply these criteria to fixed points that are
exponentially stable in their own partially synchrony subspace for $\tau=0$.
The result for the continued fixed point for $\tau>0$ small can be regarded as
an extension of Proposition 3.1 to the family under consideration: while for
$\tau=0$, the orbit is unstable with respect to perturbations that smear its
clusters, it becomes stable against the same perturbations for every positive
(and small) value of $\tau$ (Theorem 4.4). An example of application to fixed
points with two clusters with equi-distributed repressor concentrations is
given in Lemma 4.5.

###  4.1 Characteristics of partially synchronized trajectories and of fixed
points of the corresponding return map

By grouping the oscillators with equal repressor concentration into one
cluster, the population at every instant can be described by the vector
$\left(\\{n_{k},y_{k}\\}_{k=1}^{K},A\right)$ where $n_{k}\in\\{1,\cdots,N\\}$
denotes the size of the cluster $k$ and $y_{k}$ the corresponding repressor
concentration ($K\leq N$ is the total number of clusters and we have
$\sum_{k=1}^{K}n_{k}=N$).555The variables $y_{k}$ can be formally defined as
follows. For each $i\in\\{1,\cdots,N\\}$, there exists
$k_{i}\in\\{1,\cdots,K\\}$ such that $y_{k_{i}}=x_{i}$ and $x_{j}\neq x_{i}$
iff $y_{k_{j}}\neq y_{k_{i}}$. From the equations (3)-(5), it follows that the
cluster distribution $\\{n_{k}\\}_{k=1}^{K}$ remains constant in every
trajectory; hence we may consider separately the dynamics in each subspace,
called partially synchronized subspace, for which this distribution is given.
Every trajectory for which $K<N$ is called partially synchronized. Obviously,
the synchronized dynamics in Section 3 is a particular case of partially
synchronized subspace ($K=1$).

Let $N\in\mathbb{N}$, $K\in\\{1,\cdots,N-1\\}$ and a cluster distribution
$\\{n_{k}\\}_{k=1}^{K}$ be given. For the sake of the presentation, we only
consider those partially synchronized periodic orbits that are fixed points

$\left(\\{n_{k},y^{\text{FP}}_{k}\\}_{k=1}^{K},A^{\text{FP}}|_{[-\tau,0]}\right)$

of the return map to the Poincaré section $y_{K}=0$. In other words, we assume
that the period $T_{\text{FP}}>0$ coincides with the return time. However the
analysis developed below extends to arbitrary periodic orbits of such map,
without additional conceptual difficulties.

Notice also that the labelling of the clusters is irrelevant because of the
permutation symmetry. For convenience, we assume that this labelling has been
chosen so that the fixed point repressor concentrations $y_{k}^{\text{FP}}$
are ordered ie.

$y_{1}^{\text{FP}}>y_{2}^{\text{FP}}>\cdots>y_{K-1}^{\text{FP}}>0.$

Furthermore, the expression (6) and the fact that the non-negative function
$m^{\text{FP}}$ certainly does not entirely vanish over $[0,T_{\text{FP}}]$,
impose that, in any partially synchronized periodic orbit, the activator
function $A^{\text{FP}}$ must be positive.

###  4.2 Considerations about stability

As technical considerations about stability are concerned, focus will be made
on the strongest form of local Lyapunov stability, namely the exponential
stability that results when the return dynamics of small perturbations is a
contraction in an appropriate setting (space and norm). In particular, a
prerequisite for our stability criterion will be that the periodic orbit is
exponentially stable inside its proper partially synchronized space. The
criterion will then ensure exponential stability in a higher dimensional space
that contains perturbations that smear the clusters. More precisely, we shall
deal with the following notions (NB: Throughout, the symbol $\mathbf{y}$
denotes the collection $\\{y_{k}\\}_{k=1}^{K}$).

######  Definition 4.1.

(i) The fixed point
$\left(\\{n_{k},y^{\text{FP}}_{k}\\}_{k=1}^{K},A^{\text{FP}}|_{[-\tau,0]}\right)$
is said to be exponentially stable inside its proper partially synchronized
subspace if there exist $\epsilon_{0}\in\mathbb{R}_{\ast}^{+}$ and
$\gamma\in(0,1)$ such that for every $\epsilon\in(0,\epsilon_{0})$ and every
initial datum $\left(\\{n_{k},y_{k}\\}_{k=1}^{K},A|_{[-\tau,0]}\right)$ such
that

$y_{K}=0\quad\text{and}\quad\max\left\\{\|\mathbf{y}-\mathbf{y}^{\text{FP}}\|_{K-1},\|A|_{[-\tau,0]}-A^{\text{FP}}|_{[-\tau,0]}\|_{0}\right\\}\leq\epsilon,$

we have for the subsequent trajectory
$t\mapsto\left(\\{n_{k},y_{k}(t)\\}_{k=1}^{K},A(t)\right)$ of the system
(3)-(5)

$\max\left\\{\|\mathbf{y}(t_{\text{R}})-\mathbf{y}^{\text{FP}}\|_{K-1},\|A|_{[t_{\text{R}}-\tau,t_{\text{R}}]}-A^{\text{FP}}|_{[-\tau,0]}\|_{0}\right\\}\leq\gamma\epsilon,$

where the return time $t_{\text{R}}$ defined in Section 2.2.4 corresponds here
to the instant of the second firing of the oscillators in cluster $K$.

(ii) The fixed point
$\left(\\{n_{k},y^{\text{FP}}_{k}\\}_{k=1}^{K},A^{\text{FP}}|_{[-\tau,0]}\right)$
is said to be exponentially stable with respect to small perturbations that
smear cluster $K$ if there exist
$\epsilon_{0},C_{1},C_{2}\in\mathbb{R}_{\ast}^{+}$ and $\gamma\in(0,1)$ such
that for every $\epsilon\in(0,\epsilon_{0})$, every
$K^{\prime}\in\\{K+1,\cdots,N\\}$, every cluster distribution
$\\{n^{\prime}_{k}\\}_{k=1}^{K^{\prime}}$ such that $n^{\prime}_{k}=n_{k}$ for
$k\in\\{1,\cdots,K-1\\}$, and every initial datum
$\left(\\{n^{\prime}_{k},y_{k}\\}_{k=1}^{K^{\prime}},A|_{[-\tau,0]}\right)$
such that $y_{K^{\prime}}=0$ and

$\max\left\\{\|\mathbf{y}-\mathbf{y}^{\text{FP}}\|_{K-1},C_{1}\max_{k\in\\{K,\cdots,K^{\prime}-1\\}}y_{k},\|A|_{[-\tau,0]}-A^{\text{FP}}|_{[-\tau,0]}\|_{0}\right\\}\leq\epsilon,$

we have for the subsequent trajectory
$t\mapsto\left(\\{n^{\prime}_{k},y_{k}(t)\\}_{k=1}^{K^{\prime}},A(t)\right)$
of the system (3)-(5)

$\max\left\\{\|\mathbf{y}(t^{n}_{\text{R}})-\mathbf{y}^{\text{FP}}\|_{K-1},C_{1}\max_{k\in\\{K,\cdots,K^{\prime}-1\\}}y_{k}(t^{n}_{\text{R}}),\|A|_{[t^{n}_{\text{R}}-\tau,t^{n}_{\text{R}}]}-A^{\text{FP}}|_{[-\tau,0]}\|_{0}\right\\}\leq
C_{2}\gamma^{n}\epsilon,\ \forall n\in\mathbb{N}$

where the $n$th return time $t^{n}_{\text{R}}$ is the instant of the $(n+1)$th
firing of the oscillators in cluster $K^{\prime}$.

Anticipating the comment after Lemma 4.2, notice that exponential stability
inside the proper partially synchronized subspace does not depend on the phase
of the periodic orbit under consideration, viz. the fixed point
$\left(\\{n_{k},y^{\text{FP}}_{k}\\}_{k=1}^{K},A^{\text{FP}}|_{[-\tau,0]}\right)$
is exponentially stable inside its proper partially synchronized subspace iff
the fixed point

$\left(\\{n^{\prime}_{k},y^{\prime}_{k}\\}_{k=1}^{K},A^{\text{FP}}|_{[y^{\text{FP}}_{K-1}-\tau,y^{\text{FP}}_{K-1}]}\right),$

where

$n^{\prime}_{k}=\left\\{\begin{array}[]{ccl}n_{K}&\text{if}&k=1\\\
n_{k-1}&\text{if}&k\in\\{2,\cdots,K\\}\end{array}\right.\quad\text{and}\quad
y^{\prime}_{k}=\left\\{\begin{array}[]{ccl}R+\nu
A(y^{\text{FP}}_{K-1}-\tau)&\text{if}&k=1\\\
y^{\text{FP}}_{k-1}-y^{\text{FP}}_{K-1}&\text{if}&k\in\\{2,\cdots,K\\}\end{array}\right.,$

is also exponentially stable in the same sense. This is a standard consequence
of the continuity of the firing dynamics (ie. the fact that the coordinates
immediately after firing depend continuously on the coordinate immediately
after the previous firing).

###  4.3 Stability criterion for $\tau>0$

The stability criterion for partially synchronized periodic orbits, which is
given in the next statement, is reminiscent of the heuristic argument for the
stability of the fully synchronized periodic orbit given at the beginning of
section 3.3.

######  Lemma 4.2.

Assume that (11) holds. Given $N\in\mathbb{N}$, $K\in\\{1,\cdots,N-1\\}$, a
cluster distribution $\\{n_{k}\\}_{k=1}^{K}$ and $\tau>0$, assume that the
return map in the Poincaré section $y_{K}=0$ in the corresponding partially
synchronized subspace has a fixed point
$\left(\\{n_{k},y^{\text{FP}}_{k}\\}_{k=1}^{K},A^{\text{FP}}|_{[-\tau,0]}\right)$.
Assume also that the following conditions hold

  * •

the fixed point is asymptotically stable inside its proper partially
synchronized subspace,

  * •

the periodic orbit of (3)-(5) that passes through the fixed point has no
firing in any of the time intervals $[-\tau,0]\ \text{mod}\ T_{\text{FP}}$,
and the derivative of $A$ at $t=-\tau$ is negative, ie.
$\dot{A}^{\text{FP}}(-\tau)=m_{\text{FP}}(-\tau)-\beta
A^{\text{FP}}(-\tau)<0$.

Then, the fixed point is exponentially stable with respect to small
perturbations that smear cluster $K$.

That the second condition in this statement depends on the value of the
derivative $\dot{A}$ at a certain time suggests that stability with respect to
cluster smearing a priori depends on the cluster under consideration. However,
we shall see in the proof of Theorem 4.4 below that, provided that $\tau$ is
sufficiently small, this is not the case because the sign of the derivative
involved actually does not depend on the cluster under consideration.

Notice also that in the special case of full synchrony ($K=1)$, the statement
is relevant to the unique synchronized periodic orbit of Section 3. Recall
that this orbit exponentially attracts all synchronized trajectories.
Moreover, the condition $\tau<R$ ensures that no firing can happen in the
interval $[-\tau,0]$. In addition, since $m_{\text{FP}}(t)=-t$ for
$t\in[-\tau,0]$, we have $\dot{A}^{\text{FP}}(-\tau)=\tau-\beta
A^{\text{FP}}(-\tau)$. Since we showed that $A^{\text{FP}}(0)>0$, the
continuity at 0 of the map $\tau\mapsto\tau-\beta A^{\text{FP}}(-\tau)$
implies that the conditions of Lemma 4.2 certainly hold when $\tau$ is small
enough, viz. statement (ii) of Proposition 3.1 is recovered as a particular
instance where Lemma 4.2 applies.

Proof of Lemma 4.2. Given $K^{\prime}\in\\{K+1,\cdots,N\\}$, let
$\left(\\{n^{\prime}_{k},y_{k}\\}_{k=1}^{K^{\prime}},A|_{[-\tau,0]}\right)$
with $y_{K^{\prime}}=0$ be such that

$\max\left\\{\|\mathbf{y}-\mathbf{y}^{\text{FP}}\|_{K-1},C_{1}\max_{k\in\\{K,\cdots,K^{\prime}-1\\}}y_{k},\|A|_{[-\tau,0]}-A^{\text{FP}}_{[-\tau,0]}\|_{0}\right\\}\leq\epsilon$

where $\epsilon\in(0,\epsilon_{0})$ is arbitrary and $\epsilon_{0},C_{1}>0$
are to be determined. We are going to evaluate the amplitude of the
perturbation at the instants $t^{n}_{\text{R}}$ by considering separately the
coordinates $y_{k}(t^{n}_{R})$ for $k\in\\{1,\cdots K-1\\}$ to those for
$k\in\\{K,\cdots,K^{\prime}-1\\}$.

Analysis of the coordinates $\\{y_{k}(t^{n}_{\text{R}})\\}_{k=1}^{K-1}$ and
$A|_{[t^{n}_{\text{R}}-\tau,t^{n}_{\text{R}}]}$. For these coordinates, we use
the fact that the fixed point is asymptotically stable inside its proper
partially synchronized subspace. To do so, we rely on a property of Lipschitz
continuity of the return map, which is claimed and proved in Appendix A.

Let
$t\mapsto\left(\\{n_{k},y^{\text{p.sync}}_{k}(t)\\}_{k=1}^{K},A^{\text{p.sync}}(t)\right)$
be the partially synchronized solution issued from the initial datum
$\left(\\{n_{k},(1-\delta_{k,K})y_{k}\\}_{k=1}^{K},A|_{[-\tau,0]}\right)$.666$\delta_{i,j}$
is the Kronecker symbol. Notice that the value of $A(-\tau)$ in this
trajectory is the same as the value of $A(-\tau)$ in the original trajectory
$t\mapsto\left(\\{n^{\prime}_{k},y_{k}(t)\\}_{k=1}^{K^{\prime}},A(t)\right)$.
Hence, the first return time to $y_{K}=0$ for this trajectory is the same as
the first return time $t^{1}_{\text{R}}$ of the original trajectory.

Consider the decomposition

$\|\mathbf{y}(t^{1}_{\text{R}})-\mathbf{y}^{\text{FP}}\|_{K-1}\leq\|\mathbf{y}(t^{1}_{\text{R}})-\mathbf{y}^{\text{p.sync}}(t^{1}_{\text{R}})\|_{K-1}+\|\mathbf{y}^{\text{p.sync}}(t^{1}_{\text{R}})-\mathbf{y}^{\text{FP}}\|_{K-1}.$

By periodicity of the trajectory passing through
$\left(\\{n_{k},y^{\text{FP}}_{k}\\}_{k=1}^{K},A^{\text{FP}}|_{[-\tau,0]}\right)$,
we must have

$y_{1}^{\text{FP}}=\max_{k\in\\{1,\cdots,K-1\\}}y^{\text{FP}}_{k}<T_{\text{FP}}=R+\nu
A^{\text{FP}}(-\tau).$

Hence, provided that $\epsilon$ is sufficiently small, we have

$y_{1}=\max_{k\in\\{1,\cdots,K^{\prime}-1\\}}y_{k}<t^{1}_{\text{R}}=R+\nu
A(-\tau)$

Applying Lemma A.1 in Appendix A, we obtain the following inequality, when
regarding $\mathbf{y}^{\text{p.sync}}$ and
$\mathbf{y}^{\text{p.sync}}(t^{1}_{\text{R}})$ as elements of the partially
synchronized subspace defined by $\\{n^{\prime}_{k}\\}_{k=1}^{K^{\prime}}$
whose coordinates $k\in\\{K,\cdots,K^{\prime}\\}$ all vanish

$\|\mathbf{y}(t^{1}_{\text{R}})-\mathbf{y}^{\text{p.sync}}(t^{1}_{\text{R}})\|_{K-1}\leq\|\mathbf{y}(t^{1}_{\text{R}})-\mathbf{y}^{\text{p.sync}}(t^{1}_{\text{R}})\|_{K^{\prime}-1}\leq
L\|\mathbf{y}-\mathbf{y}^{\text{p.sync}}\|_{K^{\prime}-1}=L\max_{k\in\\{K,\cdots,K^{\prime}-1\\}}y_{k}\leq
L\frac{\epsilon}{C_{1}}.$

On the other hand, asymptotic stability inside the partially synchronized
subspace implies the existence of $\gamma_{1}\in(0,1)$ such that, provided
that $\epsilon$ is sufficiently small, we have

$\|\mathbf{y}^{\text{p.sync}}(t^{1}_{\text{R}})-\mathbf{y}^{\text{FP}}\|_{K-1}\leq\gamma_{1}\epsilon,$

and then

$\|\mathbf{y}(t^{1}_{\text{R}})-\mathbf{y}^{\text{FP}}\|_{K-1}\leq\left(\gamma_{1}+\frac{L}{C_{1}}\right)\epsilon.$

A similar reasoning applies to
$\|A|_{[t^{1}_{\text{R}}-\tau,t^{1}_{\text{R}}]}-A^{\text{FP}}_{[-\tau,0]}\|_{0}$,
which yields the same inequality. As a consequence, provided that $C_{1}$ is
large enough so that $\gamma_{2}:=\gamma_{1}+\frac{L}{C_{1}}<1$, we get

$\max\left\\{\|\mathbf{y}(t^{1}_{\text{R}})-\mathbf{y}^{\text{FP}}\|_{K-1},\|A|_{[t^{1}_{\text{R}}-\tau,t^{1}_{\text{R}}]}-A^{\text{FP}}_{[-\tau,0]}\|_{0}\right\\}\leq\gamma_{2}\epsilon.$

By induction, the same arguments show that for every sufficiently small
$\epsilon\in\mathbb{R}^{+}$, we have

$\max\left\\{\|\mathbf{y}(t^{n}_{\text{R}})-\mathbf{y}^{\text{FP}}\|_{K-1},\|A|_{[t^{n}_{\text{R}}-\tau,t^{n}_{\text{R}}]}-A^{\text{FP}}_{[-\tau,0]}\|_{0}\right\\}\leq\gamma_{2}^{n}\epsilon,\
\forall n\in\mathbb{N}$

provided that one can simultaneously ensure

$\max_{k\in\\{K,\cdots,K^{\prime}-1\\}}y_{k}(t^{n}_{\text{R}})\leq\frac{\gamma_{2}^{n}\epsilon}{C_{1}}.$

Analysis of the coordinates
$\\{y_{k}(t^{n}_{\text{R}})\\}_{k=K}^{K^{\prime}-1}$. For these coordinates,
we separate the cases $n=1$ and $n>1$. For $n=1$, we simply refer to Lemma A.1
to obtain

$\max_{k\in\\{K,\cdots,K^{\prime}-1\\}}y_{k}(t^{1}_{\text{R}})\leq L\epsilon,$

and the desired inequality holds for any pair $\gamma,C_{2}$ such that
$L\leq\frac{C_{2}\gamma}{C_{1}}$. For $n>1$, we proceed by induction. From the
computation of the return map in Section 2.2.4, we have in particular at the
instant $t^{2}_{\text{R}}$ of the third firing of the cluster $K^{\prime}$

$y_{k}(t^{2}_{\text{R}})=y_{k}(t^{1}_{\text{R}})+\nu\left(A(y_{k}(t^{1}_{\text{R}})+t^{1}_{\text{R}}-\tau)-A(t^{1}_{\text{R}}-\tau)\right),\
\forall k\in\\{K,\cdots,K^{\prime}-1\\}.$

The assumption that the periodic orbit does not fire in the time interval
$[T_{\text{FP}}-\tau,T_{\text{FP}}]$ is equivalent to assuming that
$T_{\text{FP}}-y_{1}^{\text{FP}}>\tau$ (NB: recall that
$y_{1}^{\text{FP}}=\max_{k}y_{k}^{\text{FP}}$). Using that
$t^{1}_{\text{R}}-T_{\text{FP}}=\nu(A(-\tau)-A^{\text{PF}}(-\tau))$, the
assumption on the initial datum
$\left(\\{n^{\prime}_{k},y_{k}\\}_{k=1}^{K^{\prime}},A|_{[-\tau,0]}\right)$
implies that we also have $t^{1}_{\text{R}}-y_{1}>\tau$ (and
$y_{1}=\max_{k}y_{k}$) when $\epsilon$ is sufficiently small, viz. the
subsequent trajectory does not fire in the time interval
$[t^{1}_{\text{R}}-\tau,t^{1}_{\text{R}}]$. Therefore, for $\epsilon$
sufficiently small, we have

$\displaystyle y_{k}(t^{2}_{\text{R}})$
$\displaystyle=y_{k}(t^{1}_{\text{R}})\left(1+\nu\dot{A}(t^{1}_{\text{R}}-\tau)\right)+{\cal
O}\left((y_{k}(t^{1}_{\text{R}}))^{2}\right)$
$\displaystyle=y_{k}(t^{1}_{\text{R}})\left(1+\nu(m(t^{1}_{\text{R}}-\tau)-\beta
A(t^{1}_{\text{R}}-\tau))\right)+{\cal
O}\left((y_{k}(t^{1}_{\text{R}}))^{2}\right),\ \forall
k\in\\{K,\cdots,K^{\prime}-1\\}$

The analysis above of
$\|A|_{[t^{1}_{\text{R}}-\tau,t^{1}_{\text{R}}]}-A^{\text{FP}}_{[-\tau,0]}\|_{0}$
showed that $|A(t^{1}_{\text{R}}-\tau)-A^{\text{FP}}(-\tau)|$ can be made
arbitrarily small by taking $\epsilon$ sufficiently small. Moreover, the fact
that no firing takes place in the time interval
$[t^{1}_{\text{R}}-\tau,t^{1}_{\text{R}}]$ implies that

$m(t^{1}_{\text{R}}-\tau)=m(t^{1}_{\text{R}})+\tau.$

Similarly, we have $m_{\text{FP}}(-\tau)=m_{\text{FP}}(0)+\tau$ for the
periodic trajectory. In addition, the arguments above showed that
$\|\mathbf{y}(t^{1}_{\text{R}})-\mathbf{y}^{\text{FP}}\|_{K^{\prime}-1}$ can
be made arbitrarily small by taking $\epsilon$ sufficiently small.
Consequently, the same property holds for
$|m(t^{1}_{\text{R}}-\tau)-m_{\text{FP}}(-\tau)|$. Altogether, using the
assumption $m_{\text{FP}}(-\tau)-\beta A^{\text{FP}}(-\tau)<0$, this proves
the existence of $\gamma_{3}\in(0,1)$ such that, provided that $\epsilon$ is
sufficiently small, we have

$y_{k}(t^{2}_{\text{R}})\leq\gamma_{3}y_{k}(t^{1}_{\text{R}}),\ \forall
k\in\\{K,\cdots,K^{\prime}-1\\},$

from where the desired induction follows with
$\gamma=\max\\{\gamma_{2},\gamma_{3}\\}$ when the results of the analysis of
the other coordinates are also taken into account. The proof of the Lemma is
complete. $\Box$

###  4.4 Instability criterion for $\tau=0$

The announced instability criterion for partially synchronized periodic orbits
is given in the next statement.

######  Lemma 4.3.

Assume that (11) holds. Given $N\in\mathbb{N}$, $K\in\\{1,\cdots,N-1\\}$ and a
cluster distribution $\\{n_{k}\\}_{k=1}^{K}$, assume that, for $\tau=0$, the
return map in the Poincaré section $y_{K}=0$ in the corresponding partially
synchronized subspace has a fixed point
$\left(\\{n_{k},y^{\text{FP}}_{k}\\}_{k=1}^{K},A^{\text{FP}}\right)$. Let
$\dot{A}^{\text{FP}}(0^{+})$ be the right limit at 0 of the derivative of the
activator concentration in the corresponding trajectory. Then, under the
condition

$\dot{A}^{\text{FP}}(0^{+})>\frac{R+\nu A^{\text{FP}}}{N}$

the fixed point is unstable with respect to some arbitrarily small
perturbations of cluster $K$.

In the case of the fully synchronized periodic orbit ($K=1$), we have
$\dot{A}^{\text{FP}}(0^{+})=R+\nu A^{\text{FP}}-\beta A^{\text{FP}}$; hence
the condition in the Proposition becomes

$\frac{N-1}{N}(R+\nu A^{\text{FP}})-\beta A^{\text{FP}}>0,$

which is exactly the identified and proved criterion in the proof of statement
(i) in Proposition 3.1.

Proof of Lemma 4.3. Given $K^{\prime}\in\\{K+1,\cdots,N\\}$, let
$\left(\\{n^{\prime}_{k},y_{k}\\}_{k=1}^{K^{\prime}},A(0)\right)$ with
$y_{K^{\prime}}=0$ be an initial condition such that

$y_{k}=\left\\{\begin{array}[]{ccl}y_{k}^{\text{FP}}&\text{if}&k\in\\{1,\cdots,K-1\\}\\\
y_{K}&\text{if}&k=K\\\
0&\text{if}&k\in\\{K+1,\cdots,K^{\prime}\\}\end{array}\right.$

where $y_{K}>0$ is sufficiently small. There is no firing in the time interval
$(0,y_{K})$. Hence, we have

$y_{K}(t_{\text{R}})=y_{K}+\nu(A(y_{K})-A(0))=y_{K}(1+\nu\dot{A}(0^{+}))+{\cal
O}(y_{K}^{2})$

Moreover, equation (5) implies that $\dot{A}(0^{+})=m(0^{+})-\beta A(0)$.
Using

$y_{k}(0^{+})=\left\\{\begin{array}[]{ccl}y_{k}^{\text{FP}}&\text{if}&k\in\\{1,\cdots,K-1\\}\\\
y_{K}&\text{if}&k=K\\\ R+\nu
A(0)&\text{if}&k\in\\{K+1,\cdots,K^{\prime}\\}\end{array}\right.\quad\text{and}\quad
y_{k}^{\text{FP}}(0^{+})=\left\\{\begin{array}[]{ccl}y_{k}^{\text{FP}}&\text{if}&k\in\\{1,\cdots,K-1\\}\\\
R+\nu A^{\text{FP}}&\text{if}&k=K\end{array}\right.$

we get

$m(0^{+})=m^{\text{FP}}(0^{+})-\frac{R+\nu
A(0)}{N}+y_{K}+\frac{n_{K}\nu(A(0)-A^{\text{FP}})}{N}.$

It follows that $m(0^{+})-\beta A(0)$ can be made arbitrarily close to
$\dot{A}^{\text{FP}}(0^{+})-\frac{R+\nu A^{\text{FP}}}{N}$ when $y_{K}$ and
$|A(0)-A^{\text{FP}}|$ are sufficiently small. The Lemma easily follows.
$\Box$

###  4.5 Delay-dependent stability in arbitrary large populations

Let $K\in\mathbb{N}$ and $\\{n_{k}\\}_{k=1}^{K}\in\mathbb{N}^{K}$ be given.
The equations (3)-(5) of the dynamics imply that if
$t\mapsto\left(\\{n_{k},y_{k}(t)\\}_{k=1}^{K},A(t)\right)$ is a partially
synchronized trajectory in the population of $\sum_{k=1}^{K}n_{k}$
oscillators, then for every $q\in\mathbb{N}$, the function
$t\mapsto\left(\\{qn_{k},y_{k}(t)\\}_{k=1}^{K},A(t)\right)$ is a partially
synchronized trajectory in the population of $q\sum_{k=1}^{K}n_{k}$
oscillators, viz. the existence and the coordinates of partially synchronized
trajectories do not depend on the (common) scaling of its cluster sizes.

The results on the stability analysis of partially synchronized periodic
orbits are summarized in the following statement, whose conclusions can be
regarded as some extension of Proposition 3.1 to fixed points of the return
map asssociated with an arbitrary cluster distribution (with given relative
cluster sizes).

######  Theorem 4.4.

Assume that (11) holds. Given $K\in\mathbb{N}$ and
$\\{n_{k}\\}_{k=1}^{K}\in\mathbb{N}^{K}$, assume that for $\tau=0$, the return
map in the partially synchronized subspace associated with
$\\{n_{k}\\}_{k=1}^{K}$, of the dynamics of a population of
$\sum_{k=1}^{K}n_{k}$ DF oscillators, has an exponentially stable fixed point,
say $\left(\\{n_{k},y^{\text{FP}}_{k}\\}_{k=1}^{K},A^{\text{FP}}\right)$.
Then, the following assertions hold.

(i) There exists $q_{0}\in\mathbb{N}$ such that for $\tau=0$ and every
$q>q_{0}$, the fixed point
$\left(\\{qn_{k},y^{\text{FP}}_{k}\\}_{k=1}^{K},A^{\text{FP}}\right)$ is
unstable with respect to some arbitrarily small perturbations of any of its
clusters.

(ii) There exists $\tau_{0}\in\mathbb{R}_{\ast}^{+}$ such that for every
$\tau\in(0,\tau_{0})$ and $q\in\mathbb{N}$, the return map of the partially
synchronized dynamics with delay $\tau$ has a fixed point which is the
continuation of
$\left(\\{qn_{k},y^{\text{FP}}_{k}\\}_{k=1}^{K},A^{\text{FP}}\right)$. This
continued fixed point is exponentially stable with respect to small
perturbations that smear any of its clusters.

As an application, we provide in the next section, an example of periodic
orbits that fit the condition of this statement for $\beta$ sufficiently
large, namely the periodic orbits with equi-distributed repressor
concentrations. Evidently, the fully synchronized fixed point also satisfies
the assumption of Theorem 4.4, as already claimed at the begining of
Proposition 3.1.777Notice that item (i) in Theorem 4.4 is slightly weaker than
the corresponding item in Proposition 3.1 (which claims that $q_{0}=1$ in this
case).

Proof of Theorem 4.4. The proof decomposes into two parts. The first part
establishes the existence of the continued fixed point for $\tau>0$
sufficiently small. The second part shows that such return map fixed points
must appropriately satisfy the conditions of Lemmas 4.2 and 4.3.

For the first part, one preliminary shows that continuity arguments ensure
that the conditions of Claim 2.1 hold for $\tau$ sufficiently small and in a
sufficiently small neighbourhood of
$\left(\\{y^{\text{FP}}_{k}\\}_{k=1}^{K},A^{\text{FP}}\right)$ in
$(\mathbb{R}^{+})^{K}$; hence the return dynamics in that neighbourhood
reduces to a mapping of $(\mathbb{R}^{+})^{K}$ into itself. Moreover, this map
can be obtained as the composition of the maps, each called a firing map, that
bring the system from the state immediately before the firing of the $(k+1)$th
cluster to the one of the $k$th cluster. The firing maps write, given a datum
$\left(\\{y_{k}\\}_{k=1}^{K},A\right)$ for which $y_{k+1}=0$ and
$y_{k}=\min_{\ell\neq k+1}y_{\ell}$

$y^{\prime}_{\ell}=\left\\{\begin{array}[]{ccl}y_{\ell}-y_{k}&\text{if}&\ell\neq
k+1\\\
R+\nu\phi_{m}^{-\tau}(A)-y_{k}&\text{if}&\ell=k\end{array}\right.\quad\text{and}\quad
A^{\prime}=\phi_{m}^{y_{k}}(A),$

where the expressions (6) and (7) are to be used with respectively

$m(t)=\frac{1}{N}\left(\sum_{k=1}^{K}n_{k}y_{k}+R+\nu\phi_{m}^{-\tau}(A)\right)-t,\
\forall t\in(0,y_{k})\quad\text{and}\quad
m(t)=\frac{1}{N}\sum_{k=1}^{K}n_{k}y_{k}-t\ \forall t\in[-\tau,0].$

Accordingly, each firing map is differentiable in $\mathbb{R}^{K}$ and,
together with its derivative, it continuously depends on $\tau$. Hence, so
does the composed return map. The assumption that
$\left(\\{n_{k},y^{\text{FP}}_{k}\\}_{k=1}^{K},A^{\text{FP}}\right)$ is an
exponentially stable fixed point implies, using the Implicit Function Theorem,
that it can be uniquely continued as an exponentially stable fixed point of
the return map for $\tau\in\mathbb{R}_{\ast}^{+}$ sufficiently small. The
first part of the proof is complete.

The main argument of the second part consists in establishing the following
signs for the left and right limits of the derivative of
$A_{\tau}^{\text{FP}}$ (where the explicit dependence on $\tau$ has been added
for clarity) at each firing

$\dot{A}_{\tau}^{\text{FP}}(\left(y_{\tau,k}^{\text{FP}}\right)^{-})<0\quad\text{and}\quad\dot{A}_{\tau}^{\text{FP}}(\left(y_{\tau,k}^{\text{FP}}\right)^{+})>0\
\forall k\in\\{1,\cdots,K\\}.$

In order to prove these signs, notice that equation (5) and the fact that
$\dot{m}(t)=-1$ between two firings imply that if it happens that
$\dot{A}(t)=0$ for some $t$ between two firings, then this can happen only
once and $\dot{A}(t)>0$ for all $t$ before (resp. $\dot{A}(t)<0$ after) that
instant. Moreover, equation (12), or the expression of $x_{i}(t_{\text{R}})$
in Section 2.2.4 implies that, for a fixed point
$\left(\\{n_{k},y^{\text{FP}}_{0,k}\\}_{k=1}^{K},A_{0}^{\text{FP}}\right)$ of
the return map for $\tau=0$, the activator concentration must be the same at
each firing, ie.

$A_{0}^{\text{FP}}(y_{0,k}^{\text{FP}})=A_{0}^{\text{FP}},\ \forall
k\in\\{1,\cdots,K-1\\}.$

Therefore, we must have $\dot{A}_{0}^{\text{FP}}(t)=0$ for some $t$ between
any two firings, and then the desired signs hold for $\tau=0$. By continuous
dependence on $\tau$ of the coordinates of the continued fixed point, the same
signs hold for $A_{\tau}^{\text{FP}}$ provided that $\tau$ is small enough.

Together with
$\dot{A}_{\tau}^{\text{FP}}(\left(y_{\tau,k}^{\text{FP}}\right)^{-})<0$, the
implicit assumption in the continuation argument that there is no firing in
any of the intervals $[y_{\tau,k}^{\text{FP}}-\tau,y_{\tau,k}^{\text{FP}}]$,
and the fact that the fixed point coordinates do not depend on $q$ imply that
the assumptions of Lemma 4.2 hold for the fixed point
$\left(\\{qn_{k},y^{\text{FP}}_{0,k}\\}_{k=1}^{K},A_{0}^{\text{FP}}\right)$,
for every $q\in\mathbb{N}$. Moreover, the same assumptions hold for the return
map fixed point of the corresponding periodic orbit in every Poincaré section
$y_{k}=0$, implying stability with respect to smearing of any cluster.

Together with
$\dot{A}_{\tau}^{\text{FP}}(\left(y_{\tau,k}^{\text{FP}}\right)^{+})>0$ and
the fact that the fixed point coordinates do not depend on $q$, we certainly
have

$\dot{A}_{\tau}^{\text{FP}}(\left(y_{\tau,k}^{\text{FP}}\right)^{+})>\frac{R+\nu
A^{\text{FP}}}{q\sum_{k=1}^{K}n_{k}},$

when $q$ is sufficiently large. In this case, the assumptions of Lemma 4.3
hold, as they do for each corresponding return map fixed point in every
Poincaré section $y_{k}=0$. The proof of Theorem is complete. $\Box$

###  4.6 Application to periodic orbits with equi-distributed repressor
concentrations

As a example of application of Theorem 4.4, we consider in this section,
return map fixed points for which the repressor coordinates are equi-
distributed among $N$ clusters, more precisely, those elements
$(\mathbf{x}^{\text{FP}},A^{\text{FP}})\in\mathbb{R}^{N+1}$ such that
$x^{\text{FP}}_{N}=0$, the function $m(s)=m^{\text{FP}}+s$ for $s\in[-\tau,0]$
and 888By letting $x^{\text{FP}}_{0}=R+\nu\phi_{m}^{-\tau}(A^{\text{FP}})$,
one could include the synchronized fixed point in this family, that would be
obtained for $N=1$.

$x^{\text{FP}}_{i}-x^{\text{FP}}_{i+1}=\frac{R+\nu\phi_{m}^{-\tau}(A^{\text{FP}})}{N},\
\forall i\in\\{1,\cdots,N-1\\}\quad\text{and}\quad
A^{\text{FP}}=\phi_{m}^{x^{\text{FP}}_{N-1}}(A^{\text{FP}}).$

It is immediate to verifiy that the coordinates of
$(\mathbf{x}^{\text{FP}},A^{\text{FP}})$ are recovered after every firing in
the trajectory, not only under the action of the return map. Moreover, the
repressor coordinates are entirely determined by $A^{\text{FP}}$. For
$\tau=0$, an equation for $A^{\text{FP}}$ can be obtained, which can be shown
to have exactly one solution, viz. for each $N>1$, there exists a unique
$N$-cluster equi-distributed fixed point. Moreover, a systematic stability
analysis can be achieved for $\tau=0$, which is rigorous for
$N\in\\{2,\cdots,5\\}$ and concludes that the fixed point is stable in its
partially synchrony subspace, provided that $\beta$ is sufficiently large. All
these results are presented in the next statement.

######  Lemma 4.5.

(i) For every $N\in\mathbb{N}$, $N>1$ and $\tau=0$, there exists a unique
fixed point with equi-distributed repressor coordinates.

(ii) Given $R,\nu$ and $N\in\\{2,\cdots,5\\}$, there exists $\beta_{N}>0$ such
that for every $\beta>\beta_{N}$, the equi-distributed fixed point is
exponentially stable in its proper partially synchronized subspace.

We believe that statement (ii) holds for every $N>1$. In order to ensure this
property, it suffices to prove that the expression of the matrix $W_{N}$
defined in the proof below, holds for every $N$. Independently, by combining
Lemma 4.5 with Theorem 4.4, one immediately obtains the following conclusion.

######  Corollary 4.6.

Given $N\in\\{2,\cdots,5\\}$, assume that (11) holds with $\beta>\beta_{N}$ so
that for $\tau=0$, the equi-distributed fixed point with $N$ clusters is
exponentially stable in its proper partially synchronized subspace.

(i) There exists $q_{N}\in\mathbb{N}$ such that for $\tau=0$ and every
$q>q_{N}$, the fixed point
$\left(\\{q,x^{\text{FP}}_{i}\\}_{i=1}^{N},A^{\text{FP}}\right)$ is unstable
with respect to some arbitrarily small perturbations of any of its clusters.

(ii) There exists $\tau_{N}\in\mathbb{R}_{\ast}^{+}$ such that for every
$\tau\in(0,\tau_{N})$ and $q\in\mathbb{N}$, the return map of the partially
synchronized dynamics with delay $\tau$ has a fixed point which is the
continuation of
$\left(\\{q,x^{\text{FP}}_{i}\\}_{i=1}^{N},A^{\text{FP}}\right)$. This
continued fixed point is exponentially stable with respect to small
perturbations that smear any of its clusters.

One can show that the continuation of $(\mathbf{x}^{\text{FP}},A^{\text{FP}})$
for $\tau>0$ has indeed equi-distributed repressor concentrations. Moreover,
recall from Claim 2.1 that $\tau_{N}$ also depends on the distance between the
repressor concentrations immediately after firings. Hence, even if we assumed
that statement (ii) in Lemma 4.5 held for all $N>1$ with
$\sup_{N>1}\beta_{N}<+\infty$, we would not be able to ensure that
$\inf_{N>1}\tau_{N}>0$. In order words, we do not know whether or not all
equidistributed fixed points can be simultanously stable for some given delay
$\tau>0$.

Proof of Lemma 4.5. (i) As usual, the return map fixed point
$(\mathbf{x}^{\text{FP}},A^{\text{FP}})$ is assumed to coincide with the state
at $t=0$ of the periodic orbit
$t\mapsto(\mathbf{x}^{\text{FP}}(t),A^{\text{FP}}(t))$ of the continuous time
system. The fixed point definition implies that there is no firing in the time
interval $(0,x^{\text{FP}}_{N-1}]$ where $x^{\text{FP}}_{N-1}=\frac{R+\nu
A^{\text{FP}}}{N}$. Hence, we have

$m^{\text{FP}}(t)=m^{\text{FP}}(0^{+})-t,\ \forall t\in(0,x_{N-1}]\
\text{where}\
m^{\text{FP}}(0^{+})=\frac{1}{N}\left(\sum_{i=1}^{N-1}x_{i}^{\text{FP}}+R+\nu
A^{\text{FP}}\right)=\frac{(N+1)(R+\nu A^{\text{FP}})}{2N}.$

Using (6), the fixed point equation
$A^{\text{FP}}=\phi_{m^{\text{FP}}}^{x^{\text{FP}}_{N-1}}(A^{\text{FP}})$ then
rewrites as $f(A^{\text{FP}})=0$ where $f$ is given by the following
expression

$f(A)=\left(\frac{(N+1)(R+\nu A)}{2N}-\beta
A+\frac{1}{\beta}\right)(1-e^{-\beta\frac{R+\nu A}{N}})-\frac{R+\nu A}{N}.$

In particular, we have $f(0)=\frac{g(\beta R)}{\beta N}$ where

$g(x)=\left(\frac{(N+1)x}{2}+N\right)(1-e^{-\frac{x}{N}})-x.$

Basic calculations yield

$g(0)=0\quad\text{and}\quad g^{\prime}(x)>0,\ \forall x>0,$

hence $f(A)>0$ for all $A>0$ sufficiently small. Moreover, recall the upper
bound $A_{\text{max}}=\frac{R}{\beta-\nu}$ of the attracting set. We have

$f(A_{\text{max}})=-\frac{\beta(N-1)A_{\text{max}}}{2N}(1-e^{-\beta^{2}\frac{A_{\text{max}}}{N}})+\frac{1}{\beta}(1-e^{-\beta^{2}\frac{A_{\text{max}}}{N}})-\frac{\beta
A_{\text{max}}}{N}<0$

because $1-e^{-x}<x$ for all $x>0$. Therefore, the function $f$ must have a
zero $A^{\text{FP}}$ in the interval $(0,A_{\text{max}})$.

In order to prove uniqueness, we successively compute the first and second
derivatives of $f$ to obtain

$f^{\prime\prime}(A)=\frac{\beta\nu}{N}\left(2(\nu-\beta)-\frac{\beta\nu}{N}\left(\frac{(N+1)(R+\nu
A)}{2N}-\beta A\right)\right)e^{-\beta\frac{R+\nu A}{N}}.$

Since $\beta>\nu>\nu\frac{N+1}{2N}$, the expression inside the parenthesis is
increasing, viz. the sign of $f^{\prime\prime}$ can only change once in
$[0,A_{\text{max}}]$, from negative to positive. Together with the facts that
$f(0)>0$ and $f(A_{\text{max}})<0$, this implies that $f$ crosses 0 only once
in this interval.

(ii) As in [4], the finite-dimensional return dynamics (that takes place in
the neighbourhod of $(\mathbf{x}^{\text{FP}},A^{\text{FP}})$) can be regarded
as the composition of the following variation of the firing maps introduced at
the begining at the proof of Theorem 4.4

$x^{\prime}_{i}=\left\\{\begin{array}[]{ccl}R+\nu A-x_{N-1}&\text{if}&i=1\\\
x_{i-1}-x_{N-1}&\text{if}&i\in\\{2,\cdots,N-1\\}\\\
\end{array}\right.\quad\text{and}\quad A^{\prime}=\phi_{m}^{x_{N-1}}(A),$ (15)

More precisely, the local stability of
$(\mathbf{x}^{\text{FP}},A^{\text{FP}})$ can be obtained from the spectrum of
the product of the derivatives of this map along the elements of the
corresponding periodic orbit. However, since
$(\mathbf{x}^{\text{FP}},A^{\text{FP}})$ is a fixed point of this firing map,
one actually has to evaluate the spectrum of $M_{N}^{N}$ where $M_{N}$ is the
derivative of (15) evaluated at $(\mathbf{x}^{\text{FP}},A^{\text{FP}})$ with

$m(t)=m(0^{+})-t=\frac{1}{N}\left(\sum_{i=1}^{N-1}x_{i}+R+\nu A\right)-t\
\text{for}\ \in(0,x_{N-1}].$

From (6) and $m(t)=m(0^{+})-t$, we have

$\phi_{m}^{t}(A)=Ae^{-\beta
t}+\left(\frac{m(0^{+})}{\beta}+\frac{1}{\beta^{2}}\right)(1-e^{-\beta
t})-\frac{t}{\beta},$

which yields

$\partial_{x_{i}}\phi_{m}^{x_{N-1}}(A)=\left\\{\begin{array}[]{ccl}\frac{1-e^{-\beta
x_{N-1}}}{\beta N}&\text{if}&i\in\\{1,\cdots,N-2\\}\\\
\left(m(0^{+})+\frac{1}{\beta}-\beta A\right)e^{-\beta
x_{N-1}}-\frac{1}{\beta}+\frac{1-e^{-\beta x_{N-1}}}{\beta
N}&\text{if}&i=N-1\end{array}\right.$

and

$\partial_{A}\phi_{m}^{x_{N-1}}(A)=e^{-\beta x_{N-1}}+\frac{\nu(1-e^{-\beta
x_{N-1}})}{\beta N}.$

Evaluating these quantities at $(\mathbf{x}^{\text{FP}},A^{\text{FP}})$ and
expanding in $\frac{1}{\beta}$, one gets that the derivative $M_{N}$ writes

$M_{N}=U_{N}+\frac{1}{\beta}V_{N}+o(\frac{1}{\beta}),$

where $U_{N}$ and $V_{N}$ are the following $N\times N$ matrices

$U_{N}=\left(\begin{array}[]{cccccc}0&\cdots&\cdots&0&-1&\nu\\\
1&0&\cdots&0&-1&0\\\ 0&1&\ddots&\vdots&\vdots&\vdots\\\
\vdots&&\ddots&\vdots&\vdots&\vdots\\\ 0&\cdots&0&1&-1&0\\\
0&\cdots&\cdots&\cdots&\cdots&0\end{array}\right)\quad\text{and}\quad
V_{N}=\left(\begin{array}[]{cccccc}0&\cdots&\cdots&\cdots&\cdots&0\\\
\vdots&&&&&\vdots\\\ \vdots&&&&&\vdots\\\ \vdots&&&&&\vdots\\\ \\\
0&\cdots&\cdots&\cdots&\cdots&0\\\
\frac{1}{N}&\cdots&\cdots&\frac{1}{N}&-\frac{N-1}{N}&\frac{\nu}{N}\end{array}\right).$

Accordingly, we have

$M_{N}^{N}=W_{N}+o(\frac{1}{\beta})\quad\text{where}\quad
W_{N}:=U_{N}^{N}+\frac{1}{\beta}\sum_{k=1}^{N-1}U_{N}^{k}V_{N}U_{N}^{N-1-k}.$

For $N\in\\{2,\cdots,5\\}$, we have checked that $W_{N}$ writes

$W_{N}=\left(\begin{array}[]{cccccc}1-\frac{(N-1)\nu}{N\beta}&\frac{\nu}{N\beta}&\cdots&\cdots&\frac{\nu}{N\beta}&-\nu(1-\frac{\nu}{N\beta})\\\
\frac{\nu}{N\beta}&1-\frac{(N-1)\nu}{N\beta}&\frac{\nu}{N\beta}&\cdots&\frac{\nu}{N\beta}&-\nu(1-\frac{\nu}{N\beta})\\\
\frac{\nu}{N\beta}&\frac{\nu}{N\beta}&\ddots&\ddots&\vdots&\vdots\\\
\vdots&&\ddots&\ddots&\frac{\nu}{N\beta}&\vdots\\\
\frac{\nu}{N\beta}&\cdots&\cdots&\frac{\nu}{N\beta}&1-\frac{(N-1)\nu}{N\beta}&-\nu(1-\frac{\nu}{N\beta})\\\
\frac{1}{N\beta}&\cdots&\cdots&\cdots&\frac{1}{N\beta}&-\frac{(N-1)\nu}{N\beta}\end{array}\right).$

The eigenvectors of $W_{N}$ and its eigenvalues can be readily obtained.
Firstly, notice that

$U_{N}(\nu,\cdots,\nu,1)^{T}=V_{N}(\nu,\cdots,\nu,1)^{T}=0,$

hence $W_{N}(\nu,\cdots,\nu,1)^{T}=0$ too. Moreover, $W_{N}$ clearly has the
following $N-2$ eigenvectors

$(1,-1,0,\cdots,0)^{T},\ (1,0,-1,0,\cdots,0)^{T},\cdots,\
(1,0,\cdots,0,-1,0,0)^{T},\ (1,0,\cdots,0,-1,0)^{T},$

(or explicitly $(v_{k})_{i}=\delta_{i,1}-\delta_{i,k}$ for
$i,k\in\\{2,\cdots,N-1\\}$) and the corresponding eigenvalue is equal to
$1-\frac{\nu}{\beta}$ in each case. Finally, the remaining eigenvector writes
$(N\beta-\nu,0,\cdots,0,1)^{T}$ and the corresponding eigenvalue is also equal
to $1-\frac{\nu}{\beta}$. Therefore, the spectrum of $W_{N}$ lies inside the
unit disk; hence so does the spectrum of $M^{N}_{N}$ when $\beta$ is
sufficiently large. Statement (ii) is proved. $\Box$

##  5 Concluding remarks

In this paper, we have presented an extended mathematical analysis of the
dynamics of the model introduced in [18] of DF oscillators coupled via a
common activator, with emphasis on the stability with respect to out-of-sync
perturbations depending on the delay in activator synthesis. After a study of
the basic properties of the flow associated with the equations (3)-(5), the
analysis has firstly considered the case of fully synchronized trajectories
and then has adressed arbitrary partially synchronized periodic orbits.

The most significant outcome of this endeavour is that the stability of
periodic orbits abruptly changes when the delay is switched on. From an
unstable solution, the orbit immediately becomes asymptotically stable with
respect to small perturbations that smear its clusters.

While this change of behaviour appears to be spectacular, it can be readily
apprehended from the general criteria of Section 4 (which themselves can be
intuited from the expression (12) of the return map) together with the profile
of the activator concentration close to firing, in the periodic orbit. More
precisely, Lemma 4.2 about stability requires that $A^{\text{FP}}(\cdot)$, as
a function of time, be decreasing immediately before firing. Lemma 4.3 about
instability needs that the derivative of this function be positive and
sufficiently large immediately after firing.

Notice that, according to the original equation (3)–(5), both properties
should be commmon features of the periodic orbits in this system, as
illustrated for the synchronized trajectory in the left panel of Fig. 2.
Indeed, the variation of the vector field acting on $A$ must be locally
minimal immediately before any firing (because $m$ decreases between two
consecutive firings) and firings trigger sudden increases of this vector
field.

Finally, we believe that these stylized features and the accompanying rapid
change in stability extend to more general, smooth models of DF oscillators
with similar ingredients, especially when firings occur on very short time
scales and the corresponding reset values are impacted by some delay in the
activator synthesis. Such extensions could be the subject of future studies.

##  Appendix A Lipschitz continuity of the return map

Let $N\in\mathbb{N}$ and $\tau>0$ be given. Consider an initial datum
$(\mathbf{x},A_{[-\tau,0]})$ in the section $x_{N}=0$,
$\|A|_{[-\tau,0]}\|_{0}<A_{\text{max}}$ and
$\text{max}_{\mathbf{x}}<x_{N}(0^{+})=t_{\text{R}}=R+\nu A(-\tau)$, so that
the trajectory lies in the attracting set $Q$ and the order in which the
oscillators fire remains the same starting from $t=0$. Assume also that
$A_{[-\tau,0]}$ is Lipschitz continuous, with maximal Lipschitz constant
$\beta A_{\text{max}}$ when in $Q$ (see the proof of Claim 2.3 above).

Given another datum $(\mathbf{x}^{\prime},A^{\prime}_{[-\tau,0]})$ with the
same constraints, let $t^{\prime}_{\text{R}}$ be the first return time to the
section $x_{N}=0$. We have the following statement

######  Lemma A.1.

In addition to the assumptions above on $(\mathbf{x},A_{[-\tau,0]})$ and
$(\mathbf{x}^{\prime},A^{\prime}_{[-\tau,0]})$, assume that

$\text{max}_{\mathbf{x}}<t^{\prime}_{\text{R}}\ \text{if}\
t^{\prime}_{\text{R}}\leq
t_{\text{R}}\quad\text{and}\quad\text{max}_{\mathbf{x}^{\prime}}<t_{\text{R}}\
\text{if}\ t_{\text{R}}\leq t^{\prime}_{\text{R}}.$

Then, there exists $L\in\mathbb{R}^{+}$ (which is independent of the data)
such that we have

$\max\left\\{\|\mathbf{x}(t_{\text{R}})-\mathbf{x}^{\prime}(t^{\prime}_{\text{R}})\|_{N-1},\|A|_{[t_{\text{R}}-\tau,t_{\text{R}}]}-A^{\prime}|_{[t^{\prime}_{\text{R}}-\tau,t^{\prime}_{\text{R}}]}\|_{0}\right\\}\leq
L\max\left\\{\|\mathbf{x}-\mathbf{x}^{\prime}\|_{N-1},\|A|_{[-\tau,0]}-A^{\prime}|_{[-\tau,0]}\|_{0}\right\\}.$

###### Proof.

Estimate of
$\|A|_{[t_{\text{R}}-\tau,t_{\text{R}}]}-A^{\prime}|_{[t^{\prime}_{\text{R}}-\tau,t^{\prime}_{\text{R}}]}\|_{0}$:
We have $t_{\text{R}}-t^{\prime}_{\text{R}}=\nu(A(-\tau)-A^{\prime}(-\tau))$.
Using (6), we consider the following decomposition similar to the one in the
proof of statement (ii) of Proposition 3.1

$\displaystyle A(t_{\text{R}}-t)-A^{\prime}(t^{\prime}_{\text{R}}-t)=$
$\displaystyle
A(t_{\text{R}}-t)(1-e^{-\beta(t^{\prime}_{\text{R}}-t_{\text{R}})})+\left(A(0)-A^{\prime}(0)+\int_{0}^{t_{\text{R}}\wedge
t^{\prime}_{\text{R}}-t}e^{\beta s}(m(s)-m^{\prime}(s))ds\right.$
$\displaystyle+\left.\int_{t_{\text{R}}\wedge
t^{\prime}_{\text{R}}-t}^{t_{\text{R}}\vee t^{\prime}_{\text{R}}-t}e^{\beta
s}\overline{m}(s)ds\right)e^{-\beta(t^{\prime}_{\text{R}}-t)}$

for $t\in[-\tau,0]$, where

$\overline{m}:=\left\\{\begin{array}[]{ccl}m&\text{if}&t_{\text{R}}>t^{\prime}_{\text{R}}\\\
-m^{\prime}&\text{if}&t_{\text{R}}<t^{\prime}_{\text{R}}\end{array}\right.$

That $A(t),A^{\prime}(t)$ are uniformly bounded implies the same property for
$t^{\prime}_{\text{R}}$ and $m,m^{\prime}$. Together with the expression of
$t_{\text{R}}-t^{\prime}_{\text{R}}$ above, this implies that the modulus of
the first, second and last terms in the above decomposition can be controlled
by $\|A|_{[-\tau,0]}-A^{\prime}|_{[-\tau,0]}\|_{0}$. For the remaining term,
we need to control $|m(s)-m^{\prime}(s)|$ for $s\in[0,t_{\text{R}}\wedge
t^{\prime}_{\text{R}}-t]$. Notice that each oscillator fires at most once in
this interval. Hence, we have

$x_{i}(s)-x^{\prime}_{i}(s)=\left\\{\begin{array}[]{ccl}x_{i}-x^{\prime}_{i}&\text{if}&s\leq
x_{i}\wedge x^{\prime}_{i}\\\
\nu(A(x_{i}-\tau)-A^{\prime}(x^{\prime}_{i}-\tau))+x_{i}-x^{\prime}_{i}&\text{if}&s>x_{i}\vee
x^{\prime}_{i}\end{array}\right.$

Accordingly, outside the intervals $[x_{i}\wedge x^{\prime}_{i},x_{i}\vee
x^{\prime}_{i}]$, the quantity $|m(s)-m^{\prime}(s)|$ is certainly bounded by

$\displaystyle\frac{\nu}{N}\sum_{i=1}^{N-1}|A(x_{i}-\tau)-A(x^{\prime}_{i}-\tau)|+\frac{\nu}{N}\sum_{i=1}^{N-1}|A(x^{\prime}_{i}-\tau)-A^{\prime}(x^{\prime}_{i}-\tau)|+\frac{1}{N}\sum_{i=1}^{N-1}|x_{i}-x^{\prime}_{i}|$
$\displaystyle\leq
L_{1}\|\mathbf{x}-\mathbf{x}^{\prime}\|_{N-1}+\nu\|A|_{[-\tau,0]}-A^{\prime}|_{[-\tau,0]}\|_{0}$

for some $L_{1}\in\mathbb{R}^{+}$ sufficiently large, where, in addition to
the fact that $A$ is Lipschitz continuous on $\mathbb{R}^{+}$, for those
$x_{i},x^{\prime}_{i}\in(0,\tau]$, the second inequality also relies on the
assumption that $A|_{[-\tau,0]}$ is Lipschitz continuous. Moreover, inside the
intervals $[x_{i}\wedge x^{\prime}_{i},x_{i}\vee x^{\prime}_{i}]$, the
quantity $|m(s)-m^{\prime}(s)|$ is uniformly bounded because the corresponding
trajectories are in $Q$. There are at most $N-1$ such intervals whose length
is bounded by $\|\mathbf{x}-\mathbf{x}^{\prime}\|_{N-1}$. Altogether, the
arguments here prove the existence of $L_{2},L_{3}\in\mathbb{R}^{+}$ (which do
not depend on $(\mathbf{x},A_{[-\tau,0]})$ and
$(\mathbf{x}^{\prime},A^{\prime}_{[-\tau,0]})$ when in $Q$) such that

$\|A|_{[t_{\text{R}}-\tau,t_{\text{R}}]}-A^{\prime}|_{[t^{\prime}_{\text{R}}-\tau,t^{\prime}_{\text{R}}]}\|_{0}\leq
L_{2}\|\mathbf{x}-\mathbf{x}^{\prime}\|_{N-1}+L_{3}\|A|_{[-\tau,0]}-A^{\prime}|_{[-\tau,0]}\|_{0}.$

Estimate of
$\|\mathbf{x}(t_{\text{R}})-\mathbf{x}^{\prime}(t^{\prime}_{\text{R}})\|_{N-1}$:
Assume that $t_{\text{R}}\leq t^{\prime}_{\text{R}}$, the other case can be
treated similarly. Consider the decomposition

$x_{i}(t_{\text{R}})-x^{\prime}_{i}(t^{\prime}_{\text{R}})=x_{i}(t_{\text{R}})-x^{\prime}_{i}(t_{\text{R}})+x^{\prime}_{i}(t_{\text{R}})-x^{\prime}_{i}(t^{\prime}_{\text{R}}).$

The assumptions
$\text{max}_{\mathbf{x}},\text{max}_{\mathbf{x}^{\prime}}<t_{\text{R}}$ imply
that all oscillators in both trajectories must have fired once when at instant
$t_{\text{R}}$. Accordingly, the first difference can be controlled using the
same estimate as above in the case $s>x_{i}\vee x^{\prime}_{i}$. Moreover, all
oscillators with $x^{\prime}_{i}>0$ have not fired a second time at instant
$t^{\prime}_{\text{R}}$. Hence, we have
$x^{\prime}_{i}(t_{\text{R}})-x^{\prime}_{i}(t^{\prime}_{\text{R}})=t^{\prime}_{\text{R}}-t_{\text{R}}$
which is also well under control. Altogether, this proves that a similar
estimate as for
$\|A|_{[t_{\text{R}}-\tau,t_{\text{R}}]}-A^{\prime}|_{[t^{\prime}_{\text{R}}-\tau,t^{\prime}_{\text{R}}]}\|_{0}$
holds for
$\|\mathbf{x}(t_{\text{R}})-\mathbf{x}^{\prime}(t^{\prime}_{\text{R}})\|_{N-1}$.
∎

## References

  * [1] T. Danino, O. Mondragón-Palomino, L. Tsimring and J. Hasty, _A synchronized quorum of genetic clocks_ , Nature 463 (2010) 326–330.
  * [2] T.S. Gardner, C.R. Cantor and J.J. Collins, _Construction of a genetic toggle switch in Escherichia coli_ , Nature 403 (2000) 339–342.
  * [3] M.B. Elowitz and S. Liebler, _A synthetic oscillatory network of transcriptional regulators_ , Nature 403 (2000) 335–338.
  * [4] B. Fernandez and L.S. Tsimring, _Corepressive interaction and clustering of degrade-and-fire oscillators_ , Phys. Rev. E 84 (2011) 051916.
  * [5] B. Fernandez and L.S. Tsimring, _Typical trajectories of coupled degrade-and-fire oscillators: from dispersed populations to massive clustering_ , J. Math. Bio. 68 (2014) 1627–1652.
  * [6] W. Gerstner and N. Kistler, _Skiping neurons models, Single neurons, populations, plasticity_ , Cambridge University Press (2002).
  * [7] D.S. Glass, X. Jin, and I.H. Riedel-Kruse, _Signaling Delays Preclude Defects in Lateral Inhibition Patterning_ , Phys. Rev. Lett. 116 (2016) 128102.
  * [8] F. Golse, _On the dynamics of large particle systems in the mean-field limit_ , in Macroscopic and large scale phenomena: coarse graining, mean field limits and ergodicity, Lect. Notes App. Math. and Mechanics 3 (2016) 1–144.
  * [9] H. Huijberts, W. Michiels and H. Nijmeijer, _Stabilizability via time-delayed feedback: An eigenvalue optimization approach_ , SIAM J. App. Dynam. Sys. 8 (2009) 1–20.
  * [10] M.H. Jensen, K. Sneppen and G. Tiana, _Sustained oscillations and time delays in gene expression of protein Hes1._ FEBS Lett. 541 (2003) 176–177.
  * [11] T. Jüngling, I. Fischer, E. Schöll and W. Just, _Synchronization of heterogeneous oscillators by noninvasive time-delayed cross coupling._ Phys. Rev. Lett. 115 (2015) 194101.
  * [12] W. Just, _Synchronization of non-identical systems by non-invasive mutual time-delayed feedback_ , Chaos 33 (2023) 033105.
  * [13] K. Yoshioka-Kobayashi, M. Matsumiya, Y. Niino et al., _Coupling delay controls synchronized oscillation in the segmentation clock_ , Nature 580 (2020) 119–123.
  * [14] R. Kageyama, A. Isomura and H. Shimojo, _Biological Significance of the Coupling Delay in Synchronized Oscillations_. Physiology 38 (2) (2023) 63–72.
  * [15] D. Lipshutz and R.J. Lipshutz, _Stability of Synchronous Slowly Oscillating Periodic Solutions for Systems of Delay Differential Equations with Coupled Nonlinearity_ , J. Dyn. Diff. Equat. 34 (2022) 2259–2314.
  * [16] A. Lezia, A. Miano and J. Hasty, _Synthetic Gene Circuits: Design, Implement, and Apply_ , Proc. of the IEEE 110 (2022) 613–630.
  * [17] W. Mather, M.R. Bennett, J. Hasty, and L.S. Tsimring, _Delay-Induced Degrade-and-Fire Oscillations in Small Genetic Circuits_ , Phys. Rev. Lett. 102 (2009) 068105.
  * [18] W. Mather, J. Hasty, and L.S. Tsimring, _Synchronization of Degrade-and-Fire Oscillations via a Common Activators_ , Phys. Rev. Lett. 113 (2014) 128102.
  * [19] H. G. Schuster, and P. Wagner, _Mutual entrainment of two limit cycle oscillators with time delayed coupling_. Prog. Theor. Phys. 81 (1989) 939–945.