Wassim Abou-Jaoudé
email: wassim.abou-jaoude@polytechnique.org

Madalena Chaves
email: madalena.chaves@inria.fr

Jean-Luc Gouzé
email: jean-luc.gouze@inria.fr

J-L G Biocore

Links between topology of the transition graph and limit cycles in a two-dimensional piecewise affine biological model

Keywords: biological regulatory networks, piecewise affine differential models, limit cycles, structural principles, transition graph Mathematics Subject Classification (2000) MSC 34K13, MSC 92B05

), has been widely used for the modelling and the analysis of biological switch-like systems, such as genetic or neural networks. Its mathematical tractability facilitates the qualitative analysis of dynamical behaviors, in particular periodic phenomena which are of prime importance in biology. Notably, a discrete qualitative description of the dynamics, called the transition graph, can be directly associated to this class of PWA systems. Here we present a study of periodic behaviours (i.e. limit cycles) in a class of two-dimensional piecewise affine biological models. Using concavity and continuity properties of Poincaré maps, we derive structural principles linking the topology of the transition graph to the existence, number and stability of limit cycles. These results notably extend previous works on the investigation of structural principles to the case of unequal and regulated decay rates for the 2-dimensional case. Some numerical examples corresponding to minimal models of biological oscillators are treated to illustrate the use of theses structural principles.

Introduction

A large range of biological phenomena present switch-like behaviors of an almost on-off nature. Genetic regulation [START_REF] Ptashne | A Genetic Switch: Phage lambda and Higher Organisms[END_REF]) and neural response [START_REF] Mcculloch | A logical calculus of the ideas immanent in nervous activity[END_REF]), for example, have been shown to follow a nonlinear, switch like process. A class of piecewise-affine (PWA) differential models, well-suited for the modeling of switch-like systems, has been proposed by Glass and Kauffman in the 70s [START_REF] Glass | The logical analysis of continuous, non-linear biochemical control networks[END_REF]; Glass (1975a); Glass (1975b); Glass (1977a); Glass (1977b)). In this formalism, switching effects are represented by thresholds on the variables involved in the equations of evolution of the PWA model. These thresholds define domains in which the evolution of the variables is continuous and linear. This formalism thus represents an intermediate method in between the "classical" continuous differential approach and purely discrete formalism such as the logical method developed by Thomas and colleagues [START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF], [START_REF] Thomas | Biological Feedback[END_REF]).

This semi-qualitative modeling approach presents several advantages. First, it represents an interesting alternative to continuous differential approaches for biochemical networks modeling as the biochemical reaction mechanisms underlying the interactions are usually incompletely or not known, and the quantitative information on kinetic parameters and molecular concentrations are generally not available. In addition, compared to purely discrete approaches, this formalism allows to integrate in a more flexible way semi-quantitative data generated in biological systems. The PWA differential formalism has thus been widely applied to model several classes of biological systems behaving in a switch-like manner, mainly genetic networks (de [START_REF] De Jong | Qualitative simulation of genetic regulatory networks using piecewise-linear models[END_REF]; [START_REF] Ropers | Qualitative simulation of the carbon starvation response in escherichia coli[END_REF]; [START_REF] Omholt | Description and analysis of switchlike regulatory networks exemplified by a model of cellular iron homeostasis[END_REF]; [START_REF] Dayarian | Shape, size, and robustness: feasible regions in the parameter space of biochemical networks[END_REF]), neural networks [START_REF] Gedeon | Global dynamics of neural nets with infinite gain[END_REF]; [START_REF] Lewis | Nonlinear and symbolic dynamics of neural networks[END_REF]) or biochemical networks [START_REF] Glass | The logical analysis of continuous, non-linear biochemical control networks[END_REF]), but also food webs (Plahte et al (1995)).

Another advantage of this class of models (relative to continuous differential models) is its mathematical tractability [START_REF] Abou-Jaoudé | A theoretical exploration of birhythmicity in the p53-mdm2 network[END_REF]; Plahte et al (1995)). Indeed, these models have mathematical properties that facilitate qualitative analysis of asymptotic and transient behavior of regulatory systems. PWA differential equations have led to extensive work on the analysis of its dynamical properties [START_REF] De Jong | Qualitative simulation of genetic regulatory networks using piecewise-linear models[END_REF]). In particular, special focus has been put on the study of oscillatory behavior in such models (Glass and Pasternack (1978b); Glass and Pasternack (1978a); [START_REF] Mestl | Periodic solutions in systems of piecewise-linear differential equations[END_REF]; [START_REF] Edwards | Analysis of continuous-time switching networks[END_REF]; [START_REF] Farcot | Geometric properties of a class of piecewise affine biological network models[END_REF]; [START_REF] Farcot | Periodic solutions of piecewise affine gene network models with non uniform decay rates: the case of a negative feedback loop[END_REF]; [START_REF] Edwards | Control design for sustained oscillation in a two-gene regulatory network[END_REF]; [START_REF] Lu | Structural principles for periodic orbits in glass networks[END_REF]; [START_REF] Lu | Structural principles for complex dynamics in glass networks[END_REF]).

Periodic phenomena are of particular importance in biology, notably in cellular regulatory networks. Cellular oscillations have been reported in various biochemical systems such as calcium signalling, circadian rhythms, cell cycle [START_REF] Goldbeter | Biochemical Oscillations and Cellular Rythms[END_REF]; [START_REF] Goldbeter | Computational approaches to cellular rhythms[END_REF]) or in the p53-Mdm2 network [START_REF] Bar-Or | Generation of oscillations by the p53-mdm2 feedback loop: a theoretical and experimental study[END_REF]). Often, biochemical oscillations are characterized by a simple pat-tern with a single oscillatory regime of stable period and amplitude. However, more complex oscillatory patterns, like birhythmicity or chaos, have been proposed to occur in biochemical networks [START_REF] Abou-Jaoudé | From structure to dynamics: frequency tuning in the p53-mdm2 network i. logical approach[END_REF]; [START_REF] Decroly | Birhythmicity, chaos, and other patterns of temporal selforganization in a multiply regulated biochemical system[END_REF]). Therefore, considering the importance of cellular rhythms in biology, predicting the conditions of emergence of oscillatory behavior in mathematical models of biological systems is of great relevance.

Since its introduction, several results have been obtained on the existence and stability of limit cycles in PWA differential models. Most of this work focused on the situation where the decay rates of the system are equal. In this case, trajectories in each domain delimited by the thresholds are straight lines and one can derive an expression of a first return map as a linear fractional map and perform an eigenvalues analysis to study the existence and stability of periodic orbits. This method for the analysis of periodic orbits was first introduced by Glass and Pasternack (Glass and Pasternack (1978b); Glass and Pasternack (1978a)) and subsequently improved by several authors [START_REF] Edwards | Analysis of continuous-time switching networks[END_REF]; [START_REF] Mestl | Periodic solutions in systems of piecewise-linear differential equations[END_REF]; [START_REF] Farcot | Geometric properties of a class of piecewise affine biological network models[END_REF]).

In particular, some of this work on the analysis of PWA models with equal decay rates led to the finding of structural principles [START_REF] Lu | Structural principles for periodic orbits in glass networks[END_REF]) linking the topology of the transition graph (i.e. the graph showing all the state domains and possible transitions between them) to the dynamics of the system, more specifically its oscillatory behavior. The first structural principles, derived by Glass and Pasternack (Glass and Pasternack (1978b)), apply to a class of configurations in the state transition graph called cyclic attractors, in which the successor of each state is unambiguous. Those authors proved that all trajectories in the regions of phase space corresponding to the cyclic attractor either converge to a unique stable limit cycle, or approach a stable equilibrium. Lately, other methods, based on an appropriate construction of focal points giving rise to the desired orbit, have been used by Lu and Edwards to derive structural principles. They notably showed that, for specific classes of cycles in the state space, there exist parameter values such that a periodic orbit can exist [START_REF] Lu | Structural principles for periodic orbits in glass networks[END_REF]).

Recently, Farcot and Gouzé proposed a new approach to analyse the existence of limit cycles in PWA differential models in the case of non-equal decay rates [START_REF] Farcot | Periodic solutions of piecewise affine gene network models with non uniform decay rates: the case of a negative feedback loop[END_REF]; [START_REF] Farcot | Limit cycles in piecewise-affine gene network models with multiple interaction loops[END_REF]), using tools from the theory of monotone systems and operators [START_REF] Smith | Cooperative systems of differential equations with concave nonlinearities[END_REF]). The study of this case is of prime relevance as it encompasses a wide range of biological contexts. Indeed, decay rates set the timescales of the consumption processes which can be very different from one biological component to the other. Their analysis was based on the monotonicity and concavity of a first return map under some specific constraints on the parameters of the system (i.e. alignment conditions of successive focal points). This method was first applied to PWA systems containing a single negative feedback loop [START_REF] Farcot | Periodic solutions of piecewise affine gene network models with non uniform decay rates: the case of a negative feedback loop[END_REF]) and successfully extended to other PWA models with multiple inter-action loops [START_REF] Farcot | Limit cycles in piecewise-affine gene network models with multiple interaction loops[END_REF]), both classes of models verifying the alignment conditions. In particular, these results allowed to extend Snoussi's theorem, stating solely existence of limit cycles, to the existence and uniqueness of limit cycles when the interaction graph consists in a single negative feedback loop [START_REF] Snoussi | Qualitative dynamics of piecewise-linear differential equations: a discrete mapping approach[END_REF]).

The aim of this paper is to propose a study, inspired by Farcot and Gouzé's approach, of the existence, number and stability of limit cycles in a class of 2dimensional PWA biological systems. This work represents, to our knowledge, the first investigation of structural principles in PWA differential models with unequal and regulated decay rates (regulated meaning that these rates may vary with the domain). The case of equal decay rates is also revisited using Farcot and Gouzé's method. Starting from a case study of a minimal PWA model for biological oscillators (section 2), general properties of monotonicity, concavity and continuity of first return maps are derived to prove structural principles on the existence, number and stability of limit cycles (sections 3 and 4). The oscillatory behavior of the case study and another example of biological oscillator is then revisited using our theoretical results (section 5). In the discussion, we also make some links with ordinary differential equations, for which in general it is not possible to obtain such detailed results.

