Population Dynamics of Globally Coupled Degrade-and-Fire Oscillators
aa r X i v : . [ m a t h . C A ] N ov Population Dynamics of Globally Coupled Degrade-and-FireOscillators
Alex Blumenthal and Bastien Fernandez ∗ July 31, 2018 Courant Institute of Mathematical SciencesNew York UniversityNew York, NY 10012, USA Laboratoire de Probabilit´es et Mod`eles Al´eatoiresCNRS - Universit´e Paris 7 Denis Diderot75205 Paris CEDEX 13 France
Abstract
This paper reports the analysis of the dynamics of a model of pulse-coupled oscillatorswith global inhibitory coupling. The model is inspired by experiments on colonies of bacteria-embedded synthetic genetic circuits. The total population can be either of finite (arbitrary)size or infinite, and is represented by a one-dimensional profile. Profiles can be discontinuous,possibly with infinitely many jumps. Their time evolution is governed by a singular differentialequation. We address the corresponding initial value problem and characterize the dynamics’main features. In particular, we prove that trajectory behaviors are asymptotically periodic,with period only depending on the profile (and on the model parameters). A criterion is obtainedfor the existence of the corresponding periodic orbits, which reveals the existence of a sharptransition as the coupling parameter is increased. The transition separates a regime where anyprofile can be obtained in the limit of large times, to a situation where only trajectories withsufficiently large groups of synchronized oscillators perdure.
To determine the amount of collective order and how its depends on parameters is a central questionin the analysis of systems of coupled oscillators [17]. A typical example is the Kuramoto model [11]where numerics have indicated that, while the oscillators evolve independently at weak coupling, assoon as the interaction strength exceeds a threshold, the overall collective behavior becomes moreand more coherent when the coupling increases [2]. This phenomenology has been identified sincethe late 1970’s. Yet, its full mathematical proof still remains to be achieved and only preliminaryresults have been obtained so far (see [20] for a summary of results and technical challenges).In the Kuramoto model, trajectories are smooth and most of the difficulties come from thepresence of heterogeneities ( i.e. the individual oscillators’ frequencies are randomly drawn), notto mention nonlinearities. However, in other circumstances, modelling results in systems with ∗ On leave from Centre de Physique Th´eorique, CNRS - Aix-Marseille Universit´e - Universit´e de Toulon, Campusde Luminy, 13288 Marseille CEDEX 09 France in vitro dynamics of assemblies of interacting self-repressor genes.The analysis here not only addresses populations of arbitrary finite size throughout the couplingparameter range, but also continua of infinite populations.As phenomenology is concerned, the critical feature of this model is the existence a sharptransition when the coupling strength increases. Below the threshold, every possible populationdistribution can be observed at large time, upon the choice of initial condition. Beyond that point,only distributions of grouped oscillators can persist - the stronger the interaction, the larger thegroups. Any initially isolated oscillator must eventually join, and must remain, in a group.These features have been identified and mathematically justified for finite size populations in[7, 8], except for the asymptotic periodicity of every trajectory (Theorem 5.3 below), whose proofin [7] turns out to be incomplete. Here, we provide a complete proof of this statement and weaddress the initial value problem. Following [1, 4, 12, 14], we also extend the analysis to thecontinuum of oscillators. In practice, this means dealing with a (singular) differential equationwhose variable is a discontinuous real function with possibly infinitely many jumps (and consistsof finitely many plateaus in the case of finite size populations). We also consider the initial valueproblem in the infinite-dimensional context and characterize the asymptotic dynamics, at least inthe weak coupling regime.The paper is organised as follows. After the definition of the model (section 2), we describe insection 3 the basic features of solutions, mainly the fact that repeated firings must occur in everycell. Thanks to the mean-field nature of the coupling, the dynamics commutes with permutations ofindividuals in the population. We analyze the consequences of this symmetry on collective behaviorsand grouping properties. In section 4, we consider the initial value problem and prove existence anduniqueness of global solutions for every piecewise constant initial condition, and also for every initialconfiguration when the coupling is sufficiently small. Section 5 focuses on asymptotic propertiesof the dynamics. We first characterize the simplest periodic orbits and establish a necessary andsufficient condition on the coupling parameter for their existence. The analysis of this conditionshows that an abrupt transition occurs as this parameter crosses a threshold. The transitionseparates a regime where all possible periodic trajectories exist, to a situation where only trajectorieswith sufficiently large plateaus perdure. In the same section, we prove that every finite-dimensionsolution must be asymptotically periodic, and usually approaches one of the previously mentionedperiodic trajectories. We also prove that the same results hold for every solution, provided againthat the coupling is not too large. Finally, we conclude the paper by a series of comments onrobustness with respect to changes in the model.2
Definitions
For simplicity, we assume that a single gene is involved in our process (instead of two genes in theoriginal experiment [5]) and that intercellular coupling is co-repressive and of mean-field type (asopposed to co-excitatory and spatially localized as in the original modelling [13]).Each oscillator represents a self-repressor gene [21] embedded in a host cell (same gene in everycell) labelled by a real number x ∈ (0 , t ∈ R + is characterisedby a real number u ( x, t ) ∈ [0 ,
1] that represents the so-called gene expression level (normalizedconcentration) in cell x [22].Intercellular coupling materializes via a repressor field. As in the original experiment, weconsider that gene products involved in the feedback loop are small enough so that they can diffusethrough membranes. Assuming in addition that strong stirring holds in the container, we alsoconsider that mixing of diffused gene levels occurs very rapidly in the medium ( i.e. on a muchshorter time scale than protein degradation processes). Accordingly, the repressor level M u ( x, t )in cell x , at time t , is given by the following linear combination M u ( x, t ) = (1 − ǫη ) u ( x, t ) + ǫη Z u ( y, t ) dy, where Z u ( y, t ) dy represents the mean field in the intercellular medium and where the diffusioncoefficient ǫη is composed by the product of the coupling parameter ǫ with the thresholdparameter η (the choice of ǫη instead of only ǫ turns out to be more convenient in our results).In principle, η can take any value in (0 , ǫ is restricted in a way that repressor levels neverbecome negative, i.e. we impose 0 < ǫ < /η . With these preliminary definitions provided, the evolution rule can be given. The dynamics ofgene expression levels is governed by the following singular differential equation ∂ t u ( x, t ) = − Sgn( u ( x, t )) if M u ( x, t ) > η (cid:26) u ( x, t ) = u ( x, t − u ( x, t + 0) = 1 if M u ( x, t ) η ∀ x ∈ (0 , , t ∈ (0 , + ∞ )and u ( x,
0) = u ( x ) ∀ x ∈ (0 , . (1)Throughout the paper, functions depending on a single variable depend on x (unless otherwisestated) and are viewed as one-dimensional profiles . As a consequence, no confusion results fromusing the symbol u ( x ) to denote the initial condition of the solution u ( x, t ).Equation (1) is inspired by the delay-differential equation that has been introduced in [13] inorder to reproduce the experimental oscillations. Both our model and the one in [13] obey similarprinciples: • A significant repressor level (
M u ( x, t ) > η ) in a cell prevents any production. In this case, thecorresponding gene expression level decays slowly (constant speed −
1) due to degradation. The dynamics can also be defined for ǫ = 0 and is obvious in this case. The assumption ǫ > We use the notation u ( x, t −
0) := lim s → t,s
M u ( x, t ) η ), repression can no longerprevent production and the gene is synthetized at fast rate (also on infinitesimal scales whencompared to degradation processes) until saturation is reached. This instantaneous produc-tion event is called a firing . (Lemma 3.1 below actually shows that we have M u ( x, t ) > η for all ( x, t ); hence firings occur exactly when M u ( x, t ) = η .) • Gene levels may locally reach and stay at negligible values ( u ( x, t ) = 0) while their repressorlevel remains above threshold. In this case, the next firing is delayed and cell x will evolvein sync with any other cell y whose expression level has reached zero in the mean time (aproperty called ’grouping process’, see section 3.2.4 below).The assumptions on the initial profile u are as follows(1) u is a real Borel measurable function (with values in (0 , M u ( x ) := M u ( x,
0) is well-defined for every x ∈ (0 , u is non-decreasing on (0 , hence the cell can always be labelledin a way that their expression level are monotonically ordered.(3) u be c`agl`ad , i.e. u is also left continuous on (0 , ,
1) follows from monotonicity.(4)
M u (0 + 0) > η and u (0 + 0) < u (1) = 1.The following comments are in order • the condition M u (0 + 0) > η ensures the existence and uniqueness of the solution u ( x, t )locally for t > M u ( x, t ) > η for all ( x, t )). • the condition u is c`agl`ad is consistent with the same feature of the solution profile at everytime (see expression (2) and Proposition 3.2 below) • the condition u (1) = 1 states that the cell(s) with highest expression level in the initialpopulation has (have) just fired at t = 0. Properties of the firing times in Lemma 3.2 belowimply that every solution satisfies this property at some moment in time (at infinitely manymoments indeed); hence that assumption amounts to a time translation. • the inequality u (0 + 0) < u (1) is also a matter of convenience. It ensures that the initialpopulation is not in full synchrony; otherwise the dynamics would be trivial (single oscillator)and does not require any elaborated investigation.We shall require more conditions below when we address the existence of global solutions ( i.e. functions u ( x, t ) that satisfy (1) for all ( x, t ) ∈ (0 , × R + , for given ǫ and η ). Namely, if the relation v ( x, t ) = u ( x, t ) , ∀ x = x , x and (cid:26) v ( x , t ) = u ( x , t ) v ( x , t ) = u ( x , t )holds for t = 0, and if u ( x, t ) is a solution of the equation (1), then v ( x, t ) is also a solution and the previous relationholds for all t ∈ (0 , + ∞ ). Basic dynamical features
Postponing the existence of global solutions of equation (1) to section 4 below, in this section, wedescribe basic but essential solution features. On these properties depend both the analysis andthe formulation of the results in the sections below. We begin with temporal features of individualoscillators before passing to the description of the collective properties of the population dynamics.
Here, focus is made on basic temporal features; in particular on the facts that (see Lemma 3.2) • every cell must fire repeatedly forever, and • between every two consecutive firings in x , there must be exactly one firing in every othercell y = x , unless y fires simultaneously with x .A first statement justifies the firing time definitions below. Lemma 3.1
For every solution ( x, t ) u ( x, t ) , the inequality M u ( x, t ) > η holds for all ( x, t ) ∈ (0 , × R + . Proof.
The proof is by contradiction. It relies on the fact that, for every x ∈ (0 , t u ( x, t ) is c`agl`ad. This property follows trivially from equation (1) when M u ( x, t ) η andis a consequence of the fact that t u ( x, t ) must be continuous when M u ( x, t ) > η . Togetherwith Lebesgue’s dominated convergence theorem, it implies that each function t M u ( x, t ) is alsoc`agl`ad.By contradiction, assume the existence of ( x , t ) such that M u ( x , t ) < η . We must have t > M u ( x, > η for all x ∈ (0 , t < t such that M u ( x , t ) < η and hence u ( x , t ) = 1 for all t ∈ ( t , t ]. Sinceevery expression level is smaller or equal to 1, it follows that M u ( x, t ) < η and hence u ( x, t ) = 1for all ( x, t ) ∈ (0 , × ( t , t ]. Using the definition of M , the latter yields M u ( x, t ) = 1 for all( x, t ) ∈ (0 , × ( t , t ], which contradicts the assumption M u ( x , t ) < η . ✷ Lemma 3.1 implies that firing events occur exactly when the repressor level reaches the threshold η . In particular, the first firing time in cell x , defined as T u ( x ) := sup { t > M u ( x, s ) > η, ∀ s < t } , can be characterized as follows T u ( x ) = inf { t > M u ( x, t ) = η } . The c`agl`ad property of t M u ( x, t ) and the initial assumption M u (0 + 0 , > η ensure that T u ( x ) > x ∈ (0 , T u ( x ) < + ∞ and that the repressor levelimmediately after firing, i.e. M u ( x, T u ( x ) + 0) > η , lies above η , the second firing time can bedefined similarly, namely T u ( x ) = inf { t > T u ( x ) : M u ( x, t ) = η } . Repeating the argument, one defines the successive firing times as follows T n +1 u ( x ) = inf { t > T n u ( x ) : M u ( x, t ) = η } , ∀ n ∈ N , T n u ( x ) n → + ∞ −−−−−→ + ∞ , one obtains the following explicitexpression of global solutions u ( x, t ) = (cid:26) ( u ( x ) − t ) + if 0 t T u ( x )(1 − t + T n u ( x )) + if T n u ( x ) < t T n +1 u ( x ) , n ∈ N ∀ x ∈ (0 , . (2)All required assumptions above are listed in the main statement of this section, which we nowformulate. Proposition 3.2
For every solution ( x, t ) u ( x, t ) , the following properties hold. • For every x ∈ (0 , and n ∈ N , the firing time T n u ( x ) is well-defined and is finite. • Every function x T n u ( x ) is non-decreasing and left continuous (and therefore c`agl`ad) andwe have T n u (1) T n +1 u (0 + 0) . • For every x ∈ (0 , we have T u ( x ) > M u ( x ) − η and T n +1 u ( x ) − T n u ( x ) > (1 − ǫη )(1 − η ) for all n ∈ N . Proof.
Most of the effort consists in proving the statement for n = 1 (and that T is well-defined;we already know that T is well-defined); the other cases will follow by induction, by applying theseconclusions to the solution at appropriate successive times.As we shall see, the proof for n = 1 only relies on the following assumptions on uu is non-decreasing and left continuous , u (1) = 1 , and M u ( x ) > η, ∀ x ∈ (0 , , (3)and does not need that u (0 + 0) < u (1), i.e. that the cells are initially out of sync. • Proof of monotonicity of the function T . Given two arbitrary points x < x , using thatmin { T u ( x ) , T u ( x ) } ∈ R + ∗ ∪{ + ∞} , let t > t min { T u ( x ) , T u ( x ) } . Expression(2), together with u ( x ) u ( x ), implies u ( x , t ) u ( x , t ) and then M u ( x , t ) M u ( x , t ) fromwhere the inequality T u ( x ) T u ( x ) follows. • Proof of the inequality T u (1) T u (0 + 0). The definition u ( x, T u ( x ) + 0) = 1 implies that M u ( x, T u ( x ) + 0) > M u ( y, T u ( x ) + 0) , ∀ y = x. Using also monotonicity of T , this implies that the second firing T u ( x ) in cell x cannot happenbefore cell 1 has first fired, i.e. T u (1) T u ( x ) for all x ∈ (0 , T u (1) T u (0+0)immediately follows. (Of note, if T u ( x ) = ∞ for some x , then there is nothing to prove. Moreover,the same argument, together with monotonicity of T , implies monotonicity of the function T .)For the proof of finiteness of T , we shall rely on the following statement. Lemma 3.3 (i) T u (0 + 0) < . (ii) Assume that T u ( x ) < + ∞ for some x ∈ (0 , . Then, there exists x = x and n ∈ { , } such that T u ( x ) < T n u ( x ) < T u ( x ) + 1 . Proof of the Lemma. (i) By contradiction, assume that T u (0 + 0) >
1. Then monotonicity of T implies T u ( x ) > x . Expression (2), together with u ( x )
1, yields u ( x,
1) = 0 for all x .This in turns gives M u ( x,
1) = 0, which is impossible by Lemma 3.1. u + := max { u, } . T u ( x ) < + ∞ , assume that we have T u ( x ) > T u ( x ) + 1 , ∀ x ∈ ( x ,
1] and T u ( x ) > T u ( x ) + 1 , ∀ x ∈ (0 , x ] . Then, as for statement (i), we conclude that
M u ( x, T u ( x ) + 1) = 0 for all x , which is impossible.Statement (ii) is nothing but a condensed formulation of the assumption negation. ✷ • Proof of finiteness of the function T . By contradiction, assume that T u ( x ) = + ∞ for some x ∈ (0 , T and Lemma 3.3 imply the existence of x ∈ (0 ,
1] such that (cid:26) T u ( x ) < + ∞ if 0 < x < x T u ( x ) = + ∞ if x x T u ( x − < + ∞ or T u ( x −
0) = + ∞ . In the first case, the property T u (1) T u (0 + 0) and monotonicity of T u imply that no firing can occur in a cell after time T u ( x − M u ( x, T u ( x −
0) + 1) = 0 for all x , which is impossible.Assume now that T u ( x −
0) = + ∞ . After time T u ( x − ǫ ), only those cells in the interval( x − ǫ , x ) can fire. Therefore, we have (cid:26) u ( x, T u ( x − ǫ ) + 1) = 0 if x ( x − ǫ , x ) u ( x, T u ( x − ǫ ) + 1) x ∈ ( x − ǫ , x )which implies M u ( x, T u ( x − ǫ ) + 1) ǫη ǫ < η, ∀ x ( x − ǫ , x )which is again impossible. • Proof of left continuity of the function T . Monotonicity obviously implies that the limit T u ( x − T u ( x − T u ( x ), for every x ∈ (0 , x ∈ (0 ,
1] and choose any y ∈ (0 , x ). Wehave M u ( x, T u ( y )) = M u ( y, T u ( y )) + (1 − ǫη ) (cid:0) ( u ( x ) − T u ( y )) + − ( u ( y ) − T u ( y )) + (cid:1) Using
M u ( y, T u ( y )) = η and left continuity of the initial profile u , we obtainlim y → x − M u ( x, T u ( y )) = η. Furthermore, that t M u ( x, t ) is c`agl`ad for every x ∈ (0 ,
1] (see proof of Lemma 3.1) implies thatthe limit here is equal to
M u ( x, T u ( x − M u ( x, T u ( x − η and the first firingtime definition implies T u ( x − > T u ( x ) from where the desired conclusion T u ( x −
0) = T u ( x )follows. • Proof of the inequality T u ( x ) − T u ( x ) > (1 − ǫη )(1 − η ) for all x ∈ (0 , T u ( x ) > M u ( x ) − η is a direct consequence of the initial assumption M u ( x ) > η together with M u ( x, t ) > M u ( x ) − t .) First, notice that the expression level in a cell that is about to fire,namely u ( x, T u ( x )), must minimise the expression levels at this time and must not lie above η .We conclude that M u ( x, T u ( x ) + 0) > M u ( x, T u ( x )) + (1 − ǫη )(1 − u ( x, T u ( x ))) > η + (1 − ǫη )(1 − η ) , from where the desired inequality immediately follows. As a by-product, this inequality also impliesthat M u ( x, T u ( x ) + 0) > η for all x ∈ (0 , T is indeedwell-defined, as claimed. 7 Induction step . At this stage, it remains to show that (an appropriate translation of) the pop-ulation profile immediately after cell 1 has fired satisfies the same assumptions (3) as the initialfunction u . There are two cases; either T u (1) < T u ( x ) for all x ∈ (0 ,
1] or T u (1) = T u ( x ) forsome x ∈ (0 , u ( x ) := u ( x, T u (1) + 0) for all x .The assumption T u (1) < T u ( x ), together with expression (2), implies (notice that all expressionlevels must be positive immediately after firing) u ( x ) = 1 − T u (1) + T u ( x ) , ∀ x ∈ (0 , , and so monotonicity and left continuity of T imply the same properties for u (and we obviouslyhave u (1) = 1). Moreover, the same assumption also implies M u ( x, T u (1)) > η , and a fortiori M u ( x, T u (1) + 0) > η , for all x such that T u ( x ) < T u (1). Therefore, the function u satisfiesall the assumptions of (3) and the induction can proceed.In the case where T u (1) = T u ( x ) for some x ∈ (0 , T u (1) in order to obtain a monotonic function on(0 , x max := sup { x ∈ (0 ,
1) : T u ( x ) = T u (1) } . Notice that the proof of Lemma 3.3 can be repeated for the function T to conclude that T ( x ) < + ∞ for all x . Then the proof of left continuity of T applies mutatis mutandis to prove that the function T must also be left continuous. Using also T u (1) > T u (1)+(1 − ǫη )(1 − η ), it follows that x max < T ( x max ) = T u (1), and we set u ( x ) = (cid:26) u ( x + x max , T u (1) + 0) = 1 − T u (1) + T u ( x + x max ) if 0 < x − x max u ( x + x max − , T u (1) + 0) = 1 if 1 − x max < x u is obviously non-decreasing, left continuous and we have u (1) = 1. Moreover, thedefinition of x max implies T u ( x ) > T u (1) for all x ∈ ( x max , M u ( x, T u (1) + 0) = M u ( x, T u (1)) > η, ∀ x ∈ ( x max ,
1] : T u ( x ) < T u (1) . As before, this easily implies
M u ( x ) > η for all x ∈ (0 ,
1] and all assumptions of (3) hold. Theinduction can proceed and the proof of Proposition 3.2 is complete. ✷ In this section we review grouping properties of the dynamics. These properties are consequencesof the label exchange symmetry. For simplicity, the properties are formulated in terms of T u . Asin the proof of Proposition 3.2 above, they extend to every finite time, by induction on profiles,after every full cycle of firings. If u ( x ) = u ( y ) then T u ( x ) = T u ( y ).Therefore, if at some time t , u ( x, t ) is constant on some interval, then it remains constant on thisinterval for all t > t . In other words, cells holding the same gene expression level define a group [15] and evolve in unison. (”Cluster” is another term for such groups [2].) Indeed, we have T u ( x ) > T u ( x ) for all x ∈ (0 , .2.2 Firing without grouping If T u ( x ) u ( x ), then T u ( x ) < T u ( y ) for all y such that u ( x ) < u ( y ).Therefore, if a cell x (or a group of cells including x ) fires before its level has hit 0, then it does sounaccompanied by any cell (or group) whose expression level differs from x at this instant. If ǫ
1, then for every initial profile u , we have T u ( x ) u ( x ) for all x ∈ (0 , ǫ ∈ (0 , ǫ T u ( x ) > u ( x ) for some x ∈ (0 , R u ( y, t ) dy
1, we get
M u ( x, u ( x )) = ǫη Z u ( y, u ( x )) dy η when ǫ , which, considering that M u ( x, > η , is incompatible with T u ( x ) > u ( x ). If u ( x ) < u ( y ) < T u ( x ), then T u ( y ) = T u ( x ).If the expression level in a cell/group hits 0 before firing, then the cell/group joins any othercell/group whose level also hits 0 before the same firing. When ǫ >
1, the maximal size of a forming/inflating group before a firing is 1 − ǫ .To see this, by contradiction again, assume there exists t such that u ( x, t ) = 0 for all x ∈ (0 , y ]with y > − ǫ . Then, using that u M u ( x, t ) ǫη (1 − y ) < ǫη ǫ = η, which is impossible.Group invariance and the grouping process imply that the total plateaus’ length ( i.e. theLebesgue measure of the set where u ( · , t ) is constant) is a non-decreasing function of time. Sincethis length cannot exceed 1, it must converge. In other words, the total length of intervals wheregrouping occurs upon firings must vanish as t → + ∞ .A by-product of the properties in section 3.2.3 and 3.2.5 above is that full grouping of cells ( i.e. complete synchrony) can never be achieved before any firing (and hence in finite time), unless allcells are initially in sync. Proposition 3.2 implies that, in order to prove existence and uniqueness global solutions, it sufficesto prove the same for the firing times T n u (see section 2). By induction, it suffices in turn to proveexistence and uniqueness of the first firing time function T u whose equation can be rewritten asfollows M u ( x, t ) > η, ∀ t ∈ [0 , T u ( x )) and M u ( x, T u ( x )) = η, ∀ x ∈ (0 , . (4)9here the quantity u ( x, t ) is given by (2) with n = 1, viz. u ( x, t ) = (cid:26) ( u ( x ) − t ) + if 0 t T u ( x )(1 − t + T u ( x )) + if T u ( x ) < t ∀ x ∈ (0 , , t ∈ [0 , T u (1)] . Existence and uniqueness of solutions to (4) can be granted in two distinct cases; either when theinitial profile is locally constant in a right neighborhood of every point in (0 , ǫ Proposition 4.1
Let η be arbitrary. There exists a unique global solution to equation (1) in thefollowing cases • ǫ is arbitrary and the profile u is such that there exists ∆ x > for every x ∈ (0 , , so that u is constant on ( x, x + ∆ x ] , • ǫ and u is arbitrary. More generally, the proof implies that T u can be uniquely determined under the following assump-tion: T u is either locally constant or it satisfies T u ( x ) u ( x ) , in the right neighborhood of everypoint . An example is given by periodic trajectories associated with profiles that are strictly increas-ing in a right neighborhood of every point, see comment after Proposition 5.1 and the inequality(8) below. However, there are cases that do not fit this setting and for which proving the existenceof global solutions remains open, especially if there exists x ∈ [0 ,
1) such that T u ( x + 0) = u ( x + 0) and T u ( x ) < T u ( y ) , ∀ y > x. Proof of Proposition 4.1.
We consider the two cases separately. • u is locally constant on right neighborhoods. Group invariance (property in section 3.2.2) impliesthat the solution T u of equation (4) must also be locally constant on right neighborhoods. Inparticular, if u ( x ) = u (∆ ) for all x ∈ (0 , ∆ ], then we must have T u ( x ) = T u (∆ ) on the sameinterval. We first aim to determine T u (∆ ) and more generally, to determine the firing time of thefirst firing plateau.Let the function S be defined by S ( x ) = Z x u ( y ) − u ( x ) dy = Z x u ( y ) dy − (1 − x ) u ( x ) , and consider separately the two cases ǫS (∆ ) ǫS (∆ ) > M u (∆ ) − u (∆ ) η . Hence, the firing time mustbe given by T u (∆ ) = M u (∆ ) − η u (∆ ) and T u (∆ ) < T u ( x ) , ∀ x > ∆ . For consistence with the second case, we set ∆ ′ = ∆ in this case.In the second case, let ∆ ′ := sup { x > ǫS ( x ) > } . Left continuity and boundedness of u imply that S is also left continuous. In addition, we claimthat it is non-increasing. (Indeed, using that the derivative of u exists and is finite a.e., it resultsthat the same property holds for S and its derivative is given by − (1 − x ) u ′ ( x ) ′ > ∆ . It is immediate to check that ∆ ′ corresponds to the size of the first firing plateau, viz. we have T u ( x ) = T u (∆ ′ ) < u (∆ ′ + 0) , ∀ x ∈ (0 , ∆ ′ ] and T u (∆ ) < T u ( x ) , ∀ x > ∆ ′ . T u (∆ ′ ) itself is uniquely defined by ǫ Z ′ u ( y ) − T u (∆ ′ ) dy = 1 . Now, by repeating the same arguments to the translated profile v ( x ) immediately after T u (∆ ′ )and defined by v ( x ) = (cid:26) u ( x + ∆ ′ , T u (∆ ′ ) + 0) if 0 < x − ∆ ′ − ∆ ′ < x T u is already defined on (0 , x ] (where x is arbitrary), then it can beuniquely extended to ( x, x + ∆ ′ x ] where ∆ ′ x > ∆ x . Moreover, recall that T u must be non-decreasingand left continuous. Hence, if it is defined on any set S , it can always be uniquely extended tothe semi-closed set S ℓ := T δ> S + [0 , δ ). The following technical statement then implies that thisprocess uniquely defines T u on (0 , Lemma 4.2
Let S ⊂ (0 , + ∞ ) be a semi-closed set such that • there exists x ∈ (0 , such that (0 , x ] ⊂ S , • for every x ∈ S , there exists ∆ x > such that ( x, x + ∆ x ] ⊂ S .Then (0 , ⊂ S . Proof of the Lemma.
The proof proceeds by transfinite induction. Starting with (0 , x ], theinduction property implies (0 , x ] ⊂ S where x > x + ∆ x , and then(0 , x n ] ⊂ S, ∀ n ∈ N . Let x ω be the limit of the increasing sequence { x n } n ∈ N . By induction, we also have (0 , x ω ) ⊂ S and then (0 , x ω ] ⊂ S ℓ = S .If x ω >
1, we are done. Otherwise, we continue the induction to successive ordinals, until weeventually reach an uncountable ordinal ω . Then we must have x ω >
1, otherwise we would havea uncountable collection of contiguous intervals whose union covers (0 , x ω ] but not (0 , , ⊂ (0 , x ω ] and the proof is complete. ✷ • ǫ
1. Firing without grouping (property in section 3.2.3) implies T u u in this case; hencesolution components u ( x, t ) remain positive for t T u (1). In particular, using the definition ofthe lower trace T u of the firing profile T u (see Appendix B), we get the following expression u ( y, T u ( x )) = − T u ( x ) + T u ( y ) if 0 < y < T u ( x ) u ( x ) − T u ( x ) if T u ( x ) = x and y = xu ( y ) − T u ( x ) if T u ( x ) < y ∀ x ∈ (0 ,
1] (5)(There is no need to define u ( T u ( x ) , T u ( x )) when T u ( x ) < x .) To proceed, we shall need thatthe profile and firing time traces must coincide in absence of grouping, as now stated. Claim 4.3 If T u ( x ) < T u ( y ) for every pair ( x, y ) ∈ (0 , such that u ( x ) < u ( y ) , then we have T u ( x ) = u ( x ) for all x ∈ (0 , . roof of the Claim. We consider the cases u ( x ) < x and u ( x ) = x separately.In the first case, the equality u ( y ) = u ( x ) for all y ∈ ( u ( x ) , x ] implies T u ( y ) = T u ( x ) for all y ∈ ( u ( x ) , x ] and thus T u ( x ) u ( x ). By contradiction, if we had T u ( x ) < u ( x ), then for any y ∈ ( T u ( x ) , u ( x )) we would have u ( y ) < u ( x ) and T u ( y ) > T u ( x ) , which contradicts the assumption of the Claim.For the second case, notice that we have u ( y ) < u ( x ) for all y < u ( x ), and then T u ( y ) < T u ( x ) forall y < u ( x ), from the Claim assumption. This implies u ( x ) T u ( x ) and the conclusion followsfrom the facts that u ( x ) = x and T u ( x ) x . ✷ The inequality T u ( x ) u ( x ) together with ’firing without grouping’ ensures that the assump-tion in Claim 4.3 holds. Using also expression (5) to manipulate equation (4), we obtain thefollowing affine functional equation(Id − L u ) T u ( x ) = (1 − ǫη ) u ( x ) − η + ǫη u ( x ) + Z u ( x ) u ( y ) dy ! , ∀ x ∈ (0 , , (6)where the linear operator L u is defined by L u v ( x ) = ǫη Z u ( x )0 v ( y ) dy, ∀ x ∈ (0 , , for every bounded Borel measurable function v defined on (0 , k · k ∞ , the assumption ǫ < /η implies k L u k ∞ <
1. Hence, the operatorId − L u is invertible with bounded inverse. Accordingly, there exists a unique bounded solution x T u ( x ) to equation (6). ✷ This section investigates the asymptotic behavior of global solutions as t → + ∞ . To that goal, weequip the set of bounded Borel measurable functions defined on (0 , L -norm k · k . Anticipating the results below on asymptotic periodicity, we present here preliminary propertiesof periodic trajectories. By a periodic trajectory , we mean a solution { u ( x, t ) } such that thereexists τ ∈ R + so that u ( x, t + τ + 0) = u ( x, t ) , ∀ x ∈ (0 , , t ∈ R + . Of note, the period τ here can only be one of the numbers { T n u (1) } n ∈ N , because the profile u ( · , t )cannot be both non decreasing and satisfy u (1 , t + 0) = 1 at other times. Here, we shall focus on T u (1)-periodic trajectories ( i.e. τ = T u (1)) because these solutions turn out to play a special rolein the asymptotic dynamics.As our next statement indicates, periodic trajectories are uniquely determined by the traces u and u of their initial profile. Recall from Appendix B that every lower trace function is en-tirely determined by a countable collection of pairwise disjoint semi-open intervals in (0 ,
1] and theknowledge of a lower trace completely determines the upper trace.12 roposition 5.1 (i)
Let η, ǫ be any parameters and let u tr be any lower trace function. Thereexists at most one non-decreasing profile u such that u = u tr and such that the trajectory issuedfrom u is periodic with period T u (1) . (ii) This periodic trajectory exists iff ǫ . R u tr ( x ) dx (cid:18) − inf x ∈ (0 ,u tr (1)) u tr ( x ) − u tr ( x )1 − u tr ( u tr ( x )+0)+ u tr ( x ) (cid:19) , (7) where the symbol . means < if the infimum is a minimum, and it means if this bound is notattained. The proof is given below. Notice that the condition (7) does not depend on η and the infimum is aminimum for every finite step trace function. In more general cases, this infimum may be attainedor not. Moreover, the condition simply reduces to ǫ . R u tr ( x ) dx for every trace u tr for which forevery δ >
0, there exists x δ ∈ (0 ,
1] such that u tr ( x δ ) − u tr ( x δ ) δ. In addition, using u tr ( x ) x and properties of the traces (Appendix B), it is easy to conclude thatthe denominator in (7) is certainly not larger than , and this bound is attained for the increasingtrace u tr ( x ) = x for all x . Moreover, the quantity in the right hand side of (7) continuously dependson u tr ( L -topology). In particular, if all existing plateaus are sufficiently small (depending on ǫ ),then the associated periodic orbit does not persist when ǫ > Corollary 5.2 (i)
There exists a T u (1) -periodic trajectory with u ( · ,
0) = u tr for every lower tracefunction u tr , iff ǫ < . (ii) For every ǫ > , there exists ℓ ǫ > such that, for every lower trace function u tr so that k u tr − u tr k < ℓ ǫ , the associated T u (1) -periodic trajectory does not exist. Proof of Proposition 5.1.
We study the equation u ( · , T u (1) + 0) = u and u = u tr , for an arbitrary trace u tr . Of note, as a profile immediately after firing, the function u must satisfy u >
0. Equation (2), together with this assumption, implies that this equation rewrites as1 − T u (1) + T u = u, which shows that the difference T u − u = T u (1) − T u by u + T u (1) −
1, and u by u tr , in equation (6) for thefirst firing time (where now ∆ = 1), one obtains after some simple algebra the following equationfor u and T u (1): ǫηu + η = 1 − T u (1) + ǫηT u (1) u tr + ǫη Z u ( y ) dy. Integrating over (0 ,
1] yields the following expression T u (1) = 1 − η − ǫη R u tr ( x ) dx u (1) = 1 in the equation for the function u above and evaluated at x = 1 in orderto eliminate the integral term, one gets that the solution writes u = 1 − T u (1)( u tr (1) − u tr ) . Clearly, this solution is unique and is left continuous. Moreover, it is non-decreasing and we have u T u (1) > i.e. iff ǫη < R u tr ( x ) dx , using the constraint η <
1. Finally, we need to make sure that T u − u is negative, viz. T u (1) < u > ǫ < R u tr ( x ) dx . (8)Altogether, we conclude that a unique initial profile with trace u tr generates a T u (1)-periodictrajectory with T u < u iff the inequality (8) holds.The analysis is similar in the second case. Here, we use the equality T u = u tr and equation(13) in Appendix B to obtain T u = u tr . Accordingly, the solution expression at time T u ( x ) isgiven by u ( y, T u ( x )) = − T u ( x ) + T u ( y ) if 0 < y u tr ( x )0 if u tr ( x ) < y u tr ( x ) u ( y ) − T u ( x ) if u tr ( x ) < y u tr ( x ) = 1.) Inserting this expression into theequation M u ( · , T u ( · )) = η , one gets after some simple algebra (using also u ( y ) = u ( x ) for all y ∈ ( u ( x ) , u ( x )]) (1 − u tr )(1 − T u (1)) + u tr − u + Z u ( y ) dy = 1 ǫ . Integrating over (0 ,
1) and using equation (16) in Appendix B yields in this case T u (1) = 2 R u tr ( x ) dx − ǫ R u tr ( x ) dx . As in the first case, by subtracting from the equation its expression for x = 1, we get that thesolution writes u = (1 − u tr )(1 − T u (1)) + 1 − u tr (1) + u tr . The existence of this solution requires T u (1) > i.e. T u − u is a non-negative function); thisinequality is equivalent to ǫ > R u tr ( x ) dx which is complementary to the existence condition (8) in the first case. We also need to impose u >
0. By monotonicity, this condition is equivalent to T u (1) − . − u tr (1)1 − u tr (0 + 0) . If u, v are two real functions, then u v (resp. u < v ) means u ( x ) v ( x ) (resp. u ( x ) < v ( x )) for all x ∈ (0 , T u = u tr , viz. we must have T u ( x ) < u ( y ) , ∀ y > u tr ( x ) , x ∈ (0 , u tr (1)) . Using T u ( x ) = u ( x ) + T u (1) −
1, the expression of u and relation (15) in Appendix B, thiscondition is equivalent to T u (1) − . inf x ∈ (0 ,u tr (1)) u tr ( x ) − u tr ( x )1 − u tr ( u tr ( x ) + 0) + u tr ( x ) . We have inf x u tr ( x ) − u tr ( x ) − u tr (1) and inf x u tr ( u tr ( x ) + 0) − u tr ( x ) u tr (0 + 0); hence the lastcondition is the strongest one. Simple algebra finally yields the inequality (7) and this completesthe proof. ✷ Proposition 4.1 implies that a global solution of equation (1) exists for every initial finite stepprofile (because every finite step function is obviously locally constant in the right neighborhoodof every point). Moreover, the grouping properties in section 3.2 imply that the number of profilesteps either remains constant or decreases before each firing; hence this number must eventuallybecome constant. Therefore, when dealing with asymptotic properties, we may assume withoutloss of generality that this number remains constant, viz. that the lower trace function is periodicafter each full cycle of firing.
Theorem 5.3
Let η, ǫ be any parameters and let u be any initial finite step function for which thenumber of clusters remains constant in the corresponding trajectory. Then, we have lim t → + ∞ k u ( · , t + T u per (1) + 0) − u ( · , t ) k = 0 , where u per is the periodic trajectory profile such that u per = u . If, in addition, the condition (7) holds with u tr = u , then we also have lim n → + ∞ k u ( · , T n u (1) + 0) − u per k = 0 . It may happen that the first limit in this statement holds although the condition (7) fails, moreprecisely when ǫ = 1 R u tr ( x ) dx (cid:18) − min x ∈ (0 ,u tr (1)) u tr ( x ) − u tr ( x )1 − u tr ( u tr ( x )+0)+ u tr ( x ) (cid:19) . In this case, the periodic cycle { u per ( · , t ) } is an example of ghost orbit mentioned in Appendix A.Besides, the proof of the Theorem below implies that when ǫ > R u tr ( x ) dx (cid:18) − min x ∈ (0 ,u tr (1)) u tr ( x ) − u tr ( x )1 − u tr ( u tr ( x )+0)+ u tr ( x ) (cid:19) , any solution issued from a finite step profile u such that u = u tr , must experience the grouping ofa least two plateaus in finite time. In this way, we have obtained a complete picture of all possibleasymptotic finite-dimensional population distributions, depending on the coupling intensity ǫ .15 roof of Theorem 5.3. The dynamics of finite step functions is entirely determined by the number N of steps, by the step lengths { ℓ n } Nn =1 and by the step expression levels { u n } Nn =1 . (Here, steplabelling follows from cell labelling on (0 , viz. n = 1 means the first group x ∈ (0 , x ] for some x > n = 2 means x ∈ ( x , x ] with x > x and so on.)Assuming that no grouping occurs in time implies that not only N but also the length distri-bution { ℓ n } Nn =1 remains constant. Only the expression levels { u n } Nn =1 may depend on time. In thissetting, we aim to prove that the time evolution of every expression level vector is asymptoticallyperiodic in R N .To that goal, it suffices to consider the discrete time dynamics that brings the expression levelsafter a firing to those after the next firing (Poincar´e return map). From thereon in this proof, by time , we mean the integer t ∈ N that labels the vector { u tn } Nn =1 after the t th firing.It turns out more convenient to combine the discrete time dynamics with the cyclic permutationof indices, so that any vector with non-decreasing components and u N = 1 is mapped after everyiteration onto a vector carrying the same properties. Ignoring systematically the last component u N = 1, this amounts to considering iterations of the ( N − { u tn } N − n =1
7→ { u t +1 n +1 } N − n =1 .Beside the original parameters η and ǫ , this firing map F ℓ is also parametrized by the steplength distribution ℓ := { ℓ n } Nn =1 (at time t ). Its action on vectors u := { u n } N − n =1 ∈ U N − , where U N − = n u = { u n } N − n =1 : 0 < u < u < · · · < u N − < o , is explicitly given by ( F ℓ u ) n = u n +1 − T ℓ u, ∀ n = 1 , · · · , N − , where T ℓ u is the time of the first firing in the trajectory starting from the finite step profile associatedwith u (at time t ). Iterations of the firing map have to incorporate permutations of the step lengthdistribution, i.e. we need to consider the composed map F Nℓ := F R N − ℓ ◦ F R N − ℓ ◦ · · · ◦ F Rℓ ◦ F ℓ , (9)where ( Rℓ ) n = ℓ n +1 mod N for all n ∈ { , · · · , N } . By identifying each { u n } N − n =1 ∈ U N − withits corresponding N -step profile, the map F Nℓ is clearly equivalent to the return map u ( · , u ( · , T u (1) + 0). Theorem 5.3 is therefore an immediate consequence of the following statement.Let k · k be any norm in R N − . Proposition 5.4
For every step length distribution ℓ = { ℓ n } Nn =1 , there exist C > and ρ ∈ [0 , such that for every pair u, v ∈ U N − for which ( F Nℓ ) k u, ( F Nℓ ) k v ∈ U N − for all k ∈ N , we have k ( F Nℓ ) k u − ( F Nℓ ) k v k Cρ k k u − v k , ∀ k ∈ N . (10) If the length distribution period N per , defined by N per = min { k : ℓ n + k mod N = ℓ n , ∀ n ∈ { , · · · , N }} , happens to be smaller than N , it actually suffices to consider the composed map F R N per − ℓ ◦ · · · ◦ F Rℓ ◦ F ℓ , becauseappropriate permutations of the profiles associated to iterates of this map are non-increasing finite step functionswith identical step length distribution as the initial profile. The same consideration suggests to consider, given anystep length distribution, the permutation that minimises the period N per . roof of the Proposition. Independently of considerations on parameters ( i.e. whether or not ǫ a priori two cases for the expression of T ℓ u , depending on whether u fires before reaching0 or after. Simple calculations yield the following expression: T ℓ u = (1 − ǫη ) u + ǫη N P n =1 ℓ n u n − η if u > T ℓ u − ℓ (cid:18) N P n =2 ℓ n u n − ǫ (cid:19) if u T ℓ u Accordingly, for every ℓ , the map F ℓ is a continuous (piecewise) affine map with at most two pieces,say F + ℓ and F − ℓ . Let DF + ℓ and DF − ℓ be the linear maps associated with each piece. We claim that,for every vector pair u, v ∈ U N − , there exists α ∈ [0 ,
1] such that F ℓ u − F ℓ v = (1 − α ) DF + ℓ ( u − v ) + αDF − ℓ ( u − v ) . (11)Equation (11) is obvious (and holds with α ∈ { , } ) if there exists s ∈ { + , −} such that F ℓ w = F sℓ w for w = u, v . Otherwise, by continuity of the firing time function w T ℓ w , there exists β ∈ (0 , F ℓ ( βu + (1 − β ) v ) = F + ℓ ( βu + (1 − β ) v ) = F − ℓ ( βu + (1 − β ) v ) . If F ℓ u = F + ℓ u and F ℓ v = F − ℓ v , one gets F ℓ u − F ℓ v = F + ℓ u − F + ℓ ( βu + (1 − β ) v ) + F − ℓ ( βu + (1 − β ) v ) − F − ℓ v = (1 − β ) DF + ℓ ( u − v ) + βDF − ℓ ( u − v ) viz. equation (11) holds with α = β . Otherwise, we must have F ℓ u = F − ℓ u and F ℓ v = F + ℓ v , and asimilar calculation shows that equation (11) holds with α = 1 − β .The assumption ( F Nℓ ) k u, ( F Nℓ ) k v ∈ U N − for all k ∈ N implies that we must have F R j ℓ ◦ · · · ◦ F ℓ ◦ ( F Nℓ ) k u and F R j ℓ ◦ · · · ◦ F ℓ ◦ ( F Nℓ ) k v ∈ U N − , ∀ j ∈ { , · · · , N − } , k ∈ N . By induction, the composed map F Nℓ defined by (9) is also a continuous (piecewise) affine mapwith a priori N − pieces; each piece writes F s N − R N − ℓ ◦ F s N − R N − ℓ ◦ · · · ◦ F s Rℓ ◦ F s ℓ for some { s n } N − n =0 ∈ {− , + } N . Moreover, equation (11) implies that, for every vector pair u, v ,there exists α { s n } N − n =0 ∈ [0 ,
1] for every piece, such that P { s n } N − n =0 ∈{− , + } N α { s n } N − n =0 = 1 and F Nℓ u − F Nℓ v = X { s n } N − n =0 ∈{− , + } N α { s n } N − n =0 DF s N − R N − ℓ · · · DF s ℓ ( u − v ) . (12)Consider now the joint spectral radius (see e.g. [18]) of the finite collection A = { DF sR j ℓ } j ∈{ , ··· ,N − } ,s ∈{− , + } defined byJoinSpecRad( A ) := lim sup k → + ∞ max ((cid:13)(cid:13)(cid:13)(cid:13)(cid:13) k Y i =1 A i (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) : A i ∈ A , ∀ i ∈ { , · · · , k } )! k .
17y applying relation (12) repeatedly to the iterates ( F Nℓ ) k u and ( F Nℓ ) k v , a straightforward argu-ment shows that the inequality (10) holds under the condition JoinSpecRad( A ) <
1. Furthermore,since A is a finite set, we have [3]JoinSpecRad( A ) = lim sup k → + ∞ max ( SpecRad k Y i =1 A i ! : A i ∈ A , ∀ i ∈ { , · · · , k } )! k , where SpecRad is the standard spectral radius of a matrix. Therefore, in order to prove theProposition, it suffices to show that the right hand side here is smaller than 1.To proceed, we observe that, by applying the change of variable { u n } N − n =1
7→ { v n } N − n =1 , where v n = (cid:26) u n − u n +1 if n ∈ { , · · · , N − } u N − if n = N − , the matrix DF sℓ transforms into A a ( ℓ,s ) , where the vectors a ( ℓ, s ) = { a n ( ℓ, s ) } N − n =1 for s = − , + arerespectively defined by a n ( ℓ, +) = 1 − ǫη N X m = n +1 ℓ m and a n ( ℓ, − ) = 11 − ℓ n +1 X m =2 ℓ m , ∀ n = 1 , · · · N − , and where, given an arbitrary vector a = { a k } Kk =1 ( K ∈ N ), the matrix A a is the following K × K companion matrix A a = · · · · · · · · · − a − a · · · · · · − a K . Using that the coefficients a n ( ℓ, +) and a n ( ℓ, − ) are positive and decreasing, a classical result innumerical analysis [9] implies that SpecRad( A a ( ℓ,s ) ) < s ∈ {− , + } , for every length distribution ℓ = { ℓ n } N − n =1 . However, the products k Q i =1 A a ( ℓ i ,s i ) of such matrices usually do not commute (and theassociated semi-group is usually infinite) and do not appear to have any special form that wouldallow one to immediately conclude about their spectral radius. Instead, we shall need the followingstatement that takes advantage of the matrix characteristics. Lemma 5.5 [10]
Let { a ( j ) } Jj =1 be a finite collection of K -dimensional vectors a ( j ) = { a ( j ) k } Kk =1 whose components satisfy the following constraint min k ∈{ , ··· ,K } a ( j ) k > , ∀ j ∈ { , · · · , J } . Then, for any collection { j i } ki =1 with j i ∈ { , · · · , J } for all i , we have SpecRad k Y i =1 A a ( ji ) ! max j ∈{ , ··· ,J } max k ∈{ , ··· ,K } a ( j ) k a ( j ) k +1 ! k where we have set a ( j ) K +1 := 1 for all j . { a ( R j ℓ, s ) } j ∈{ , ··· ,N − } ,s ∈{− , + } and letting a N ( R j ℓ, s ) = 1,we immediately concludeJoinSpecRad( A ) = max j ∈{ , ··· ,N − } ,s ∈{− , + } ,n ∈{ , ··· ,N − } a n ( R j ℓ, s ) a n +1 ( R j ℓ, s ) < , as required. ✷ In this section, we go back to arbitrary non-decreasing left continuous initial profiles. Let µ c ∈ (0 . , .
47) be the positive solution of the equation e µ + µ − µ = 2 . Proposition 5.6
Let η and ǫ be such that ǫη < µ c . There exists ρ ∈ [0 , so that for every pairof initial profiles u, v satisfying u = v and T u u and T u (1) < T u ( resp. T v v and T v (1) < T v ) , we have k u ( · , T u (1) + 0) − v ( · , T v (1) + 0) k ρ k u − v k . By choosing v = u ( · , T u (1) + 0), the previous condition obviously holds for every trajectory forwhich no grouping ever happens. As the next statement says, this is the case of every trajectory inthe weak coupling regime and the conclusion holds provided that η < µ c (recall that this thresholdis assumed to be small from the experiment under modeling). Corollary 5.7
Let η < µ c and ǫ < . For every initial profile, the subsequent trajectory is asymp-totically periodic. More generally, Proposition 5.6 implies asymptotic stability of periodic trajectories (associatedwith infinite dimensional profiles). Indeed, the periodic trajectory associated with the infinitelymany step profile with trace u tr is known to exist iff ǫ Z u tr . Moreover, the analysis in the proof of Proposition 5.1 shows this condition implies the one inProposition 5.6. Accordingly for ǫη < µ c , this trajectory attracts all solutions in its neighborhoodwhose profiles after every cycle of firing have trace u tr . Proof of the Proposition.
The assumption T u (1) < T u implies that the solution u ( · , T u (1) + 0)after a full cycle of firing is given by u ( x, T u (1) + 0) = 1 − T u (1) + T u ( x ) , ∀ x ∈ (0 , . Similarly, we have v ( x, T v (1) + 0) = 1 − T v (1) + T v ( x ) for all x . The assumption T u u impliesthat the first firing time T u is defined by expression (6) in Section 4. A similar expression holdsfor T v . In addition, the assumption u = v implies L u = L v .19y inverting the operator Id − L u by means of a Neumann series, and letting µ = ǫη , we obtainafter some simple algebra u ( x, T u (1) + 0) =(1 − µ ) u ( x )+ µ + ∞ X k =1 L ku (cid:18)Z u ( y ) dy − u (cid:19) ( x ) − L ku (cid:18)Z u ( y ) dy − u (cid:19) (1)+ C ( x ) , where C ( x ) does not depend on u . A similar expression holds for v ( x, T v (1) + 0).To obtain an upper bound on k u ( · , T u (1) + 0) − v ( · , T v (1) + 0) k , we need to estimate thenorm k L ku w − L ku w (1) k for every measurable function w defined on (0 ,
1] and every k ∈ N . First,notice that for such functions, by reversing the order of integration, we obtain k L u w k µ Z Z u ( x )0 | w ( y ) | dydx µ Z u (1)0 Z u ( y ) dx ! | w ( y ) | dy µ k w k . Moreover, a similar argument implies the following inequality | L k +1 u w ( x ) − L k +1 u w (1) | µ k L ku w k , ∀ x ∈ (0 , , k ∈ N , and then the estimate k L ku w − L ku w (1) k µ k +1 k w k for all k ∈ N , by using induction. For constantfunctions w ( x ) = w , this estimate can be improved by using a better estimate on k L u w k . Indeed,using u ( x ) x in a induction implies k L ku w k µ k ( k + 1)! k w k . Altogether, we obtain the inequality k u ( · , T u (1) + 0) − v ( · , T v (1) + 0) k (cid:18) − µ + µ ( e µ − µ −
1) + µ − µ (cid:19) k u − v k , and the Proposition follows from the fact that 1 − µ + µ ( e µ − µ −
1) + µ − µ < µ < µ c . ✷ The model equation (1) has been inspired by the delayed-differential equation in [13] that mimicsexperimental synchronized oscillations reported in [5]. Its features have been selected in orderto combine biological relevance, mathematical rigour and non-trivial phenomenology dependingon parameters. While some of these features may appear specific, e.g. permutation symmetry ofthe interaction term in the repressor definition, possible temporary vanishing of expression levels,instantaneous firings, not to mention the absence of heterogeneities or noise fluctuations, we wouldlike to stress that the phenomenology here does not depend on these features and can be extendedby continuity.Indeed, the contraction property of (sufficiently high iterates of) the return maps after full firingcycles is robust to smooth perturbations of the dynamics. It means for instance that any sufficientlysmall C perturbation of the map F Nℓ , where ℓ is an arbitrary finite dimensional step lengthdistribution, satisfies the contraction property (10). Hence, if the associated periodic trajectory20xists ( i.e. condition (7) holds for the increasing profile associated with ℓ ), then this trajectorycan be continued as a stable periodic orbits of the C -perturbed map. As a consequence, out ofbifurcation points, all attractive periodic trajectories of finite populations persist for sufficientlysmall perturbation of the original dynamics, independently of the mentioned specific properties ofequation (1).Naturally, this does not precludes a separate investigation into more general degrade-and-firedynamics. This will be the subject of future studies. Acknowledgments
We are grateful to E. Key and H. Volkmer for promptly providing us a proof of Lemma 5.5.
A Discontinuous dependence on initial conditions
Instantaneous resetting simplifies the analysis of equation (1). However, it makes the proof ofglobal existence of solutions rather delicate (and apparently unaccessible by standard approachessuch as the Picard operator). It also implies that the solution dependence on initial profiles hasdiscontinuities. In particular, this is the case for the first firing time function T u .Indeed, there are examples of sequences { u n } of profiles that uniformly converge to a limit profile u ∞ and for which we have lim n → + ∞ T u n ( x ) = T u ∞ ( x ) for some x ∈ (0 , ǫ
1, let u ∞ be a profile with a left plateau, i.e. u ∞ ( x ) = u ∞ ( x ) for all x ∈ (0 , x ]( x > u n ( x ) = (cid:26) u ∞ ( x ) − n if 0 < x x u ∞ ( x ) if x < x T u ∞ ( x ) = T u ∞ ( x ) for all x ∈ (0 , x ] and direct calculations yield the followingresult lim n → + ∞ T u n ( x ) = (cid:26) T u ∞ ( x ) , ∀ x ∈ (0 , x ] T u ∞ ( x ) + ǫη x ( T u ∞ ( x ) − u ∞ ( x ) + 1) , ∀ x ∈ ( x , x ]and the inequalities 0 u ∞ ( x ) − T u ∞ ( x ) < n → + ∞ T u n ( x ) > T u ∞ ( x ) for all x ∈ ( x , x ].In addition, discontinuities may also result in the existence of attracting ghost orbits , de-pending on parameters. Ghost orbits are periodic cycles of profiles, viz. { u ( x, t ) } ( x,t ) ∈ (0 , × R + with u ( · , t + τ + 0) = u ( · , t ) for some τ >
0, which, while they do not satisfy equation (1), attract alltrajectories in their neighborhood (uniform topology). As shown after Theorem 5.3, ghost orbitsexist at bifurcation points in the parameter space, when a periodic orbit collapses.
B The lower and upper traces of a non-decreasing function
Let u : (0 , → (0 ,
1] be a left continuous and non-decreasing function. Its lower trace u andrespectively upper trace u are defined as follows u ( x ) = inf { y ∈ (0 ,
1] : u ( y ) > u ( x ) } , ∀ x ∈ (0 , , and u ( x ) = sup { y ∈ (0 ,
1] : u ( y ) u ( x ) } , ∀ x ∈ (0 , . These functions satisfy the following basic properties.21 u ( x ) x u ( x ) x ∈ (0 , • either u ( x ) < x or x < u ( x ) implies u ( y ) = u ( x ) for all y ∈ ( u ( x ) , u ( x )]. • If u is strictly increasing, then u ( x ) = u ( x ) = x for all x ∈ (0 , • u ◦ u ( x ) = u ( x ) iff u is continuous at x . • u ◦ u = u . • Both functions u and u are left continuous and non-decreasing. (We prove the property for u here; the proof for u is similar and is left to the reader. Monotonicity is obvious and implies u ( x − u ( x ). Left continuity is also evident in the case u ( x ) < x . If, otherwise u ( x ) = x ,there must be a sequence { x n } n ∈ N such that u ( x n ) < u ( x n +1 ) and lim n → + ∞ x n = x . The formercondition implies u ( x n ) > x n − . Together with the latter, we obtain u ( x − > x = u ( x ) asdesired.)In our context, the traces provide information about the group structure of a population at time t : u ( x, t ) = u ( x, t ) = x means that cell x is isolated, while u ( x, t ) < u ( x, t ) means that all cells y ∈ ( u ( x, t ) , u ( x, t )] belong to the same group.The properties of the lower trace above imply that this function can be entirely determined byits plateaus; namely by considering the following decomposition(0 ,
1] = C < ∪ C = , where C < = { x ∈ (0 ,
1] : u ( x ) < x } and C = = { x ∈ (0 ,
1] : u ( x ) = x } , the second item above imposes the existence of a countable (possibly empty) set D such that C < = S i ∈D ( x − i , x + i ] where x − i < x + i x − i +1 for all i . (Notice that C = is empty when u (or u ) isa step function.) In other words, every countable (possibly empty) collection of pairwise disjointsemi-open intervals in (0 ,
1] uniquely defines a lower trace function.The upper trace function depends only on the lower trace, i.e. u = u ◦ u (and vice-versa , we have u = u ◦ u ). One can prove this fact using the sets C < and C = and the analogous decomposition for theupper trace. However, for our purpose, it is more convenient to use the following characterization u ( x ) = inf { u ( y ) : x < u ( y ) } , ∀ x ∈ (0 , , (13)with the convention that inf ∅ = 1 in this expression. To prove this relation, notice first that wemust have u ( x ) inf { u ( y ) : x < u ( y ) } . Indeed, otherwise there existed y such that x < u ( y ) and u ( y ) < u ( x ). Using that the former inequality is equivalent to u ( x ) < y , it results that we musthave u ( y ) < u ( x ) < y, (14)which is clearly incompatible with the definition of the traces. Secondly, still by using contradiction,assume that u ( x ) < inf { u ( y ) : x < u ( y ) } . This implies the existence of z such that u ( x ) < z < Notice that the lower trace can be alternatively defined as u − ◦ u where the generalized inverse (also called thequantile function in Probability Theory) u − can be defined as u − = inf { y ∈ (0 ,
1] : u ( y ) > x } , ∀ x ∈ (0 , . In this viewpoint, the property in this item reads u − ◦ u = Id for every strictly increasing function u . A similarcomment applies to the upper trace. { u ( y ) : x < u ( y ) } . However, the first inequality here implies u ( x ) < u ( z ) and then x < u ( z )which contradicts the second inequality.Similar arguments prove the following relation u ( x ) = u ( u ( x ) + 0) = inf { u ( y ) : u ( x ) < y } , ∀ x ∈ (0 , u (1)) . (15)Indeed, as before, we must have u ( x ) inf { u ( y ) : u ( x ) < y } because the converse would yield tothe double inequality (14) otherwise. Now if there existed z such that u ( x ) < z < inf { u ( y ) : u ( x ) < y } , the right inequality would imply z < u ( y ) y for all u ( x ) < y , hence z u ( x ) holds, whichcontradicts the left inequality.In the main text, we also refer to the following relation Z u ( x ) dx + Z u ( x ) dx = 1 . (16)In order to prove this relation, consider again the decomposition C < ∪ C = with C < = S i ∈D ( x − i , x + i ].A moment’s reflexion yields Z u ( x ) dx = X i ∈D x − i ( x + i − x − i ) + Z C = xdx and Z u ( x ) dx = X i ∈D x + i ( x + i − x − i ) + Z C = xdx, and then Z u ( x ) dx + Z u ( x ) dx = X i ∈D ( x + i ) − ( x − i ) + 2 Z C = xdx. Equation (16) then directly follows from the fact that( x + i ) − ( x − i ) = 2 Z x + i x − i xdx. References [1]
L.F. Abbott and C. van Vreeswijk , Asynchronous states in networks of pulse-coupledoscillators , Phys. Rev. E (1993), 1483–1490.[2] J.A. Acebron, L.L. Bonilla, C.J. Perez-Vicente, F. Ritort, and R. Spigler , TheKuramoto Model: A simple paradigm for synchronization phenomena , Rev. Mod. Phys. (2005), 137–185.[3] M.A. Berger and Y. Wang , Bounded semigroups of matrices , Linear Algebra Appl. (1992), 21–27.[4]
P.C. Bressloff , Mean-field theory of globally coupled integrate-and-fire neural oscillatorswith dynamic synapses , Phys. Rev. E (1999), 2180–2190.[5] T. Danino, O. Mondragon-Palomino, L.S. Tsimring, and J. Hasty , A synchronizedquorum of genetic clocks , Nature (2010), 326–330.236]
U. Ernst, K. Pawelzik, and T. Geisel , Synchronization induced by temporal delays inpulse-coupled oscillators , Phys. Rev. Lett. (1995), 1570.[7] B. Fernandez and L.S. Tsimring , Corepressive interaction and clustering of degrade-and-fire oscillators , Phys. Rev. E (2011), 051916.[8] B. Fernandez and L.S. Tsimring , Typical trajectories of coupled degrade-and-fire oscilla-tors: from dispersed populations to massive clustering , J. Math. Bio. (2014) 1627–1652.[9] E.I. Jury , Theory and applications of the z-transform method , Wiley, New York, 1964.[10]
E.S. Key and H. Volkmer , Private Communication.[11]
Y. Kuramoto , Chemical oscillations, waves and turbulence , Springer, 1984.[12]
Y. Kuramoto , Collective synchronization of pulse-coupled oscillators and excitable units ,Physica D (1991), 15–30.[13] W. Mather, M.R. Bennet, J. Hasty, and L.S. Tsimring , Delay-induced degrade-and-fireoscillations in small genetic circuits , Phys. Rev. Lett. (2009), 068105.[14]
A. Mauroy and R. Sepulchre , Global analysis of a continuum model for monotone pulse-coupled oscillators , IEEE Trans. Autom. Control (2013), 1154–1166.[15] R.E. Mirollo and S.H. Strogatz , Synchronization of pulse-coupled biological oscillators ,SIAM J. Appl. Math. (1990), 1645–1662.[16] C.S. Peskin , Mathematical aspects of heart physiology , Courant Institute of MathematicalScience Publication, New York, 1975.[17]
A. Pikovsky, M. Rosenblum, and J. Kurths , Synchronization: A universal concept innonlinear sciences , Cambridge University Press, 2001.[18]
G.-C. Rota and G. Strang , A note on the joint spectral radius , Proc. Netherlands Academy (1960), 379–381.[19] W. Seen and R. Urbanczik , Similar nonleaky integrate-and-fire neurons with instantaneouscouplings always synchronize , SIAM J. Appl. Math. (2000), 1143–1155.[20] S.H. Strogatz , From Kuramoto to Crawford: exploring the onset of synchronization in pop-ulations of coupled oscillators , Physica D (2000), 1–20.[21]
R. Thomas and R. D’Ari , Biological feedbacks , CRC Press, 1990.[22]
J.J. Tyson, K. Chen, and B. Novak , Network dynamics and cell physiology , Nat. Rev.Mol. Cell Bio.2