Hafiz Ahmed
email: hafiz.ahmed@inria.fr

Rosane Ushirobira

Denis Efimov

On Robustness of Phase Resetting to Cell Division under Entrainment

Keywords: Oscillation control, Phase resetting, Cell division

The problem of phase synchronization for a population of genetic oscillators (circadian clocks, synthetic oscillators, etc.) is considered in this paper, taking into account a cell division process and a common entrainment input in the population. The proposed analysis approach is based on the Phase Response Curve (PRC) model of an oscillator (the first order reduced model obtained for the linearized system and inputs with infinitesimal amplitude). The occurrence of cell division introduces state resetting in the model, placing it in the class of hybrid systems. It is shown that without common entraining input in all oscillators, the cell division acts as a disturbance causing phase drift, while the presence of entrainment guarantees boundedness of synchronization phase errors in the population. The performance of the obtained solutions is demonstrated via computer experiments for two different models of circadian/genetic oscillators (Neurospora's circadian oscillation model and the repressilator).

I. INTRODUCTION

The interest in the analysis and synthesis of genetic oscillators is continuously growing these last decades [START_REF] Guevara | Phase locking, period-doubling bifurcations, and irregular dynamics in periodically stimulated cardiac cells[END_REF], [START_REF] Kuramoto | Chemical Oscillations, Waves and Turbulence[END_REF], [START_REF] Tass | Phase Resetting in Medicine and Biology. Stochastic Modeling and Data Analysis[END_REF], [START_REF] Winfree | The Geometry of Biological Time[END_REF]. Any periodic oscillation is characterized by its frequency (or frequency spectrum), phase and amplitude. The amplitude and frequency are mainly governed by external stimulus applied to oscillators, a phenomenon called entrainment [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting[END_REF], [START_REF] Pikovsky | Synchronization. A Universal Concept in Nonlinear Sciences[END_REF], while the phase value is dependent on properties of the oscillator and characteristics of entrainment. This phase feature has attracted the attention of many researchers and in particular, the phase synchronization phenomenon studies are very popular [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting[END_REF], [START_REF] Pikovsky | Synchronization. A Universal Concept in Nonlinear Sciences[END_REF]. Phase synchronization is frequently observed in networks of oscillators, like a colony of the smallest free-living eukaryotes [START_REF] Thommen | Robustness of circadian clocks to daylight fluctuations: Hints from the picoeucaryote ostreococcus tauri[END_REF], the mammalian circadian pacemaker neural network [START_REF] Antle | Gates and oscillators II: zeitgebers and the network model of the brain clock[END_REF], [START_REF] Zhao | Phase organization of circadian oscillators in extended gate and oscillator models[END_REF] or networks of neural oscillators [START_REF] Smeal | Phase-response curves and synchronized neural networks[END_REF], [START_REF] Tass | Phase Resetting in Medicine and Biology. Stochastic Modeling and Data Analysis[END_REF], [START_REF] Canavier | Pulse coupled oscillators and the phase resetting curve[END_REF], to mention a few. Controlled phase resetting has been studied in [START_REF] Bagheri | Circadian phase entrainment via nonlinear model predictive control[END_REF], [START_REF] Danzl | Spike timing control of oscillatory neuron models using impulsive and quasi-impulsive charge-balanced inputs[END_REF], [START_REF] Efimov | Controlling the phase of an oscillator: a phase response curve approach[END_REF], [START_REF] Efimov | Phase resetting control based on direct phase response curve[END_REF] and for a population of oscillators in [START_REF] Efimov | Phase resetting for a network of oscillators via phase response curve approach[END_REF].

A simple but effective approach for analysis of phase resetting and dynamics for a single oscillator is based on PRC [START_REF] Glass | Predicting the entrainment of reentrant cardiac waves using phase resetting curves[END_REF], [START_REF] Govaerts | Computation of the phase response curve: A direct numerical approach[END_REF], [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting[END_REF], [START_REF] Sacré | Sensitivity analysis of oscillator models in the space of phase-response curves: oscillators as open systems[END_REF]. The PRC map is calculated for the system linearized around the limit cycle and inputs with small amplitudes. If the entraining input is a series of pulses, then a Poincaré phase map based on PRC can be calculated to predict the phase behavior [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting[END_REF]. Such a reduced phase model has been used in [START_REF] Efimov | Controlling the phase of an oscillator: a phase response curve approach[END_REF], [START_REF] Efimov | Phase resetting for a network of oscillators via phase response curve approach[END_REF] for pulse amplitude and timing calculation for a controlled phase resetting.

Another interesting problem that has emerged recently in [START_REF] Gonze | Modeling the effect of cell division on genetic oscillators[END_REF], deals with the influence of cell division on the behavior of genetic oscillators. It has been observed that oscillations persist across cell divisions in Repressilator [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF], similarly for circadian oscillations in cyanobacteria cells [START_REF] Mihalcescu | Resilient circadian oscillator revealed in individual cyanobacteria[END_REF]. In [START_REF] Nagoshi | Circadian gene expression in individual fibroblasts: Cell-autonomous and self-sustained oscillators pass time to daughter cells[END_REF], the persistence of circadian oscillations in culture fibroblasts under cell division has been demonstrated, and it has been noted that cell division can shift the phase in circadian cycle. A rapid phase decorrelation between daughter cells has been remarked in [START_REF] Geva-Zatorsky | Oscillations and variability in the p53 system[END_REF] for oscillations in the p53/Mdm2 system. Since cell division introduces a discontinuity in the oscillator dynamics (that is usually described by a system of nonlinear differential equations), then the analysis of division influence leads to the study of a hybrid or impulsive nonlinear oscillating system, which is a rather complicated problem [START_REF] Churilov | Periodical solutions in a pulse-modulated model of endocrine regulation with time-delay[END_REF], [START_REF] Efimov | On existence of oscillations in hybrid systems[END_REF]. In [START_REF] Gonze | Modeling the effect of cell division on genetic oscillators[END_REF], this problem has been investigated using a stochastic simulation approach, and in [START_REF] Tourigny | Geometric phase shifts in biological oscillators[END_REF], the geometric phase approach has been adopted from quantum mechanics.

The goal of the present work is to analyze the phase behavior and synchronization under cell division in genetic oscillators using PRC formalism. A motivating example given by a simple biological model of circadian oscillations in Neurospora, is studied in Section II. The analysis of cell division influence on the phase dynamics is presented in Section III. An illustration by simulations of the obtained results is given in Section IV. General results about phase dynamics are summarized in the Appendix.

II. MOTIVATING EXAMPLE

Let us consider a simple biological model of circadian oscillations in Neurospora in the following form [START_REF] Leloup | Limit cycle models for circadian rhythms based on transcriptional regulation in drosophila and neurospora[END_REF]:

Ṁ (t) = (v s + u(t)) K n I K n I + F n N (t) -v m M (t) K m + M (t)
,

ḞC (t) = k s M (t) -v d F c (t) K d + F c (t) -k 1 F c (t) + k 2 F n (t), ḞN (t) = k 1 F c (t) -k 2 F N (t), (1)
where M (t), F c (t) and F N (t) are the concentrations (defined with respect to the total cell volume) of the frqmRNA, the cytosolic and nuclear forms of FRQ, respectively. The parameter v s defines the rate of frq transcription (this parameter increases in the light phase) while the influence of light (the external entraining input in the model ( 1)) is denoted by u(t) ≥ 0. A description of the other parameters appearing in these equations can be found in [START_REF] Leloup | Limit cycle models for circadian rhythms based on transcriptional regulation in drosophila and neurospora[END_REF]. The following values of parameters are proposed there:

v m = 0.505, v d = 1.4, k s = 0.5, k 1 = 0.5, k 2 = 0.6, K m = 0.5, K I = 1, K d = 0.13, n = 4 and 1 ≤ v s + u(t) ≤ 2.5.
For all these values, the system (1) for u(t) = 0 has single unstable equilibrium and globally attractive limit cycle that represents a rhythmic behavior of the circadian rhythm in Neurospora with a period T > 0. It is a continuous-time dynamical system that for any initial conditions M (0) > 0, F C (0) > 0 and F N (0) > 0 has a continuous positive solution for all t ≥ 0. To model the cell division in [START_REF] Guevara | Phase locking, period-doubling bifurcations, and irregular dynamics in periodically stimulated cardiac cells[END_REF], it is necessary to introduce an increasing series of time instants t k > 0, k = 1, 2, . . . with a division at each t k . During the division, the state variables are resetted [START_REF] Gonze | Modeling the effect of cell division on genetic oscillators[END_REF], i.e.

M (t + k ) = λ M k M (t k ), F C (t + k ) = λ F C k F C (t k ) and F N (t + k ) = λ F N k F N (t k ), where M (t + k ) is the value of the concentration M after division at instant t k ; λ M k > 0, λ F C k > 0 and λ F N k
> 0 are parameters. The cell division cycle can be larger than the period of oscillations T [START_REF] Tourigny | Geometric phase shifts in biological oscillators[END_REF] or similar, as in proliferating human cells [START_REF] Bernard | Why do cells cycle with a 24 hour period?[END_REF] (where the circadian clock is a major synchronizing factor, which orchestrates daily rhythms regulating the cell division cycle), or two times faster as in cyanobacteria [START_REF] Mori | Circadian gating of cell division in cyanobacteria growing with average doubling times of less than 24 hours[END_REF].

The values λ M k , λ F C k , λ F N k
have been selected around 0.5 in [START_REF] Gonze | Modeling the effect of cell division on genetic oscillators[END_REF] (for the Goodwin model), but in [START_REF] Cookson | Cell cycledependent variations in protein concentration[END_REF] it has been observed in vivo that concentrations do not jump significantly after cell division. In the present work, we will accept the latter hypothesis taking

λ M k , λ F C k , λ F N k close to 1.
The modeling of such a hybrid oscillator corresponds to a mother cell in the population, then after each division the daughter cells have a similar dynamics and forthcoming divisions augmenting the population. It is assumed that division instants t k for each cell are different, then the phase synchronization behavior in a population (assuming that there is no interconnection between cells) can be analyzed using [START_REF] Guevara | Phase locking, period-doubling bifurcations, and irregular dynamics in periodically stimulated cardiac cells[END_REF]. If the phase converges to a steady-state in this hybrid system under some conditions, then the population will be phase synchronized in some sense.

Taking the previously mentioned parameter values and v s = 1.11, the period of the autonomous oscillation of ( 1) is obtained as T = 19.25 min. For these values of parameters and for the case u(t) = 0 and t k = kT -υ k , k ≥ 1, where υ k ∈ [0.15T, 0.30T ] is a uniformly distributed random variable, the results of the Neurospora's circadian oscillation model simulation for the same initial conditions and different realizations of υ k for 4 different cells undergoing divisions can be seen in Fig. 2. As we can conclude from these results the phase is diverging as it has been noted in [START_REF] Nagoshi | Circadian gene expression in individual fibroblasts: Cell-autonomous and self-sustained oscillators pass time to daughter cells[END_REF], [START_REF] Geva-Zatorsky | Oscillations and variability in the p53 system[END_REF] and in some experiments of [START_REF] Gonze | Modeling the effect of cell division on genetic oscillators[END_REF]. Next, by taking u(t) = max {0, 0.2 sin(ωt)} as the common external entraining input and repeating the same experiments, the results are given in Fig. 1. From Fig. 1, it is evident that the oscillations converge to a common entrained mode.

In this paper, we will try to find conditions providing both these two types of phase behavior using the PRC phase model for small inputs. This Section begins with the introduction of the formalized problem statement. Next, the reduced PRC model is introduced and the phase synchronization analysis is presented.

Details of the standard procedure for a phase model derivation for an oscillator can be found in [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting[END_REF], [START_REF] Pikovsky | Synchronization. A Universal Concept in Nonlinear Sciences[END_REF], [START_REF] Efimov | Phase resetting for a network of oscillators via phase response curve approach[END_REF] and they are briefly summarized in the Appendix.

A. Problem statement

Let us consider a population of N > 0 cells (genetic oscillators) with s = 1, 2, . . . , N :

ẋs (t) = f s (x s (t), u(t)) t ∈ [t s,k , t s,k+1 ), k ≥ 0; (2) x s (t + s,k+1 ) = Λ s,k,n x s (t s,k+1 ) k ≥ 1, where x s (t) ∈ R n is the state (concentrations of different products) of the s th cell and the input u(t) is a periodical train of pulses u(t) = +∞ i=0 w(t -iT),
with a pulse w(t), w(t) = 0 for all t ≥ T or t ≤ 0, sup 0≤t≤T |w(t)| ≤ W < +∞ and T < T where T > 0 is the period of u; t s,0 = 0 and t s,k , k ≥ 0 is a strictly increasing sequence of impulses (discontinuous jumps in ( 2)) for all

s = 1, 2, . . . , N , Λ s,k,n = diag[λ s,k,1 , . . . , λ s,k,n ] ∈ R n×n with λ s,k,1 ∈ [1 -ε, 1]
for some ε > 0 sufficiently small. The periodical input u(t) models the common entrainment for all cells and the discontinuities at instants t s,k , k ≥ 1 represent the cell division, the diagonal matrix Λ s,k,n determines changes in the state vector (in concentrations) after division. The instants of division t s,k and the concentration changes Λ s,k,n may be different for each cell. Remark 1. Note that, formally, at each t s,k , k ≥ 1, the population should be augmented by a daughter cell, that has dynamics similar to mother one. Then, the number N is continuously growing. In the present work, we will consider a fixed size of the population N , since as it will be shown below (and due to similarity of dynamics for newborn cells and initial conditions), the problem of phase synchronization can be analyzed using the model even for a single cell.

It is necessary to establish conditions (restrictions on f s , u and t s,k or Λ s,k,n ) under which there exists a synchronization phenomenon in the cell population (2).

B. Reduced phase model under cell division

The presence of divisions can be alternatively modeled by an additional impulsion input:

ẋs (t) = f s (x s (t), u(t)) + +∞ k=1 χ s,k δ(t -t s,k ), (3)
where

δ(t) is a delta-impulse function, χ s,k = (Λ s,k,ns -I ns ) x s (t s,k+1
) and I ns is the identity matrix of dimension n s , s = 1, 2, . . . , N . Assumption 1. For each s = 1, 2, . . . , N , the s th subsystem in (2), with u = 0 and Λ s,k,ns = I ns for all k ≥ 0, has a limit cycle Γ s ⊂ A s with an open set of attraction A s ⊂ R ns , and with period T s > 0, ω s = 2πT -1 s . This assumption says if there is no entrainment u and cell division, then each cell in the population is an oscillator with the limit cycle in Γ s and period T s . Under Assumption 1 and using the theory presented in Appendix, for each cell in (2) it is possible to define its asymptotic phase θ s ∈ [0, 2π). Under additional restrictions that ε and W are sufficiently small, we can design a phase dynamical model of (2) in some vicinity of Γ s as in Appendix. Since the model derived in [START_REF] Canavier | Pulse coupled oscillators and the phase resetting curve[END_REF] is based on the first order approximation and in the system (3) there are two inputs (u and the train of impulses), by superposition principle, eq. ( 11) takes the form in this case for s = 1, 2, . . . , N :

θs = ω s + Q s (t + θ s,0 ω -1 s ) b s (t + θ s,0 ω -1 s ) u(t) +Q s (t + θ s,0 ω -1 s ) +∞ k=1 χ s,k δ(t -t s,k ), (4)
where

b s (t + θ s,0 ω -1 s ) = ∂f s (x s , u) ∂u xs=γs(t+θs,0ω -1 s )
and γ s (t + θ s,0 ω -1 s ) is a trajectory of the s th cell in (2) for u = 0 and Λ s,k,ns = I ns for all k ≥ 0 with initial conditions in Γ s with the initial phase θ s,0 ∈ [0, 2π], Q s (t) is the infinitesimal PRC derived in the Appendix. This model is constructed around the base trajectory γ s (t + θ s,0 ω -1 s ) under the assumption that the perturbed trajectory with u = 0 and Λ s,k,ns = I ns stays close to that one [START_REF] Efimov | Phase resetting for a network of oscillators via phase response curve approach[END_REF]. Since such a closeness assumption is rather restrictive and may be invalid on a sufficiently long time interval (the excited trajectory can belong to a small vicinity of Γ s for sufficiently small ε and W , but moving away from γ s (t+θ s,0 ω -1 s ) due to a phase shift induced by external inputs), then it is better to recalculate the phase of base trajectory γ s (t + θ s,0 ω -1 s ) after a period T, for example (that is the idea of Poincaré phase map approach [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting[END_REF]). In this case, by recurrent integration of (4) (similarly as for ( 12)) the phase shift over the interval [iT, (i + 1)T] can be evaluated as follows:

θ s,i+1 = ω s (T -T s ) + θ s,i + PRC s (θ s,i ) + ∆ s,i , (5) PRC s (θ) = ˆT 0 Q s (t + θω -1 s ) b s (t + θω -1 s ) w(t)dt, ∆ s,i = k∈Ks,i χ s,k Q s (t s,k + θ s,i ω -1 s )
for all s = 1, 2, . . . , N , where

K s,i = {k ≥ 1 | t s,k ∈ [iT, (i + 1)T]} is the set of indexes whose impulses happen in the interval [iT, (i + 1)T], PRC s : [0, 2π) → [0, 2π
) is the PRC of the s th oscillator for the pulse w and ∆ s,i ∈ R is the phase perturbation imported by cell division on the interval [iT, (i + 1)T].

Remark 2. Formally, the set K s,i can be decomposed on two parts:

K s,i = K 1 s,i ∪ K 2 s,i , K 1 s,i = {k ≥ 1 | t s,k ∈ [iT, iT + T ]} , K 2 s,i = [iT + T , (i + 1)T],
where K 1 s,i characterizes the impulses arrived for u(t) = 0 and K 2 s,i is for u(t) = 0 on the interval [iT, (i + 1)T]. Then the model ( 12) can be rewritten as follows:

θ s,i+1 = ω s (T -T s ) + θ s,i + PRC s (θ s,i + ∆ 1 s,i ) + ∆ s,i , ∆ s,i = ∆ 1 s,i + ∆ 2 s,i , ∆ j s,i = k∈K j s,i χ s,k Q s (t s,k + θ s,i ω -1 s ), j = 1, 2.
The difference with respect to [START_REF] Bagheri | Circadian phase entrainment via nonlinear model predictive control[END_REF] is that the perturbation caused by cell division appears nonlinearly in the last model. For brevity of consideration only the case of ( 12) is studied below.

C. Phase synchronization

The model [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting[END_REF] for each s = 1, 2, . . . , N is a scalar nonlinear integrator-like discrete-time system (that is a considerable advantage with respect to (2)) with the state θ s,i and external input ω s (T-T s )+∆ s,i , where the constant part represents the influence of entrainment and ∆ s,i is the perturbation originated by cell division.

Assume that there is no common entrainment and u(t) = 0, then the model ( 5) can be simplified to a pure integrator of the input T∆ -1 s,i , which is an average influence of pulses other interval [iT, (i + 1)T]:

θ s,i+1 = θ s,i + ∆ s,i .
If ∆ s,i are different for each s = 1, 2, . . . , N and have not a zero mean, then the phase θ s,i will be drifting in a unique manner for each s = 1, 2, . . . , N . Thus, there is no phase synchronization. This is the case presented in Fig. 2 of Section II and also observed in [START_REF] Geva-Zatorsky | Oscillations and variability in the p53 system[END_REF], [START_REF] Nagoshi | Circadian gene expression in individual fibroblasts: Cell-autonomous and self-sustained oscillators pass time to daughter cells[END_REF].

Therefore, the synchronous properties of ( 5) depend critically on the nonlinear function PRC s . In this work, as in [START_REF] Danzl | Spike timing control of oscillatory neuron models using impulsive and quasi-impulsive charge-balanced inputs[END_REF], [START_REF] Efimov | Controlling the phase of an oscillator: a phase response curve approach[END_REF], [START_REF] Efimov | Phase resetting for a network of oscillators via phase response curve approach[END_REF], we assume that the PRC map has particular properties (it is similar to type II PRC from [START_REF] Hansel | Synchrony in excitatory neural networks[END_REF]).

Assumption 2. For all s = 1, 2, . . . , N , the map PRC s is continuously differentiable and there exist 0 < β s ≤ 1 and

Θ s ∈ [0, 2π) such that the equation PRC s (θ 0 s ) = ω s (T s -T) has a solution θ 0 s ∈ [0, 2π) with -2 + β s ≤ ∂ PRC s (θ) ∂θ ≤ -β s ∀θ ∈ [θ 0 s -Θ s , θ 0 s + Θ s ].
Obviously, for ∆ s,i = 0 (no cell division), θ 0 s corresponds to a stable equilibrium of the system (5) for given s with the domain of attraction

[θ 0 s -Θ s , θ 0 s + Θ s ] [5]. Lemma 3. For each s = 1, 2, . . . , N , under the Assumption 2, if |θ s,0 -θ 0 s | ≤ Θ s -β -1 s ∆ s where ∆ s = sup i≥0 |∆ s,i | < +∞, then |θ s,i -θ 0 s | ≤ Θ s ∀i ≥ 0, lim i→+∞ |θ s,i -θ 0 s | ≤ β -1 s ∆ s . (6)
Proof. Taking into account definition of θ 0 s , the model ( 5) can be rewritten as follows:

θ s,i+1 = θ s,i + PRC s (θ s,i ) -PRC s (θ 0 s ) + ∆ s,i .
Using the Mean value theorem (PRC s is continuously differentiable by Assumption 2), we obtain:

θ s,i+1 = θ s,i + ∂ PRC s (θ s,i ) ∂θ s,i (θ s,i -θ 0 s ) + ∆ s,i ,
where θ s,i = µθ s,i + (1 -µ)θ 0 s for some µ ∈ [0, 1]. Define the phase error e s,i = θ s,i -θ 0 s , then

e s,i+1 = 1 + ∂ PRC s (θ s,i ) ∂θ s,i e s,i + ∆ s,i . By Assumption 2, -2 + β s ≤ ∂ PRCs(θ s,i ) ∂θ s,i ≤ -β s provided that |e s,i | ≤ Θ s . Then taking Lyapunov function V (e) = |e| we have: V (e s,i+1 ) -V (e s,i ) = 1 + ∂ PRC s (θ s,i ) ∂θ s,i e s,i + ∆ s,i -|e s,i | ≤ 1 + ∂ PRC s (θ s,i ) ∂θ s,i -1 |e s,i | + |∆ s,i | ≤ -β s V (e s,i ) + |∆ s,i |, that implies |e s,i | ≤ (1 -β s ) i |e s,0 | + β -1 s ∆ s i ≥ 0
under assumption that |e s,i | ≤ Θ s for all i ≥ 0. However, if |e s,0 | ≤ Θ s -β -1 s ∆ s then |e s,i | ≤ Θ s for all i ≥ 0 as needed, and

lim i→+∞ |e s,i | ≤ β -1 s ∆ s .
Consequently, if the influence of cell division quantified by ∆ s is sufficiently small and the initial phase θ s,0 lies sufficiently close to θ 0 s , then the phase θ s,i stays in the domain of attraction of θ 0 s and asymptotically converges to a vicinity of that equilibrium. Since all cells in population (2) yield this kind of behavior, then under these conditions, the phases are asymptotically synchronized with the error of synchronization proportional to superposition of max{β -1 s1 ∆ s1 , β -1 s2 ∆ s2 } and |θ 0 s1 -θ 0 s2 | for any 1 ≤ s 1 = s 2 ≤ N . Theorem 4. Let assumptions 1, 2 be satisfied and ε, W be sufficiently small in (2).

If |θ s,0 -θ 0 s | ≤ Θ s -β -1 s ∆ s for all s = 1, 2, . . . , N , then for any 1 ≤ s 1 = s 2 ≤ N |θ s1,i -θ s2,i | ≤ |θ 0 s1 -θ 0 s2 | + Θ s1 + Θ s2 ∀i ≥ 0, lim i→+∞ |θ s1,i -θ s2,i | ≤ |θ 0 s1 -θ 0 s2 | + β -1 s1 ∆ s1 + β -1 s2 ∆ s2 .
Proof. If Assumption 1 holds and ε, W are sufficiently small, then the results presented in Section V imply that a first order approximation of (2) can be used for analysis of the population behavior, and the reduced PRC model ( 5) can be derived for each s = 1, 2, . . . , N . Next, since all conditions of Lemma 3 are satisfied, then the relations ( 6) are valid for all s = 1, 2, . . . , N . Consider the phase difference θ s1,i -θ s2,i of any two oscillators with

1 ≤ s 1 = s 2 ≤ N . Since θ s1,i -θ s2,i = θ s1,i -θ 0 s1 -θ s2,i + θ 0 s2 + θ 0 s1 -θ 0 s2 from (6)
:

|θ s1,i -θ s2,i | ≤ |e s1,i -e s2,i | + |θ 0 s1 -θ 0 s2 | ≤ |θ 0 s1 -θ 0 s2 | + Θ s1 + Θ s2 and lim i→+∞ |θ s1,i -θ s2,i | ≤ lim i→+∞ |e s1,i -e s2,i | + |θ 0 s1 -θ 0 s2 | ≤ |θ 0 s1 -θ 0 s2 | + β -1 s1 ∆ s1 + β -1 s2 ∆ s2
as required.

This theorem establishes phase-lock behavior in the population (2), which may be composed by different cells. If all cells are identical, then the following synchronization conditions can be obtained.

Corollary 5. Let all conditions of Theorem 4 be satisfied and PRC s (θ) = PRC(θ) for all s = 1, 2 . . . N and all θ ∈ [0, 2π) (then also θ 0 s = θ 0 , Θ s = Θ and

β s = β). If |θ s,0 -θ 0 | ≤ Θ -β -1 ∆ s for all s = 1, 2 . . . N , then |θ i -θ 0 | ≤ Θ, ∀i ≥ 0, lim i→+∞ |θ i -θ 0 | ≤ 1 N β N s=1 ∆ s ,
where θ i = N -1 N s=1 θ s,i is the average phase of the population.

Proof. The result follows from Lemma 3 and Theorem 4 under assumption that all PRC s are identical. If there is no cell division, then

1 N β N s=1 ∆ s = 0
and we recover a well-known result on phase synchronization under a periodical entrainment [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting[END_REF], [START_REF] Efimov | Phase resetting for a network of oscillators via phase response curve approach[END_REF].

IV. EXAMPLE

Let us illustrate the theoretical findings obtained in the previous Section.

A. Circadian oscillations in Neurospora

Consider a population of circadian oscillators in Neurospora [START_REF] Guevara | Phase locking, period-doubling bifurcations, and irregular dynamics in periodically stimulated cardiac cells[END_REF]. Take all cells in the population identical with the values of parameters given in Section II. Then the Assumption 1 is satisfied for T s = T = 19.25 minutes (ω = 2πT -1 ). Select

w(t) = max{0, W sin(ωt)} if 0 ≤ t ≤ T 0 otherwise
with W = 0.02. Setup the same Λ s,k,n for all s and n as Λ k defined by: Λ k = diag[0.99 0.98 0.98]

with t k = iT -υ k , k ≥ 1, where υ k ∈ [0.15T, 0.
30T ] is a uniformly distributed random variable as before. The values ε, W are chosen sufficiently small. For this pulse w(t), the obtained PRC(θ) and PRC (θ) = ∂ PRC(θ) ∂θ are shown in Fig. 3. From these plots, θ 0 = 2.34, Θ = 1.04 and β = 0.05, and the Assumption 2 is also satisfied. Thus, all conditions of Theorem 4 are verified.

Simulated phase behavior of ( 1) is shown in Fig. 4 by the blue curve (the phase value was computed by finding the closest point on the limit cycle at instants iT for i ≥ 0 and by assigning the phase of that point as θ i ). The values of phase obtained by the model [START_REF] Bagheri | Circadian phase entrainment via nonlinear model predictive control[END_REF] are presented in the same figure by the red curve. As we can see, both curves are very close and that confirms all theoretical developments presented in this work, and phase asymptotically converges to a vicinity of θ 0 , then synchronization of phase would be observed for a population of circadian oscillators of Neurospora (2) as in [START_REF] Gonze | Modeling the effect of cell division on genetic oscillators[END_REF] for the Godwin model.

B. The Repressilator

The repressilator [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF] is a very simple genetic oscillator consisting of three genes, which can be modeled as below [START_REF] Gonze | Modeling the effect of cell division on genetic oscillators[END_REF]:

Ṁ1 (t) = α 0 + α 1 K n K n + P n 3 (t)
-δM 1 (t) + u(t)

Ṁ2 (t) = α 0 + α 1 K n K n + P n 1 (t) -δM 2 (t) Ṁ3 (t) = α 0 + α 1 K n K n + P n 2 (t) -δM 3 (t) (7)
Ṗ1 (t) = βM 1 (t) -γP 1 (t) Ṗ2 (t) = βM 2 (t) -γP 2 (t) Ṗ3 (t) = βM 3 (t) -γP 3 (t)
where the variables M i (t) and P i (t) (with i = 1, 2, 3) represent the concentrations of mRNA and protein of the three components of the repressilator respectively and u(t) ≥ 0 represents the external entraining input. Details about other parameters of the model can be found in [START_REF] Gonze | Modeling the effect of cell division on genetic oscillators[END_REF]. We will consider the following values for [START_REF] Thommen | Robustness of circadian clocks to daylight fluctuations: Hints from the picoeucaryote ostreococcus tauri[END_REF]:

α 1 = 1, α 0 = 0.01, K = 1, n = 2 and δ = β = γ = 0.1.
With these values, this model has a single equilibrium and one limit cycle and the period of autonomous oscillation is obtained as T = 116.6 minutes. Now, let us consider a population of identical repressilators with the previously mentioned parameters. The Assumption 1 is satisfied for

T s = T = 116.6 min ω = 2πT -1 . Select w(t) = max{0, W sin(ωt)} if 0 ≤ t ≤ T 0 otherwise
with W = 0.002. As seen previously, set:

Λ k = diag[0.99 0.98 0.98] with t k = iT -υ k , k ≥ 1, where υ k ∈ [0.15T, 0.
30T ] is a uniformly distributed random variable as before. The values ε, W are chosen sufficiently small. For this pulse w(t) the obtained PRC(θ) and PRC (θ) = ∂ PRC(θ) ∂θ are shown in Fig. 5, from these plots θ 0 = 3.9 ( which is -2.4 in a scale between -π to π), Θ = 1.6 and β = 0.01, and Assumption 2 is also satisfied. Thus, all conditions of Theorem 4 are verified.

Simulated phase behavior of ( 7) is shown in Fig. 6 by the blue curve (the phase value was computed by finding the closest point on the limit cycle at instants iT for i ≥ 0 and by assigning the phase of that point as θ i ). The values of phase obtained by the model ( 12) are presented in the same figure by the red curve. As we can see, both curves are very close and that confirms all theoretical developments presented in this work, and phase asymptotically converges to a vicinity of θ 0 , then synchronization of phase would be observed for a population of repressilator model (2) as in [START_REF] Gonze | Modeling the effect of cell division on genetic oscillators[END_REF].

V. CONCLUSION

The influence of cell division on the dynamics of a population of genetic oscillators is analyzed. As it has been observed in vivo [START_REF] Geva-Zatorsky | Oscillations and variability in the p53 system[END_REF], [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF], [START_REF] Mihalcescu | Resilient circadian oscillator revealed in individual cyanobacteria[END_REF], [START_REF] Nagoshi | Circadian gene expression in individual fibroblasts: Cell-autonomous and self-sustained oscillators pass time to daughter cells[END_REF], oscillations in cells are frequently quite resilient to cell division. Recently, this phenomenon has been analyzed by a stochastic simulation in [START_REF] Gonze | Modeling the effect of cell division on genetic oscillators[END_REF], where phase synchronization in the population has been observed. In the present work (modeling cell division by impulses places the dynamics of population in the class of hybrid systems), analytical conditions are established of phase synchronization applying PRC model approach for small inputs. The results are illustrated by numerical experiments with two different circadian/genetic oscillator models.

APPENDIX: PHASE MODEL IN VICINITY OF A LIMIT CYCLE

A. Linearized model Consider a (smooth) dynamical system

ẋ = f (x, u), x ∈ R n , u ∈ [-U, U ] ⊂ R, U > 0. (8)
Denote by x(t, x 0 , u) a solution of ( 8) with the initial condition x 0 and input u and assume that for u(t) ≡ 0, t ≥ 0 and some x 0 ∈ R n the system ( 8) has (non-constant) T -periodic solution x(t, x 0 , 0) = γ(t) = γ(t + T ) ∈ R n , t ≥ 0. Then the corresponding limit cycle, described by the set Γ = {x ∈ R n | x = γ(t), 0 ≤ t < T }, attracts a non-empty open bounded set of initial conditions A ⊂ R n , Γ ⊂ A, and the linearized system δ ẋ(t) = A(t) δx(t) + b(t)u(t) + d (δx(t), γ(t), u(t)) , ( 9)

A(t) = ∂f (x, u) ∂x x=γ(t) , b(t) = ∂f (x, u) ∂u x=γ(t)
has n -1 multipliers strictly inside the unit cycle and one multiplier equals to 1 [START_REF] Andronov | Theory of oscillators[END_REF], [START_REF] Yakubovich | Linear differential equations with periodic coefficients[END_REF], where δx(t) = x(t)-γ(t), the matrix function A and the vector function b are T -periodic due to properties of γ; the function d (δx(t), γ(t), u(t)) represents the higher order terms with respect to δx(t) in the system (8) linearization and for all x ∈ A and |u| ≤ U there exist d 1 > 0, d 2 > 0 such that (the function d contains products of δx and u with power 2 and higher):

|d(δx, γ, u)| ≤ d 1 |δx| 2 + d 2 u 2 .
Multipliers are the eigenvalues of the monodromy matrix M = Φ(T ) defined via the fundamental matrix function Φ of the system (9) and the solution of adjoint system Ψ:

Φ(t) = A(t)Φ(t), Φ(0) = I; Ψ(t) = -A(t) T Ψ(t), Ψ(0) = I,
where I is the identity matrix and Φ(t) T Ψ(t) = I.

B. Phase variables

Any point x 0 ∈ Γ can be characterized by a scalar phase ϕ 0 ∈ [0, 2π), that uniquely determines the position of the point x 0 on the limit cycle Γ (Γ is a one-dimensional closed curve in R n ) [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting[END_REF], [START_REF] Pikovsky | Synchronization. A Universal Concept in Nonlinear Sciences[END_REF]. The smooth bijective phase map ϑ : Γ → [0, 2π) assigns to each point x 0 ∈ Γ the corresponding phase ϕ 0 = ϑ(x 0 ). Any solution of the system (8) x(t, x 0 , 0) with x 0 ∈ Γ satisfies x(t, x 0 , 0) = γ(t + ϕ 0 ω -1 ), where ω = 2πT -1 is the system frequency, provided we choose the convention γ(t) = x(t, ϑ -1 (0), 0), then we can define ϑ -1 (ϕ) = γ(ϕω -1 ). The phase variable ϕ : R + → [0, 2π) is defined for the trajectories x(t, x 0 , 0), x 0 ∈ Γ as ϕ(t) = ϑ (x(t, x 0 , 0)) = ϑ γ(t + ϕ 0 ω -1 ) . Due to the periodic nature of γ(t), the function ϕ(t) is also periodic. Moreover the function ϑ can be defined providing ϕ(t) = ωt + ϕ 0 and φ(t) = ω [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting[END_REF], [START_REF] Pikovsky | Synchronization. A Universal Concept in Nonlinear Sciences[END_REF].

The notion of phase can be extended to any solution x(t, x 0 , 0) starting in the attraction set A of the limit cycle. By definition, for all x 0 ∈ A there exists an asymptotic phase θ 0 ∈ [0, 2π) such that lim t→+∞ |x(t, x 0 , 0) -γ(t + θ 0 ω -1 )| = 0.

Then there exists the asymptotic phase map υ : A → [0, 2π) connecting a point x 0 ∈ A and the corresponding phase θ 0 , i.e. θ 0 = υ(x 0 ) and by construction υ(x 0 ) = ϑ(x 0 ) for all x 0 ∈ Γ. The asymptotic phase variable θ : R + → [0, 2π) is derived as θ(t) = υ (x(t, x 0 , 0)), t ≥ 0. In the case ϕ(t) = ωt + ϕ 0 we have θ(t) = ωt + θ 0 and θ(t) = ω, which implies time invariance of this map: if υ(x 1 ) = υ(x 2 ), then υ (x(t, x 1 , 0)) = υ (x(t, x 2 , 0)) for all t ≥ 0 and x 1 , x 2 ∈ A [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting[END_REF]. The initial conditions x 1 , x 2 ∈ A having the same asymptotic phase determine the isochrone curves [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting[END_REF].

The notion of asymptotic phase variable can be extended to a generic u(t) = 0, t ≥ 0 provided that the corresponding trajectory x(t, x 0 , u) stays in the set A for all t ≥ 0. In this case the asymptotic phase variable can be defined in a trivial way as θ(t) = υ (x(t, x 0 , u)), t ≥ 0. Then the variable θ(t ) at an instant t ≥ 0 evaluates the asymptotic phase of the point x(t , x 0 , u) if one would pose u(t) = 0 for t ≥ t . Dynamics of the asymptotic phase variable θ(t) in the generic case for u(t) = 0, t ≥ 0 is hard to derive. A local model obtained in a small neighborhood of the limit cycle for infinitesimal inputs is presented below [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting[END_REF], [START_REF] Efimov | Phase resetting for a network of oscillators via phase response curve approach[END_REF].

C. Infinitesimal PRC

Consider the case u(t) = 0 for t ≥ 0, then by definition γ(t) = f (γ(t), 0), γ(t) = A(t) γ(t) and γ(t) = Φ(t) γ(0) for all t ≥ 0. Therefore, γ(0) = f (γ(0), 0) is the left eigenvector of the matrix M for the eigenvalue λ 1 (M ) = 1. There exists a right eigenvector m ∈ R n such that m T M = m T and m T γ(0) = ω. Finally, define

Q(t) = m T Ψ(t) T then Q(t)f (γ(t), 0) = m T Ψ(t) T f (γ(t), 0) = m T Ψ(t) T Φ(t) γ(0) = m T γ(0) = ω.
From another side, θ(t) = υ (γ(t)) = ωt + θ(0) and

ω = θ = ∂υ(x) ∂x x=γ(t) f (γ(t), 0). Therefore Q(t) = ∂υ(x) ∂x x=γ(t) + ζ(t)
, where ζ(t) is a rowvector orthogonal to f (γ(t), 0). Since m is the eigenvector corresponding to λ 1 (M ) = 1 (or movement on the limit cycle), then Q(t) = m T Ψ(t) T is independent of perturbations orthogonal to the limit cycle flow f (γ(t), 0) and the convention

Q(t) = ∂υ(x) ∂x x=γ(t) = m T Ψ(t) T (10)
is adopted. The first equality in [START_REF] Smeal | Phase-response curves and synchronized neural networks[END_REF] explains the physical meaning of Q(t), while the last equality in (10) is used for numerical calculation. The function Q(t) is T -periodic by construction.

The function Q(φω -1 ) for phase ϕ ∈ [0, 2π) is called infinitesimal PRC [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting[END_REF], it serves as the phase response characteristics for a delta-impulse input.

D. Phase dynamics

Consider the case u(t) = 0 (we assume that x(t, x 0 , u) ∈ A for all t ≥ 0), then θ(x(t)) = θ(γ(t) + δx(t)) This model has been derived around the solution γ(t), due to the periodicity of the solution γ(t + φω -1 ), φ ∈ [0, 2π) and u, the model for γ(t + φω -1 ) has a similar form [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting[END_REF], [START_REF] Efimov | Phase resetting for a network of oscillators via phase response curve approach[END_REF]:

=
θ = ω + Q(t + φω -1 ) b(t + φω -1 ) u(t)
+g γ(t + φω -1 ), δx(t), u(t) .

Skipping the residual function g we obtain the first order approximation of the phase model:

θ = ω + Q(t + φ ω -1 ) b(t + φω -1 ) u(t). (11)
Since the property |g(γ, δx, u)| ≤ g 1 u 2 + g 2 |δx| 2 holds for all x ∈ A and |u| ≤ U , such an approximation is rather accurate for a sufficiently small U . Assume that the input u(t) = w(t), where w(t) has a pulselike form, i.e. |w(t)| ≤ U for all 0 < t < T < T and w(t) = 0 for all t ≥ T or t ≤ 0. Then integration of (11) yields for

PRC(θ) = ˆT 0 Q(τ + θω -1 ) b(τ + θω -1 ) u(τ ) dτ.

The map PRC(θ), θ ∈ [0, 2π) is defined for the particular pulse w (by definition -π ≤ PRC(θ) < π for all θ ∈ [0, 2π)), it tabulates the phase shift by the pulse w [5], [START_REF] Pikovsky | Synchronization. A Universal Concept in Nonlinear Sciences[END_REF]. For the w(t) = δ(t -η) with η ∈ (0, T ), where δ(t) is impulse input we obtain the infinitesimal PRC iPRC(θ) = Q(η + θω -1 ) b(η + θω -1 ), which defines the phase shift under an impulse input.

Figure 1 .

 1 Figure 1. Oscillations of different single cells with cell divisions and common external entraining input

Figure 2 .

 2 Figure 2. Oscillations of different single cells with cell divisions and without any common input

Figure 3 .

 3 Figure 3. PRC(θ) and PRC'(θ) for the input w(t) for Neurospora model.

Figure 4 .

 4 Figure 4. Phase behavior of (1)

Figure 5 .

 5 Figure 5. PRC(θ) and PRC'(θ) for the input w(t) for the repressilator model.

Figure 6 .

 6 Figure 6. Phase behavior of[START_REF] Thommen | Robustness of circadian clocks to daylight fluctuations: Hints from the picoeucaryote ostreococcus tauri[END_REF]

t

  ≥ T :θ(t) = ωt + θ(0) + ˆt 0 Q(τ + θ(0)ω -1 ) b(τ + θ(0)ω -1 ) u(τ ) dτ = ωt + θ(0) + PRC[θ(0)],

  where the term r 1 (γ(t), δx(t), u(t))T δx(t) corresponds to the powers of δx(t) higher than one in the Taylor series of the function θ(γ(t) + δx(t)) with respect to the variable δx(t). From above, the quantity θ(t) = ω should be satisfied for u(t) = 0, therefore ∂υ(x) (γ(t), δx(t), 0) T δx(t) = ω, which implies the propertyr 1 (γ(t), δx(t), 0) = 0. Next, T δx + r 2 (γ, u) u 2 ,where r 2 (γ, u) u 2 represents the terms with powers two and higher for the Taylor series of the function

		∂υ(x) ∂x x=γ(t)+δx(t)	f (γ(t) + δx(t), u(t))
	=	∂υ(x) ∂x x=γ(t)	f (γ(t), u(t))
		+r 1 (γ(t), δx(t), u(t))	T δx(t),
				∂x x=γ(t)	f (γ(t), 0) +
	r 1 θ(x(t)) =	∂υ(x) ∂x x=γ(t)	f (γ(t), 0)
			+	∂υ(x) ∂x x=γ(t)	∂f (γ(t), u) ∂u
	∂υ(x)
	∂x x=γ(t)

u=0 u(t) +g (γ(t), δx(t), u(t)) , g(γ, δx, u) = r 1 (γ, δx, u) f (γ(t), u(t)) with respect to the control u. For all x ∈ A and |u| ≤ U there are g 1 > 0 and g 2 > 0 such that |g(γ, δx, u)| ≤ g 1 u 2 + g 2 |δx| 2 . Recalling the previously introduced designations, we obtain θ = ω + Q(t) b(t) u(t) + g (γ(t), δx(t), u(t)) .

This work was partially supported by the Government of Russian Federation (Grant 074-U01) and the Ministry of Education and Science of Russian Federation (Project 14.Z50.31.0031).

Hafiz Ahmed is partly supported by the regional council of Nord-Pas de Calais.


Franck Charmantray

Philippe Dellis

Virgil Hélaine

Soth Samreth

Laurence Hecquet

come

Chemoenzymatic synthesis of 5-thio-D-xylopyranose.