A minimal piecewise affine model for biological oscillators

Cellular regulatory networks contain multiple positive and negative feedback loops to ensure an appropriate regulation of its behaviour. Whereas positive circuits are involved in differentiation and memory processes, negative ones are crucial to maintain homeostasis and set biological rhythms [START_REF] Thomas | Biological Feedback[END_REF]). Cellular rhytms are present in important biological phenomena like calcium signalling, circadian rhythms, cell cycle [START_REF] Goldbeter | Biochemical Oscillations and Cellular Rythms[END_REF]; [START_REF] Goldbeter | Computational approaches to cellular rhythms[END_REF]) or in the p53-Mdm2 network [START_REF] Bar-Or | Generation of oscillations by the p53-mdm2 feedback loop: a theoretical and experimental study[END_REF]). Interestingly, these two types of circuits can act in concert to form building blocks of cellular regulatory networks and ensure robust cellular rhythms [START_REF] Kim | Coupled positive and negative feedback circuits form an essential building block of cellular signaling pathways[END_REF]; [START_REF] Tsai | Robust, tunable biological oscillations from interlinked positive and negative feedback loops[END_REF]).

A broad class of biological systems, among which genetic and neural networks [START_REF] Ptashne | A Genetic Switch: Phage lambda and Higher Organisms[END_REF]; [START_REF] Mcculloch | A logical calculus of the ideas immanent in nervous activity[END_REF]) and cell signaling pathways [START_REF] Ferrell | Tripping the switch fantastic: how a protein kinase cascade can convert graded inputs into switch-like outputs[END_REF]), are characterized by switch-like behaviors. A formalism well suited to model such systems, initially proposed by Glass and Kauffman [START_REF] Glass | The logical analysis of continuous, non-linear biochemical control networks[END_REF]), is the piecewise affine (PWA) differential framework where the biological processes are represented by step functions (Fig. 1a):

{ s + (x, θ) = 0 if x < θ s + (x, θ) = 1 if x > θ
for activation processes and:

{ s -(x, θ) = 1 if x < θ s -(x, θ) = 0 if x > θ
for inhibition processes, where θ is the process threshold and x the level of the biological component regulating the process.

We now consider the following class of 2-dimensional PWA model for biological oscillators:

     dx dt = k 1x • s + (y, θ 2 y ) -d x • x dy dt = k 1y • s -(x, θ x ) + k 2y • s + (y, θ 1 y ) + k 3y • s + (y, θ 2 y ) -d y • y (1)
This model is composed of two auto-regulatory feedback loops on y which encompass the two step processes, k 2y • s + (y, θ 1 y ) and k 3y • s + (y, θ 2 y )), and one two-element negative feedback loop between the two components x and y of the model (Fig. 1b). d x • x and d y • y represent the linear decay processes of x and y respectively. This class of model is a general minimal model which can be used to represent biological oscillators composed of intertwined positive and negative feedback loops (start of cell cycle in budding yeast, calcium oscillations,...) [START_REF] Kim | Coupled positive and negative feedback circuits form an essential building block of cellular signaling pathways[END_REF]). For the calcium oscillations model of Keizer et al. (Kim et al (2007); [START_REF] Keizer | Insp3-induced Ca 2+ excitability of the endoplasmic reticulum[END_REF]), x represents the SERCA ATPases pumping calcium out into the endoplasmic reticulum lumen and y the level of cytoplasmic calcium (cytCa). The two step processes which form the positive feedback loops of the model represent the IP3R-cytCa and RYR-cytCa circuits which are activated to increase cytoplasmic calcium. The negative feedback loop models the regulation between the SERCA ATPases and cytoplasmic calcium.

Although minimal, this 2-dimensional class of models can already account for a rich variety of dynamical phenomena. For appropriate parameter settings, this model displays different types of oscillatory behavior, from simple to more complex patterns, among which: damped oscillations towards a single equilibrium point (simulation not shown), a single oscillatory regime of stable period and amplitude (Fig. 2a), birhythmicity with the coexistence of two stable oscillatory regime (Fig. 2b) separated by an unstable limit cycle (not shown). As precise quantitative information on kinetic parameters are generally missing, predicting the conditions of emergence of oscillatory patterns from a qualitative description of biological models, inferred from qualitative information on the parameter values, is of particular interest. In the following, we derive theoretical results which set constraints on the oscillatory patterns which can emerge in a general class of 2-D PWA biological models according to topological properties of its transition graph, a discrete qualitative description of this type of models. These structural principles will be derived using concavity and continuity properties of Poincaré maps associated to cycles of the transition graph (section 4). As a preliminary step, we first define the class of PWA biological models on which our analysis will be applied and the theoretical objects and tools which will be used to state our results (section 3).

Piecewise affine differential models

The general model

This paper focuses on the study of a general class of 2-dimensional piecewise affine (PWA) systems, with positive variables, (x, y)

∈ ℜ + × ℜ + , ℜ + = [0, ∞[:      dx dt = f x (x, y) -g x (x, y) • x dy dt = f y (x, y) -g y (x, y) • y (2)
where f i : ℜ + × ℜ + → ℜ + and g i : ℜ + × ℜ + → ℜ * + (i ∈ {x, y}) are piecewise constant functions, representing the interactions between the various components of the system. The degradation of each component is assumed to be a linear process regulated by the components of the system.

The analysis of the system can be restricted to the phase space region: [0, M x ]× [0, M y ], with M x = max{ fx(x,y) gx(x,y) : (x, y) ≥ 0} and M y = max{ fy(x,y) gy(x,y) : (x, y) ≥ 0}, which defines a compact set that all trajectories enter and never leave (de Jong et al ( 2004)).

The variables x and y of the system are each associated with thresholds which set the switching values of the vector fields. We assume that x and y have n x and n y thresholds respectively:

0 < θ 1 x < . . . < θ nx x < M x 0 < θ 1 y < . . . < θ ny y < M y with θ 0 x = θ 0 y = 0.
For convenience of notation, M x and M y will be renamed θ nx+1

x and θ ny+1 y respectively.

These thresholds partition the state space into (n x + 1) • (n y + 1) regular domains in which the vector fields are given by an affine function. These domains will be labelled using the following notation:

D ij : (i, j) ∈ {1, 2, . . . , n x + 1} × {1, 2, . . . , n y + 1}, θ i-1 x < x < θ i x and θ i-1 y < y < θ i y
The segments and threshold intersections defining the borders of the regular domains are called switching domains. Such segments will be called switching segments and threshold intersections will be renamed switching points.

In each regular domain D ij , the functions f x , f y , g x and g y take constant values and the system can be rewritten as follows:

     dx dt = k ij x -d ij x x dy dt = k ij y -d ij y y (3)
where

f x = k ij x , f y = k ij y , g x = d ij x and g y = d ij y for (x, y) ∈ D ij .
From Equations 3, we can define the so called focal points which are the points, ϕ ij , towards which the system tends monotically from each domain D ij (Glass and Pasternack (1978b)):

ϕ ij = ( k ij x d ij x , k ij y d ij y )
.

The solutions of this system can be explicitly written:

{ x(t) = ( x(0) -ϕ ij x ) • exp -d ij x t +ϕ ij x y(t) = ( y(0) -ϕ ij y ) • exp -d ij y t +ϕ ij y (4)
The equation of the trajectory in D ij can furthermore be derived by eliminating time t. From Equations 4, we obtain:

y(t) = ( y(0) -ϕ ij y ) • ( x(t) -ϕ ij x x(0) -ϕ ij x ) d ij y d ij x + ϕ ij y (5)
which defines the equation of the trajectory of the system in domain D ij . In the case of equal decay rates, d ij x = d ij y = d ij for all (i, j) and the trajectory in D ij is reduced to a straight line.

Throughout the paper, we will make the following assumptions: Assumption 1. The focal points do not belong to the switching domains of the phase space. Therefore, if a regular domain D ij does not contain its focal point, the system will eventually escape this domain.

Assumption 2. Switching segments reached by a trajectory from a regular domain are transparent walls i.e. the flow in these segments is well defined. Following Assumption 2, a trajectory reaching a switching segment can evolve into a contiguous regulatory domain, thus originating a transition between two regular domains. Stable solutions on switching segments (i.e called stable sliding modes in control theory [START_REF] Casey | Piecewise-linear models of genetic regulatory networks: equilibria and their stability[END_REF]) are thus excluded by this assumption. However, it may not be possible to continuously extend solutions reaching switching points. An approach originally proposed by [START_REF] Filippov | Differential equations with discontinuous righthand sides[END_REF]) and more recently applied to PWA systems [START_REF] Gouzé | A class of piecewise linear differential equations arising in biological models[END_REF]) can then be used to define the solutions on switching points (see proof of Theorem 4 in Appendix A). An important consequence of Assumption 2 is the following lemma:

Lemma 1 For any initial condition, a solution of (2) which does not cross switching points is unique.

The proof of this Lemma can be found in Appendix A. This property will be notably used in the analysis of the continuity of first return maps (see section 4.2).

D ij D ij D ij D ij D ij D ij D ij D ij

Transition graph

A discrete, qualitative description of the dynamics of a PWA system, initially proposed by Glass (Glass (1975a)), is called the transition graph: it is a directed graph whose vertices are the regular domains of the system and whose edges represent the possible transitions between these domains. The transition graph is obtained from the positions of the focal points.

Due to Assumption 2, each D ij will have either zero, one or two successors depending on the position of its focal point ϕ ij : D ij has no successor if ϕ ij belongs to D ij , and one successor (resp. two successors) if ϕ ij belongs to a contiguous regular domain which shares a switching segment (resp. a switching point) boundary in common with D ij . The different types of escaping transition configurations from a vertex of the transition graph are summed up in Fig. 3. A vertex which has 2 successors is called branching vertex and the corresponding domain will be called branching domain. Each branching vertex gives rise to a curve (called separatrix ) emerging from the switching point defined by the intersection of the two switching segments crossed by the transitions leaving the branching vertex. The separatrix curve thus corresponds to the trajectory of the phase space which reaches this switching point. It delimits the branching domain into two subsets from which the system will either enter one successor domain or the other.

Combining two of these transitions gives rise to 3-vertex paths in the transition graph which will be of special interest in this work. These paths can be classified in two categories: 3-vertex parallel paths, whose three regular domains are adjacent along two parallel switching segments, and 3-vertex perpendicular paths, whose three regular domains are adjacent along two perpendicular switching segments. These two classes of paths will be called parallel motifs and perpendicular motifs respectively. The different types of parallel and per-pendicular motifs are listed in Fig. 4 and 5. A trajectory passing through the domains composing a parallel (resp. perpendicular) motif will thus enter and escape the second domain through two parallel (resp. perpendicular) switching segments (Fig. 7). Finally, perpendicular motifs can be further classified into two subtypes: clockwise perpendicular motifs and counterclockwise perpendicular motifs (Fig. 5).

Transition cycles and n-cyclic attractors

We now introduce the following object. A transition cycle C of length n is defined as a periodic sequence of n vertices and n transitions in the transition graph:

D r1s1 → D r2s2 → . . . → D rnsn → D r1s1 with (r i , s i ) ∈ {1, . . . , n x + 1} × {1, . . . , n y + 1}
, each vertex of the sequence being connected to its successor by a transition and no vertex appearing more than once in the sequence (see also Glass and Pasternack (1978b)). Note that the existence of a transition cycle does not imply the existence of a limit cycle for trajectories.

As we are in dimension two, we can moreover define the inside and the outside of a transition cycle C. The inside (resp. outside) of a transition cycle is the set of domains (regular and/or switching) which are located inside (resp. outside) the transition cycle. Therefore, if a transition cycle contains a branching vertex, the transition from which the system can escape the cycle crosses a switching segment located either in the inside or the outside of the cycle. The former (resp. latter) type of transition will be called inside (resp. outside) branching transition of transition cycle C (see Fig. 10 for examples of inside and outside branching transitions). Transition cycles which do not contain one or the other type of branching transition will be considered when dealing with structural principles (section 4.3).

Transition cycles can be further classified in two broad categories according to the type of perpendicular motif composing the transition cycle: transition cycles which contain both clockwise and counterclockwise perpendicular motifs (which will be called transition cycles with turn change) and transition cycles which contain only clockwise or only counterclockwise perpendicular motif (which will be called transition cycle with no turn change, see Fig. 10 for examples). The relevance of this classification will appear in section 4.2.

An important class of cycles called cyclic attractors has been proposed by Glass and Pasternack (Glass and Pasternack (1978b)), which are cycles that do not contain branching vertices. We extend this notion of cyclic attractor to that of n-cyclic attractor. A n-cyclic attractor C n is defined as the union of n transition cycles C i of the transition graph: C n = ∪C i which does not contain vertices from which the system can escape the union of cycles. Thus, once a trajectory enters the union of the domains composing C n , it will remain in these domains. Cyclic attractors defined by Glass and Pasternack are therefore 1-cyclic attractors.

D ij D i(j+1) D i(j-1) D ij D i(j+1) D i(j-1) D ij D (i-1)j D (i+1)j D ij D (i-1)j D (i+1)j Fig. 4 Parallel motifs D ij D i(j+1) D (i+1)j D ij D i(j+1) D (i+1)j D ij D (i-1)j D i(j-1) D ij D (i-1)j D i(j-1) D i(j-1) D ij D (i+1)j D i(j-1) D ij D (i+1)j D ij D i(j+1) D (i-1)j D ij D i(j+1) D (i-1)j
In this paper, we will limit the scope of our work to transition-connected ncyclic attractors, which are n-cyclic attractors whose transition cycles share at least one transition in common (see Fig. 10 and 12 for examples). For sake of simplicity, transition-connected n-cyclic attractors will be renamed n-cyclic attractors.

Elementary maps and first return maps

Elementary maps

Since we assume there is no solution along switching segments (Assumption 2), a trajectory reaching the boundary of a regular domain will evolve into an adjacent domain by simply crossing the switching segments that separate the two domains. Therefore, to each 3-vertex path D k-1 → D k → D k+1 contained in the transition graph, we can define an elementary map

F k (F k : ℜ + × ℜ + → ℜ + × ℜ + ) which maps the entering switching segment D k-1 s
to the escaping switching segment D k s which border D k (see also [START_REF] Edwards | Analysis of continuous-time switching networks[END_REF]). The image of a point of D k-1 s by the map F k is defined as the intersection of the trajectory starting from this point with D k s . Note that (in this 2-dimensional framework), each map F k has a fixed coordinate at both the entering and escaping segments (see Fig. 6).

In order to compute elementary maps, we define a scalar function f k (f k : ℜ + → ℜ + ) associated to the elementary map F k by setting the direction and the origin of the axes supporting the entering and escaping switching segment, D k-1 s and D k s , where f k is computed. For each f k , two possible orientations can be chosen for the axis such that the origin coincides with either of the switching point extremities, while the other extremity is assumed to be positive (see Fig. 6). Once the origin and direction of the axis are set, the switching segment is said to be oriented. For convenience of notation, f k will also be called elementary map of the 3-vertex path D k-1 → D k → D k+1 and the axes where f k is computed will be called entering and escaping axes of f k .

Let I k be the interval of definition of f k and l k the length of the entering switching segment of F k . First, f k can be straightforwardly extended to the boundaries of I k which are switching points. Then, if D k is a not a branching domain, the escaping switching segment from D k is unambiguous: all the points of D k-1 s map D k s . In this case,

I k = [0, l k ]. If D k is a

First return map of a transition cycle

We now assume that the transition graph of a PWA system contains a transition cycle C of length n:

D t1u1 → D t2u2 → . . . → D tnun → D t1u1
Fig. 6 Two choices of orientation for the axes of an elementary map of a 3-vertex parallel path D (i-1)j → D ij → D (i+1)j . Two possible orientations can be chosen for each axis such that the origin coincides with either of the switching point extremities, while the other extremity is assumed to be positive. Red arrows represent the transition of the 3-vertex path. Blue ones are the entering and escaping axes where the elementary map of the 3-vertex path is computed.

with (t i , u i ) ∈ {1, 2, . . . , n x + 1} × {1, 2, . . . , n y + 1}.
Given a switching segment crossed by C, the associated first return map F (or Poincaré map) (F :

ℜ + × ℜ + → ℜ + × ℜ +
) is a mapping from this segment to itself, computed from two consecutive crossings of a trajectory of the system with this segment [START_REF] Edwards | Analysis of continuous-time switching networks[END_REF]). Let D i s be the switching segment crossed by the transition D ti-1ui-1 → D tiui for i ∈ {2, . . . , n} and D 1 s the switching segment crossed by the transition from D tnun to D t1u1 .

If F i is the elementary map associated with the 3-vertex path:

D tnun → D t1u1 → D t2u2 for i = 1, D ti-1ui-1 → D tiui → D ti+1ui+1 for i ∈ {2, . . . , n -1}
and D tn-1un-1 → D tnun → D t1u1 for i = n, the first return map F of the transition cycle C from and to the switching segment D 1 s is the composite of the elementary maps F i for i ∈ {1, 2, . . . , n}:

F = F n • F n-1 • . . . F 1
The domain of definition of F is called the returning cone of F [START_REF] Edwards | Analysis of continuous-time switching networks[END_REF]).

As for an elementary map, in order to compute the first return map of C, we can define a scalar function f (f : ℜ + → ℜ + ) associated with the first return map F by setting the direction and the origin of the axis supporting the switching segment D 1 s where f is computed. The origin of the axis is set in the same manner as the origin of the axes of elementary maps (see previous section). f will also be called the first return map from and to D 1 s of transition cycle C.

Let the escaping axis of D tiui and the entering axis of D ti+1ui+1 have the same orientation. Let f i be the elementary map associated with the 3-vertex path:

D tnun → D t1u1 → D t2u2 for i = 1, D ti-1ui-1 → D tiui → D ti+1ui+1 for i ∈ {2, . . . , n -1} and D tn-1un-1 → D tnun → D t1u1 for i = n. f is then the composite of the elementary maps f i for i ∈ {1, 2, . . . , n}: f = f n • f n-1 • . . . f 1

First return map of a n-cyclic attractor

We assume that the transition graph contains an n-cyclic attractor C n composed of the n transition cycles C i for i ∈ {1, 2, . . . n}. Let D s be an oriented switching segment crossed by a transition common to all C i and let l be the length of the segment D s .

Let D bv1 , D bv2 , . . . , D bvm be the branching vertices of C n . Each branching vertex D bvi gives rise to a separatrix curve. Let x 1 s < x 2 s < . . . < x p s (p m) be the coordinates of the last intersection (if it exists) between the trajectories lying in the separatrix curves and D s before these trajectories reach a switching point. The x i s partition D s into (p + 1) segments I i :

I 1 = [ 0, x 1 s [ , I i = ] x i-1 s , x i s [ for 2 ≤ i ≤ p and I p+1 = ]x p s , l].
Each trajectory starting from a point whose coordinate belongs to I i will thus follow a distinct transition cycle C ui in the transition graph before a first return in D s . Note that all the cycles composing C n are not necessarily followed by a trajectory starting from D s .

I i then defines the interval of definition of the first return map f i associated to C ui from and to the oriented switching segment D s . We can then define a map f on the union of intervals I = ∪I i as follows:

for i ∈ {1, 2, . . . , n}: f (x) = f i (x) for x ∈ I i f will be called n-cycle first return map of the n-cyclic attractor C n from and to the oriented switching segment D s .

In the following, we will limit our work to the case n ≤ 2. In the case of a 1-cyclic attractor (which corresponds to a cyclic attractor as defined by Glass and Pasternack (1978b)), there is no branching vertex and I = [0, l]. In the case of a 2-cyclic attractor, we state the following lemmas (see Appendix A for the derivation of the proofs). Let C 2 be a 2-cyclic attractor and D C 2 the union of the domains composing C 2 .

Lemma 2 C 2 contains a single branching vertex.

There is therefore a single separatrix curve in the domains composing a 2cyclic attractor which emerges from the unique branching vertex. There is also a unique vertex where the transition cycles composing C 2 merge. The properties of the subgraph composed of the two distinct pathways linking the branching and the merging vertices will be of special interest when dealing with the continuity of the first return map of 2-cyclic attractors (section 4.2).

Lemma 3 Two different trajectories cannot intersect in D C 2 .
This Lemma will be notably used in the proof of Theorem 4.

Results

Monotonicity and concavity properties of an elementary map

In this section, the monotonicity and the concavity properties of the elementary maps f k of the different types of parallel and perpendicular motif (listed in Fig. 4 and5) are studied. The proofs of the following statements can be found in Appendix A.

We first state a lemma which will be used to derive the monotonicity and concavity properties of elementary maps established in Propositions 1 and 2.

Lemma 4 Let f k be an elementary map. Assume f k is monotone and of constant concavity. Then, changing the orientation of the entering axis of f k changes the monotonicity but does not change the concavity of f k . Changing the orientation of the escaping axis of f k changes the monotonicity and concavity of f k .

Proposition 1 Let f k be the elementary map of a parallel motif. Then f k is an affine function. If the entering and escaping axes of f k have the same orientation, f k is an increasing function. Otherwise, f k is a decreasing function.

Proposition 2 Let f k be the elementary map of a perpendicular motif and S k the switching point at the intersection of the entering and escaping segment of f k .

(1) If the entering and escaping axes of f k are both oriented either towards or away from S k , f k is increasing. Otherwise f k is decreasing.

(2) If the origin of the escaping axis is S k , f k is strictly concave. Otherwise, it is strictly convex.

Note that elementary maps are continuous functions (see analytical expression of f k in the proofs of Propositions 1 and 2). These results on the monotonicity and concavity of elementary maps will be used in the next section to derive monotonicity and concavity properties of first return maps. We now state general properties of first return maps concerning their monotonicity, concavity and continuity. Theorems 1, 2 and 3 concern the properties of monotonicity and concavity of the first return map of a transition cycle. Theorem 4 deals with continuity properties of 2-cycle first return maps. These properties will be used in the next section to derive structural principles linking the topology of the cycles in the transition graph to the number and stability of limit cycles. The proofs of the theorems can be found in Appendix A.

For Theorems 1, 2 and 3, we assume that the transition graph contains a transition cycle C. We now assume that the transition graph contains a 2-cyclic attractor C 2 , composed of two transition cycles C 1 and C 2 . According to Lemma 2, C 2 contains a single branching and a single merging vertex. Let D S and D T be the branching and merging vertices respectively of C 2 and SG S→T the subgraph composed of the two distinct pathways linking D S to D T .

Let S be the switching point from which emerges the separatrix curve in D S (see Fig. 9). Let D C be an oriented switching segment crossed by a common transition of C 1 and C 2 . Assume that the separatrix curve emerging from the branching vertex of C 2 intersects D C and let α be the coordinate of the last intersection between D C and the trajectory lying in the separatrix curve before it reaches S. Let f be the 2-cycle first return map of C 2 calculated on the oriented segment

D C : { f (x) = f 1 (x) for x ∈ [0, α[ f (x) = f 2 (x) for x ∈ ]α, l] (6)
where f 1 and f 2 are the first return maps of C 1 and C 2 respectively calculated on D C , and l the length of D C . The following theorem states the topological conditions on a 2-cyclic attractor for its 2-cycle first return map to be extended to a continuous function. x -θ x > 0 and ϕ 1 y -θ y < 0). In the left case, D 2 and D 3 communicate with D 4 (ϕ 2

x -θ x > 0, ϕ 3 y -θ y < 0) leading to a 4-node SG S→T , whereas in the right case D 4 does not communicate with D 2 (ϕ 2

x -θ x < 0 and ϕ 4

x -θ x > 0) leading to a SG S→T composed of more than 4 nodes. According to Theorem 4, the left case gives a continuous first return map, which is discontinuous in the right case.

Structural principles

From the previous theorems, we can derive structural principles, i.e. properties which emerge from the structure of the transition graph, on the number and stability of the limit cycles [START_REF] Lu | Structural principles for periodic orbits in glass networks[END_REF]). These structural principles are limited to the cases where either C has no turn change or the decay rates are equal, which are cases for which information about the concavity of the first return map have been obtained (see Theorems 2 and 3). The analysis of the number of limit cycles (i.e. fixed points of the first return maps) is only considered in the generic case, when the intersections between the first return map and the identity are transverse (see Fig. 14). The first structural principle is a general principle which states the maximum number of stable and unstable limit cycles a system can have in the domains of the phase space crossed by C. The second and third structural principles state constraints on the number and stability of limit cycles for specific topological properties of the transition graph. Finally, the last theorem concerns the number of unstable limit cycles for a 2-cyclic attractor when specific topological conditions on the structure of the transition graph are fulfilled (see Theorem 4). The proofs of the following statements can be found in Appendix A.

For the first three structural principles (Theorems 5, 6 and 7), we assume that the transition graph contains a transition cycle C and call D C the union of the domains composing C. Let f be a first return map of C, I the interval of definition of f and l the length of the switching segment where f is calculated.

Theorem 5 (First structural principle) Assume that either (1) C is a transition cycle with no turn change or (2) the decay rates are equal. Then there are at most two limit cycles in D C . If there are two limit cycles in D C , one is stable and the other is unstable.

The following lemma will be used to prove the next two structural principles.

Lemma 5 We assume that the axis where f is calculated is oriented towards the outside of C. Theorem 6 (Second structural principle) Assume that C is a transition cycle with no turn change. If C contains no inside branching transition, the system does not admit an unstable limit cycle in D C but either (1) has no limit cycle or (2) has one single stable limit cycle in D C . Theorem 7 (Third structural principle) Assume that either (1) C has no turn change or (2) the decay rates are equal. If C contains no outside branching transition and if the system admits an unstable limit cycle in D C then it admits a single stable limit cycle in D C .

Note that information about the concavity of the first return map is required for Theorem 6 whereas only a constant concavity is needed for Theorem 5 and 7, which explains why Theorem 6 only applies to the case where C has no turn change and not the case where the decay rates are equal.

For the last structural principle, we assume that the transition graph contains a 2-cyclic attractor C 2 , composed of transition cycles C 1 and C 2 . Let D C be a switching segment crossed by a common transition to C 1 and C 2 , and D C 2 the union of the domains composing C 2 . Let SG S→T be the subgraph composed of the two distinct pathways linking the branching vertex of C 2 to the vertex where C 1 and C 2 merge.

Theorem 8 (Fourth structural principle) Assume that either (1) C 1 and C 2 have no turn change or (2) the decay rates are equal. Assume also that SG S→T is composed of 4 vertices. If the system admits two stable limit cycles in D C 2 , there exists a unique unstable limit cycle in D C 2 which delimits the basin of attraction of the two stables limit cycles.

Applications

Two applications of the previous theoretical results will be described in this section. Both applications represent PWA models of biological oscillators whose transition graph is composed of a 2-cyclic attractor. All the possible dynamical configurations in terms of the number and the stability of the limit cycles are derived using the structural principles stated in the previous section. In the first application, which corresponds to a case where the 2-cycle first return map is continuous, the oscillatory behavior of the case study presented in section 2 is revisited in the light of our theoretical results. The second application concerns a more complex biological oscillator whose transition graph reproduces the one of a reduced version of the p53-Mdm2 network model (Abou-Jaoudé et al ( 2009); Abou-Jaoudé et al ( 2011)) and for which the 2-cycle first return map is discontinuous.

A continuous first return map

We consider the class of 2D-piecewise affine differential models introduced in section 2 which represents a minimal model for biological oscillators (see Equations 1). The interaction graph, shown in Fig. 1b, is composed of one 2-element negative feedback loop and two 1-element positive feedback loops. The step functions partition the phase space into 6 domains delimited by threshold θ x along x, and thresholds θ 1 y and θ 2 y along y. The parameter value sets have been chosen such that transition graph is the 2-cyclic attractor indicated in Let D C1 (resp. D C2 ) be the union of the domains composing C 1 (resp. C 2 ). C 1 and C 2 have no turn change. Moreover, cycle C 1 has no inside branching transition. Therefore, according to the second structural principle, the system admits no unstable limit cycle and at most one stable limit cycle in D C1 . Cycle C 2 has no outside branching transition. Thus, according to third structural principle, if the system admits an unstable limit cycle in D C2 , it admits a stable limit cycle in D C2 . At most two stable limit cycles and one unstable limit cycle is therefore possible in the union of domains D C1 and D C2 .

Finally, the subgraph composed of the two paths linking the branching vertex to the vertex where both cycles merge contains 4 vertices: D 22 , D 21 , D 11 and D 12 . Therefore, according to the fourth structural principle, if the system admits two stable limit cycles, there exists a unique unstable limit cycle which delimits the basin of attraction of the two stable limit cycles.

Thus, the dynamical configurations of this system in terms of the number and the stability of the limit cycles are limited to at most the following five cases: (a) one large stable limit cycle located in D C2 ; (b) one small stable limit cycle located in D C1 ; (c) one large stable limit cycle and one unstable limit cycle both located in D C2 ; (d) one small stable limit cycle located in D C1 , and one large stable limit cycle and one unstable limit cycle located D C2 . The unstable limit cycle delimits the basins of attraction of the stable limit cycles; (e) no limit cycle.

We could actually find parameter value sets corresponding to each of the five dynamical configurations stated above. Fig. 11 shows numerical simulations in the phase space corresponding to cases (a),(b),(c) and (d) (cases for which there is at least one limit cycle). The parameter sets for the two oscillatory patterns displayed by the case study in Fig. 2 correspond to cases (a) and (d). The corresponding 2-cycle return maps are described in Appendix B.

Fig. 11

Numerical simulations for the continuous case corresponding to the possible dynamical configurations derived from the structural principles: case with one large stable limit cycle (top left); one small stable limit cycle (top right); one large stable limit cycle and one unstable limit cycle (bottom left); one large and one small stable limit cycles and one unstable limit cycle (bottom right). The stable (resp. unstable) limit cycles are drawn in blue (resp. in red).

A discontinuous first return map

We now consider the following class of 2D-piecewise affine differential models:

           dx dt = [ k 1x + k 2x • s + (y, θ 1 y ) ] • s -(y, θ 2 y ) • s + (x, θ x ) + s + (y, θ 3 y )) • [ k 3x + k 4x • s + (x, θ x ) ] -d x • x dy dt = k 1y • s -(x, θ x ) -d y • y (7)
where s + and s -are step functions as previously defined in section 2.

The interaction graph, shown in Fig. 12, is composed of two 2-element negative feedback loops and two positive feedback loops (one 1-element and one 2-element). The step functions partition the phase space into 8 domains delimited by threshold θ x along x, and thresholds θ 1 y , θ 2 y and θ 3 y along y. The parameter value sets have been chosen such that transition graph is the 2-cyclic attractor indicated in Let D C1 (resp. D C2 ) be the union of the domains composing C 1 (resp. C 2 ). As for the continuous case, both cycles have no turn change. Cycle C 1 has no inside branching transition while cycle C 2 has no outside branching transition. The second and third structural principles can thus both be applied for C 1 and C 2 respectively. However, the subgraph composed of the two paths linking the branching vertex to the vertex where both cycles merge contains now 6 vertices: D 23 , D 22 , D 21 , D 11 , D 12 and D 13 . Therefore, the fourth structural principle cannot be applied. One supplementary dynamical configurations in terms of the number and the stability of the limit cycles could thus arise in addition to the five dynamical configurations which appear in the previous example: (f) two stable limit cycles, one located in D C1 , the other in D C2 . In this additional configuration, the two stable limit cycles are not separated by an unstable limit cycle but by the separatrix curve emerging in the branching vertex D 23 (Fig. 12) from the threshold intersection (θ x , θ 2 y ).

We could actually find parameter value sets corresponding to each of the six possible dynamical configurations. Fig. 13 shows numerical simulations corresponding to cases (a),(b),(c), (d) and (f) (cases for which there is at least one limit cycle). The corresponding 2-cycle return maps are listed in Appendix B. Fig. 10B of Abou-Jaoudé et al (2011) gives another numerical example of the case where the system has two stable limit cycles and no unstable limit cycle (case (f)).

Discussion

In this work, we derived structural principles linking the topology of the transition graph of a class of 2-dimensional piecewise affine biological models to the number and stability of limit cycles. To do so, we analyzed the continuity, monotonicity and concavity properties of Poincaré maps associated to the transition cycles of the transition graph. In the case of nonequal decay rates, structural principles linking the topology of a transition cycle to the dynamics of the system in terms of number and stability of limit cycles have been determined when the transition graph contains no turn change. For 2-cyclic attractors in the transition graph, structural principles have been derived on the number of unstable limit cycles from the continuity of the first return map associated to the attractor. The results of our work have then been applied to two biological cases whose transition graph are 2-cyclic attractors: a case for which the first return map is continuous, a case for which the first return map is discontinuous.

Fig. 13

Numerical simulations for the discontinuous case corresponding to the possible dynamical configurations derived from the structural principles: case with one large stable limit cycle (top left); one small stable limit cycle (top right); one large stable limit cycle and one unstable limit cycle (middle left); one small and one large stable limit cycle and one unstable limit cycle (middle right); one small and one large stable limit cycles (bottom left). Stable (resp. The mathematical tractability of the class of PWA biological models analyzed in our work allows to derive the stated structural conditions on the number and stability of limit cycles. Such results on the number of limit cycles cannot be derived in ODE systems. For planar ODE systems, Poincaré-Bendixson theorem gives mathematical conditions for the existence of limit cycles but there is no general theoretical results on the number of limit cycles. More specific results concern the number or stability of limit cycles, but for more specific systems (Lienard systems, polynomial systems...) and may be rather complex (see (Perko, 1991, p. 234) for references). To investigate whether the structural results obtained for the class of PWA models considered in our work are conserved in the continuous framework, we translated the PWA system of the discontinuous case (section 5.2) into an ODE model by replacing step functions with Hill functions. Interestingly, for appropriate parameter values, a first numerical analysis of this model suggests that the discontinuity observed in the PWA model has some counterpart in the continuous framework (Appendix C).

The extension of our results to higher dimension PWA systems seems to be difficult to achieve since the conclusions stated in our work strongly rely on the topological constraint imposed by the 2-dimensional space. However, several developments of this work can be considered. First, for unequal decay rates, the structural principles derived here are restricted to the case where transition cycles have no turn change. Under this condition, properties on the concavity of first return maps associated with transition cycles can be derived. These properties are then used to determine the structural principles on the number and the stability of limit cycles in the domain of the phase plane covered by a transition cycle of this type. An extension of the concavity results when transition cycles contain a turn change still has to be investigated. First results on this issue show that we can conclude on the monotonicity and concavity properties for some specific topological structures of the transition graph (results not shown). Secondly, this work focused on the continuity properties of 2-cycle first return maps. An extension of these results to n-cycle first return maps associated with n-cyclic attractors for n ≥ 3 can also be investigated.

Finally, applications of these structural principles contribute to gain intuition on the dynamical behavior of biological networks and provide guidance on parameters or experimental conditions that generate a given behavior.

Acknowledgements

We would like to thank Denis Thieffry for critical reading of the manuscript, and also on of the reviewers for his careful reading and many useful suggestions. W. Abou-Jaoudé was supported in part by the LabEx MemoLife (www.memolife.biologie.ens.fr). M. Chaves and J.L. Gouzé were supported in part by the projects GeMCo (ANR 2010 BLANC020101), ColAge (Inria-Inserm Large Scale Initiative Action), RESET (Investissements dAvenir, Bioinformatique), and also by the LABEX SIGNALIFE (ANR-11-LABX-0028-01).

A Proofs

Lemma 1.

Proof In regular domains, the evolution of the system is described by continuous affine differential equations. Therefore, for a given an initial condition, the solution of system (2) in each regulatory domain is unique. Moreover, the switching segments are all transparent. Thus, when a trajectory reaches a switching segment, it can be continued into its contiguous regular domain. Therefore, a solution of the system which does not cross a switching point is unique, which ends the proof.

Lemma 2.

Proof The proof of this lemma can be straightforwardly derived from the definition of a 2-cyclic attractor, which is the union of two transition cycles sharing at least a transition in common and which does not contain vertices from which the system can escape the 2-cyclic attractor.

Lemma 3.

Proof First, two different trajectories cannot intersect in a point which is not a switching point according to Lemma 1. Then assume that they intersect on a switching point. The two trajectories would then pass through two different domains before reaching the switching point. This implies that C 2 would contain more than 1 branching vertex which is forbidden by Lemma 2.

Lemma 4.

Proof let l be the length of the entering switching segment of f k . Changing the direction of the entering axis is equivalent to transform z → l -z. Moreover:

d [f k (l -z)] dz = - df k dz (l -z) and d 2 [f k (l -z)] dz 2 = d 2 f k dz 2 (l -z)
Thus changing the direction of the entering axis changes the monotonicity but does not change the concavity of f k .

Changing the direction of the escaping axis is equivalent to transform f k (z) → l -f k (z). Moreover:

d [l -f k (z)] dz = - df k dz (z) and d 2 [l -f k (z)] dz 2 = - d 2 f k dz 2 (z)
Therefore, changing the direction of the escaping axis changes the monotonicity and concavity of f k which ends the proof.

Proposition 1.

Proof let f 1 k , f 2 k , f 3 k and f 4 k be the elementary maps of the four parallel motifs (listed in Fig. 4) :

D i(j-1) → D ij → D i(j+1) , D i(j+1) → D ij → D i(j-1) , D (i+1)j → D ij → D (i-1)j
and D (i-1)j → D ij → D (i+1)j respectively, with the entering and escaping axes oriented along axis [0,x), for f 1 k and f 2 k , and [0,y), for f 3 k and f 4 k .

The analytical expression of f 1 k and f 2 k are derived by replacing in Equation 5:

-(x(0), y(0)) = (z + θ i-1 x , θ i-1 y ) and (x(t), y(t)) = (f k (z) + θ i-1 x , θ i y ) for f 1 k ; -(x(0), y(0)) = (z + θ i-1 x , θ i y ) and (x(t), y(t)) = (f k (z) + θ i-1 x , θ i-1 y ) for f 2 k : f 1 k (z) = ( z + θ i-1 x -ϕ ij x ) • ( θ i y -ϕ ij y θ i-1 y -ϕ ij y ) d ij x d ij y + ϕ ij x -θ i-1 x with ϕ ij y > θ i y . f 2 k (z) = ( z + θ i-1 x -ϕ ij x ) • ( θ i-1 y -ϕ ij y θ i y -ϕ ij y ) d ij x d ij y + ϕ ij x -θ i-1 x with ϕ ij y < θ i-1 y . Let a = θ i-1 x -ϕ ij x , b = θ i y -ϕ ij y θ i-1 y -ϕ ij y and c = d ij x d ij y
.

Then we have:

f 1 k (z) = b c • (z + a) -a f 2 k (z) = (1/b) c • (z + a) -a Therefore, f 1 k and f 2 k are affine functions. As b > 0, f 1 k and f 2 k are increasing affine functions.
The analytical expression of f 3 k and f 4 k are derived by exchanging x and y in the expression of f 1 k and f 2 k respectively. Thus, f 3 k and f 4 k are also increasing affine functions.

Finally, according to Lemma 4, changing the orientation of either the entering or escaping axis changes the monotonicity of an elementary map, which ends the proof.

Proposition 2. Proof Let f 1 k , f 2 k , f 3 k , f 4 k , f 5 k , f 6 k , f 7 k , f 8
k be the elementary maps of the 8 perpendicular motifs (listed in Fig. 5): j+1) respectively, and S i k the switching point at the intersection of the entering and escaping segment of f i k . We assume that the origin of the entering and escaping axes of

D (i+1)j → D ij → D i(j+1) , D (i-1)j → D ij → D i(j-1) , D i(j+1) → D ij → D (i-1)j , D i(j-1) → D ij → D (i+1)j , D i(j+1) → D ij → D (i+1)j , D i(j-1) → D ij → D (i-1)j , D (i+1)j → D ij → D i(j-1) and D (i-1)j → D ij → D i(
f i k is S i k .
The analytical expression of f 1 k , f 2 k , f 3 k and f 4 k are derived by replacing in Equation 5: -(x(0), y(0)) = (θ i

x , θ i y -z) and (x(t),

y(t)) = (θ i x -f k (z), θ i y ) for f 1 k ; -(x(0), y(0)) = (θ i-1 x , z + θ i-1 y ) and (x(t), y(t)) = (f k (z) + θ i-1 x , θ i-1 y ) for f 2 k ; -(x(0), y(0)) = (θ i-1 x + z, θ i y ) and (x(t), y(t)) = (θ i-1 x , θ i y -f k (z)) for f 3 k ; -(x(0), y(0)) = (θ i x -z, θ i-1 y
) and (x(t), y(t)) = (θ i x , θ i-1 y + f k (z)) for f 4 k : The analytical expression of f 5 k , f 6 k , f 7 k , f 8 k are derived by exchanging x and y in the expression of f 1 k , f 2 k , f 3 k , f 4 k respectively. Thus, f 5 k , f 6 k , f 7 k , f 8 k are also increasing and concave.

f 1 k (z) = - ( θ i y -ϕ ij y θ i y -ϕ ij y -z ) d ij x d ij y • ( θ i x -ϕ ij x ) + θ i x -ϕ ij x with ϕ ij x < θ i x and ϕ ij y > θ i y f 2 k (z) = - ( θ i-1 y -ϕ ij y θ i-1 y -ϕ ij y +z ) d ij x d ij y • ( ϕ ij x -θ i-1 x ) + ϕ ij x -θ i-1 x with ϕ ij x > θ i-1 x and ϕ ij y < θ i-1 y f 3 k (z) = - ( θ i-1 x -ϕ ij x θ i-1 x -ϕ ij x +z ) d ij y d ij x • ( θ i y -ϕ ij y ) + θ i y -ϕ ij y with ϕ ij y < θ i y and ϕ ij x < θ i-1 x f 4 k (z) = - ( θ i x -ϕ ij x θ i x -ϕ ij x -z ) d ij y d ij x • ( ϕ ij y -θ i-1 y ) + ϕ ij y -θ i-1 y with ϕ ij y > θ i-1 y and ϕ ij x > θ i x f i k (z) for i = {1,
Finally, according to Lemma 4, changing the orientation of the entering axis or the escaping axis of an elementary map changes its monotonicity, while concavity only changes with orientation of the escaping axis, which ends the proof.

Theorem 1.

Proof Let the entering and escaping switching segments of all the motifs composing C be oriented towards the outside of C. Let f k be the elementary maps of C.

In this case, the entering and escaping axis of the parallel motifs composing C have the same orientation. Therefore, according to Proposition 1, the elementary maps of the parallel motifs of C are increasing functions.

The entering and escaping axis of the perpendicular motifs composing C both point either towards or away from the switching point at the intersection of the entering and escaping segment of these motifs. Therefore, according to Proposition 2, the elementary maps of the perpendicular motifs of C are increasing functions.

The first return map f of C is thus an increasing and continuous function as the composite of increasing and continuous functions.

Moreover, if f (x) is increasing, -f (-x) is also increasing. Therefore, changing the orientation of the axis where f is calculated does not change the monotonicity of f , which ends the proof.

Theorem 2.

Proof By taking d ij x = d ij y in the analytical expression of the elementary maps determined in the proofs of Propositions 1 and 2, it is straightforward to deduce that the elementary maps of C are linear fractional functions (i.e. homographic functions). The first return map of C is thus homographic as composite of homographic functions, and has, therefore, a constant and strict concavity.

Theorem 3.

Proof Let the entering and escaping axes of the elementary maps of C be oriented towards the outside of C. By definition of a transition cycle with no turn change, the perpendicular motifs composing C are either all clockwise or all counterclockwise.

Thus, the origins of the entering and escaping axes of the perpendicular motifs are the switching point at the intersection of the entering and escaping segments of these motifs. Then, according to Proposition 2, the elementary maps of the perpendicular motifs of C are strictly concave and increasing functions.

Moreover, the elementary maps of parallel motifs of a transition cycle are increasing affine functions (Proposition 1).

From the expression of the 2nd derivative of the composite of two functions:

d 2 (f •g) dx 2 (x) = ( df dx • g)(x) • d 2 g dx 2 (x) + ( d 2 f dx 2 • g)(x) • ( dg dx (x)
) 2 , we deduce that the composite of two increasing concave (resp. an increasing concave, and an increasing and strictly concave) functions is concave (resp. strictly concave) (and obviously increasing).

Therefore, the first return map f of C is strictly concave. Finally, it is straightforward to see that if f (x) is strictly concave, -f (-x) is strictly convex. Thus changing the direction of the axis of f changes the concavity of f , which ends the proof.

Theorem 4.

Proof f is continuous in each interval [0, α[ and ]α, l] (Theorem 1). To study if f can be continuously extended at x = α, one approach is to analyze the uniqueness of the solution of initial condition x = α until the first return in D C . D S is the unique branching domain of C 2 (Lemma 2). Thus the trajectory will not cross a switching point before it enters D S . The trajectory is therefore unique before it reaches the switching point S (Lemma 1).

Let t 1 be the time when the trajectory reaches S (Fig. 9). For t = t 1 , the PWA differential equations are not defined. In order to define the solution on S, we use the Filippov approach [START_REF] Gouzé | A class of piecewise linear differential equations arising in biological models[END_REF]) as we will see later in the proof. We now make use of the assumption on the property of SG S→T .

Part (a). Assume SG S→T is composed of 4 vertices (trajectories merge within 4 domains).

Let the regular domains composing SG S→T be labeled D 1 , D 2 , D 3 and D 4 (as in Fig. 9 left). We can rewrite the vector field f i (x, y) = (f i x (x), f i y (y)) in D i as follows:

f i (X) = Γ i (ϕ i -X) with X = (x, y), ϕ i = (ϕ i x , ϕ i y ) and Γ i = [ d i x 0 0 d i y ]
(1) Without loss of generality, we can consider the case where the separatrix is in D 1 with:

ϕ 1 x -θx > 0 and ϕ 1 y -θy < 0.
(2) To avoid stable sliding modes (excluded by Assumption 2), we have to exclude the following four cases:

(2.a) ϕ 3 x -θx < 0 from (1) (2.b) ϕ 2 y -θy > 0 from (1) (2.c) ϕ 4 y -θy > 0 and ϕ 3 y -θy < 0 (2.d) ϕ 2

x -θx > 0 and ϕ 4 x -θx < 0

(3) We also need to satisfy the assumption that trajectories merge within D 2 , D 3 or D 4 , that is: either (3.a) ϕ 3 y -θy < 0 and ϕ 4 x -θx < 0 (merging in D 2 ) or (3.b) ϕ 2

x -θx > 0 and ϕ 4 y -θy > 0 (merging in D 3 ) or (3.c) ϕ 3 y -θy < 0 and ϕ 2 x -θx > 0 (merging in D 4 ) (trajectories cannot merge in more than 1 domain due to (2)).

(4) Following Fillipov approach, the vector field on S = (θx, θy) is a vector in the convex hull of the adjacent vector fields, computed at S. That is:

f (S) ∈ co { f 1 , f 2 , f 3 , f 4 } = { a 1 Γ 1 (ϕ 1 -S) + a 2 Γ 2 (ϕ 2 -S) + a 3 Γ 3 (ϕ 3 -S) + a 4 Γ 4 (ϕ 4 -S) : a i > 0, ∑ 4 i=1 a i = 1
} Consider the case (3.a). From (1), (2.b) and (2.c), we have: ϕ i y -θy < 0 for all i. Hence the vertical component of f (S) is negative as well so a trajectory starting at S may evolve towards D 2 or D 4 . Moreover, from (2.d) we have: ϕ 2

x -θx < 0. Therefore the horizontal component on D 2 and D 4 is oriented towards D 2 . It follows that there exists only one absolutely continuous trajectory from S (evolving towards D 2 ) that satisfies the differential inclusion in S and the piecewise affine system. The case (3.b) is similar, but on the domains D 3 and D 4 (ϕ i

x -θx > 0 for all i using (1), (2.a) and (2.d), and ϕ 3 y -θy > 0 from (2.c)).

For the case (3.c), we have from (1) an (2): ϕ i y -θy < 0 and ϕ i x -θx > 0 for all i. Thus, the vertical and horizontal component of f (S) are negative and there exists only one absolutely continuous trajectory from S (evolving towards D 4 ) that satisfies the differential inclusion in S and the piecewise affine system. Therefore, if SG S→T is composed of 4 vertices, the trajectory can be extended continuously in S (t = t 1 ). As D S is the unique branching domain of C 2 , the trajectory will not cross a switching point until its first return in D C (at t = t 2 ). The solution of initial condition x = α is thus unique for t 1 < t ≤ t 2 (Lemma 1). This solution is therefore unique for 0 ≤ t ≤ t 2 . f can then be extended by continuity at x = α.

Part (b). Assume SG S→T is composed of more than 4 vertices (trajectories merge within more than 4 domains) (see Fig. 9 right).

(1) Let's again assume without loss of generality that the separatrix is in D 1 (D 1 is the unique branching domain in C 2 ), that is:

ϕ 1
x -θx > 0 and ϕ 1 y -θy < 0.

(2) To avoid sliding modes (excluded by Assumption 2), one needs:

ϕ 3 x -θx > 0 and ϕ 2 y -θy < 0 from (1)
(3) For a merging of the two trajectories within more than 4 domains, either D 2 or D 3 cannot communicate with D 4 . This means that:

(3.a) either ϕ 4

x -θx > 0 and ϕ 2 x -θx < 0 (3.b) or ϕ 4 y -θy < 0 and ϕ 3 y -θy > 0 (4) All together the information on the focal points yields:

(4.a) ϕ 1 x -θx > 0, ϕ 3 x -θx > 0, ϕ 4

x -θx > 0 and ϕ 2 x -θx < 0 (4.b) or ϕ 1 y -θy < 0, ϕ 2 y -θy < 0, ϕ 4 y -θy < 0 and ϕ 3 y -θy > 0

Given the convex hull at S, in each case there are at least two absolutely continuous solutions to the system that satisfy the differential inclusion in S and the piecewise affine differential system. In case (4.a), a solution from S evolves to D 2 ; another solution from S to D 3 or D 4 (depending on the vertical coordinates). In case (4.b), a solution from S evolves to D 3 ; another solution from S to D 2 or D 4 (depending on the horizontal coordinates).

Therefore, if SG S→T is composed of more than 4 vertices, there are at least two distinct trajectories emerging from S. According to Lemma 3, two different trajectories cannot intersect. Each will thus cross D C in two distinct points after a first return. f then cannot be extended by continuity at x = α, which ends the proof.

Theorem 5.

Proof First, we deduce from Theorems 2 and 3 that f has a constant and strict concavity, and from Theorem 1 that f is continuous.

Let g(x) = f (x) -x. Thus g is continuous, and dg dx is strictly monotone and has at most one zero. Therefore, by applying the intermediate value theorem on each interval where g is strictly monotone (one interval if dg dx has no zeros, two if dg dx has one zero), we deduce that g has at most 2 zeros i.e. f has at most two fixed points. Then, there are at most two limit cycles in D C . Moreover, if f admits two fixed points at x = x 1 and x = x 2 , dg dx (x 1 ) • dg dx (x 2 ) < 0. In addition, df dx > 0 (Theorem 1). Thus one of the two limit cycles is stable and the other one is unstable [START_REF] Strogatz | Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, And engineering[END_REF], p.281), which ends the proof.

Lemma 5.

Proof Let all the axes of the elementary maps f k composing f be oriented towards the outside of C. We have:

f = fn • f n-1 • . . . f 1 . Let I k be the interval of definition of f k and l k the length of the entering switching segment where f k is calculated. Let g k = f k • . . . • f 1 for k ∈ {1, 2, . . . , n}.
(1) To prove statement (a), assume that C has no inside branching transition. This means that each I j is of form:

I j = [0, α j ] (for some α j ≥ 0) or I j = ∅ (2) Assume that I ̸ = ∅, that is, ∃a ≥ 0 : a ∈ I. This statement is equivalent to: ∀j ∈ {1, 2, . . . n -1} , g j (a) ∈ I j+1
We also have that I j ̸ = ∅ for all j (otherwise I = ∅).

(3) To get a contradiction, assume that 0 / ∈ I. Then:

∃p: ∀i ∈ {1, . . . , p -2} , g i (0) ∈ I i+1 and g p-1 (0) / ∈ Ip
Now since Ip = [0, αp] (for some αp ≥ 0), we have that:

g p-1 (a) αp from (2), and g p-1 (0) > αp from (3).

Since g p-1 is increasing (because all axes are oriented outside of C, see Propositions 1 and 2), we have:

αp < g p-1 (0) g p-1 (a) αp
which is a contradiction (αp < αp). Therefore 0 ∈ I which ends the proof of statement (a). Statement (b) can be straightforwardly proved using the same reasoning.

Theorem 6.

Proof f is continuous according to Theorem 1. Assume that the axis where f is calculated is oriented to the outside of C, and that C is a transition cycle with no turn change. Then f is strictly concave (Theorem 3). Assume moreover that C contains no inside branching transition. Then we deduce from Lemma 5a that I contains 0 and f (0) 0.

Let g(x) = f (x) -x. g is continuous, strictly concave and g(0) 0. Therefore, by applying the intermediate value theorem on each interval where g is strictly monotone, we show that g has either no zero, a single zero, or two zeros one of which being x = 0. Therefore D C contains either no limit cycle or a single limit cycle.

Moreover if f admits a non-zero fixed point for x = x 0 , dg dx (x 0 ) < 0. Therefore, it corresponds to a stable limit cycle, which ends the proof. Note that the case f (0) = 0 corresponds to an equilibrium point of the system but also belongs to a switching domain. By Assumption 1, it cannot be a focal point, so it is a Filippov-type equilibrium. This case can only happen when C consists of 4 consecutive perpendicular motifs with the same orientation (see transition cycle C 1 in Fig. 10 for an example).

Theorem 7.

Proof Assume that the axis where f is calculated is oriented to the outside of C, and that C contains no outside branching transition. Then we deduce from Lemma 5b that I contains l and f (l) ≤ l. From Assumption 1, we further deduce that f (l) < l.

Assume the system admits an unstable limit cycle in D C , that is ∃x 0 : f (x 0 ) = x 0 and | df dx (x 0 )| > 1 (Strogatz (2001), p.281). Since f is an increasing function, df dx (x 0 ) > 1. Therefore, ∃x 1 > x 0 : f (x 1 ) > x 1 .

Assume either C has no turn change or the decay rates are equal. Then, f has a constant and strict concavity (Theorems 2 and 3). Moreover, f is a continuous and increasing function (Theorem 1).

Let g(x) = f (x) -x. We have g(l) < 0, g(x 1 ) > 0 and g is continuous, and has a constant and strict concavity. Therefore, by applying the intermediate value theorem, we show that f admits a single fixed point x = x 2 . Moreover, dg dx (x 2 ) < 0. Thus, x = x 2 corresponds to a stable limit cycle [START_REF] Strogatz | Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, And engineering[END_REF], p.281) and the system admits a single stable limit cycle in D C , which ends the proof.

Theorem 8.

Proof Let f be the 2-cycle first return map of C 2 calculated on the oriented segment D C , f 1 and f 2 the first return maps of C 1 and C 2 respectively, and I 1 = [0, a[ and I 2 =]a, l] the intervals of definition of f 1 and f 2 .

Let g(x) = f (x) -x. Assume either C 1 and C 2 have no turn change, or the decay rates are equal. Thus, according to Theorems 2 and 3, g has a constant and strict concavity in each interval I 1 and I 2 . Assume also that SG S→T is composed of 4 vertices. Thus f and therefore g is continuous (Theorem 4). Assume moreover the system admits two stable limit cycles in D C 2 , that is there exists x 1 < a and x 2 > a such that g(x 1 ) = g(x 2 ) = 0.

By applying the intermediate value theorem, we show that g admits a unique fixed point x 3 ∈ ]x 1 , x 2 [ which corresponds to an unstable limit cycle separating the basins of attraction of the two stable limit cycles, which ends the proof. Fig. 14 First return maps for the continuous case. Positive fixed points of the maps (e.g. positive intersections between the maps f (in blue) and line of equation f (z) = z (in red)) correspond to limit cycles: (a) one large stable limit cycle; (b) one small stable limit cycle; (c) one large stable limit cycle, one unstable limit cycle and one stable equilibrium point at z = 0; (d) one large and one small stable limit cycles and one unstable limit cycle; (e) one stable equilibrium point at z = 0. The parameter values for each case are indicated in the caption of Fig. Fig. 15 First return maps for the discontinuous case. Positive fixed points of the maps (e.g. positive intersections between the maps f (in blue) and line of equation f (z) = z (in red)) correspond to limit cycles: (a) one large stable limit cycle; (b) one small stable limit cycle; (c) one large stable limit cycle, one unstable limit cycle and one stable equilibrium point at z = 0; (d) one small and one large stable limit cycle, and one unstable limit cycle; (e) one small and one large stable limit cycles; (f) one stable equilibrium point; (g) one large stable limit cycle and one stable equilibrium point at z = 0. For cases (a) to (f), the parameter values are indicated in the caption of Fig. 13. For (g), the parameter values are the same as the other cases except: θ x = 2, k 4x = 0

B First return maps for the applications considered

C Translation to the continuous differential framework

The class of PWA biological model corresponding to the discontinuous case (section 5.2) is translated into a continuous ordinary differential model by replacing step functions s + (x, θ) and s -(x, θ) by Hill function of the form x n

x n +θ n and θ n x n +θ n , respectively. The continuous differential system obtained is:

                   dx dt = [ k 1x + k 2x • y n y n + (θ 1 y ) n ] • (θ 2 y ) n y n + (θ 2 y ) n • x n x n + θ n x + y n y n + (θ 3 y ) n • [ k 3x + k 4x • x n x n + θ n x ] -dx • x dy dt = k 1y • θ n x x n + θ n x -dy • y (8)
Simulations of this model for parameter setting corresponding to the situation where the discontinuous case admits a large stable limit cycle and an equilibrium point (case (g), Fig. 15) and the computation of a first return map are presented in Fig. 16. Interestingly, the simulations show an abrupt change of the derivative of the first return map (at x = 91.9, Fig. 16b) which corresponds to the discontinuity observed in the corresponding PWA model (Fig. 15).

Fig. 1

 1 Fig. 1 (a) Example of step function modeling a switch-like activation process. The step function is defined as: s + (x, θ) = 0 if x < θ and s + (x, θ) = 1 if x > θ, with θ = 1. (b) Regulatory network of the class of biological model represented by Equations 1. Each arrow represents a step function involved in the system. Normal arrows correspond to step functions related to activation processes, blunt arrows to step functions related to inhibition processes

Fig. 2

 2 Fig.2Temporal simulations of x and y levels whose evolution is described by Equations 1: (a) monorhythmic case with one large stable limit cycle; (b) birhythmic case with one large and one small stable limit cycles. A pulse of y is applied at t = 5 a.u. This pulse induces a shift from a small amplitude to a large amplitude oscillatory regime. Simulation of the level of x (resp. y) is shown in blue (resp. red). The initial conditions are x = 5.1,y = 54.3 for (a) and x = 10.6,y = 53.6 for (b). The parameter values are: θ 1 y = 50, θ 2 y = 53, k 2y = 15, d x = 1, d y = 1, k 1y = 58, k 1x = 20, θ x = 10 and: k 3y = 10 for (a); k 3y = 5 for (b).

Fig. 3

 3 Fig. 3 Transition configurations arising from one vertex according to its focal point position. The four configurations of the bottom are branching vertices

Fig. 5

 5 Fig. 5 Perpendicular motifs. Top: clockwise perpendicular motifs. Bottom: counterclockwise perpendicular motifs.

  branching domain, two cases can be distinguished according to the position of the separatrix curve which emerges in D k : (1) the separatrix curve does not intersect D k-1 s . In this case, either all the points of D k-1 s map D k s and I k = [0, l k ], or none of the points of D k-1 s map D k s and I k = ∅ (Fig. 8, left). (2) The separatrix curve intersects D k , splitting the entering switching segment of D k in two segments. One of these segments will map D k s and the other not and I k ⊂ [0, l k ] (Fig. 8, right).

Fig. 7 Fig. 8

 78 Fig.7An example of oriented perpendicular motif (left) and parallel motif (right) with entering axis [0, z 1 ) and escaping axis [0, z 2 ). The red arrows are the transitions compos ing each motif. The elementary map associated to the perpendicular (resp. parallel) motif is increasing and strictly concave (resp. affine). A trajectory of the system from the entering to the escaping axis is indicated in blue arrow

Theorem 1 (

 1 Monotonicity)The first return map of a transition cycle C is an increasing and continuous function.Theorem 2 (Concavity for equal decay rates) Assume the decay rates are equal in each domainD ij , i.e. d ij x = d ij y for (i, j) ∈ {1, 2, . . . , n x + 1} × {1, 2, . . . , n y + 1}).Then the first return map of C has a constant and strict concavity.Theorem 3 (Concavity for transition cycles with no turn change)Assume C is a transition cycle with no turn change. If the axis of definition of the first return map f of C is oriented towards the outside of C, f is strictly concave. Otherwise, f is strictly convex.

Theorem 4 (Fig. 9

 49 Fig. 9 Two examples of transition graph containing a 2-cyclic attractor. (left) SG S→T is a 4-node subgraph. (right) SG S→T contains more than 4 nodes (here 6 nodes). The red arrows are the transitions composing the subgraph SG S→T . The black arrows are the transitions connecting SG S→T to the rest of the transition graph (not shown). The trajectory lying in the separatrix reaches the switching point S and either does not split (left) or split into two distinct trajectories (right). In both cases, D 1 corresponds to a branching vertex (ϕ 1x -θ x > 0 and ϕ 1 y -θ y < 0). In the left case, D 2 and D 3 communicate with D 4 (ϕ 2x -θ x > 0, ϕ 3 y -θ y < 0) leading to a 4-node SG S→T , whereas in the right case D 4 does not communicate with D 2 (ϕ 2x -θ x < 0 and ϕ 4x -θ x > 0) leading to a SG S→T composed of more than 4 nodes. According to Theorem 4, the left case gives a continuous first return map, which is discontinuous in the right case.

  (a) Assume C contains no inside branching transition. Then if I ̸ = ∅, 0 ∈ I. (b) Assume C contains no outside branching transition. Then if I ̸ = ∅, l ∈ I.

Fig. 10 .

 10 This attractor is composed of two embedded transition cycles: one 4-elements cycle D 22 → D 12 → D 13 → D 23 → D 22 (transition cycle C 1 ) and one 6-elements cycle D 11 → D 12 → D 13 → D 23 → D 22 → D 21 → D 11 (transition cycle C 2 ). The branching vertex is D 22 and the vertex where both cycles merge is D 12 .

Fig. 10

 10 Fig. 10 Transition graph for the continuous case. The transition graph is composed of two transition cycles which form a 2-cyclic attractor: cycle D 22 → D 12 → D 13 → D 23 → D 22 (transition cycle C 1 ) and cycle D 11 → D 12 → D 13 → D 23 → D 22 → D 21 → D 11 (transition cycle C 2 ). C 1 contains one outside branching transition (D 22 → D 21 ) whereas C 2 contains one inside branching transition (D 22 → D 12 ). C 1 and C 2 have no turn change

  Fig. 11Numerical simulations for the continuous case corresponding to the possible dynamical configurations derived from the structural principles: case with one large stable limit cycle (top left); one small stable limit cycle (top right); one large stable limit cycle and one unstable limit cycle (bottom left); one large and one small stable limit cycles and one unstable limit cycle (bottom right). The stable (resp. unstable) limit cycles are drawn in blue (resp. in red). The parameter values for the top left and bottom right figures are indicated in the legend of Fig. 2. The parameter values for the others figures are: θ 1 y = 50, θ 2 y = 53, k 2y = 15, d x = 1, d y = 1 and: k 1y = 58, k 1x = 20, k 3y = 5, θ x = 20 (top right); k 1y = 54, k 1x = 30, k 3y = 0, θ x = 10 (bottom left). The system admits one stable equilibrium point at (θ x , θ 2 y ) and no limit cycle for k 1y = 58, k 1x = 20, k 3y = 0, θ x = 10 (simulation not shown).

Fig. 12

 12 Fig. 12 Discontinuous case: graph of interactions (left) and transition graph (right). The transition graph is composed of two transition cycles which form a 2-cyclic attractor: cycle D 23 → D 13 → D 14 → D 24 → D 23 (transition cycle C 1 ) and cycle D 11 → D 12 → D 13 → D 14 → D 24 → D 23 → D 22 → D 21 → D 11 (transition cycle C 2 ). C 1 contains one outside branching transition (D 23 → D 22 ) whereas C 2 contains one inside branching transition (D 23 → D 13 )

  Fig. 12. This graph reproduces the transition graph of a reduced version of the p53-Mdm2 network model (Abou-Jaoudé et al (2009)) analyzed in Abou-Jaoudé et al (2011). The transition graph is composed of two embedded transition cycles: one 4element cycle D 23 → D 13 → D 14 → D 24 → D 23 (transition cycle C 1 ) and one 8-element cycle D 11 → D 12 → D 13 → D 14 → D 24 → D 23 → D 22 → D 21 → D 11 (transition cycle C 2 ). The branching vertex is D 23 and the vertex where both cycles merge is D 13 .

  Fig. 13Numerical simulations for the discontinuous case corresponding to the possible dynamical configurations derived from the structural principles: case with one large stable limit cycle (top left); one small stable limit cycle (top right); one large stable limit cycle and one unstable limit cycle (middle left); one small and one large stable limit cycle and one unstable limit cycle (middle right); one small and one large stable limit cycles (bottom left). Stable (resp. unstable) limit cycles are drawn in blue (resp. in red). The parameter values are: θ 1 y = 10.1, θ 2 y = 10.2, θ 3 y = 10.5, k 1x = 29, k 2x = 2, k 3x = 35, k 1y = 12.5, d x = 1, d y = 1 and:θ x = 5, k 4x = 4 (top left), θ x = 50, k 4x = 4 (top right), θ x = 8, k 4x = 0 (middle left), θ x = 30, k 4x = 4 (middle right), θ x = 20, k 4x = 4 (bottom left). The system admits no limit cycle (one stable equilibrium point at (θ x , θ 3 y )) for θ x = 30, k 4x = 0 (simulation not shown)

  b, c) > 0. The first and second derivative of this function are respectively positive and negative. Therefore, f 1 k , f 2 k , f 3 k and f 4 k are increasing and strictly concave.

  Fig.14First return maps for the continuous case. Positive fixed points of the maps (e.g. positive intersections between the maps f (in blue) and line of equation f (z) = z (in red)) correspond to limit cycles: (a) one large stable limit cycle; (b) one small stable limit cycle; (c) one large stable limit cycle, one unstable limit cycle and one stable equilibrium point at z = 0; (d) one large and one small stable limit cycles and one unstable limit cycle; (e) one stable equilibrium point at z = 0. The parameter values for each case are indicated in the caption of Fig.11

Fig. 16

 16 Fig. 16 Simulations of the continuous differential equations translated from the discontinuous example studied in section 5.2. (a) Numerical simulations of the continuous model in the phase plane. The initial conditions of the trajectories (shown in blue) which converge to a stable point and a limit cycle are (25.7, 90) and (25.7, 100) respectively. The initial conditions of the trajectories simulated in reverse-time (shown in red) are: (26.5, 64.475), (26.5, 64.465), (27.1, 62.4), (27.1, 62.41). The trajectories simulated in reverse-time give an approximation of the separatrix curve which delimits the basins of the two attractors. Two additional unstable points appear: one at (27.35, 55) and the other at (26.8, 63.5). (b) Numerical simulation of the first return map from and to the half line of equation: {x = 25.7, y ≥ 0} (black dashed line in Figure (a)). The fixed point x 1 = 97.8 of the first return map corresponds to the stable limit cycle towards which the trajectory starting from (25.7, 100) converges (Figure (a)). The parameter values are: n = 10 , θ x = 25, θ 1 y = 50, θ 2 y = 70, θ 3 y = 90, k 1x = 5, k 2x = 50, k 3x = 55, k 4x = 30, k 1y = 190, d x = 1 and d y = 1