, a necessary and sufficient condition for the existence of a limit cycle of a specific 2-dimensional HGRN has been discussed. There is no general analysis method of Ndimensional HGRNs. So the study of such analysis methods is of great interest, which is the main research direction of this thesis.

Among different kinds of dynamical properties to describe gene regulatory networks, we mostly focus in this thesis on long-term dynamical properties, in other words, the behaviors of trajectories when time approaches infinity. For instance, a trajectory can converge to a limit cycle when time approaches infinity. Such limit cycles represent stable oscillations which can be found in many gene regulatory networks, such as the circadian clock. In discrete models, long-term dynamical properties correspond to discrete attrac-

Organization of the Manuscript

Chapter 2 introduces different modeling frameworks that have been used to model gene regulatory networks, including HGRNs. The notations used in this chapter are also used in the rest of this thesis.

Chapters 3 and 4 introduce our general analysis method of HGRNs: In Chapter 3 we present our limit cycle analysis method and in Chapter 4 we present our reachability analysis method. Some new concepts of HGRNs are introduced in Chapter 3, and are also used in the work of Chapter 4.

In Chapters 5, 6 and 7, we focus on some specific networks. Chapter 5 introduces our analysis of a HGRN of 3-dimensional repressilator, where we find conditions for the existence of sustained oscillations in this HGRN. This work relies on some methods of Chapter 3. In Chapter 6, we present a parameter identification method of HGRNs using the HGRN of repressilator in Chapter 5 as a case study. The work of Chapter 7 is a continuation of the work of Chapter 5: In Chapter 7, we present conditions for the existence of a discrete periodic attractor in discrete models of 4-dimensional repressilators.

Finally, Chapter 8 concludes the whole thesis.

Notations

d s or d s i : A discrete state, introduced in Section 2.2 c s or c s i : The celerity vector of a discrete state, introduced in Section 2.3.2 h or h i : A hybrid state, introduced in Section 2.3.2 e or e i : A boundary of HGRN, introduced in Section 2.3.2 τ or τ i : A hybrid trajectory of HGRN, introduced in Section 2.3.2 D or D i : A discrete domain, introduced in Section 3.2.1

Chapter 1

INTRODUCTION

Context and Motivations

Gene Regulatory Networks

Gene expression has two major steps: Transcription, which generates mRNA from genes, and translation, which generates proteins from mRNA. In general, each mRNA molecule produces a specific protein (or set of proteins). This protein can be structural (giving it a particular structural property) or be an enzyme (catalyzing a certain reaction). Additionally, it can also serve only to activate or inhibit other genes. These regulations (activation or inhibition) between genes constitute gene regulatory networks, which play a key role in various cellular processes and pathways.

One major problematic of gene regulatory networks is analyzing their dynamical properties. In the literature, different kinds of dynamical properties have been used to describe the behaviors of gene regulatory networks, such as the existence of stable states or oscillations, the reachability between state(s) (for example, whether the trajectory/trajectories from a state can reach certain state/states), or some properties described by temporal logic (for example, the sequence of states on a trajectory always satisfies certain property, or will eventually satisfy certain property, or will satisfy certain property until another property is satisfied). Knowing the possible dynamical properties of gene regulatory networks can help us understand better the underlying nature of these biological systems, and can further guide us to develop control methods to change the behaviors of certain systems for medical purposes. Actually, the complexity of gene regulatory networks can increase dramatically when the number of genes increases, which makes it hard to analyze these properties.

Modeling of Gene Regulatory Networks

In order to analyze these dynamical properties, one solution is to model mathematically the gene regulatory networks. Such approaches can be more efficient and cheaper compared to biological experiments, and can also allow us to understand the dynamics of the networks at system level. First, a gene regulatory network is abstracted by a directed graph called influence graph, where the vertices represent the genes in the system and the arcs represent the regulations between genes. In fact, the initial understanding of gene regulatory networks normally only contains the relations between genes, so it is natural that they are abstracted by directed graphs as a first step. As stated before, a gene can generate different products (mRNA and protein), while in an influence graph we usually abstract all products of one gene by one entity (one vertex). The structure of the influence graphs can give some insights about the dynamical properties of gene regulatory networks, but to determine the dynamical properties, influence graphs are not enough, dynamical models are also required.

In the literature, different modeling frameworks have been applied to model gene regulatory networks, mainly continuous models [START_REF] Almeida | Control of synchronization ratios in clock/cell cycle coupling by growth factors and glucocorticoids[END_REF][START_REF] Barik | A model of yeast cell-cycle regulation based on multisite phosphorylation[END_REF][3], discrete models [START_REF] Kauffman | Metabolic stability and epigenesis in randomly constructed genetic nets[END_REF][START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF][START_REF] Thomas | Regulatory networks seen as asynchronous automata: a logical description[END_REF][START_REF] Khalis | The SMBioNet method for discovering models of gene regulatory networks[END_REF][START_REF] Bernot | Gene regulatory networks with multiplexes in European simulation and modelling conference proceedings[END_REF][START_REF] Boyenval | What is a cell cycle checkpoint? the TotemBioNet answer in Computational Methods in Systems Biology[END_REF] and hybrid models [START_REF] Behaegel | A hybrid model of cell cycle in mammals[END_REF][START_REF] Cornillon | Hybrid gene networks: a new framework and a software environment[END_REF][START_REF] Comet | A formal model for gene regulatory networks with time delays in International Conference on Computational Systems[END_REF]. In continuous models, a state is a real vector representing the continuous expression of each gene. In discrete models, the continuous state space is divided into discrete regions and a state, also called a discrete state, is a vector of integers representing a discrete region. Hybrid models have both continuous and discrete components. In hybrid models, like in discrete models, the continuous state space is divided into discrete regions, but contrary to discrete models where we do not consider the dynamics inside these discrete regions, the dynamics of hybrid models inside each discrete region is described by a simplified continuous model. More details about these modeling frameworks are introduced in Chapter 2.

One major difference between discrete models and models having continuous behaviors (continuous or hybrid models) is that: Discrete models can be applied on systems of large scale, but, for continuous or hybrid models, it is generally difficult to study dynamical properties of models in high dimension, mainly due to their dynamical complexities.

In fact, the dynamics of a discrete model can be described by all the possible transitions between discrete states, which is a directed graph with a finite number of nodes, so theoretically we can determine all dynamical properties, the only constraint is that, when the dimension (the number of genes) increases, the computational cost can be too high. Moreover, there are less parameters in discrete models (there are only parameters that describe the qualitative properties of the system), so it requires less data or biological knowledge to construct the models.

For continuous models like ordinary differential equations, there are many open questions related to dynamical analysis, for instance, there is no general method to prove the existence of a stable limit cycle for models that have more than three dimensions. Some hybrid models can be considered as simplifications of ordinary differential equations, and their dynamical properties are easier to analyze, but still, the dimension of the systems that can be analyzed by these hybrid models is limited. Additionally, parameter identification of continuous or hybrid models requires more data.

Scope of the Study and Motivations

A recent proposed hybrid modeling framework: Hybrid gene regulatory network (HGRN) [START_REF] Cornillon | Hybrid gene networks: a new framework and a software environment[END_REF], which is an extension of Thomas' discrete modeling framework [START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF][START_REF] Thomas | Regulatory networks seen as asynchronous automata: a logical description[END_REF], has lower dynamical complexity compared to other hybrid models. In HGRNs, the time derivative is a constant vector in each discrete region. So it is easier to analyze and could be potentially applied to reveal continuous dynamical properties of larger (compared to networks that are modeled by continuous models or other hybrid models) gene regulatory networks. Moreover, HGRN has a special dynamical behavior: The existence of sliding modes, which means that if a trajectory reaches a boundary of a discrete region but can not cross this boundary, then it slides along this boundary (such boundaries are also called "black walls" in the literature). The physical meaning of sliding modes is that some genes reach temporal stable stables, so their values will not change for some time, while the values of other genes keep changing. tors. One limit is that, based on discrete models, it is hard to derive continuous behaviors inside discrete attractors; for example, we can not distinguish sustained oscillations and damped oscillations, while such continuous behaviors can be derived from HGRNs. The study of such long-term behaviors of gene regulatory networks is of great concern. Indeed, the stable limit cycles correspond to robust oscillations (oscillations will always recover after small perturbations) which exist in some genetic systems, for example circadian clock, and the analysis of limit cycles can help us verify whether the networks are correctly constructed.

Contributions

The major contributions of this thesis lie in both general analysis methods and the analysis of some specific networks. The contributions for the general analysis methods include a limit cycle analysis method presented at CMSB 2022 as regular paper [START_REF] Sun | Limit cycle analysis of a class of hybrid gene regulatory networks[END_REF] and a reachability analysis method presented at RP 2023 as regular paper [START_REF] Sun | Reachability Analysis of a Class of Hybrid Gene Regulatory Networks in Reachability Problems[END_REF].

General Analysis Methods

Limit Cycle Analysis

Our limit cycle analysis method of HGRNs can automatically identify limit cycles of HGRNs in any dimension and analyze their stability. This method can find all limit cycles that do not cross the same discrete region for more than once in one period. It is indeed important to analyze limit cycles, because many gene regulatory networks have stable oscillations which correspond to stable limit cycles. In this work, we also propose new concepts associated to HGRNs (for instance, in order to facilitate the analysis of HGRNs, we propose the new concept "discrete domain", based on which the dynamics of the original system can be abstracted by transitions on a directed graph), which can be useful for the future studies of HGRNs.

In the literature, even though few works exist about limit cycles analysis in HGRNs (only a 2-dimensional limit cycle was analyzed in [START_REF] Cornillon | Hybrid gene networks: a new framework and a software environment[END_REF]), limit cycles were studied in other hybrid models of gene regulatory networks. Most of these works are based on the Poincaré map, which is also the main idea of our approach. In [START_REF] Belgacem | Control of negative feedback loops in genetic networks in[END_REF][START_REF] Firippi | Topology-induced dynamics in a network of synthetic oscillators with piecewise affine approximation[END_REF], the Poincaré map is used to study the limit cycle of simple piecewise affine systems in two dimensions. In these works, since the system is planar, it is easy to compute and analyze the Poincaré map. In [START_REF] Mestl | Chaos in high-dimensional neural and gene networks[END_REF][START_REF] Edwards | Analysis of continuous-time switching networks[END_REF][START_REF] Edwards | A calculus for relating the dynamics and structure of complex biological networks[END_REF], methods are proposed to find and analyze limit cycles in higher dimensions of piecewise affine system with a uniform decay rate. The hypothesis of a uniform decay rate in these works makes it always possible to calculate a Poincaré map because they have a simple shape. However, for a general piecewise affine system, it is difficult to prove theoretically the existence of limit cycle except for some particular examples such as negative loops [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] Chaves | Hierarchy of models: From qualitative to quantitative analysis of circadian rhythms in cyanobacteria[END_REF]. Compared to previous works about limit cycles in hybrid systems, our work has two major novelties: (1) We consider limit cycles with sliding modes, and [START_REF] Barik | A model of yeast cell-cycle regulation based on multisite phosphorylation[END_REF] We use an abstraction method in order to find automatically cycles of discrete regions, which might contain limit cycles (in other words, cycles of discrete regions which contain at least one continuous trajectory).

Reachability Analysis

The second proposed general analysis method is a reachability analysis algorithm of HGRNs, which can verify automatically if the trajectory from a singular state can eventually reach certain region (a set of states). This algorithm always stops in finite time, and in most cases, it gives the correct answer of this reachability problem. In fact, only in some particular cases, for example in the presence of chaos, the algorithm remains inconclusive. Such cases exist for general HGRNs, but they are not identified in any preexisting HGRN of gene regulatory network. In fact, we have only identified one hybrid system which is equivalent to a HGRN and which has chaotic attractor, and this is a hybrid system of an electric circuit.

Actually, this work is related to the decidability problem, which is whether we can find an algorithm that always stops in finite time and answers this reachability problem. The decidability problem has not been investigated for HGRNs, but it has been studied among hybrid systems that are close to HGRNs. The hybrid systems that are the closest to HGRNs are piecewise-constant derivative systems (PCD systems), which are decidable in 2 dimensions [START_REF] Maler | Reachability analysis of planar multi-linear systems in Computer Aided Verification[END_REF] and become undecidable in 3 dimensions [START_REF] Asarin | Reachability analysis of dynamical systems having piecewise-constant derivatives[END_REF]. There exists also a class of decidable hybrid systems called initialized rectangular automata which does not include HGRNs. So whether HGRNs are decidable or not is still an open question. In our work, we formalize a set of trajectories for which this reachability problem is decidable. Meanwhile, for other trajectories, it is unknown.

Analysis of Specific Networks

Besides general analysis methods, we also analyze some specific networks of interest, called the repressilators, which led to two publications: The first one is about the 3dimensional repressilator and has been presented at BIOINFORMATICS 2023 as regular paper [START_REF] Sun | Condition for Sustained Oscillations in Repressilator Based on a Hybrid Modeling of Gene Regulatory Networks[END_REF]; the second one is about 4-dimensional repressilators and has been presented at CMSB 2023 as regular paper [START_REF] Sun | Condition for Periodic Attractor in 4-Dimensional Repressilators in Computational Methods in Systems Biology[END_REF].

3-dimensional Repressilator

The 3-dimensional canonical repressilator is a widely studied network in the field of synthetic biology. In this work, we find a sufficient and necessary condition for the existence of sustained oscillations in a HGRN of 3-dimensional canonical repressilator, and we propose a method to compute simple constraints on parameters that satisfy this condition. This result could potentially contribute to the biological construction of synthetic repressilators. In the literature, most studies of the canonical repressilator use ordinary differential equations. This work is the first one based on HGRNs.

4-dimensional Repressilators

In our study of 4-dimensional repressilators in the form of discrete models, we provide a sufficient and necessary condition, which is described by the topological features of the influence graphs, for the existence of a discrete periodic attractor in 4-dimensional repressilators. This condition corresponds to a sufficient condition for the existence of oscillations in associated HGRNs. Other theoretical works investigate the relation between the topology of influence graphs and the dynamical properties. The major novelty of this work is that we discover a new topological feature.

Parameter Identification

Finally we also investigate parameter identification problems, which is the first step toward the understanding of systems. Indeed our analysis methods can be applied only after that the parameters are identified. There are different kinds of parameter identification problems that are related to gene regulatory networks, for example, there are works aiming to identify the regulations between genes, and there are other works that want to identify the parameters of dynamical models assuming that the influence graph is known. Our work belongs to the second case and is based on time series data.

Collaborations

The analysis work of the 3-dimensional repressilator and a part of the analysis work of the 4-dimensional repressilators come from a collaboration with the team of Jean-Paul Comet at the I3S laboratory (Laboratoire d'Informatique, Signaux et Systèmes de Sophia Antipolis) in Nice. For this collaboration, I stayed 3 weeks at I3S in June 2022. T = (D i , D i+1 , ..., D j-1 , D j ): A sequence of discrete domains, introduced in Section 3.2.2 M T : The transition matrix of the sequence of discrete domains T , introduced in Section 3.2.2 S: A compatible zone, introduced in Section 3.2.2 N d (h, r): The neighborhood in the same discrete state of the hybrid state h with the radius r, introduced in Section 3.3 N D (h, r): The neighborhood in the same discrete domain of the hybrid state h with the radius r, introduced in Section 3.

STATE OF THE ART OF GENE REGULATORY NETWORKS MODELING

A gene regulatory network is abstracted mathematically by a directed graph IG = (V, A) called influence graph, where V is the set of vertices representing the genes in the system and A is the set of arcs representing the regulations between genes. To understand the dynamical properties of gene regulatory networks, influence graphs are not enough, we also need dynamical models. Different classes of dynamical models have been used in the literature to model gene regulatory networks. The most used ones are ordinary differential equations (ODE) and discrete models, which are presented briefly in the first two sections of this chapter. Another modeling approach is hybrid modeling which is not yet widely used for now but is attracting more attention in recent years. In the third section of this chapter, the piecewise affine systems are presented at first, which are the most used hybrid systems to model gene regulatory networks; then the hybrid systems studied in this thesis, Hybrid gene regulatory networks, are presented.

Ordinary Differential Equations (ODE)

Ordinary differential equations (ODE) are mostly used in the literature to model gene regulatory networks [START_REF] Almeida | Control of synchronization ratios in clock/cell cycle coupling by growth factors and glucocorticoids[END_REF][START_REF] Barik | A model of yeast cell-cycle regulation based on multisite phosphorylation[END_REF][3][START_REF] Chen | Biocomputing'99 29[END_REF][START_REF] Goldbeter | Computational approaches to cellular rhythms[END_REF][START_REF] Tyson | Regulation of the eukaryotic cell cycle: molecular antagonism, hysteresis, and irreversible transitions[END_REF]. The use of ODEs is based on the assumption that the numbers of molecules describing the expressions of genes are very big so that real numbers can be used to represent the expressions of genes. A widely used form to describe the dynamics of one gene by ODEs, for example G 0 , is given in Eq (2.1).

dx 0 dt = k 0 + f (x 1 , x 2 , ..., x m ) -γx 0 (2.1)
where x i is a real number describing the expression value of gene G i ; G 1 , G 2 , ..., G m are the genes that influence gene G 0 ; k 0 is the production rate of gene G 0 when it is not regulated by other genes; γ is the degradation rate which describes, for example, the efficacy of protease if the expression of the gene is measured by the amount of generated protein; function f describes the regulations of G 1 , G 2 , ..., G m on G 0 , and is usually based on Hill function, for example f can be of the following form.

f (x 1 , x 2 , ..., x m ) = f 1 (x 1 ) + f 2 (x 2 ) + ... + f m (x m ) (2.2)
where

f i (x i ) = k i x n i θ n i +x n i if G i activates G 0 and f i (x i ) = k i θ n i θ n i +x n i if G i inhibits G 0 .
There are two major limits of ODEs: 1) Their dynamics are hard to analyze; 2) There are many parameters to be identified. For now it is still difficult (or impossible) to analyze large networks based on ODEs.

Discrete Modeling

Discrete models are also widely used to model gene regulatory networks, which in fact attract more and more attention in the last 20 years. In discrete models [START_REF] Kauffman | Metabolic stability and epigenesis in randomly constructed genetic nets[END_REF][START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF][START_REF] Thomas | Regulatory networks seen as asynchronous automata: a logical description[END_REF][START_REF] Khalis | The SMBioNet method for discovering models of gene regulatory networks[END_REF][START_REF] Bernot | Gene regulatory networks with multiplexes in European simulation and modelling conference proceedings[END_REF][START_REF] Boyenval | What is a cell cycle checkpoint? the TotemBioNet answer in Computational Methods in Systems Biology[END_REF], the continuous expression of a gene is abstracted by an integer (e.g. 0, 1, 2, ...), called discrete level, which describes the discrete expression level of a gene. We note a i as the discrete level of gene G i . A discrete model is a logic program which is a set of logic rules. The form of a logic rule is shown as follows:

a i = k ← ϕ i (2.3)
where k is a possible discrete level of G i , and ϕ i is a logic formula. The form of a logic formula ϕ is given as follows:

ϕ :== ⊤ | ⊥ | a i ∼ k | ϕ 1 ∧ ϕ 2 | ϕ 1 ∨ ϕ 2 (2.4)
where k is a possible discrete level of gene G i , ∼ is one of the relations {>, <, =, ≥, ≤}, ϕ 1 and ϕ 2 are also logic formulas.

A logic rule (see Eq (2.3)) indicates that if at discrete time t (t is an integer) ϕ i is satisfied, then at time t + 1 the discrete level of G i can be updated to k. For example, a

a 0 = 1 ← ⊤ a 1 = 1 ← (a 0 = 1)
a 1 = 0 ← (a 0 = 0) a 2 = 1 ← (a 0 = 0) ∧ (a 1 = 1)
a 2 = 0 ← (a 0 = 1) ∨ (a 1 = 0) (2.5)
In this example, the second line indicates that, if at time t the discrete level of G 0 is 1, then at time t + 1 the discrete level of G 1 can be updated to 1; this logic rule corresponds to the activation from G 0 to G 1 in the influence graph. The first line indicates that for any moment t, the discrete level of G 0 can be updated to 1 at time t + 1 (in fact, once the discrete level of G 0 reaches 1, it will remain at 1).

The simulation of a discrete model is not solely dependent on these logic rules, but also relies on the semantics of the model. Intuitively, the semantics dictates the number of genes that can be updated simultaneously. Various semantics, such as synchronous, asynchronous, general and most permissive, have been proposed in the literature [START_REF] Ribeiro | Learning any semantics for dynamical systems represented by logic programs[END_REF][START_REF] Paulevé | Reconciling qualitative, abstract, and scalable modeling of biological networks[END_REF]. For example, the synchronous semantics means that all genes are updated simultaneously and the asynchronous semantics means that only one gene can be updated at a time.

Based on the choice of semantics, we can obtain the transition graph of discrete states of a discrete model. A discrete state is an integer vector of length N (if there are N genes in the system), noted by d s , which assigns the discrete level d i s to gene G i , where i ∈ {1, 2, 3, ..., N } and d i s is the i th component of d s . The transition graph of discrete states is a directed graph containing all possible transitions between discrete states from t to t + 1. The transition graph of discrete states describes the dynamics of a discrete model. Note that depending on the semantics, the dynamics can be non-deterministic, as is the case for the asynchronous semantics.

For example, Fig 2 .1-Middle presents the transition graph of discrete states, derived by asynchronous semantics, of the discrete model described by the logic program in Eq (2.5). Consider the discrete state 101 (representing the assignment a 0 = 1, a 1 = 0, a 2 = 1), according to the logic program in Eq (2.5), G 1 can be updated from 0 to 1 and G 2 can be updated from 1 to 0. Since we use asynchronous semantics, only one gene can be updated, so it can reach 100 or 111, but it cannot make two updates at the same time to reach 110. Fig 2 .1-Right presents the transition graph of discrete states, derived by synchronous semantics, of the same discrete model. It can be seen that the dynamics changes when the semantics changes.

Compared to ODEs, discrete models are easier to analyze and implement, while they can still capture the qualitative dynamical properties of gene regulatory networks, which is sometimes enough in certain biological studies. Due to the dynamical simplicity of discrete models, they can be used to study large networks. A part of this thesis focuses on discrete models (Chapter 7).

Hybrid Modeling

Hybrid models are dynamical systems that contain at the same time discrete and continuous components, which, following ODEs and discrete models, have been used to model gene regulatory networks [START_REF] Belgacem | Control of negative feedback loops in genetic networks in[END_REF][START_REF] Firippi | Topology-induced dynamics in a network of synthetic oscillators with piecewise affine approximation[END_REF][START_REF] Mestl | Chaos in high-dimensional neural and gene networks[END_REF][START_REF] Edwards | Analysis of continuous-time switching networks[END_REF][START_REF] Edwards | A calculus for relating the dynamics and structure of complex biological networks[END_REF]. Like discrete models, in hybrid models, the continuous state space is divided into several discrete regions, and the dynamics in each discrete region is described by simplified ODE. Compared to ODE, generally the dynamical properties of hybrid models are easier to analyze, and compared to discrete models, we can know the continuous behaviors of systems using hybrid models.

In the following subsections, firstly the piecewise affine systems are presented which are the most studied hybrid models of gene regulatory networks, then hybrid gene regulatory networks are presented which is the major formalism studied in this thesis.

Piecewise Affine Systems (PWA)

The piecewise affine systems (PWA) [START_REF] Belgacem | Control of negative feedback loops in genetic networks in[END_REF][START_REF] Firippi | Topology-induced dynamics in a network of synthetic oscillators with piecewise affine approximation[END_REF][START_REF] Mestl | Chaos in high-dimensional neural and gene networks[END_REF][START_REF] Edwards | Analysis of continuous-time switching networks[END_REF][START_REF] Edwards | A calculus for relating the dynamics and structure of complex biological networks[END_REF] are simplifications of ODEs. One major reason behind these simplifications is that the widely used Hill functions in ODEs can be approximated by step functions. A general form of PWAs to model gene regulatory networks is given in Eq (2.6).

dx 0 dt = k 0 + g(x 1 , x 2 , ..., x m ) -γx 0 (2.6)
which is similar to Eq (2.1), but here the function g is based on step functions, for example it can be of the following form:

g(x 1 , x 2 , ..., x m ) = g 1 (x 1 ) + g 2 (x 2 ) + ... + g m (x m ) (2.7)
where

g i (x i ) =    k + i if x i ≥ θ i k - i if x i < θ i (2.8)
with k + i and k - i real numbers. In fact, these θ i divide the continuous state space into several discrete regions and in each discrete region the dynamics is described by affine functions.

Hybrid Gene Regulatory Networks (HGRN)

In this thesis, we mainly focus on a class of hybrid models called Hybrid gene regulatory networks (HGRN). The difference between HGRNs and PWAs mainly lies in the behaviors of trajectories on boundaries. HGRN and its basic notions are defined formally as follows.

Consider a gene regulatory network with N genes. For i ∈ {1, 2, 3, ..., N }, the i th gene, noted G i , has n i + 1 discrete levels which are represented by integers: {0, 1, 2, ..., n i }. The set of all discrete states is

E d = d s = (d i s ) i∈[1;N ] ∈ N N | ∀i ∈ {1, 2, ..., N } , d i s ∈ {0, 1, ..., n i } .
In fact, the notion of discrete state in HGRNs is similar as the one in discrete models.

Definition 1 (Hybrid gene regulatory network (HGRN))

A Hybrid gene regula- tory network (HGRN) is noted H = (E d , c). E d is the set of all discrete states. c is a function from E d to R N . For each d s ∈ E d , c(s)
, also noted c s , is called the celerity of discrete state d s and describes the time derivative of the system in d s .

A B

C A C B 0 0 0.6 -0.7 0 1 -0.7 -0.9 1 0 0.7 0.8 1 1 -0.6 0.9

Definition 3 (Hybrid trajectory)

A hybrid trajectory τ is a function from a time interval [0, t 0 ] to E τ = E h ∪ E sh , where t 0 ∈ R + ∪ {∞}, E h is the set of all hybrid states, and E sh is the set of all finite or infinite sequences of hybrid states:

E sh = {(h 0 , h 1 , ..., h m ) ∈ (E h ) m+1 | m ∈ N ∪ {∞}}.
A hybrid trajectory τ is called a periodic hybrid trajectory if it is defined on [0, ∞[ and

∃T > 0, ∀t ∈ [0, ∞[, τ (t) = τ (t + T ). A 2-dimensional HGRN is shown in Fig 2.2.
In this example, there are two genes, that we rename for simplicity: A = G 1 and B = G 2 . Gene A activates gene B, so we can consider that there exists a threshold of gene A such that, when the expression value of gene A is above this threshold, the derivative of gene B is positive, and when it is below this threshold, the derivative of gene B is negative. Similarly for gene B, since B inhibits gene A, when the expression value of B is above a threshold, the derivative of gene A is negative, and when it is below this threshold, the derivative of gene A is positive. Since each gene has one threshold, they have two discrete levels: 0 and 1. Therefore, this gene regulatory network has four discrete states: 00, 01, 10 and 11, where, for example, 01 is a shorthand representing that the discrete level of gene A is 0 and the discrete level of gene B is 1.

In the illustration of the state space (Fig 2 .2-Right), black arrows represent the celerities (time derivatives) of each discrete state and red arrows represent a possible hybrid trajectory of this system, which happens, in this particular case, to be a periodic hybrid trajectory.

Definition 4 (Boundary

) A boundary in a discrete state d s is a set of hybrid states defined by e(G i , π 0 , d s ) = {(π, d s ) ∈ E h | π i = π 0 , }, where i ∈ {1, 2, ..., N } , d s ∈ E d and π 0 ∈ {0, 1}.
In the rest of this thesis, we simply use e to represent a boundary.

Still in the example of Fig 2.2, the hybrid state h

M = ((π 1 M , 1), (1, 1)) of point M (π M is the fractional part of h M and π 1
M is the first value of the vector π M ) belongs to e 1 = (B, 1, (1, 1)), that is, the upper boundary in the second dimension (the dimension of gene B) of the discrete state 11. Since there is no other discrete state on the other side of e 1 , the hybrid trajectory from h M cannot cross e 1 and has to slide along e 1 (e 1 can be called a black wall). Boundaries like e 1 , which can be reached by hybrid trajectories but cannot be crossed, are defined as attractive boundaries. Let us mention that the notion of black wall has been used in other hybrid systems. Generally, the behavior on black wall is not easy to define because the derivatives might be different on the different sides of a black wall. In HGRNs, by using hybrid states, a black wall is separated into two boundaries, therefore the system can have different derivatives on the different sides of the wall. There exist other methods to define behaviors of the system on a black wall [START_REF] Gouzé | A class of piecewise linear differential equations arising in biological models[END_REF][START_REF] Plahte | Analysis and generic properties of gene regulatory networks with graded response functions[END_REF].

The hybrid state h P = ((π 1 P , 0), (0, 1)) of point P belongs to e 2 = (B, 0, (0, 1)), the lower boundary in the second dimension of the discrete state 01. The hybrid trajectory from h P reaches instantly h Q = ((π 1 Q , 1), (0, 0)), which belongs to e 3 = (B, 1, (0, 0)), the upper boundary in the second dimension of discrete state 00, because the celerities on both sides allow this (instant) discrete transition. e 2 is called an output boundary of 01 and e 3 is called an input boundary of 00. ■ If c i s < 0 (c i s is the i th value of the vector c s ), then the hybrid trajectory from the current hybrid state enters the interior of the current discrete state. e is called an input boundary of d s .

Although the hybrid trajectory in

dτ (t + ) dt = c s , dτ (t -) dt = c s ′ and τ (t) = ((π ′ , d s ′ ), (π, d s )),
where π ′k = π k for all k ̸ = i, and π ′i = 0, which means that there is an instant transition from (π ′ , d s ′ ) to (π, d s ) at t. For instance, the hybrid state Q in Fig 2 .2 belongs to an input boundary.

■ If c i s = 0, then the hybrid trajectory from the current hybrid state will slide along the boundary e, which is then called a neutral boundary of d s .

dτ (t + ) dt = dτ (t -) dt = c s and τ (t) = (π, d s ).
■ If c i s > 0, and either d i s is the maximal discrete level of the i th gene, or d i s is not the maximal discrete level of the i th gene but the i th component of c s ′ is negative, then the hybrid trajectory from the current hybrid state will slide along the boundary e, which is called an attractive boundary of d s . If τ reaches e at t, then:

dτ (t + ) dt k = c k s (where dτ (t + ) dt k is the k th component of the vector dτ (t + ) dt ) for all k ̸ = i, dτ (t + ) dt i = 0, dτ (t -) dt = c s and τ (t) = (π, d s ). If τ reaches e at t 0 < t, then: dτ (t) dt k = c k s for all k ̸ = i,
= c s ′ , dτ (t -) dt = c s and τ (t) = ((π, d s ), (π ′ , d s ′ ))
, where π ′k = π k for all k ̸ = i, and π ′i = 0. For instance, the hybrid state P in Fig 2 .2 belongs to an output boundary.

-If h belongs to several boundaries, then the previous cases can be mixed: ■ If in these boundaries there is no output boundary, then the hybrid trajectory from the current hybrid state will exit all input boundaries and slide along all attractive or neutral boundaries.

■ If in these boundaries there is only one output boundary, then the hybrid trajectory from the current hybrid state will cross this output boundary.

■ If in these boundaries there are several output boundaries, then the hybrid trajectory from the current hybrid state can cross any of them, but can only cross one boundary at one time, which causes non-deterministic behavior. HGRNs are not yet widely used to model gene regulatory networks, but they could potentially solve some problems that are difficult to solve by using ODEs or PWAs. For example, the analysis of limit cycles in ODEs and PWAs is difficult and sometimes it is not possible to prove theoretically the existence or the stability of limit cycle, while it is feasible in HGRNs; and HGRNs can be potentially applied on large networks.

This thesis majorly focuses on HGRNs. In the following chapters, we present our works on different aspects of HGRNs, including analysis of limit cycles (Chapter 3), reachability analysis (Chapter 4), condition for sustained oscillations in the canonical repressilator (Chapter 5) and parameter identification (Chapter 6).

Chapter 3

SEARCH AND STABILITY ANALYSIS OF LIMIT CYCLES

Many gene regulatory networks have periodic behavior, for instance the cell cycle or the circadian clock. Therefore, the study of formal methods to analyze limit cycles in mathematical models of gene regulatory networks is of interest. In this chapter, we study the limit cycles of HGRNs, and we present a new formal method to find all limit cycles that are simple and deterministic, and analyze their stability, that is, the ability of the model to converge back to the cycle after a small perturbation. Before this work, only 2-dimensional limit cycles of HGRNs have been studied [START_REF] Cornillon | Modèles qualitatifs de réseaux génétiques: réduction de modèles et introduction d'un temps continu[END_REF]; our work fills this gap by proposing a generic approach applicable in higher dimensions. This method is based on two major ideas: Firstly, the hybrid states are abstracted to consider only their borders, in order to enumerate all simple abstract cycles containing possible concrete trajectories; secondly, a Poincaré map is used, based on the notion of transition matrix of the concrete continuous dynamics inside these abstract paths. The lower dynamical complexity of HGRN allows the application of this method on HGRN in any dimension. We apply the method on three HGRNs of negative feedback loops with 3 components, and a HGRN of the cell cycle with 5 components. This work has been presented at CMSB 2022 as regular paper [START_REF] Sun | Limit cycle analysis of a class of hybrid gene regulatory networks[END_REF].

Limit Cycles of HGRNs

We define that a periodic hybrid trajectory of HGRN is a limit cycle if there exists another hybrid trajectory that converges to (when time approaches positive or negative infinity) (see Fig 3.1) or reaches this periodic hybrid trajectory (see Fig 3 .2).

The original definition of limit cycles, which comes from the field of non-linear dynamical systems, does not consider periodic trajectories that can be reached by other trajectories. In fact, there is no such periodic trajectories in non-linear dynamical sys-Figure 3.1 -Left: A limit cycle (the blue periodic hybrid trajectory) with another hybrid trajectory (the red hybrid trajectory) converging to it when time approaches positive infinity. Right: A limit cycle (the blue periodic hybrid trajectory) with another hybrid trajectory (the red hybrid trajectory) converging to it when time approaches negative infinity. tems, while the existence of sliding modes in HGRNs makes such behaviors exist. In our work, we choose to consider these periodic hybrid trajectories as limit cycles because they also show the asymptotic behaviors of the systems.

As there exist limit cycles reached by other hybrid trajectories, there exist also limit cycles such that from a hybrid state on the limit cycle, there exists a hybrid trajectory which leaves the limit cycle. Such limit cycles must contain non-deterministic hybrid

Limit Cycles of HGRNs

state(s). In our work, we do not consider this type of limit cycles because they are not realistic when considering real-life systems.

There exist also limit cycles of HGRN such that another hybrid trajectory converges to it with an infinite number of transitions. For example, the periodic hybrid trajectory in Fig 3 .3 with hybrid states A, B, C, D, which is a limit cycle with only instant transition.

We define that a transition from h 1 to h 2 is a non-instant transition if h 1 and h 2 are two different hybrid states in the same discrete state and the hybrid trajectory from h 1 reaches a new boundary for the first time at h 2 . For example the transition from N to M in Fig 2 .2 is a non-instant transition. And an instant transition is a transition from a hybrid state in a given discrete state to another hybrid in a different discrete state which takes no time; for example, in Fig 2 .2, the transition from P to Q is an instant transition.

We see that in Fig 3 .3, there is hybrid trajectory τ converging to the limit cycle, and the intersection between this hybrid trajectory and the boundary e 1 (the lower boundary in the second dimension of the discrete state 01) are hybrid states denoted by the sequence P 1 , P 2 , P 3 , .... From P 1 , this hybrid trajectory can be divided into shorter hybrid trajectories τ P i →P i+1 (hybrid trajectory from P i to P i+1 ) by this sequence of intersection hybrid states. The duration of each shorter hybrid trajectory τ P i →P i+1 is denoted by t i . Based on the above notations, the whole duration of τ is the sum t ϵ + t 1 + t 2 + t 3 + ... where t ϵ is the duration of τ before reaching P 1 . By the similarity between each τ P i →P i+1 , we can derive that this sum can be expressed by t ϵ + t 1 + αt 1 + α 2 t 1 + α 3 t 1 ... where α is a positive real number that is strictly less than 1. In fact, we can see that t (by the property of similarity), so α = LD ED . Since t 1 + αt 1 + α 2 t 1 + α 3 t 1 ... is the sum of a geometric sequence of which the ratio (α) is between 0 and 1, the whole duration of τ is a finite value, which means that τ converges to this limit cycle in finite time.

In fact, only limit cycles that have only instant transitions can have hybrid trajectories converging to them in finite time. All hybrid states in such limit cycle are related to the same point in Euclidean space. We say that the hybrid state (π, d s ) and the point x ∈ R N in Euclidean space are related if x = π + d s . Such periodic hybrid trajectory could be called a Zeno fixed point. In this chapter, we do not consider such limit cycles because they are not associated to limit cycles in Euclidean space.

There are also some periodic hybrid trajectories that are not limit cycles (see Fig 3.4).

Search of Periodic Hybrid Trajectories

This section presents our methods to find periodic hybrid trajectories which are potential limit cycles. In this chapter, we make two assumptions about periodic hybrid trajectories in HGRNs.

Assumption 1 Any non-instant transition on a periodic hybrid trajectory does not reach more than one new boundary at the same time.

Assumption 2 For any hybrid state h on a periodic hybrid trajectory, there is at most one output boundary to which h belongs. In other words, any periodic hybrid trajectory does not have non-deterministic hybrid state.

For Assumption 1, in real-life systems, it is indeed very unlikely for parameters to be that constrained due to measurement noise. An example of a non-instant transition reaching two new boundaries at the same time is given in Fig 2 .3.

If Assumption 1 is satisfied, then there is no non-instant transition that reaches a non-deterministic hybrid state. If we further assume that a threshold of one gene only influences at most one another gene, which is a reasonable assumption for the modeling of gene regulatory networks, see for example [START_REF] Abou-Jaoudé | From structure to dynamics: frequency tuning in the p53-mdm2 network: I. logical approach[END_REF], then Assumption 2 is satisfied.

The major reason why we add these two assumptions is that the cases that do not satisfy these two assumptions are very rare and more complicated to analyze. So it is not interesting to consider them in practice. However, our approaches can be generalized to these cases.

Our method to find periodic hybrid trajectories has three steps which are described in the following.

Abstraction of HGRN with Discrete Domains

First, the HGRN is transformed into a graph of discrete domains. A discrete domain is a new concept proposed in our work which is defined as follows.

Definition 5 (Discrete domain) A discrete domain D(d s , S -, S + ) is a set of hybrid states inside one discrete state d s , defined by:

D(d s , S -, S + ) = {(π, d s ) | ∀i ∈ {1, 2, ..., N }, π i ∈          {1} if i ∈ S + {0} if i ∈ S - ]0, 1[ if i ̸ ∈ S -∪ S + }
where S + and S -are subsets of {1, 2, ..., N } such that S + ∩ S -= ∅ and S + ∪ S -̸ = ∅. In fact, S + (S -) represents the dimensions in which the upper (lower) boundaries are reached by any hybrid state h ∈ D(d s , S -, S + ). In the rest of this thesis, we simply use D to represent a discrete domain when there is no ambiguity.

In the rest of this thesis, as a notation, we add exponents to the vector representation of a discrete state to indicate which upper (lower) boundaries are reached. For instance, 11 + denotes the discrete domain inside discrete state 11 where the upper boundary is It is then possible to build the graph of discrete domains of a HGRN, such as in Fig 3.6, where the nodes are the discrete domains and the edges are computed by considering only the signs of the celerities. In this graph, a discrete domain D j is a successor of discrete domain D i if:

-There exists an instant transition from D i to D j , which means that at least one hybrid trajectory from D i crosses a boundary and instantly reaches D j ; see for example 0 + 0 and 1 -0 in is a hybrid trajectory from 00 + which reaches 0 + 0. We exclude cases where two new boundaries are reached at the same time; for instance, there is no edge between 00 + and 0 + 0 -.

Need to mention that this idea of discrete abstraction is commonly used for the analysis of hybrid dynamical systems [START_REF] Alur | Discrete abstractions of hybrid systems[END_REF][START_REF] Alur | Counter-example guided predicate abstraction of hybrid systems[END_REF].

Search of Closed Discrete Trajectories

Before giving the definition of closed discrete trajectories and the method to find them, some concepts are introduced at first. A walk in the graph of discrete domains gives a sequence of discrete domains T = (D 0 , D 1 , ...D p ). A hybrid trajectory is said to be inside such a sequence of discrete domains if it begins from the first discrete domain and reaches by order all discrete domains in the sequence. Based on this, we define two new notions on such a sequence: The transition matrix, which allows to compute the final hybrid state of a hybrid trajectory inside a given sequence of discrete domains, when it exists, and the compatible zone, which is the set of initial hybrid states so that such a hybrid trajectory exists. These two notions are presented formally as follows. Note that, there exist hybrid trajectories which are inside sequences of discrete domains that are not in the graph of discrete domains introduced in the previous subsection because the cases where two new boundaries are reached at the same time are excluded for the construction of the graph. In fact, we do not need to consider such sequences of discrete domains because Assumption 1 is made. Definition 6 (Transition matrix) Consider two different discrete domains D i and D j such that there exists a sequence of discrete domains T from D i to D j . If there exists a hybrid state h i = (π i , d s i ) in D i such that from h i there is a hybrid trajectory τ (defined on [0, t 0 ]) which is inside T and reaches D j on h j = (π j , d s j ) at t 0 , then there exists a transition matrix M which describes the relation between π i and π j , that is:

π j = s -1 (M s(π i ))
, where s is a function that adds an extra dimension and the value in the extra dimension is always 1: s((a 1 , a 2 , ..., a N )) = (a 1 , a 2 , ..., a N , 1). The transition matrix M only depends on T . Consider Fig 3 .5, the transition matrix of the sequence of discrete domains (00

+ , 0 + 0) is M (00 + ,0 + 0) =      0 0 1 - c 2 00 c 1 00 1 c 2 00 c 1 00 0 0 1     
where (c 1 00 , c 2 00 ) is the celerity vector of discrete state 00.

M (00 + ,0 + 0) can describe the relation between hybrid state h B = ((π 1 B , π 2 B = 1), (0, 1)) and hybrid state

h C = ((π 1 C = 1, π 2 C ), (0, 0)):      π 1 C π 2 C 1      = M (00 + ,0 + 0)      π 1 B π 2 B 1      . ( 3.1)
Similarly, for the sequence of discrete domains (0 + 0, 1 -0), the transition matrix is

M (0 + 0,1 -0) =      1 0 -1 0 1 0 0 0 1     
which describes the relation between

h C = ((π 1 C = 1, π 2 C ), (0, 0)) and h D = ((π 1 D = 0, π 2 D ), (1, 0)):      π 1 D π 2 D 1      = M (0 + 0,1 -0)      π 1 C π 2 C 1      . ( 3.2)
The transition matrix can be computed by the multiplication of the transition matrices of all elementary sequences (sequences of length 2), for example M (00 + ,0 + 0,1 -0) = M (0 + 0,1 -0) M (00 + ,0 + 0) .

To illustrate the method to compute the transition matrix of an elementary sequence of discrete domains, the following examples are given.

First, we consider an elementary sequence in a 4-dimensional HGRN that is not instant (meaning that all hybrid trajectories inside it are not instant): (a + bc -d, a + bcd + ) where abcd is a discrete state. The celerity vector of this discrete state is denoted by (c

t 0 = 1 -π 4 i c 4 s (3.3)
In other dimensions, the relations between (1, π 2 i , 0, π 4 i ) and (1, π 2 j , π 3 j , 1) can be described based on t 0 :

1 = 1 + 0 × t 0 (3.4) π 2 j = π 2 i + c 2 s × t 0 (3.5)
π 3 j = 0 + c 3 s × t 0 (3.6)
The temporal derivative in the first dimension is 0 (Eq (3.4)) because of the existence of the sliding mode. By reformulating all these above relations and eliminating t 0 , the transition matrix M (a + bc -d,a + bcd + ) is obtained (Eq (3.7)).

           1 π 2 j π 3 j 1 1            =            1 0 0 0 0 0 1 0 - c 2 S c 4 s c 2 s c 4 s 0 0 1 -c 3 s c 4 s c 3 s c 4 S 0 0 0 0 1 0 0 0 0 1                       1 π 2 i 0 π 4 i 1            (3.7)
Second, we consider another elementary sequence in a 4-dimensional HGRN that is instant (meaning that all hybrid trajectories inside it are instant): (a + bcd + , a + bc(d + 1) -). Any hybrid trajectory inside (a + bcd + , a + bc(d + 1) -) can be represented by ((1, π 2 j , π 3 j , 1), abcd) -→ ((1, π 2 j , π 3 j , 0), abc(d + 1)). Since the only difference between (1, π 2 j , π 3 j , 1) and (1, π 2 j , π 3 j , 0) is the value in the fourth dimension, the transition matrix M (a + bcd + ,a + bc(d+1) -)) can be obtained by modifying only one line of the identity matrix to switch the value in the fourth dimension from 1 to 0 (Eq (3.8)).

           1 π 2 j π 3 j 0 1            =            1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 -1 0 0 0 0 1                       1 π 2 j π 3 j 1 1            (3.8)
The general method to compute the transition matrix M (D i ,D i+1 ) of any elementary sequence of discrete domains (D i , D i+1 ) is given as follows.

-Case 1: The transition between D i and D i+1 is instant (e.g. (0 + 0, 1 -0)). In this case, by assuming that the boundary in dimension n 0 is crossed,

M (D i ,D i+1 ) = (a kl ) ∈ R (N +1)×(N +1)
, where -

a kl = 1, if k = l. -a kl = x, if k = n 0 ∧ l = N + 1, where -x = -1,
, M (D i ,D i+1 ) = (a kl ) ∈ R (N +1)×(N +1)
, where -

a kl = 1, if k = l ̸ = n 0 . -a kl = -c k c n 0 , if k ̸ = n 0 ∧ l = n 0
and there is no sliding mode in dimension k (c is the celerity vector of the discrete state containing D i and D i+1 , and c k is the -

k th component of c). -a kl = c k c n 0 x, if k ̸ = n 0 ∧ l = N + 1
a kl = x, if k = n 0 ∧ l = N + 1.
(The value of x is same as the above case.) a kl = 0, otherwise.

In case 2 of the general method to compute the transition matrix, a priori, several new boundaries can be reached, in such case, we can choose the dimension of any of them as n 0 . While in this chapter, since we only consider sequences on the graph of discrete domains, there is always a unique n 0 for any elementary sequence that is not instant.

The second proposed notion associated to sequences of discrete domains is the compatible zone. Definition 7 (Compatible zone) Consider a sequence of discrete domains T = (D 0 , D 1 , D 2 , ..., D m ). The compatible zone S is the maximal subset of D 0 such that from any hybrid state of S, there is a hybrid trajectory that contains a sub-trajectory inside T . More formally, for any hybrid state h ∈ S, there is a hybrid trajectory τ from h and there exists t 0 such that the restriction of τ on [0, t 0 ] is a hybrid trajectory inside T .

For a general sequence of discrete domains (a sequence of discrete domains that is not necessarily in the graph of discrete domains), the compatible zone is of the form:

S = {(π, d s 0 ) ∈ D 0 | W π > c ∧ W ′ π = c ′ }
where W, W ′ are matrices and c, c ′ are vectors. However, the compatible zone S of any sequence of discrete domains in the graph of discrete domains can be expressed only with linear inequalities: S = {(π, d s 0 ) ∈ D 0 | W π > c} because, inside these sequences of discrete domains, there is no hybrid trajectory that has a non-instant transition reaching two new boundaries at the same time. The idea to compute the compatible zone is based on Theorem 1.

., m}).

Proof: Proof of sufficient condition: We can easily see that if h = (π, d s 0 ) belongs to the compatible zone S of T = (D 0 , D 1 , D 2 , ..., D m ), then ∀i ∈ {1, 2, ..., m}, h also belongs to the compatible zone of (D

0 , D 1 , D 2 , ..., D i ), so (s -1 (M (D 0 ,D 1 ,...,D i ) s(π)), d s i ) ∈ D i .
Proof of necessary condition: By induction. Consider a sequence of discrete domains of length 2: (D 0 , D 1 ), it is evident that h belongs to the compatible zone of (D 0 , D 1 ) if (s -1 (M (D 0 ,D 1 ) s(π)), d s 1 ) ∈ D 1 . Now suppose that it is true for any sequence of discrete domains of length k + 1, and consider a sequence of discrete domains of length k + 2: (D 0 , D 1 , D 2 , ..., D k+1 ). Since it is true for a sequence of discrete domains of length k + 1, h belongs to the compatible zone of (D 0 , D 1 , D 2 , ..., D k ), so the hybrid trajectory from h will stay inside (D 0 , D 1 , D 2 , ..., D k ) and will reach

D k at h k = (s -1 (M (D 0 ,D 1 ,...,D k ) s(π)), d s k ). Let h k = (π k , d s k ). We can easily see that s -1 (M (D 0 ,D 1 ,...,D k+1 ) s(π)) = s -1 (M (D k ,D k+1 ) s(π k )). Since (s -1 (M (D 0 ,D 1 ,...,D k+1 ) s(π)), d s k+1 ) ∈ D k+1 , we have (s -1 (M (D k ,D k+1 ) s(π k )), d s k+1 ) ∈ D k+1
. So h k belongs to the compatible zone of (D k , D k+1 ). Therefore, h belongs to the compatible zone of (D 0 , D 1 , D 2 , ..., D k+1 ).

□

Thus, all hybrid states that satisfy the constraints in Theorem 1 constitute the compatible zone. And we can see that the computation of compatible zone replies on the computation of transition matrix. Now we give the definition of closed discrete discrete trajectories. A sequence of discrete domains T is called a discrete trajectory if the compatible zone of T is not empty, in other words, there exists at least one hybrid trajectory inside T . A discrete trajectory

T = (D 1 , D 2 , ...D m ) is said closed if D 1 = D m .
For example, the sequence of discrete domains (01 -, 00 + , 0 + 0, 1 -0, 10 + ) in Fig 3 .5 and Fig 3 .6 is a discrete trajectory because its compatible zone is not empty as there exists at least one hybrid trajectory inside this sequence (e.g. the hybrid trajectory from A to E in Fig 3 .5).

In order to find closed discrete trajectories, we use a depth first algorithm. For this, we rely on the notion of Poincaré section (more explanations about Poincaré section will be given in the next subsection, here we can simply consider that a Poincaré section is a boundary of a discrete state). We first choose one or several input boundaries of discrete states as Poincaré sections, more precisely an input boundary is chosen if it is crossed by a cycle of discrete states, and then on each discrete domain on these Poincaré sections, we apply this depth first algorithm. In each step of this depth first algorithm, the compatible zone is computed and the search continues if the compatible zone is not empty and the current path does not reach the same discrete domain for more than one time and does not return to the initial discrete state. This algorithm finds all discrete trajectories that begin from a discrete domain and return to the initial discrete state without crossing the same discrete state more than once. An execution of this algorithm on discrete domain 00 + of the HGRN in Fig 2 .2 is illustrated in Fig 3 .7. Among these discrete trajectories, we can easily find the closed ones. Note that in Fig 3 .7, all discrete trajectories from 00 + that return to 00 happen to be closed but generally we can find discrete trajectories that are not closed.

Consider the HGRN in Fig 2 .2, we can easily see that there is only one cycle of discrete states in this system, which is 00 -→ 10 -→ 11 -→ 01 -→ 00. Therefore, for this system, we only need one Poincaré section and any boundary in this cycle can take this role. Let us choose for instance the input boundary of discrete state 00 (the boundary between 01 and 00) as Poincaré section, that is, the union of the three discrete domains 0 -0 + , 00 + and 0 + 0 + . We thus apply the depth first algorithm on each of these three discrete domains. As a result, we can find 5 discrete trajectories that begin from the Poincaré section and return to the initial discrete state:

1 : 0 -0 + - → 00 -- → 0 + 0 -- → 1 -0 -- → 10 + - → 11 -- → 11 + - → 1 -1 + - → 0 + 1 + - → 01 -- → 00 + 2 : 00 + - → 0 + 0 - → 1 -0 - → 10 + - → 11 -- → 1 -1 - → 0 + 1 - → 01 -- → 00 + 3 : 00 + - → 0 + 0 - → 1 -0 - → 10 + - → 11 -- → 11 + - → 1 -1 + - → 0 + 1 + - → 01 -- → 00 + 4 : 00 + - → 00 -- → 0 + 0 -- → 1 -0 -- → 10 + - → 11 -- → 11 + - → 1 -1 + - → 0 + 1 + - → 01 -- → 00 + 5 : 0 + 0 + - → 1 -0 + - → 1 -1 -- → 0 + 1 -- → 0 + 0 +
Examples of hybrid trajectories inside each of these 5 discrete trajectories are shown in Fig 3 .8. We note that there always exists at least one hybrid trajectory inside a discrete trajectory since, by definition, its compatible zone is not empty. Among the 5 discrete trajectories above, only the first one is not closed.

Identification of Periodic Hybrid Trajectories with

Poincaré Map

Consider a closed discrete trajectory T = (D 0 , D 1 , ...D m , D 0 ). We define that a periodic hybrid trajectory τ is inside T if τ begins from a hybrid state h ∈ D 0 , reaches by order all discrete domains of T and finally reaches back the hybrid state h. To check if there -∃(π 0 , d s 0 ) ∈ D 0 such that s -1 (M T s(π 0 )) = π 0 , and -(π 0 , d s 0 ) belongs to the compatible zone of T .

Then (π 0 , d s 0 ) is called a fixed point of T , the function s -1 (M T s()) is the Poincaré map of the periodic hybrid trajectory that corresponds to this fixed point and D 0 (or the boundary that contains D 0 ) is the associated Poincaré section.

Under Assumption 1, any periodic hybrid trajectory of the system must be inside one of the closed discrete trajectories found by the depth first algorithm. Meanwhile, if a periodic hybrid trajectory reaches more than one new boundary at the same time, then it is not inside any closed discrete trajectories found by the algorithm.

Among the five closed discrete trajectories in the HGRN in Fig 2 .2, we find only one periodic hybrid trajectory of interest, inside the third closed discrete trajectory; it is the hybrid trajectory labeled "3" in Fig 3 .8. For this periodic hybrid trajectory, π 0 = (0.222, 1),

M T =      0 0 0.222 0 0 1 0 0 1     
, and the compatible zone is Actually, there is another periodic hybrid trajectory inside the fifth closed discrete trajectory; it is the periodic hybrid trajectory labeled "5" in Fig 3 .8. This periodic hybrid trajectory only contains instant transitions. As stated before, our analysis method of limit cycles does not consider this type of periodic hybrid trajectories.

The Poincaré map has been widely used in the literature to study limit cycles of hybrid systems [START_REF] Belgacem | Control of negative feedback loops in genetic networks in[END_REF][START_REF] Firippi | Topology-induced dynamics in a network of synthetic oscillators with piecewise affine approximation[END_REF][START_REF] Mestl | Chaos in high-dimensional neural and gene networks[END_REF][START_REF] Edwards | Analysis of continuous-time switching networks[END_REF][START_REF] Edwards | A calculus for relating the dynamics and structure of complex biological networks[END_REF][START_REF] Znegui | Design of an explicit expression of the Poincaré map for the passive dynamic walking of the compass-gait biped model[END_REF][START_REF] Flieller | Computation and stability of limit cycles in hybrid systems[END_REF][START_REF] Girard | Computation and stability analysis of limit cycles in piecewise linear hybrid systems[END_REF][START_REF] Hiskens | Stability of hybrid system limit cycles: Application to the compass gait biped robot[END_REF]. The dynamical simplicity of HGRNs and the existence of sliding modes make the form of Poincaré maps of HGRNs particular: A Poincaré map is a piecewise affine function on a Poincaré section.

Stability Analysis of Limit Cycles

In order to define the stability of limit cycles, we give the definition of the neighborhood in the same discrete state at first.

Definition 8 (Neighborhood in the same discrete state)

The neighborhood in the same discrete state of a hybrid state h = (π 0 , d s ) is a set of hybrid states defined as:

N d (h, r) = (π, d s ) | d(π, π 0 ) < r, π ∈ [0, 1] N ,
with r > 0 the radius of this neighborhood, and d the maximum norm between vectors:

d(π, π 0 ) = max i∈{1,2,...,N } | π i -π i 0 |.
Based on the definition of neighborhood in the same discrete state, the stability of limit cycle is defined as follows.

Definition 9 (Stability of limit cycles in HGRNs)

A limit cycle τ of HGRN is stable if, for any hybrid state h on τ , there exists a neighborhood in the same discrete state of radius r such that any hybrid trajectory τ 0 that begins from this neighborhood N d (h, r) satisfies: lim t→∞ (Dis min (τ 0 (t), τ )) = 0 where Dis min (h ′ , τ ) is defined as

Dis min (h ′ , τ ) = min h 0 ∈τ d(x(h ′ ), x(h 0 )), with h ′ ∈ E h , x(.

) the sum (dimension by dimension) of the fractional part and the discrete state of a hybrid state, and

d(v 1 , v 2 ) = max 1≤i≤n |v i 1 -v i 2 | (v 1 , v 2 are two vectors of length n), otherwise it is unstable.
It is noteworthy that in most cases, a value of t high enough is sufficient to obtain Dis min (τ 0 (t), τ ) = 0, without needing a limit computation. This definition of stability is illustrated in Fig 3 .10 by a 2-dimensional stable limit cycle.

Continuity in the Neighborhood of Periodic Hybrid Trajectories

We call neighborhood of a hybrid trajectory a union of neighborhoods in the same discrete state of all the hybrid states in this hybrid trajectory. A periodic hybrid trajectory is said to respect the continuity of neighborhood if there exists a neighborhood of this

(τ ) = (h 1 , h 2 , ..., h 9 ).
hybrid trajectory that is small enough so that all hybrid trajectories starting from this neighborhood remain in this neighborhood. The continuity of neighborhood is formally defined after the following definitions.

Definition 10 (Neighborhood in the same discrete domain)

The neighborhood in the same discrete domain of a hybrid state h = (π 0 , d s ) is a set of hybrid states defined as:

N D (h, r) = {(π, d s ) ∈ D i | d(π, π 0 ) < r}, where D i is the discrete domain which includes h.
Definition 11 (Turning states of a periodic hybrid trajectory) Consider a periodic hybrid trajectory τ , we can find a sequence of hybrid states on τ : (h 1 , h 2 , ..., h m ), such that each h i represents a hybrid state when τ reaches a new discrete domain and there is a transition from h i to h i+1 where i ∈ {1, 2, ..., m -1} and from h m to h 1 . All hybrid states in this sequence are the turning states of τ and this sequence is called the sequence of turning states of τ denoted by T s (τ ).

The sequence of turning states of a periodic hybrid trajectory is given in Fig 3 .11.

Definition 12 (Continuity of neighborhood of a periodic hybrid trajectory)

The neighborhood of a periodic hybrid trajectory τ is continuous if for any neighborhood in the same discrete domain of a turning state h 0 of τ : N D (h 0 , r 0 ) and for any hybrid state h 1 on τ , there exists a small neighborhood in the same discrete state of h 1 : N d (h 1 , r 1 ), such that all hybrid trajectories from N d (h 1 , r 1 ) will reach N D (h 0 , r 0 ) in finite time.

When a periodic hybrid trajectory does not have continuity of neighborhood, some hybrid trajectories in the neighborhood may undergo a "disruption" by touching another boundary and thus follow another sequence of discrete states.

Without Assumption 1 and Assumption 2, some neighborhoods of a periodic hybrid trajectory might not respect this continuity, no matter how small they are. For example, consider a periodic hybrid trajectory that contains a transition ((1, 1), (a, b)) → ((0, 1), (a+1, b)), where the upper boundaries in the first and second dimensions of discrete state ab are both output boundaries. In the neighborhood of ((1, 1), (a, b)), no matter how small it is, we can always find a hybrid state from which the hybrid trajectory reaches the boundary in the second dimension at first, and as it will reach a different discrete state, it might never return to the neighborhood of this periodic hybrid trajectory.

We claim that Assumption 1 and Assumption 2 together constitute a sufficient condition for the continuity of neighborhood of any periodic hybrid trajectories in HGRNs. To prove this, we give some definitions and propositions as follows.

Proposition 1 Consider a periodic hybrid trajectory τ with the sequence of turning states

T s (τ ) = (h 1 , h 2 , .

.., h m ). For any two adjacent turning states h i and h

j (j = i + 1 if i ∈ {1, 2, ..., m -1}, j = 1 if i = m)
, for any neighborhood in the same discrete domain of h j : N D (h j , r j ), if Assumption 1 and Assumption 2 are satisfied, there exists a neighborhood in the same discrete domain of h i : N D (h i , r i ), such that any hybrid trajectory from N D (h i , r i ) reaches directly N D (h j , r j ) (we say that a hybrid trajectory reaches directly N D (h j , r j ) if it reaches a new boundary for the first time at a hybrid state that belongs to N D (h j , r j )).

Proposition 1 is illustrated in Fig 3 .12.

Proof: If the transition between h i and h j is an instant transition, then h i and h j are related to the same point in Euclidean space, with Assumption 2, we can easily prove this proposition by choosing the same N D (h i , r i ) as N D (h j , r j ). In fact, without Assumption 2, a priori, h i can be non-deterministic, in such case we do not have Proposition 1.

If the transition is not instant, without loss of generality, consider a 4-dimensional example:

h i = ((1, 0, π 3 i , π 4 i ), d s ) and h j = ((1, π 2 j , 1, π 4 j ), d s ).
We suppose that the celerity Figure 3.12 -Illustration of Proposition 1. h i and h j are two adjacent turning states on a periodic hybrid trajectory. Black boxes represent neighborhoods in the same discrete domain of h i /h j that correspond to Proposition 1.

related to this transition is (c 1 s , c 2 s , c 3 s , c 4 s ) and we can see that there is a sliding mode in the first dimension. The duration of this transition is noted as t 0 , we have:

1 = 1 + t 0 × 0 π 2 j = 0 + t 0 × c 2 s π 4 j = π 4 i + t 0 × c 4 s (3.9)
where

t 0 = 1 -π 3 i c 3 s . (3.10)
Any hybrid state from the neighborhood in the same discrete domain of h i can be represented by

h ′ i = ((1, 0, π 3 i + δ 3 , π 4 i + δ 4 ), d s )
. Suppose that the hybrid trajectory from h ′ i still reaches directly the discrete domain which includes h j at hybrid state

h ′ j = ((1, π ′2 j , 1, π ′4 j ), d s )
. The duration of this new transition becomes:

t 1 = 1 -π 3 i -δ 3 c 3 s = 1 -π 3 i c 3 s - δ 3 c 3 s = t 0 - δ 3 c 3 s . ( 3.11)
And

1 = 1 + t 1 × 0 π ′2 j = 0 + t 1 × c 2 s π ′4 j = π 4 i + δ 4 + t 1 × c 4 s (3.12)
By combining Eq (3.9), Eq (3.11) and Eq (3.12), we have:

π ′2 j = π 2 j - δ 3 c 3 s × c 2 s π ′4 j = π 4 i + δ 4 - δ 3 c 3 s × c 4 s (3.13)
To ensure that h ′ j is included in N D (h j , r j ), firstly these constraints must be satisfied:

| δ 3 c 3 s × c 2 s |< r j | δ 4 - δ 3 c 3 s × c 4 s |< r j (3.14)
Secondly, we need to ensure that the hybrid trajectory from h ′ i reaches the boundary in the third dimension at first, so these constraints are also needed to be satisfied:

1 -π 3 i -δ 3 c 3 s < 1 -π 4 i -δ 4 c 4 s 1 -π 3 i -δ 3 c 3 s < 1 -0 c 2 s (3.15)
Here we suppose that c 4 s is positive. We can derive from Eq (3.14) and Eq (3.15) by using triangle inequality that a sufficient condition, for the constraints in Eq (3.14) and Eq (3.15) to be satisfied, is:

| δ 3 |< r j | c 3 s | | c 2 s | | δ 3 | + | δ 4 |< r j max(1, |c 4 s | |c 3 s | ) | δ 3 | + | δ 4 |< 1-π 4 i c 4 s - 1-π 3 i c 3 s max( 1 |c 3 s | , 1 |c 4 s | ) | δ 3 |<| c 3 s | ( 1 c 2 s - 1 -π 3 i c 3 s ) (3.16)
From any hybrid state h ′ i = ((1, 0, π 3 i + δ 3 , π 4 i + δ 4 ), d s ) that satisfies the above constraints (Eq (3.16)), the hybrid trajectory reaches directly N D (h j , r j ). So any hybrid trajectories from N D (h i , r i ) reaches directly N D (h j , r j ) if r i satisfies:

r i < min(r j | c 3 s | | c 2 s | , 1 2
r j max(1, |c 4 s | |c 3 s | ) , 1 2
1-π 4 i c 4 s - 1-π 3 i c 3 s max( 1 |c 3 s | , 1 |c 4 s | ) , | c 3 s | ( 1 c 2 s - 1 -π 3 i c 3 s )) (3.17)
In the above demonstration, we consider a transition of a HGRN in four dimensions and propose a method to find a N D (h i , r i ) that satisfies the Proposition 1. The same method can be used to find such a N D (h i , r i ) in any HGRN in N dimension (the only difference is that there will be more inequalities in Eq (3.14) and Eq (3.15)), if Assumption 1 and Assumption 2 are satisfied. In fact, Assumption 1 ensures that the constraints on these δ can be represented by some inequalities of the sum of absolute values of these δ (like Eq (3.16)). If Assumption 1 is not satisfied, then we need stronger constraints on these δ represented by linear equations. □ Proposition 2 Consider a periodic hybrid trajectory τ with the sequence of turning states

T s (τ ) = (h 1 , h 2 , ..., h m ).
For any two turning states h i and h j , for any neighborhood in the same discrete domain of h j : N D (h j , r j ), if Assumption 1 and Assumption 2 are satisfied, there exists a neighborhood in the same discrete domain of h i : N D (h i , r i ), such that all hybrid trajectories from N D (h i , r i ) reach by order all discrete domains between h i and h j (the discrete domains containing respectively h i ,h i+1 ,h i+2 ,...,h j-1 ,h j ) and finally reach N D (h j , r j ).

Proposition 2 is illustrated in Fig 3 .13.

Proof: Consider two turning states h i and h j , since τ is a periodic hybrid trajectory, we can find turning states between h i and h j , noted as h i+1 , h i+2 , ..., h j-1 , such that there are non-instant or instant transitions between h k and h k+1 , ∀k ∈ {i, i + 1, ..., j -1}. (If j < i, we can take h i as the first turning state in the series of turning states of τ and renumber h j ). Consider a neighborhood in the same discrete domain of h j : N D (h j , r j ), according to Proposition 1 we can find a neighborhood in the same discrete domain of

h j-1 : N D (h j-1 , r j-1 )
, such that all hybrid trajectories from N D (h j-1 , r j-1 ) reach directly N D (h j , r j ). Then, in the same way, we can find a neighborhood in the same discrete domain of h j-2 : N D (h j-2 , r j-2 ), such that all hybrid trajectories from N D (h j-2 , r j-2 ) reach directly

N D (h j-1 , r j-1
). Using this backward method, we can finally find a neighborhood in the same discrete domain of h i : N D (h i , r i ), such that all hybrid trajectories from N D (h i , r i ) firstly reach N D (h i+1 , r i+1 ), then reach N D (h i+2 , r i+2 ), ..., and finally reach N D (h j , r j ). □ Proposition 3 Consider a periodic hybrid trajectory τ with the sequence of turning states

T s (τ ) = (h 1 , h 2 , .

.., h m ). For any two adjacent turning states h i and h

j (j = i + 1 if i ∈ {1, 2, ..., m -1}, j = 1 if i = m
) such that the transition from h i to h j is a non-instant transition (the hybrid trajectory from h i to h j is denoted by τ i ), for any neighborhood in the same discrete domain of h j : N D (h j , r j ) and for any hybrid state h k on τ i (including

h i but not including h j ), if Assumption 1 is satisfied, there exists a neighborhood in the same discrete state of h k : N d (h k , r k ), such that all hybrid trajectories from N d (h k , r k ) reach N D (h j , r j ).
Proposition 3 is illustrated in Fig 3 .14.

Proof: Consider two adjacent turning states:

h i = ((1, 0, π 3 i , π 4 i ), d s ) and h j = ((1, π 2 j , 1, π 4 j ), d s )
(same example used in the proof of Proposition 1) and the neighborhood in the same discrete domain of h j : N D (h j , r j ). Firstly, we consider h i . Any hybrid state in a neighborhood of the same discrete state of

h i can be represented by h ′ i = ((1 + δ 1 , 0 + δ 2 , π 3 i + δ 3 , π 4 i + δ 4 ), d s )
(we can see that δ 1 must be negative and δ 2 must be positive). We suppose that the hybrid trajectory from h ′ i firstly reaches the upper boundary in the first dimension (an attractive boundary) at the hybrid state

h ′′ i = ((1, 0 + δ 2 -δ 1 c 2 s c 1 s , π 3 i + δ 3 -δ 1 c 3 s c 1 s , π 4 i + δ 4 -δ 1 c 4 s c 1 s ), d s ).
To ensure that the hybrid trajectory from h ′′ i reaches N D (h j , r j ), we can use the same method as in the proof of Proposition 1 to get some inequalities of the sum of absolute values of these δ (like Eq (3.16)) and we call these linear inequalities Constraint 1. To ensure that the hybrid trajectory from h ′ i firstly reaches the upper boundary in the first dimension, the following inequalities must be satisfied:

-δ 1 c 1 s < 1 -δ 2 c 2 s -δ 1 c 1 s < 1 -δ 3 -π 3 i c 3 s -δ 1 c 1 s < 1 -δ 4 -π 4 i c 4 s (3.18)
Figure 3.14 -Illustration of Proposition 3. h i and h j are two adjacent turning states on a periodic hybrid trajectory. h k is a hybrid state on the hybrid trajectory from h i to h j . Black box represents a neighborhood in the same discrete domain of h j and blue box represents a neighborhood in the same discrete state of h k . These neighborhoods correspond to Proposition 3.

By using triangle inequality, we can get a sufficient condition of Equation 3.18, described by inequalities of the sum of absolute values of these δ:

| δ 1 | + | δ 2 |< 1 c 2 s max( 1 |c 1 s | , 1 |c 2 s | ) | δ 1 | + | δ 3 |< 1-π 3 i c 3 s max( 1 |c 1 s | , 1 |c 3 s | ) | δ 1 | + | δ 4 |< 1-π 4 i c 4 s max( 1 |c 1 s | , 1 |c 4 s | ) (3.19)
These inequalities are called Constraint 2.

By combing Constraint 1 and Constraint 2, same as in the proof of Proposition 1 ,we can get a N d (h i , r i ) such that all hybrid trajectories from N d (h i , r i ) reach N D (h j , r j ). In this example, there is only one sliding boundary (the first dimension) between h i and h j . If there are several sliding boundaries, to compute the constraints of δ, we need to ensure that any of these sliding boundaries is reached before any non sliding boundary, so there will be more inequalities in Eq (3.18). But we can always find a sufficient condition of these inequalities which is described by inequalities of the sum of absolute values of these δ.

For any other hybrid state h k in τ i (not including h j ), we can use exactly the same method to find such a N d (h k , r k ).

□ Proposition 4 Consider a periodic hybrid trajectory τ with the sequence of turning states

T s (τ ) = (h 1 , h 2 , ..., h m ). h i 0 , h i 0 +1 , ..., h j 0 is a sequence in T s (τ ) such that ∀i ∈ {i 0 , i 0 + 1, .
.., j 0 -2} the transition from h i to h i+1 is an instant transition and the transition from h j 0 -1 to h j 0 is a non-instant transition (such sequence exists because in this chapter we do not consider periodic hybrid trajectories with only instant transitions). For any neighborhood in the same discrete domain of h j 0 : N D (h j 0 , r j 0 ) and for any hybrid state h i , i ∈ {i 0 , i 0 + 1, ..., j 0 -2}, if Assumption 1 and Assumption 2 are satisfied, then there exists a neighborhood in the same discrete state of

h i : N d (h i , r i ), such that all hybrid trajectories from N d (h i , r i ) reach N D (h j 0 , r j 0 ).
Proof: Without loss of generality, we consider four turning states h 1 , h 2 , h 3 , h 4 , where

h 1 = (π 1 = (1, 1, 1, π 4 1 ), d s 1 = (a, b, c, d)), h 2 = (π 2 = (0, 1, 1, π 4 1 ), d s 2 = (a + 1, b, c, d)), h 3 = (π 3 = (0, 0, 1, π 4 1 ), d s 3 = (a + 1, b + 1, c, d)), h 4 = (π 4 = (π 1 4 , π 2 4 , π 3 4 , 1), d s 3 = (a + 1, b + 1, c, d))
. We can see that there is an instant transition from h 1 to h 2 which crosses the upper boundary in the first dimension, an instant transition from h 2 to h 3 which crosses the upper boundary in the second dimension and a non-instant transition from h 3 to h 4 . The celerity vector of

d s 1 is c s 1 = (c 1 s , c 2 s , c 3 s , c 4 s
) and we suppose that all values in c s 1 are positive. The celerity vector of

d s 2 is c s 1 = (c 1′ s , c 2′ s , c 3 s , c
= (π ′ 1 = (1 + δ 1 , 1 + δ 2 , 1 + δ 3 , π 4 1 + δ 4 ), d s 1 = (a, b, c, d))
. We suppose that in d s 1 , the upper boundaries in the second and third dimension are attractive boundaries and the upper boundary in the fourth dimension is an output boundary (there is only one reached output boundary because of Assumption 2). After the upper boundary e 1 in the first dimension in d s 1 is crossed, in d s 2 the upper boundary e 2 in the second dimension becomes an output boundary (e 1 is related to a threshold of first gene that influences the second gene), so e 2 is crossed by the hybrid trajectory.

We want that the hybrid trajectory from h ′ 1 will firstly cross the upper boundary in the first dimension and then cross the upper boundary in the second dimension. To ensure this, we only need to ensure that e 1 and e 2 are reached before the upper boundary e 4 in the fourth dimension is reached. We do not need to ensure that e 1 is reached before e 2 because e 2 is an attractive boundary. We need to satisfy the following constraints:

- δ 1 c 1 s < 1 -δ 4 -π 4 i c 4 s - δ 2 min(c 2 s , c 2′ s ) < 1 -δ 4 -π 4 i c 4 s (3.20)
A sufficient condition, for these constraints, represented by inequalities of the sum of absolute values of δ, to be satisfied, is:

| δ 1 | + | δ 4 |< 1-π 4 i c 4 s max( 1 |c 1 s | , 1 |c 4 s | ) | δ 2 | + | δ 4 |< 1-π 4 i c 4 s max( 1 |min(c 2 s ,c 2′ s )| , 1 |c 4 s | ) (3.21)
This condition is noted Constraint 1.

After the hybrid trajectory from h ′ 1 crosses e 2 , the hybrid state can be represented by

h ′ 3 = (π ′ 3 = (0 + ∆ 1 , 0, 1 + ∆ 3 , π 4 1 + ∆ 4 ), d s3 = (a + 1, b + 1, c, d))
where ∆ 1 , ∆ 3 and ∆ 4 are linear combinations of δ 1 , δ 2 , δ 3 and δ 4 . To ensure that the hybrid trajectory from h ′ 3 reaches N D (h 4 , r 4 ), using a similar method as in the proof of Proposition 3, we can get constraints represented by inequalities of the sum of absolute values of these δ (like Eq (3.16)). By combining these constraints with Constraint 1, we can get N d (h 1 , r 1 ) such that all hybrid trajectories from N d (h 1 , r 1 ) reach N D (h 4 , r 4 ).

The above demonstration only considers one specific example, but this method can be used to find such N d (h i , r i ) in any case if Assumption 1 and Assumption 2 are satisfied. In this example, we consider a sequence of two instant transitions. If there is a sequence of N instant transitions, we only need to ensure that any boundary, that is crossed in the sequence, is reached before any other boundary is reached, so we have more inequalities in Eq (3.20). □ Now based on the above propositions, we have the following theorem.

Theorem 2 The neighborhood of a periodic hybrid trajectory is continuous if Assumption 1 and Assumption 2 are satisfied.

Proof: Consider any neighborhood in the same discrete domain of a turning state h 0 of τ : N D (h 0 , r 0 ) and consider any h 1 on τ .

If h 1 is not a turning state of τ , then we can find two adjacent turning states h i and h j such that there is a non-instant transition between h i and h j (a straight hybrid trajectory) and h 1 belongs to this straight hybrid trajectory. According to Proposition 2, we can find a neighborhood in the same discrete domain of h j : N D (h j , r j ), such that all hybrid trajectories from N D (h j , r j ) will finally reach N D (h 0 , r 0 ). Then according to Proposition 3, we can find a neighborhood in the same discrete state of h 1 : N d (h 1 , r 1 ), such that all hybrid trajectories from N d (h 1 , r 1 ) will reach N D (h j , r j ). So all hybrid trajectories from N d (h 1 , r 1 ) will firstly reach N D (h j , r j ) and finally reach N D (h 0 , r 0 ).

If h 1 is a turning state of τ and the transition following h 1 is a non-instant transition, that means that there is a straight hybrid trajectory from h 1 that reaches the next turning state h 2 . Then, same as the previous case, we can find a neighborhood in the same discrete state of h 1 : N d (h 1 , r 1 ), such that all hybrid trajectories from N d (h 1 , r 1 ) will finally reach

N D (h 0 , r 0 ).
If h 1 is a turning state of C τ and the transition following h 1 is an instant transition, that means that the hybrid trajectory from h 1 firstly crosses a boundary and reaches the next turning state h 2 . We note h j as the first turning state after h 2 such that the transition before h j is a non-instant transition. According to Proposition 2, we can find a neighborhood in the same discrete domain of h j : N D (h j , r j ), such that all hybrid trajectories from N D (h j , r j ) will finally reach N D (h 0 , r 0 ). Then according to Proposition 4, we can find a neighborhood in the same discrete state of h 1 : N d (h 1 , r 1 ), such that all hybrid trajectories from N d (h 1 , r 1 ) will reach N D (h j , r j ). So all hybrid trajectories from N d (h 1 , r 1 ) will firstly reach N D (h j , r j ) and finally reach N D (h 0 , r 0 ). Now we have proved that for any h 1 in τ , we can find a neighborhood in the same discrete state of h 1 : N d (h 1 , r 1 ), such that all hybrid trajectories from N d (h 1 , r 1 ) will finally reach N D (h 0 , r 0 ). So the neighborhood of τ is continuous. □

Eigenanalysis

In this subsection, we present our method to analyze the stability of limit cycles. Consider a periodic hybrid trajectory τ inside the closed discrete trajectory T = (D 1 , D 2 , ...D m , D 1 ). Suppose that τ reaches D 1 at h = (π, d s 1 ), thus we have:

π = s -1 (M T s(π)) (3.22)
For π, there might be some dimensions in which the values are 0 or 1 because in these dimensions the upper or lower boundaries are reached. If we only consider the dimensions in which the boundaries are not reached, Eq (3.22) becomes:

x = Ax + b (3.23)
where x, called the reduction vector of h (denoted by x = r(h)), is a short version of π which only contains the dimensions in which the boundaries are not reached. The matrix A is called the reduction matrix of T and vector b is called the constant vector of T . Based on the concept of reduction vector, the compatible zone of a sequence of discrete domains

T = (D 0 , D 1 , D 2 , ..., D m ) can be described by S = {(π, d s 0 ) ∈ D 0 | r(π) ∈ S r }
where S r , called the reduction compatible zone, is a set of reduction vectors of hybrid states in D 0 .

If the length of x is not 0, meaning that D 1 does not contain only a singular hybrid state (a counter example is the discrete domain 1 + 1 + in Fig 3 .5 which contains only a singular hybrid state), then the stability analysis method of the limit cycles is based on -

If max i∈{1,2,...,p} |λ i | < 1 then τ is a stable limit cycle. -If max i∈{1,2,...,p} |λ i | > 1 then τ is an unstable limit cycle. -If max i∈{1,2,...,p} |λ i | = 1 and ∃i 0 ∈ {1, 2, ..., p} , |λ i 0 | < 1, then τ is an unstable limit cycle. -If ∀i ∈ {1, 2, ..., p} , |λ i | = 1, then τ is not a limit cycle.
Proof: Consider a periodic hybrid trajectory τ that exists inside a closed discrete trajectory T . The intersection of τ with the chosen Poincaré section e is h 0 = (π 0 , d s 0 ). The Poincaré map is noted as x k+1 = Ax k + b, where x is the reduction vector of the hybrid states in the compatible zone of T (x 0 = r(h 0 )). The stability of the fixed point(s) of the system x k+1 = Ax k + b depends on the eigenvalues of A.

If the absolute values of all eigenvalues of A are less than 1, then x 0 is asymptotically stable for the system x k+1 = Ax k + b. And since the neighborhood of τ is continuous, we can find a small neighborhood in the same discrete domain of h 0 : N D (h 0 , r 0 ), such that any hybrid trajectory from N D (h 0 , r 0 ) stays inside the neighborhood of τ and converges asymptotically to or reaches τ . Also, based on the fact that the neighborhood of τ is continuous, for any hybrid state h ′ on τ , we can find a neighborhood in the same discrete state of h ′ : N d (h ′ , r), such that any hybrid trajectory from N d (h ′ , r) reaches N D (h 0 , r 0 ). Thus, for any hybrid trajectory τ ′ from N d (h ′ , r), we have: lim t→∞ Dis min (τ ′ (t), τ ) = 0, which proves that τ is a stable limit cycle.

If the maximum absolute value of all eigenvalues of A is greater than 1, then x 0 is unstable for system x k+1 = Ax k + b, so we can always find a hybrid trajectory from any small neighborhood in the same discrete domain of h 0 that does not converge to or reach τ . Therefore, τ is unstable.

If the maximum absolute value of all eigenvalues of A is 1 and there is at least one eigenvalue that differs from 1, then in any small neighborhood in the same discrete domain of h 0 , we can find hybrid trajectories that do not converge to or reach τ (which might converge to or reach another limit cycle), so τ is unstable.

If the absolute values of all eigenvalues are 1, then all hybrid trajectories from any small neighborhood in the same discrete domain of h 0 are periodic hybrid trajectories (like Fig 3 .4), so τ is not a limit cycle.

□

In case that the length of x is 0, in other words D 1 contains only a singular hybrid state, if there exists i ∈ {1, 2, ..., m} such that D i does not contain only a singular hybrid state, then D i can be chosen as the new D 1 to apply Theorem 3. If any discrete domain in this sequence contains only a singular hybrid state (see the example in Fig 3 .15), then we have directly that τ is a stable limit cycle because the neighborhood of τ is continuous.

Applications

In this section, we apply the above limit cycle analysis method on three HGRNs of negative feedback loop in 3 dimensions and one HGRN of cell cycle in 5 dimensions. The negative feedback loop in 3 dimensions can be used to describe real biological oscillators, for example the p53 system [START_REF] Geva-Zatorsky | Oscillations and variability in the p53 system[END_REF]. The signs of the celerities in these three HGRNs are determined by the influence graph (positive for an activation and negative for an inhibition) and their absolute values of celerities are randomly selected. The parameters of the HGRN in 5 dimensions are generated randomly respecting the constraints in Table 3 of [START_REF] Behaegel | A hybrid model of cell cycle in mammals[END_REF]. The influence graphs of both systems can be found in Fig 3 .16. Details about the implementation can be found at https://doi.org/10.5281/zenodo.6524936.

A B C C A C B C C 0 0 0 1 -0.
6 -0.7 0 0 1 1 0.7 -0.9 0 1 0 -0.8 -0.8 -0.7 0 1 1 -0.8 0.6 -0.9 1 0 0 0.7 -0.6 0.6 1 0 1 0.7 0.7 0.5 1 1 0 -0.9 -0.8 0.6 1 1 1 -0.9 0.6 0.5

A B C C A C B C C 0 0 0 3 -0.
6 -0.7 0 0 1 3 0.7 -2.9 0 1 0 -2.8 -0.8 -0.7 0 1 1 -2.8 0.6 -2.9 1 0 0 2.7 -0.6 2.6 1 0 1 2.7 0.7 0.5 1 1 0 -2.9 -0.8 2.6 1 1 1 -2.9 0.6 0.5

A B C C A C B C C 0 0 0 3 -0.
6 -0.7 0 0 1 3 0.7 -2.9 0 1 0 -0.8 -0.8 -0.7 0 1 1 -0.8 0.6 -2.9 1 0 0 0.7 -0.6 2.6 1 0 1 0.7 0.7 0.5 1 1 0 -2.9 -0.8 2.6 1 1 1 -2.9 0.6 0.5

3-dimensional Models of Negative Feedback Loop

The parameters of these three HGRNs of negative feedback loop in 3 dimensions are shown in Table 3.1. The signs of the celerities in these three models are the same so they have the same graph of discrete states. There is only one cycle of discrete states in each of these systems, which is:

111 - → 011 - → 010 - → 000 - → 100 - → 101 - → 111
Therefore, for these three models, we choose the input boundary e of 000 in the cycle as the Poincaré section. Simulations depicting the convergence to the stable cycle or to the fixed point (see below) in these three HGRNs are shown in Fig 3 .17.

In the first HGRN, by using our limit cycle analysis method, we find one sta-ble limit cycle and one periodic hybrid trajectory which only contains instant transitions, which is considered here as a fixed point of this system. Regarding the stable limit cycle, the fixed point of this limit cycle in discrete domain 0 -0 + 0 is ((0, 1, 0.125), (0, 0, 0)), the transition matrix is

        0 0 0 0 0 0 0 1 0 0 0 0.125 0 0 0 1        
, the compatible zone is {(π, (0, 0, 0)) | π 1 = 0, π 2 = 1, π 3 ∈]0, 0.7[} and the reduction matrix is 0 , therefore hybrid trajectories from the neighborhood of this limit cycle will reach this limit cycle very quickly (less than one turn if the neighborhood is small enough).

In this HGRN we can also prove that all hybrid trajectories will reach this limit cycle except hybrid trajectories which can reach the fixed point of the system. All discrete trajectories which begin from the Poincaré section and return to the Poincaré section in this HGRN are shown in Fig 3 .18 A. Since in this HGRN there is only one cycle of discrete states which is also a global attractor, any hybrid trajectory from the Poincaré section must return to the Poincaré section and it must begin from the compatible zone or the boundary of the compatible zone of one of the discrete trajectories in Fig 3 .18 A. We see that all discrete trajectories which are not closed will finally reach closed discrete trajectories [START_REF] Gouzé | A class of piecewise linear differential equations arising in biological models[END_REF][START_REF] Plahte | Analysis and generic properties of gene regulatory networks with graded response functions[END_REF][START_REF] Cornillon | Modèles qualitatifs de réseaux génétiques: réduction de modèles et introduction d'un temps continu[END_REF][START_REF] Boyenval | What is a cell cycle checkpoint? the TotemBioNet answer in Computational Methods in Systems Biology[END_REF][START_REF] Behaegel | A hybrid model of cell cycle in mammals[END_REF]. Discrete trajectories 31, 32 and 33 have the same transition matrix and their reduction matrix is 0 so any hybrid trajectory from 0 -0 + 0 will reach the limit cycle. For the discrete trajectory 10, the two eigenvalues of the reduction matrix are 7.0306 and 0.0368, so hybrid trajectories inside discrete trajectory 10 will finally leave the compatible zone and reach 0 -0 + 0. An example of such hybrid trajectory is illustrated in Fig 3 .19, this hybrid trajectory starts at the hybrid state h 1 on the Poincaré section and returns to the Poincaré section for the first time at hybrid state h 2 . The compatible zone of the discrete trajectory 10 is illustrated by the red triangle (not including the boundaries) and v 1 and v 2 represent the two eigenvectors. Since the absolute value of the first eigenvalue is bigger than 1 and the absolute value of the second eigenvalue is smaller than 1, h 2 is farther from the up-right corner (which is a fixed point on the Poincaré section) in the direction of v 1 and is closer to the up-right corner in the direction of v 2 , so this hybrid trajectory leaves the compatible zone of the closed discrete trajectory 10 and reaches the compatible zone of the closed discrete trajectory 11. From here, we can see that any hybrid trajectories from the Poincaré section will reach the limit cycle except the hybrid trajectory inside discrete trajectory 9 which is related to a fixed point. As any hybrid trajectory in this system will finally reach this Poincaré section, all hybrid trajectories will reach this limit cycle except hybrid trajectories which can reach the fixed point.

For the second HGRN, by using our method, we can also find one stable limit cycle and one fixed point. Unlike the first HGRN, hybrid trajectories from the neighborhood of the limit cycle converge asymptotically to the limit cycle: The limit cycle is inside the discrete trajectory which begins from 0 -0 + 0 and the reduction matrix of the limit cycle is 0.0298 . We can also prove that all hybrid trajectories will converge to this limit cycle except hybrid trajectories which can reach the fixed point, by using the same method as for the first HGRN.

Contrary to the first and the second HGRN, we cannot find a limit cycle in the third HGRN but only a fixed point which is related to the discrete trajectory 2 in Fig 3 .18 C. For the discrete trajectories 3 and 4, the fixed point of the Poincaré map of discrete trajectories 4 is strictly inside the compatible zone of discrete trajectory 3, the fixed point of the Poincaré map of the discrete trajectory 3 is the fixed point of this HGRN and the absolute values of all eigenvalues of the reduction matrix of discrete trajectories 3 and 4 are strictly inferior than 1 (the two eigenvalues of the reduction matrix of discrete trajectory 3 are 0.3091 and 0.0062, and the two eigenvalues of the reduction matrix of the discrete trajectory 4 are 0 and 0.1269), so we can prove that all hybrid trajectories in this system will converge to the fixed point. An illustration of this fixed point is given in

5-dimensional Model of Cell Cycle

For the HGRN in 5 dimensions, the transition graph of discrete states is more complex (Fig 3.20). By using a depth first algorithm, we find that there are 1104 cycles of discrete states in which 930 cycles contain the discrete transition 01010 -→ 01011, 94 cycles contain 00011 -→ 01011 and all the rest contain 00110 -→ 00010. So for this system, we use 3 Poincaré sections which are the input boundary of 01011 between 01010 and 01011, the input boundary of 01011 between 00011 and 01011 and the input boundary of 00010 between 00110 and 00010. Any limit cycle of this system must cross one of these Poincaré sections. By using the method above, we find 184 closed discrete trajectories which cross the first Poincaré section, among these 184 closed discrete trajectories, there is only one stable limit cycle; we find 50 closed discrete trajectories which cross the second Poincaré section, these 50 closed discrete trajectories do not contain any limit cycle; and we find 111 closed discrete trajectories which cross the third Poincaré section, among these 111 closed discrete trajectories, there is one stable limit cycle and one unstable limit cycle. In fact, these two stable limit cycles are the same which crosses both the first and the third Poincaré sections. This stable limit cycle is the same one studied in [START_REF] Behaegel | A hybrid model of cell cycle in mammals[END_REF] to calculate the constraints of parameters. The simulations of both cycles are shown in Fig 3 .21 A andB. For now we have not identified any biological behavior related to this unstable limit cycle yet. For the stable limit cycle of the cell cycle model, the fixed point of this limit cycle in the discrete domain 01 + 01 + 1 -is ((0.3714, 1, 0.8581, 1, 0), (0, 1, 0, 1, 1)), and the reduction matrix is

  0 0 0 0   .
For the unstable limit cycle of cell cycle model, the fixed point of this limit cycle in the discrete domain 000 + 10 -is ((0.6375, 0.2552, 1, 0.3472, 0), (0, 0, 0, 1, 0)), and the reduction matrix A is

     4.95359512 • 10 3 0 1.37489884 • 10 -13 -5.25996267 • 10 2 0 -1.45993292 • 10 -14 -7.15619779 • 10 2 0 -1.98624389 • 10 -14     
. The eigenvalues of A are 0, 4.95359512 • 10 3 and 3.15544362 • 10 -30 , making it unstable.

The current implementation in Python reaches its limits with respect to execution time when the size of the system increases: Finding the limit cycles above takes less than one minute for the HGRNs in 3 dimensions, and 8 hours for the HGRN in 5 dimensions 1 . Figure 3.21 -Simulation of the two limit cycles found in the HGRN of 5 dimensions. A: Stable limit cycle. B: Unstable limit cycle.

Summary

This chapter presented a method to find all limit cycles of HGRNs with some minor restrictions, mainly to remove non-deterministic behaviors and complex loops, and to analyze their stability. This method is the first one to find and analyze limit cycles of HGRNs in N dimensions. We showed the merits of this method on randomly generated HGRNs of a negative feedback loop with 3 components and a HGRN of the cell cycle with 5 components taken from the literature.

As stated above, a first limitation of this method is that we do not handle complex periodic hybrid trajectories which are periodic hybrid trajectories crossing a discrete state for more than once in one period. Such periodic hybrid trajectories are not considered because they are hard to enumerate.

Another limitation in the application of this method is that we first need to construct a HGRN of a specific gene regulatory network; however, the observation of real biological systems is limited and it is not always possible to determine all parameters. In some cases, some parameters can only be described by constraints, or remain unknown. Thus, considering extensions of this method that are parameterized or that take into account a set of constraints on parameters is of interest. Another work about analysis of oscillations in a parameterized HGRN of a specific network is presented in chapter 5.

In this chapter, we focused on the existence of limit cycles and their stability, which is one way to describe the dynamical properties of a system. They are of particular interest because they represent the long-term behaviors of the system. In the next chapter, we Chapter 4

REACHABILITY ANALYSIS

In this chapter, a new reachability method on HGRNs is presented. The reachability problem concerned in this chapter is, given a singular hybrid state and a region (a set of hybrid states), to determine whether the hybrid trajectory from this singular hybrid state will reach this region. This problem is undecidable for general hybrid automata, and is known to be decidable only for a restricted class of hybrid automata, but this restricted class does not include HGRNs. So, a priori, this reachability problem in HGRNs is not decidable; however, we find that it is decidable in some cases. Based on this fact, the main idea of this work is that if the decidable cases can be determined automatically, then the reachability problem can be solved partially. The method of this chapter relies on the concepts introduced in the previous chapter, e.g., transition matrix, discrete domain, compatible zone, etc. The two major contributions are the following: firstly, we classify hybrid trajectories into different classes and provide theoretical results about decidability; then based on these theoretical results, we propose a reachability analysis algorithm which always stops in finite time and answers partially the reachability problem, meaning that it gives the correct answer if it is not inconclusive. Finally, by applying our implementation of this method on a 5-dimensional HGRN of cell cycle, we show a potential application of this method to estimate the basin of attraction. This work has been presented at RP 2023 as regular paper [START_REF] Sun | Reachability Analysis of a Class of Hybrid Gene Regulatory Networks in Reachability Problems[END_REF].

Problem Statement

One important question concerning the analysis of dynamical system is reachability: Whether a certain state (or set of states) is reachable from an initial state (or set of states), which is widely used to describe the dynamical properties of gene regulatory networks [START_REF] Batt | Symbolic reachability analysis of genetic regulatory networks using discrete abstractions[END_REF][START_REF] Xue | Reachable set estimation for genetic regulatory networks with time-varying delays and bounded disturbances[END_REF].

Reachability analysis methods have been studied on different formalisms, mainly on discrete systems [START_REF] Folschette | Sufficient conditions for reachability in automata networks with priorities[END_REF][START_REF] Paulevé | Reduction of qualitative models of biological networks for transient dynamics analysis[END_REF][START_REF] Chai | Static analysis and stochastic search for reachability problem[END_REF] and hybrid systems [START_REF] Asarin | Reachability analysis of dynamical systems having piecewise-constant derivatives[END_REF][START_REF] Batt | Symbolic reachability analysis of genetic regulatory networks using discrete abstractions[END_REF][START_REF] Henzinger | What's decidable about hybrid automata[END_REF][START_REF] Asarin | Low dimensional hybrid systems-decidable, undecidable, don't know[END_REF][START_REF] Frehse | SpaceEx: Scalable verification of hybrid systems in Computer Aided Verification[END_REF][START_REF] Dang | Reachability Analysis for Polynomial Dynamical Systems Using the Bernstein Expansion[END_REF][START_REF] Sandler | Deciding Reachability for Piecewise Constant Derivative Systems on Orientable Manifolds in Reachability Problems: 13th International Conference[END_REF]. In this chapter, we study a There are generally two kinds of reachability problems in hybrid systems: The first one is, given two states (or sets of states) A and B, to determine whether trajectories (or a trajectory) from A reach(es) B [START_REF] Asarin | Reachability analysis of dynamical systems having piecewise-constant derivatives[END_REF][START_REF] Asarin | Low dimensional hybrid systems-decidable, undecidable, don't know[END_REF][START_REF] Sandler | Deciding Reachability for Piecewise Constant Derivative Systems on Orientable Manifolds in Reachability Problems: 13th International Conference[END_REF][START_REF] Oliveira Oliveira | Mortality and Edge-to-Edge Reachability are Decidable on Surfaces in 25th[END_REF]; the second one is to compute all states reachable from certain initial states [START_REF] Frehse | SpaceEx: Scalable verification of hybrid systems in Computer Aided Verification[END_REF][START_REF] Dang | Reachability Analysis for Polynomial Dynamical Systems Using the Bernstein Expansion[END_REF][START_REF] Alur | Reachability analysis of hybrid systems via predicate abstraction in Hybrid Systems: Computation and Control: 5th International Workshop[END_REF][START_REF] Asarin | Approximate reachability analysis of piecewise-linear dynamical systems in Hybrid Systems: Computation and Control: Third International Workshop[END_REF]. The reachability problem concerned in this chapter belongs to the first case, which is to determine whether the hybrid trajectory from certain hybrid state can reach a certain region (a set of hybrid states) with the restriction that we only consider single hybrid trajectories, that is, starting from a single hybrid state and avoiding any non-determinism. The reachability problem of this work is defined formally as follows, where we assume that the system has N genes (N dimensions). The following assumptions are made in this chapter.

Problem 1 (Reachability) Consider a hybrid state

h ini = (π ini , d s ini ) and a region R target = (π, d starget ) | π i ∈ [a i , b i ], i ∈ {1, 2, ..., N } , where ∀i ∈ {1, 2, ..., N } , a i , b i ∈ R and 0 ≤ a i ≤ b i ≤ 1.
Assumption 3 For any sequence of discrete domains T of which the compatible zone is not empty, we assume that all eigenvalues of the reduction matrix of T are real.

For now, we have not found a reduction matrix with complex eigenvalues. For the definition of reduction matrix, see Section 3.3.2.

Assumption 4

The hybrid trajectory from h ini has no non-deterministic behavior.

Generally, hybrid trajectories with non-deterministic behaviors exist, but among state-ofthe-art HGRNs of gene regulatory networks, the probability of a randomly chosen initial hybrid state that leads to non-deterministic behaviors is almost 0. Therefore, we ignore this kind of hybrid trajectory in this chapter. Actually, the method of this work could also be adapted for non-deterministic hybrid trajectories (each time when a non-deterministic hybrid state is reached, the current hybrid trajectory splits into two or several hybrid trajectories, and the same method is applied on each of these new hybrid trajectories).

Assumption 5 Any non-instant transition on a limit cycle does not reach more than one new boundary at the same time.

This assumption is very similar to Assumption 1 in Chapter 3. In fact, Assumption 5 is less restrictive than Assumption 1, because there exist periodic hybrid trajectories that are not limit cycles and have non-instant transition reaching more than one new boundary at the same time; HGRNs with such periodic hybrid trajectories are ignored in Chapter 3 but are not ignored in this chapter.

In this chapter, we mainly focus on the decidability problem, that is, whether we can find an algorithm to determine the reachability problem such that this algorithm always stops in finite time and gives a correct answer.

The decidability problem among some hybrid systems that are close to HGRNs has been investigated in the literature. It has been proved that, for PCD systems (Piecewise-Constant Derivative systems), it is decidable in 2 dimensions [START_REF] Maler | Reachability analysis of planar multi-linear systems in Computer Aided Verification[END_REF] but it is undecidable in 3 dimensions [START_REF] Asarin | Reachability analysis of dynamical systems having piecewise-constant derivatives[END_REF]. For general hybrid automata, there exists a restricted class called initialized rectangular automata which is decidable in any dimension [START_REF] Henzinger | What's decidable about hybrid automata[END_REF], but this class does not include HGRNs.

Up to now, there is no theoretical result of the decidability of the reachability problem on HGRNs. A priori, we can expect that it is not decidable because of the existence of chaos. However, if we can show that it is decidable in certain cases, for example, when the hybrid trajectory considered in a reachability problem converges asymptotically to a n-dimensional limit cycle, and if these cases can be identified automatically, then the reachability problem can be answered partially. Based on this idea, we classified the hybrid trajectories of HGRNs into different classes which are introduced in the next section.

Different Classes of Hybrid Trajectories

We classify hybrid trajectories of HGRNs into three classes: Hybrid trajectories halting in finite time, hybrid trajectories attracted by cycles of discrete domains and chaotic hybrid trajectories, which are introduced as follows.

Hybrid Trajectories Halting in Finite Time

A hybrid trajectory τ is a hybrid trajectory halting in finite time if ∃t 0 such that the time derivative of τ (t 0 ) is 0 in any dimension; in other words, τ (t 0 ) is a fixed point. Therefore, the discrete domains reached by a hybrid trajectory halting in finite time form a finite sequence like (D 1 , D 2 , ..., D m ). The hybrid trajectory in Fig 4 .1 is a hybrid trajectory halting in finite time.

Theorem 4 Problem 1 is decidable if the hybrid trajectory from h ini is a hybrid trajectory halting in finite time.

Proof: In this case, the hybrid trajectory is a composition of a finite number of ndimensional "straight lines"; to verify if this hybrid trajectory reaches R target , we only need to verify if any of these "straight lines" crosses R target , which can be verified in finite time. □

Hybrid Trajectories Attracted by Cycles of Discrete Domains

A hybrid trajectory τ is a hybrid trajectory attracted by a cycle of discrete domains if ∃t 0 such that after t 0 , τ always stays inside a cycle of discrete domains

C T = (D 0 , D 1 , D 2 , ..., D p , D 0 )
, meaning that τ crosses this cycle for an infinite number of times without leaving it. Intuitively, if a hybrid trajectory τ is attracted by a cycle of discrete domains, then τ converges to or reaches a periodic hybrid trajectory. We To prove the decidability of hybrid trajectories attracted by cycles of discrete domains, we introduce the notion of predecessor in the same discrete state.

Definition 13 (Predecessor in the same discrete state) For any set of hybrid states in the same discrete state defined by R = {(π, d s ) | π ∈ E} where E ⊆ [0, 1] N is a closed set, the predecessor of R in the same discrete state, noted by P re ds (R), is the union of sets of hybrid states: P re ds (R) = i∈{1,2,...,q} Z i , such that:

-∀i ∈ {1, 2, ..., q}, ∃D i (a discrete domain) on an input boundary of d s , such that Proof: According to the Definition 13, we can easily see that if τ reaches P re ds target (R target ) then it must reach R target . In case that R target is reached by τ , the only way that P re ds target (R target ) is not reached by τ is that h ini is also in the discrete state d starget and τ reaches directly R target inside d starget , however since we add the condition "τ from h ini has already crossed at least one discrete state without reaching the region R target ", this is not possible. So if R target is reached by τ , then P re ds target (R target ) is also reached by τ . □ Based on Proposition 5, we give the following theorem.

Z i ⊆ D i . -∀i, j ∈ {1, 2, ..., q} , i ̸ = j, D i and D j are different discrete domains. -∀i ∈ {1, 2, ..., q}, ∀h ∈ Z i , the hybrid trajectory τ from h (that is, τ (0) = h) first reaches R at t 0 (that is, τ (t 0 ) ∈ R and ∀t ∈ [0, t 0 [, τ (t) / ∈ R)

Theorem 5 Problem 1 is decidable if the hybrid trajectory from h ini is a hybrid trajectory attracted by a cycle of discrete domains.

The idea of Theorem 5 can be explained intuitively by In Fig 4 .2 left, the hybrid trajectory which reaches hybrid state A, noted by τ , can be considered as two hybrid trajectories: The first one is the part of τ before reaching A and the second one is the part of τ after reaching A. This first one can be considered as a hybrid trajectory halting in finite time so whether it reaches R target is decidable, and in this example it does not reach R target . For the second one, these two following statements can be verified automatically in finite time: 1. The intersection points between this hybrid trajectory and the "right" boundary of discrete state 01 must be located in the line segment AB (because it converges to the limit cycle with only instant transitions); 2. The line segment AB does not intersect with the predecessor of R target in the same discrete state (the predecessor is illustrated by blue boxes). Based on these two statements, we can prove that this second part cannot reach R target either. In this way, we prove theoretically that R target is not reached by τ , and since this process can be done automatically in finite time, the problem is decidable. Note that in the general case, this "line segment AB" is a (n-1)-dimensional region such that the hybrid trajectory always returns to this region and this region does not intersect the predecessor of R target in the same discrete state.

In Fig 4 .2 right, it can be verified automatically in finite time that the limit cycle with only instant transitions (at the center) reaches R target , and that τ converges to this limit cycle, so we can prove that τ finally reaches R target , and this case is thus decidable too.

The formal proof of Theorem 5 is given as follows.

Proof: For a hybrid trajectory τ attracted by a cycle of discrete domains, ∃t 0 such that after t 0 , τ always stays in a cycle of discrete domains C T = (D 0 , D 1 , D 2 , ..., D p , D 0 ) and before t 0 , τ has already crossed at least one discrete state. τ can be separated on two hybrid trajectories by t 0 : The restriction of τ on time interval [0, t 0 ], noted by τ [0,t 0 ] , and the restriction of τ on time interval [t 0 , ∞], noted by τ [t 0 ,∞] .

For the hybrid trajectory τ [0,t 0 ] , whether it reaches R target is decidable, the proof is the same as the proof of Theorem 4. If it reaches R target , then it is reachable. And if it does not reach R target , then we need to consider τ [t 0 ,∞] .

For the hybrid trajectory τ [t 0 ,∞] , in case that τ [0,t 0 ] does not reach R target , according to Proposition 5, whether τ [t 0 ,∞] reaches R target is equivalent to whether it reaches P re ds target (R target ) = i∈{1,2,...,q} Z i . If ∀i ∈ {0, 1, 2, ..., p}, ∀j ∈ {1, 2, 3, ..., q}, Z j ̸ ⊆ D i , then P re ds target (R target ) is not reachable. Otherwise, if ∃i 0 ∈ {0, 1, 2, ..., p} and ∃j 0 ∈ {1, 2, 3, ..., q}, such that Z j 0 ⊆ D i 0 , then to verify if Z j 0 is reachable, we need to consider the relation between the intersection points of τ [t 0 ,∞] with D i 0 noted by an infinite sequence of hybrid states (h p1 , h p2 , ...) and Z j 0 as follows. When a hybrid trajectory reaches a discrete domain, it might slide inside this discrete domain, in this case, based on the general meaning of "intersection", the intersection points of this hybrid trajectory and this discrete domain are not singular points. In this chapter, we define that the intersection points between a hybrid trajectory and a discrete domain are the points where the hybrid trajectory just reaches the discrete domain.

Since τ always stays inside this cycle of discrete domains, the sequence (h p1 , h p2 , ...) converges asymptotically to a hybrid state or reaches a hybrid state which belongs to the closure of the compatible zone of C T . For both cases, this particular hybrid state is noted by h ∞ .

If h ∞ ∈ Z j 0 , then P re ds target (R target ) is reachable. In case that h ∞ is on the boundary of Z j 0 and the sequence (h p1 , h p2 , ...) converges asymptotically to h ∞ without reaching reaching directly Z j 0 , Z j 0 is also considered as reachable, which in fact is reached after infinite times of intersections.

If h ∞ / ∈ Z j 0 , then we can find a small neighborhood of h ∞ inside D i 0 , noted by N h∞ , such that ∀h ∈ N h∞ , h / ∈ Z j 0 . Since the sequence (h p1 , h p2 , ...) converges to h ∞ or finally reaches h ∞ , we can find n 0 such that ∀n > n 0 , h pn ∈ N h∞ . This means that after that τ [t 0 ,∞] reaches h pn 0 , it will never reach Z j 0 . So, in this case, to verify if Z j 0 is reached, we only need to verify if the finite sequence (h p1 , h p2 , ..., h pn 0 ) ever reaches Z j 0 , which is decidable. Now we see that, for a certain Z i , whether τ [t 0 ,∞] reaches Z i or not is decidable. Since P re ds target (R target ) has a finite number of such Z i , whether τ [t 0 ,∞] reaches P re ds target (R target ) is also decidable. □

From the above proof of Theorem 5, we can derive a method to check the reachability of hybrid trajectories attracted by cycles of discrete domains, which is introduced later in Section 4.3.

We also develop the following theorem to determine automatically whether a hybrid trajectory is attracted by a cycle of discrete domains. In order to simplify this theorem, for the cycle of discrete domains C T = (D 0 , D 1 , D 2 , ..., D p , D 0 ) and the hybrid state h 0 considered in this theorem, we note that:

-The reduction matrix and the constant vector of C T are A and b respectively (the reduction matrix and the constant vector are defined in Section 3.3.2). -The reduction compatible zone of C T is described by linear constraints {x | W x > c} where c is a vector and W is a matrix. W is of size n 0 × n 1 , where n 1 is the number of dimensions of r(h 0 ) (r(h 0 ) is the reduction vector of h 0 which is defined in Section 3.3.2). W i is the

i th line of matrix W (W i is of size 1 × n 1 ) and c i is the i th component of vector c. -r ∞ = lim n→∞ f n (r(h 0 )) where f (x) = Ax + b.
-The eigenvalues and eigenvectors of A are {λ i | i ∈ {1, 2, ..., n 1 }} and {v i | i ∈ {1, 2, ..., n 1 }} respectively. λ 1 is chosen as the eigenvalue with the maximum absolute value among the eigenvalues that differ from 1.

-The decomposition of r(h 0 ) -r ∞ in the directions of eigenvectors of the reduction matrix A is noted as

r(h 0 ) -r ∞ = n 1 i=1 α i v i .
Theorem 6 A hybrid trajectory τ is attracted by a cycle of discrete domains if and only if τ reaches h 0 which belongs to the compatible zone of a cycle of discrete domains

C T = (D 0 , D 1 , D 2 , ..., D p , D 0
) such that D 0 has no free dimension (meaning that, in D 0 , boundaries are reached in all dimensions) or the following conditions are satisfied.

-D 0 has at least one free dimension.

-∀i ∈ {1, 2, ..., n 1 }, |λ i | ≤ 1 ∧ λ i ̸ = -1.
-∀i ∈ {1, 2, ..., n 0 }, we have either

W i r ∞ = c i or W i r ∞ > c i .
We use I e to represent the maximum set of integers such that ∀i ∈ I e , W i r ∞ = c i and we use I n to represent the maximum set of integers such that ∀i ∈

I n , W i r ∞ > c i . -If λ 1 ̸ = 0 (we assume that λ 1 is unique if λ 1 ̸ = 0) and I e is not empty, then λ 1 is positive. -If λ 1 ̸ = 0, then ∀i ∈ I e , ∀j ∈ {2, ..., n 1 } , |W i v 1 α 1 | > n 1 |W i v j α j | (we ignore the case where ∃i ∈ I e , W i v 1 = 0). -If λ 1 ̸ = 0, then ∀i ∈ I n , max β∈{-1,1} n 1 ∥ n 1 j=1 β j α j v j ∥ 2 < W i r∞-c i ∥W i ∥ 2 .
The main idea of Theorem 6 is illustrated in Fig 4 .3 where the huge rectangle represents a discrete domain D which has two free dimensions and the zone surrounded by dashed lines represents the compatible zone S (which is a open set) of a certain cycle of discrete domains C T . Each dashed line l ci represents a linear constraint of the form w T x > c where w, x are vectors and c is a real number. The fact that a hybrid trajectory τ is attracted by C T is equivalent to the fact that the intersection points between τ and D, noted by the sequence (h 1 , h 2 , ...), always stay inside S and converge to (λ 1 ̸ = 0) or reach (λ 1 = 0) h ∞ , which belongs to the closure of S. This idea of using the intersection points between a hybrid trajectory and a hyperplan to study the properties of this hybrid trajectory is, like the previous chapter, also based on Poincaré map.

Whether h ∞ belongs to the closure of S or not can be easily verified by using these linear constraints. A necessary condition for this sequence to always satisfy these linear constraints is that the absolute values of all eigenvalues of the reduction matrix of C T are less than or equal to 1. In case that these eigenvalues satisfy this necessary condition, to verify if this sequence always satisfies these linear constraints, we separate these constraints on two classes: The first class contains all constraints which are not reached by h ∞ : l c2 , l c3 , l c4 , the second class contains all constraints which are reached by h ∞ : l c1 , l c5 . To verify if l c2 , l c3 , l c4 are always satisfied, we can verify if this sequence enters and stays in a circle centered by h ∞ which only contains hybrid states satisfying constraints l c2 , l c3 , l c4 (this is related to the condition:

If λ 1 ̸ = 0, then ∀i ∈ I n , max β∈{-1,1} n 1 ∥ n 1 j=1 β j α j v j ∥ 2 < W i r∞-c i ∥W i ∥ 2
). Such circle can always be found if it is sufficiently small, for example: the circle in Fig 4 .3. To verify if l c1 , l c5 are always satisfied, we can verify if this sequence is sufficiently "close" to v 1 which is the eigenvector related to the eigenvalue with the maximum absolute value among the eigenvalues that differ from 1 and which also "points into" S (this is related to the condition:

If λ 1 ̸ = 0, then ∀i ∈ I e , ∀j ∈ {2, ..., n 1 } , |W i v 1 α 1 | > n 1 |W i v j α j |).
Here, sufficiently "close" to v 1 means intuitively that the angle between ---→ h ∞ h i and v 1 is sufficiently small. The formal proof of Theorem 6 is given as follows.

Proof: The proof of the sufficiency part:

In the case that a hybrid trajectory τ is attracted by a cycle of discrete domains, then after certain moment, τ will always stay in the same cycle of discrete domains, noted by

C T = (D 0 , D 1 , D 2 , ..., D p , D 0 ).
We consider, at first, the case that there is one D i of C T such that D i has at least one free dimension (meaning that there exists one dimension such that the boundaries are not reached in this dimension), without loss of generality we assume that D 0 has at least one free dimension, so the first condition is satisfied.

A priori, we cannot ensure the second condition, because the decomposition of r(h 0 )r ∞ in the directions of certain eigenvectors of the reduction matrix of C T can be zero and these eigenvectors happen to be associated to eigenvalues of which the absolute values are greater than 1, however this is extremely rare. So in this chapter we ignore this case, and based on this assumption, we must have that for any eigenvalue λ i , |λ i | ≤ 1, otherwise we cannot ensure that τ always stays inside this cycle, because if there exists an eigenvalue such that |λ i | > 1, then r ∞ is infinitely far away from r(h 0 ). If there exists λ i = -1, we can double this cycle of discrete domains (the doubled cycle of discrete domains is (D 0 , D 1 , D 2 , ..., D p , D 0 , D 1 , ..., D p , D 0 )) to eliminate this eigenvalue. So the second condition is satisfied.

The intersection points of τ with D 0 will converge to or reach a certain hybrid state h ∞ = (π ∞ , d s 0 ) which belongs to the closure of the compatible zone. We have that r ∞ is the short version of π ∞ by only considering the dimensions where boundaries are not reached in D 0 . Based on Assumption 5, the reduction compatible zone can be described by {x | W x > c}. If h ∞ belongs to the compatible zone, then we have ∀i ∈ {1, 2, ..., n 0 }, W i r ∞ > c i . Else if h ∞ does not belong to the compatible zone but belongs to the closure of the compatible zone, then we have that for some i ∈ {1, 2, ..., n 0 }, W i r ∞ = c i . This means that some "boundaries" of the reduction compatible zone are reached, and for others i ∈ {1, 2, ..., n 0 }, we have W i r ∞ > c i . So the third condition is satisfied.

We note the intersection hybrid states of τ and D 0 as the sequence (h 0 , h 1 , h 2 , ...). We have that r(h n ) = r ∞ + n 1 i=1 λ n i α i v i where v i are eigenvectors of the reduction matrix of C T and r(h 0 ) -r ∞ = n 1 i=1 α i v i . By choosing λ 1 as the eigenvalue with maximum absolute eigenvalue among the eigenvalues that differ from 1, if λ 1 ̸ = 0, then we have:

lim n→∞ |λ n 1 α 1 | > > |λ n i α i | where i ̸ = 1, and lim n→∞ n 1 i=1 λ n i α i v i → 0.
This proves that the fifth and the sixth conditions are satisfied.

Suppose that λ 1 ̸ = 0. Since the sequence (h 0 , h 1 , h 2 , ...) is inside the compatible zone, if I e is not empty, then for any j ∈ I e , we should have ∀n ∈ N,

W j n 1 i=1 λ n i α i v i > 0.
When n is sufficiently big, for any j ∈ I e , the sign of

W j n 1 i=1 λ n i α i v i is dominated only by W j λ n 1 α 1 v 1 , meaning that W j n 1 i=1 λ n i α i v i and W j λ n 1 α 1 v 1 have the same sign, in fact, when n is sufficiently big, we have |W j λ n 1 α 1 v 1 | > n 1 i=2 |W j λ n i α i v i |, this means that adding n 1 i=2 W j λ n i α i v i to W j λ n 1 α 1 v 1 does not change the sign of W j λ n 1 α 1 v 1 . So the sign of λ 1 is
positive, otherwise we can not ensure that ∀j ∈ I e , ∀n ∈ N, W j n 1 i=1 λ n i α i v i > 0, which is equivalent to ∀j ∈ I e , ∀n ∈ N, W j λ n 1 α 1 v 1 > 0 when n is sufficiently big. This proves that the fourth condition is satisfied.

In the case that any D i of C T has no free dimension, we can easily prove the sufficiency part.

The proof of the necessity part:

If D 0 has no free dimension, then τ reaches a periodic hybrid trajectory which is a special case of τ being attracted by a cycle of discrete domains.

If these six conditions are satisfied, to prove that τ is attracted by a cycle of discrete domains, we only need to prove that τ always returns to the compatible zone of C T . This can be easily proved if λ 1 = 0. From now on, we suppose λ 1 ̸ = 0. We note that the hybrid trajectory from h 0 returns to D 0 at h 1 . We have that r(h 1 )-r ∞ = n 1 i=1 λ i α i v i . A sufficient condition for "h 0 belongs to the compatible zone" is (∀j

∈ I e , W j n 1 i=1 α i v i > 0) ∧ (∀j ∈ I n , max β i ∈{-1,1} ∥ n 1 i=1 β i α i v i ∥ 2 < W j r∞-c i ∥W j ∥ 2
) which is satisfied in this case. In fact, the first part of this sufficient condition is a necessary condition for h 0 to belong to the compatible zone and the second part is ensured by the sixth condition. Now we want to prove that this sufficient condition is also satisfied for h 1 . The fifth condition indicates that for any j ∈ I e the sign of W j n 1 i=1 α i v i is dominated by the sign of W j α 1 v 1 , so W j α 1 v 1 is also positive. Since ∀i ̸ = 1, we have either

|λ 1 | > |λ i | or |λ 1 | < |λ i | ∧ α i = 0, so we have ∀j ∈ I e , ∀i ∈ {2, ..., n 1 } , |W j v 1 λ 1 α 1 | > n 1 |W j v i λ i α i |.
This means that for any j ∈ I e the sign of W j n 1 i=1 α i λ i v i is also dominated by the sign of W j α 1 λ 1 v 1 . Since λ 1 is positive if I e is not empty, we have ∀j ∈ I e , W j n 1 i=1 λ i α i v i > 0, because the sign of W j n 1 i=1 λ i α i v i is same as the sign of W j α 1 λ 1 v 1 which is same as the sign of W j α 1 v 1 of which the sign is positive. Since ∀λ i , |λ i | ≤ 1, we have max

β i ∈{-1,1} ∥ n 1 i=1 β i λ i α i v i ∥ 2 ≤ max β i ∈{-1,1} ∥ n 1 i=1 β i α i v i ∥ 2 . So we also have ∀j ∈ I n , max β i ∈{-1,1} ∥ n 1 i=1 β i λ i α i v i ∥ 2 < W j r∞-c i ∥W j ∥ 2
. By now we proved that h 1 also belongs to the compatible zone. By mathematical induction, we can prove that τ always returns to the compatible zone of C T . So τ is attracted by a cycle of discrete domains.

□ Theorem 6 is used later in the algorithm to analyze reachability.

Chaotic Hybrid Trajectories

In this thesis, a hybrid trajectory of HGRN is called a chaotic hybrid trajectory if it does not reach a fixed point and it is not attracted by a cycle of discrete domains. So all hybrid trajectories which are not included in the previous two classes are considered chaotic hybrid trajectories. It is worth mentioning that the dynamics of chaotic hybrid trajectories, a priori, can be different from the chaotic dynamics of classic nonlinear dynamical systems. The reason why we still use the terminology "chaotic" is that a similar concept of chaos has been used in some pre-existing works of other hybrid systems [START_REF] Edwards | A calculus for relating the dynamics and structure of complex biological networks[END_REF][START_REF] Hamatani | Analysis of a 3-Dimensional Piecewise-Constant Chaos Generator Without Constraint[END_REF].

As stated before, the sequence of discrete domains reached by a hybrid trajectory of the previous two classes is of the form: (D 1 , D 2 , ..., D m ) or (D 1 , D 2 , ..., D k , (D k+1 , D k+2 , ..., D k+p ) n ). While a chaotic hybrid trajectory corresponds to an irregular sequence of discrete domains, like an irrational number with each digit replaced by a discrete domain.

To prove such chaotic hybrid trajectories exist, we have constructed a HGRN with chaotic hybrid trajectories based on a pre-existing model of electrical circuit with a chaotic attractor [START_REF] Hamatani | Analysis of a 3-Dimensional Piecewise-Constant Chaos Generator Without Constraint[END_REF].

In our work, we have not yet found a method to check reachability for chaotic hybrid trajectories, which, a priori, can be undecidable. So, in this subsection, we only introduce a method to predict whether a hybrid trajectory is chaotic, which is based on a necessary condition.

We prove that a chaotic hybrid trajectory has the following properties.

Property 1 For a chaotic hybrid trajectory τ , there exist t 0 and a finite set of discrete domains L D , such that after t 0 , τ cannot reach any discrete domain which does not belong to L D , and for any discrete domain D 0 ∈ L D , D 0 is reached by τ an infinite number of times.

Proof: This can be proved by the fact that the number of discrete domains is finite and the hybrid trajectory does not stay in a particular discrete domain. In fact, this property is true for any hybrid trajectory which is not a hybrid trajectory halting in finite time. □ Following Property 1, for any D 0 ∈ L D , τ will return to D 0 an infinite number of times, and each time it stays inside a sequence of discrete domains of the form (D 0 , ..., D 0 ) (only the first and the last discrete domain are D 0 ). The set of all such sequences of discrete domains is noted by L D 0 , and we make the following assumption.

Assumption 6 L D 0 is a finite set.
This Assumption is based on the fact that the number of discrete domains is limited and the dynamics in the discrete states is simple (a constant vector).

Property 2 Following the notions of Property 1, for any D 0 ∈ L D , we can find t 1 > t 0 such that τ reaches D 0 at t 1 and there exists a finite set of sequences of discrete domains L T which satisfies the following properties:

-∀T ∈ L T , T is of the form (D 0 , ..., D 0 ).

-From t 1 , τ will return to D 0 an infinite number of times, and each time it stays in a sequence of discrete domains which belongs to L T . More formally, if τ reaches D 0 at t α , then from t α , it will return to D 0 at t β , and the restriction of τ on [t α , t β ] is inside a sequence of discrete domains of L T . -From t 1 , ∀T ∈ L T , T is crossed by τ an infinite number of times.

-L T contains at least two sequences of discrete domains.

Proof: For the first two properties, according to Property 1, from t 1 , τ will return to D 0 an infinite number of times, and each time it stays in a sequence of discrete domains of the form (D 0 , ..., D 0 ), all these sequence constitute L T . And according to Assumption 6, L T contains a finite number of sequences.

For the third property, for any T ∈ L T , if it is crossed by τ a finite number of times, then we can increase the value of t 1 such that from t 1 , T is not crossed by τ any more. This result is also based on the fact that L T is a finite set which is a result of Assumption 6.

For the fourth property, if there is only one sequence of discrete domains, then after certain moment, τ will always stay in the same cycle of discrete domains. In this case, τ is attracted by a cycle of discrete domains and is not a chaotic hybrid trajectory.

□ Following the notations of Property 2, the sequence of discrete domains crossed by τ from t 1 can be described by the sequence (T 1 , T 2 , T 3 , ...), where ∀i ∈ N, T i ∈ L T . This sequence has the following properties.

Property 3 ∃i ∈ N, such that T i ̸ = T i+1 .
Proof: This is a direct result of the fact that all elements of L T must appear in the sequence (T 1 , T 2 , T 3 , ...) and L T has at least two elements. □ Property 4 ∀i ∈ N, ∃k ∈ N, such that T i = T i+k .

Proof: If it is not true, then it means that T i is crossed by τ for finite times, which contradicts with Property 2. □ By combining Property 3 and Property 4, we get a necessary condition for τ to be a chaotic hybrid trajectory:

Property 5 ∃i ∈ N, ∃k ∈ N, k ̸ = 1, such that T i ̸ = T i+1 and T i = T i+k .
This condition is used in the following section to predict whether a hybrid trajectory is a chaotic hybrid trajectory.

Reachability Analysis Algorithm

In this section, we present our reachability analysis algorithm. The main algorithm is presented in Algorithm 1. Its main idea is that we continue simulating the hybrid trajectory τ from h ini until reaching one of the following situations:

τ reaches R target . In this case, the algorithm stops and returns "R target is reached".

-At some time t 0 , τ is already attracted by a cycle of discrete domains and the periodic hybrid trajectory, which attracts τ , reaches R target . In this case, the algorithm stops and returns "R target is reached". -At some time t 0 , τ is already attracted by a cycle of discrete domains and after t 0 there is no more chance to reach R target . In this case, the algorithm stops and returns "R target is not reached". -Property 5 is satisfied, therefore τ is probably chaotic. In this case, the algorithm stops and returns "unknown result". More details about Algorithm 1 are discussed as follows.

To

check if a transition h → h ′ reaches R target = (π, d starget ) | π i ∈ [a i , b i ], i ∈ {1, 2
, ..., N } , we need to consider two cases: Firstly, if h → h ′ is an instant transition, then we only need to verify if h ′ belongs to R target ; secondly, if h → h ′ is a non-instant transition, meaning that h → h ′ is an N -dimensional straight line, then we only need to verify if h → h ′ crosses a boundary of R target . To give an example of the second case, if a boundary of

R target is (π, d starget ) | π i 0 = a i 0 , π i ∈ [a i , b i ], i ∈ {1, 2, ..., N } \ {i 0 } , we only need to ver- ify if h → h ′ contains a hybrid state h t = (π t , d starget ) such that π i 0 t = a i 0 and ∀i ∈ {1, 2, ..., N } \ {i 0 } , π i t ∈ [a i , b i ]. To determine if
(D 0 , D 1 , D 2 , ..., D p , D 0 ) Input 3: A region R target = (π, d starget ) | π i ∈ [a i , b i ], i ∈ {1, 2, ..., N }
Output (three possibilities): "Yes", "No", and "Reached" return "No" 18: end if cessor in the same discrete state of R target as P re ds target (R target ) = i∈{1,2,...,q} Z i . We note r(Z i ) as the set of the reduction vectors of all hybrid states of Z i . This Non-intersection condition (its formal definition will be given later) is a sufficient condition for that ∀j ∈ {1, 2, ..., q} , if Z j ⊂ D i (D i is a discrete domain that belongs to the cycle of discrete domains attracting the hybrid trajectory), then

1: Compute the hybrid states h ∞ 0 , h ∞ 1 , h ∞ 2 , ...,
x ∈ S r (D i ) | ∥x -r(h ∞ i )∥ 2 ≤ max β∈{-1,1} n 1 ∥ n 1 k=1 β k α k v k ∥ 2 ∩ r(Z j ) = ∅ where S r (D i ) is the reduction compatible zone of the cycle of discrete domains (D i , D i+1 , ..., D p , D 0 , ..., D i ), n 1 is the dimension of S r (D i ), r(h i ) -r(h ∞ i ) = n 1 k=1 α k v k , v
., q} , if Z j ⊂ D i , then x ∈ S r (D i ) | ∥x -r(h ∞ i )∥ 2 ≤ max β∈{-1,1} n 1 ∥ n 1 k=1 β k α k v k ∥ 2 ∩ r(Z j ) = ∅
. This is equivalent to say that the black box does not intersect the blue boxes in Fig 4 .4 left. The reason why we want to verify this condition, noted here as Condition1, is that once it is satisfied, there is no more chance that a hybrid trajectory can reach R target directly from D i (here "directly" means that it does not reach another discrete state before reaching R target ). However, it is complicated to compute directly the intersection of these sets. So, firstly we overestimate the set of all hybrid states in this discrete state which can be reached by hybrid trajectories from h ∈ S(D The Non-intersection condition is defined formally as follows. The hybrid trajectory from h ∞ i can be represented by

i ) | ∥r(h) -r(h ∞ i )∥ 2 ≤ max β∈{-1,1} n 1 ∥ n 1 k=1 β k α k v k ∥
h ∞ i = h ∞ i,0 → h ∞ i,1 → ... → h ∞ i,K where h ∞ i,j , j ∈ {1, 2
, ..., K} represents all hybrid states at which this hybrid trajectory reaches new boundaries (similar to the turning states in Chapter 3) and h ∞ i,K represents the hybrid state at which this hybrid trajectory reaches an output boundary of d starget for the first time. We note u j , j ∈ {1, 2, ..., K} as the left derivative of hybrid state h ∞ i,j . We use w j ∈ {0, 1} N , j ∈ {1, 2, ..., K} to describe the dimension where the new boundary is reached in each transition, for the transition h ∞ i,j-1 → h ∞ i,j , j ∈ {1, 2, ..., K}: if the new boundary is reached in the m j th dimension, then w m j j = 1 and w k j = 0, k ∈ {1, 2, ..., N } \ {m j }. Based on the notations of u j and w j , we compute a sequence l j , j ∈ {0, 1, 2, ..., K -1, K} as follows:

l 0 = max β∈{-1,1} n 1 ∥ n 1 k=1 β k α k v k ∥ 2 and ∀j ∈ {1, 2, ..., K -1, K} , l j = l j-1 | u j •w j ∥u j ∥ 2 |
where u j • w j is the inner product of u j and w j . The Non-intersection condition is satisfied if the hybrid trajectory from h ∞ i,0 does not reach directly the region

R ′ target = (π, d starget ) | π i ∈ [max(a i -l K , 0), min(b i + l K , 1)
], i ∈ {1, 2, ..., N } . The reason why we compute this sequence l j is illustrated in Fig 4 .5 by a 2-dimensional example. Intuitively, when a bunch of hybrid trajectories with radius l j-1 reaches a new boundary, the radius can increase. We can see that the sequence l j is increasing. To simplify the computation, we choose the maximum value l K to overestimate the hybrid states that can be reached by these hybrid trajectories. This choice of the value l K is only for overestimation, it is possible that a smaller value also works.

Here we prove that the Non-intersection condition is a sufficient condition so that ∀j ∈ {1, 2, ..., q} , if Z j ⊂ D i , then x ∈ S r (D i ) | ∥x -r(h ∞ i )∥ 2 ≤ max β∈{-1,1} n 1 ∥ n 1 k=1 β k α k v k ∥ 2 ∩ r(Z j ) = ∅. We note this hybrid trajectory h ∞ i = h ∞ i,0 → h ∞ i,1 → ... → h ∞ i,K which begins from h ∞
i,0 and finally reaches h ∞ i,K as τ i . If the Non-intersection condition is satisfied, then for any hybrid state h x which belongs to R target and for any hybrid state h y which belongs to τ i , the distance between h x and h y is greater than l K . On the other hand, for any hybrid trajectory τ x which is inside this discrete state and begins from the

region h ∈ S(D i ) | ∥r(h) -r(h ∞ i )∥ 2 ≤ max β∈{-1,1} n 1 ∥ n 1 k=1 β k α k v k ∥ 2
where S(D i ) is the compatible zone of this cycle, for any hybrid state h ′

x on τ x , we can find a hybrid state h ′ y on τ i , such that the distance between h ′ x and h ′ y is smaller than l K . If τ x crosses R target , then there exists a hybrid state h ′′ y which is on τ x and also belongs to R target , such that ∃h ′′

x on τ i , the distance between h ′′ x and h ′′ y is smaller than l K which is not possible as h ′′ y belongs to R target . This proves that ∀j ∈ {1, 2, ..., q} , if

Z j ⊂ D i , then x ∈ S r (D i ) | ∥x -r(h ∞ i )∥ 2 ≤ max β∈{-1,1} n 1 ∥ n 1 k=1 β k α k v k ∥ 2 ∩ r(Z j ) = ∅.
In fact, if the hybrid trajectory from h ini does not reach R target and h i is sufficiently close to h ∞ i , then the Non-intersection condition must be satisfied, because L K is sufficiently close to 0 in this case. Based on this fact, it can be proved that if the hybrid trajectory from h ini is attracted by a cycle of discrete domains and does not reach R target , then the Non-intersection condition must be verified in finite time for all discrete domains which belong to an input boundary of the discrete state of R target and belong to the attractive cycle of discrete domains (the cycle of discrete domains which attracts the hybrid trajectory from h ini ). In this case, this algorithm returns "Yes", meaning that the hybrid trajectory from h ini has no more chance to reach R target .

It can be proved that Algorithm 1 always stops in finite time. Firstly, if the hybrid trajectory from h ini is a hybrid trajectory halting in finite time, then the algorithm stops after a finite number of transitions. Secondly, if the hybrid trajectory is a chaotic hybrid trajectory, then Property 5 will be satisfied after a finite number of transitions, and once it is satisfied, the algorithm stops. Thirdly, if the hybrid trajectory is attracted by a cycle of discrete domains, then there are three cases: 1. The hybrid trajectory reaches R target in finite time; 2. The hybrid trajectory reaches R target after an infinite number of transitions; 3. The hybrid trajectory does not reach R target . We assume here that Property 5 is not satisfied before the hybrid trajectory reaching the attractive cycle of discrete domains. For case 1, the algorithm stops in finite time, as the hybrid trajectory will eventually reach R target . For case 2, the function Stop_condition returns "Reached" in finite time. For case 3, the function Stop_condition returns "Yes" in finite time. We need to mention that, since Property 5 is a necessary condition for the hybrid trajectory from h ini to be a chaotic hybrid trajectory, a priori, it could be true even if the hybrid trajectory belongs to other classes; however, empirically, we never identified such case, which should be rare. Yet, an "unknown result" should be taken as a possible missed formal positive or negative case.

Application: Estimation of Basins of Attraction

In this section, we present a potential application of this reachability analysis method which is to estimate the basin of attraction of attractors. An attractor can be a stable fixed point, a stable limit cycle or a chaotic attractor. The basin of attraction of an attractor is the set of all hybrid states such that hybrid trajectories from these hybrid states converge to or reach this attractor. Here, estimating the basin of attraction means estimating the percentage of hybrid states which belong to the basin of attraction in a given region. A concrete example of application is presented in follows.

In this example, we consider a pre-existing 5-dimensional HGRN of cell cycle from [START_REF] Behaegel | A hybrid model of cell cycle in mammals[END_REF][START_REF] Sun | Limit cycle analysis of a class of hybrid gene regulatory networks[END_REF]. This HGRN has a stable limit cycle and a stable fixed point. We apply our reachability analysis method on this HGRN to estimate, for each discrete state d s , the percentage of hybrid states which belong to the basin of attraction of this limit cycle, noted by P LC (d s ), and the percentage of hybrid states which belong to the basin of attraction of this stable fixed point, P F P (d s ).

To estimate these, we define two regions: R LC target and R F P target , such that for any hybrid trajectory, whether it reaches or converges to the stable limit cycle or the stable fixed point is equivalent to whether it reaches R LC target or R F P target respectively. These regions are found manually and defined as:

R LC target = {(π, (0, 0, 0, 0, 1)) | π 1 ∈ [0.82, 0.84], π 2,3,4 ∈ [0, 0.01], π 5 ∈ [0.99, 1]} R F P target = {(π, (2, 1, 0, 1, 0)) | π 1,2,3,4,5 ∈ [0, 1]}. For R LC
target , it can be easily verified that it contains hybrid states from the stable limit cycle and we can also prove that any hybrid trajectory which reaches this region must reach this stable limit cycle. R F P target covers the whole discrete state which contains the stable fixed point, so reaching the stable fixed point is equivalent to reaching R F P target . Then, in each discrete state d s , we randomly generate n s hybrid states by uniform distribution and we check whether these hybrid states reach R LC target or R F P target by using our reachability analysis method. We note n LC (d s ) the number of generated hybrid states that reach R LC target and we note n F P (d s ) the number of generated hybrid states that reach R F P target . P LC (d s ) and P F P (d s ) can be estimated by P LC (d s ) = n LC ns and P F P (d s ) = n F P ns . In our experimentation, we choose n s = 20. A higher value of n s can increase the precision of estimation, but it also takes more time. This setting takes nearly 30 minutes on a standard laptop computer 1 . The estimation results of P LC and P F P are illustrated in the transition graph of discrete states in Fig 4 .6. For each node in Fig 4 .6, the two floating point numbers represent estimated P LC and P F P . For example, for discrete state 20000, P LC (20000) = 0.75 and P F P (20000) = 0.25. The cycle of discrete states which contains the stable limit cycle is illustrated by the red arrows and the discrete state which contains the stable fixed point is illustrated by a red node. These results could provide useful information if we want to control a gene regulatory network such that the hybrid trajectory moves from one attractor to another.

Summary

In this chapter, we propose a reachability analysis method for HGRNs. In the first part of this work, we classify hybrid trajectories of HGRNs into different classes: Hybrid trajectories halting in finite time, hybrid trajectories attracted by cycles of discrete domains and chaotic hybrid trajectories, and provide some theoretical results about these hybrid trajectories regarding the reachability problem. Then, based on these theoretical results, we provide the first reachability analysis algorithm for HGRNs. Based on a first implementation, we apply this method on a 5-dimensional HGRN of cell cycle to estimate the basin of attraction of a stable limit cycle and a stable fixed point.

Since we use Property 5 (a necessary condition for chaotic hybrid trajectories) on line 14 in Algorithm 1, the algorithm might return inconclusive results ("unknown result") even in the cases that are decidable (non-chaotic hybrid trajectories). In fact, among HGRNs of gene regulatory networks, the cases that satisfy Property 5 are likely to be 0. very rare: There is no identified HGRN of gene regulatory network with either chaos or non-chaotic hybrid trajectory that satisfies Property 5. The example with chaos in this chapter is not a biological model but comes from physics. For now, this algorithm can be considered sufficient for checking reachability in practice. However, it could be interesting to further investigate the decidability of chaotic hybrid trajectories and the condition for a hybrid trajectory to be chaotic.

This chapter and the previous one focus on the general analysis methods of HGRNs. These methods are very adapted to the study of oscillations in some specific networks. The following chapter focuses on this topic.

Chapter 5

CONDITION FOR SUSTAINED OSCILLATION IN CANONICAL REPRESSILATOR

The better control of synthetic oscillations in repressilator is an important question in synthetic biology. Most previous works use ordinary differential equations to study repressilators. In this chapter, we explore the possibility to use HGRN to study oscillations in the canonical repressialtor and aim to find conditions for the existence of sustained oscillations described by separable constraints on parameters. The methods in this chapter rely on the work in Chapter 3. Two major results are the following: 1) we develop, by using the Poincaré map, a sufficient and necessary condition for the existence of sustained oscillations (Section 5.4); 2) based on this condition, we give a method using the range enclosure property of Bernstein coefficients to compute compatible separable constraints (Section 5.5). The work of this chapter follows a 3-week research stay at the I3S laboratory (Laboratoire d'Informatique, Signaux et Systèmes de Sophia Antipolis) in Nice in June 2022 and has been presented at conference BIOINFORMATICS 2023 as regular paper [START_REF] Sun | Condition for Sustained Oscillations in Repressilator Based on a Hybrid Modeling of Gene Regulatory Networks[END_REF].

Introduction of Canonical Repressilator

Repressilators represent gene regulatory networks consisting of at least one feedback loop, in which the expression of each gene inhibits the expression of the next gene. In this chapter, we focus on a specific repressilator called canonical repressilator, which is a network with three genes having a unique feedback loop with only inhibitions and not having any other regulation. The influence graph of the canonical repressilator is illustrated in Fig 5 .1.

The canonical repressilator is widely studied in synthetic biology because it could generate synthetic sustained oscillations, which could be the treatment of diseases related to circadian rhythms, for instance, by allowing drug delivery at a particular pace. The first biological implementation of the canonical repressilator with sustained oscillations is realized in [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF] by using three natural repressor proteins, the TetR, LacI and CI repressors and the stability of these oscillations is improved later in other works [START_REF] Potvin-Trottier | Synchronous longterm oscillations in a synthetic gene circuit[END_REF].

In fact, how to biologically implement a canonical repressilator with sustained oscillations is still an open question, particularly in eukaryotic cells, therefore exploring new models to search for conditions for sustained oscillations is of high interest.

Previous works about mathematical analysis of oscillations in the canonical repressilator are mainly based on differential equations. Many models of the canonical repressilator with three components are developed from a differential equation model using 6 variables [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF], with 3 variables for repressor-protein concentrations and 3 variables for corresponding mRNA concentrations. This 6-variable model can also be reduced to a 3-variable model with only repressor-protein variables under certain assumptions [START_REF] Buşe | Dynamical properties of the repressilator model[END_REF]. These models and their variations are extensively studied in the literature [START_REF] Dukarić | On three genetic repressilator topologies[END_REF][START_REF] Dilão | The regulation of gene expression in eukaryotes: bistability and oscillations in repressilator models[END_REF][START_REF] Kuznetsov | Heteroclinic cycles in the repressilator model[END_REF][START_REF] Buşe | Existence of limit cycles in the repressilator equations[END_REF][START_REF] Müller | A generalized model of the repressilator[END_REF][START_REF] El Samad | Repressilators and promotilators: Loop dynamics in synthetic gene networks in Proceedings of the 2005[END_REF]. One major limit about differential equations is that some dynamical properties are hard to analyze. So in this chapter, we use HGRNs to study oscillations in the canonical repressilator, as the dynamical properties of HGRNs are easier to analyse.

A HGRN of Canonical Repressilator

In this chapter, we assume that each gene has one threshold when it influences another gene. Although, a priori, one gene can have multiple thresholds for another gene, here we only consider the simplest case. Based on this assumption, since each gene only influences one other gene in the canonical repressilator, it has only two discrete levels separated by only one threshold.

The parameters (celerities) of this HGRN of the canonical repressilator are shown symbolically in Table 5.1. Each parameter in Table 5.1 is strictly positive and is denoted by C xyizj , which represents the absolute value of the celerity of variable x when the discrete level of variable y is i and the discrete level of variable z is j. Consider the influence of A on B: When the discrete level of A is 1, meaning that the expression of A is above the threshold to inhibit B, then the time derivative of B is always negative, no matter the discrete level of B (0 or 1) which corresponds to the negative values -C ba1b0 and -C ba1b1 . On the other hand, when the expression of A is below the threshold to inhibit B, the time derivative of B is always positive, corresponding to parameters C ba0b0 and C ba0b1 . From the parameters in this table, we can also see that the number of different parameters [START_REF] Comet | A formal model for gene regulatory networks with time delays in International Conference on Computational Systems[END_REF] is smaller than the multiplication of the number of dimensions by the number of discrete states [START_REF] Sun | Condition for Sustained Oscillations in Repressilator Based on a Hybrid Modeling of Gene Regulatory Networks[END_REF], because some discrete states have celerities in common (same regulation on some variables).

In addition to the threshold which separates the discrete levels 0 and 1, each gene also has a maximal value and a minimal value. For example, when A inhibits B (see Figure 5.1), B will continue to decrease until it reaches the minimal value (most of the time this minimal value is 0) which is related to the lower boundary in the second dimension (the dimension of gene B) of discrete state 10 * , where * can be 0 or 1. Similarly, when A does not inhibit B, B will continue to increase until it reaches the upper boundary in the second dimension of 01 * .

Table 5.1 -Parameters of the HGRN of the canonical repressilator. A B C C

A C B C C 0 0 0 C ac0a0 C ba0b0 C cb0c0 0 0 1 -C ac1a0 C ba0b0 C cb0c1 0 1 0 C ac0a0 C ba0b1 -C cb1c0 0 1 1 -C ac1a0 C ba0b1 -C cb1c1 1 0 0 C ac0a1 -C ba1b0 C cb0c0 1 0 1 -C ac1a1 -C ba1b0 C cb0c1 1 1 0 C ac0a1 -C ba1b1 -C cb1c0 1 1 1 -C ac1a1 -C ba1b1 -C cb1c1

Qualitative Behaviors in this HGRN of Canonical

Repressialtor

In this section, we discuss different qualitative properties of this HGRN. To analyze dynamical properties of a HGRN, we need to firstly analyze the transition graph of discrete states, which is determined by the signs of celerities. The signs of celerities of this HGRN of the repressilator can be found in Table 5.1, based on which the transition graph of discrete states can be constructed using the classical discrete asycnhronous semantics, see Figure 5.3.

From the transition graph, we can see that there is a unique cycle of discrete states, which is 001 -→ 011 -→ 010 -→ 110 -→ 100 -→ 101 -→ 001. This cycle is a global attractor of discrete states (the unique terminal strongly connected component), which means that any hybrid trajectory in this HGRN will finally enter this cycle. Following the same sequence of discrete states, there is a special hybrid trajectory: ((1, 1, 0), (0, 0, 1)) -→ ((1, 0, 0), (0, 1, 1)) -→ ((1, 0, 1), (0, 1, 0)) -→ ((0, 0, 1), (1, 1, 0)) -→ ((0, 1, 1), (1, 0, 0)) -→ ((0, 1, 0), (1, 0, 1)) -→ ((1, 1, 0), (0, 0, 1)). It is illustrated in Figure 5.4 by green arrows. This hybrid trajectory contains 6 different hybrid states (the last and the first one are identical). From each hybrid state in this hybrid trajectory, there is an instant transition (transition which crosses boundaries between discrete states and takes no time), which reaches the next hybrid state. This hybrid trajectory represents also a closed hybrid trajectory, meaning that starting from each of these 6 hybrid states, the hybrid trajectory will return to the initial hybrid state. When representing hybrid trajectories, we often use an embedding of hybrid states in R N : A hybrid state (π, d s ) is represented in R N by summing its discrete and fractional parts: π + d s . We say that the hybrid state (π, d s ) and the point π +d s ∈ R N are related. When doing so, the six previous hybrid states are embedded in the same point H 0 in R 3 : H 0 = (1.0, 1.0, 1.0). H 0 is called a characteristic state of this HGRN. The characteristic state is formally defined as follows: Definition 14 (Characteristic state) A characteristic state of a HGRN is a state H in Euclidean space such that: For any hybrid state h 0 related to H, all hybrid trajectories from h 0 will never reach a hybrid state which is not related to H, and there exist oscillations in any small neighborhood of H.

In this chapter, we say that a hybrid trajectory of a HGRN is an oscillation if it is related to an oscillation in Euclidean space. Likewise, the nature of an oscillation of a HGRN (damped or sustained) and the relation between an oscillation and a characteristic state (converging to the characteristic state, moving away from it, etc.) are determined by the related oscillation in Euclidean space. A neighborhood of a characteristic state H is defined as a set N (H) = {(π, d s ) ∈ E h | ∥π + d s -H∥ < r} where r ∈ R is the radius of the neighborhood. In Euclidean space, a characteristic state is a fixed point, because from any hybrid state h related to a characteristic state, all hybrid trajectories can only reach hybrid states related to this characteristic state.

We can easily prove that H 0 is a characteristic state: Apart from the 6 hybrid states on this closed hybrid trajectory, there are two other hybrid states which are related to H 0 : ((1, 1, 1), (0, 0, 0)) and ((0, 0, 0), (1, 1, 1)) From each of these two other hybrid states, all hybrid trajectories reach directly the closed hybrid trajectory, and once they do, they can never leave it. Finally, we can find oscillations which follow the unique cycle of discrete states in any neighborhood of H 0 . In this HGRN, there is only one characteristic state. This can be proved by verifying all "corners" of discrete states.

All hybrid trajectories in this HGRN will oscillate in this unique cycle of discrete states, except some special hybrid trajectories from discrete state 000 or 111 which can reach directly a hybrid state related to the characteristic state, see for example the black arrows in Figure 5.4. These oscillations can have several different dynamical properties.

To better illustrate the possible dynamical properties in this repressilator, different qualitative behaviors of a HGRN of negative feedback loop in 2 dimensions are illustrated in Figure 5.5, where black arrows represent celerities in each discrete state and red arrows represent some hybrid trajectories. The three figures are obtained by choosing three different parameterisations, which represent the three qualitative behaviors of this HGRN. The closed hybrid trajectory with hybrid states P, M, N, Q contains only instant transitions. The state in Euclidean space that is related to states P, M, N, Q is a characteristic state, and is the only characteristic state in this HGRN. On the left, the characteristic state is stable as all hybrid trajectories tend to converge to it. In this case, there is no sustained oscillation. On the right, the characteristic state is unstable as all hybrid trajectories from a small neighborhood of the characteristic state will move away from it and will finally reach a stable limit cycle containing at least one sliding mode. Between these two cases, in the middle, there is a third possibility in which all hybrid trajectories circle around the characteristic state without getting closer or moving away, which we call parallel cycles and can be considered as a special case of sustained oscillations.

The HGRN of the repressilator in 3 dimensions (Figure 5.1) is more complicated than the example HGRN of negative feedback loop in 2 dimensions (Figure 5.5), but we think that it also has three similar qualitative behaviors, since each discrete state in the only cycle of discrete states has only one successor (Figure 5.3). However, we have no proof for the non-existence of other possibilities yet, for example chaos or the co-existence of sustained oscillations and a stable fixed point. So in this chapter, we make Hypothesis 1. 5.1), either the characteristic state is stable and all oscillations are damped, or the characteristic state is unstable and all oscillations are sustained.

Hypothesis 1 In this HGRN of canonical repressilator (see Figure 5.1 and Table

In this HGRN of the canonical repressilator, the characteristic state is said stable if we can find a small neighborhood around the characteristic state such that all oscillations which begin from this neighborhood converge to the characteristic state, and the characteristic state is said unstable if there is no damped oscillation which converges to it. Now, based on Hypothesis 1, we can use the stability of the characteristic state to determine the existence of sustained oscillations in this HGRN of the canonical repressilator: All oscillations are sustained if and only if the characteristic state is unstable. This is the main idea used in this chapter to find conditions for the existence of sustained oscillations. A similar idea was used in other works with differential models, see for example [START_REF] El Samad | Repressilators and promotilators: Loop dynamics in synthetic gene networks in Proceedings of the 2005[END_REF][START_REF] Page | Degradation rate uniformity determines success of oscillations in repressive feedback regulatory networks[END_REF][START_REF] Wang | Construction of genetic oscillators with interlocked feedback networks[END_REF].

The problem now is how to analyze the stability of the characteristic state. To do so, we apply a method based on the Poincaré map.

A Sufficient and Necessary Condition for Sustained Oscillation

In this section, firstly we present in details how we compute symbolically a Poincaré map that determines the stability of the characteristic state and then we develop a sufficient and necessary condition for sustained oscillation by symbolic eigenanalysis of this Poincaré map. The idea of the method to compute the Poincaré map in this chapter is same as the computation of the transition matrix in Chapter 3.

Symbolic Computation of Poincaré Map

In order to analyze the stability of the characteristic state, we only need to consider hybrid trajectories without sliding mode around this state. Indeed, the characteristic state is stable if we can find a small neighborhood around the characteristic state such that all oscillations which begin from this neighborhood converge to the characteristic state. Such a hybrid trajectory exists because otherwise, the celerities would prevent the characteristic state from existing. Without loss of generality, we choose the lower boundary in the second dimension of discrete state 011 as Poincaré section; see the blue boundary in Figure 5.4. Now, we consider any hybrid trajectory τ which begins from a hybrid state s 1 = ((x 1 , 0, z 1 ), (0, 1, 1)) on the Poincaré section and returns to the Poincaré section for the first time at s 7 = ((x 7 , 0, z 7 ), (0, 1, 1)) without sliding mode; such a hybrid trajectory is illustrated in red in Figure 5.4. Thus, the Poincaré map is an affine application describing the relation between (x 1 , 0, z 1 ) and (x 7 , 0, z 7 ). Since s 1 is on an input boundary of the discrete state 011, from s 1 , τ will first (continuously) cross the discrete state 011 and reach a hybrid state ((x 2 , y 2 , 0), (0, 1, 1)) on the lower boundary in the third dimension of 011 which is the output boundary of this discrete state towards 010. We name this output boundary e 1 . Then, it crosses instantly e 1 and reaches an input boundary of 010 in hybrid state ((x 2 , y 2 , 1), (0, 1, 0)). The duration of crossing in discrete state 011 is:

t 1 = 0 -z 1 -C cb1c1 (5.1)
It should be noted that, to ensure that this hybrid trajectory τ has no sliding mode in 011, we also need to ensure that the lower boundary in the first dimension and the upper boundary in the second dimension of 011 are not reached before e 1 , which gives us two additional inequalities:

t 1 < 0 -x 1 -C ac1a0 (5.2) t 1 < 1 -0 C ba0b1 (5.3)
In fact, these inequalities can always be satisfied if τ is sufficiently close to the characteristic state, which is the case we consider here. Therefore, in the rest of this section, we do not consider these additional constraints.

Based on the duration t 1 and the fact that there is no sliding mode, we can get the relation between (x 1 , 0, z 1 ) and (x 2 , y 2 , 1):

x 2 = x 1 -C ac1a0 × t 1
(5.4)

y 2 = 0 + C ba0b1 × t 1 (5.5)
Following the same process, we can get the duration of τ in each discrete state and the relations between hybrid states from one input boundary to another input boundary:

t 2 = 1 -x 2 C ac0a0 t 3 = 0 -y 3 -C ba1b1 t 4 = 1 -z 4 C cb0c0 (5.6) t 5 = 0 -x 5 -C ac1a1 t 6 = 1 -y 6 C ba0b0
(5.7)

y 3 = y 2 + C ba0b1 × t 2 z 3 = 1 -C cb1c0 × t 2 (5.8
)

x 4 = 0 + C ac0a1 × t 3 z 4 = z 3 -C cb1c0 × t 3
(5.9)

x 5 = x 4 + C ac0a1 × t 4 y 5 = 1 -C ba1b0 × t 4
(5.10)

y 6 = y 5 -C ba1b0 × t 5 z 6 = 0 + C cb0c1 × t 5
(5.11)

x 7 = 1 -C ac1a0 × t 6 z 7 = z 6 + C cb0c1 × t 6 (5.12)
where t 2 , t 3 , t 4 , t 5 , t 6 are the durations of τ in discrete states 010, 110, 100, 101, 001 respectively, and (0, y 3 , z 3 ), (x 4 , 1, z 4 ), (x 5 , y 5 , 0), (1, y 6 , z 6 ), (x 7 , 0, z 7 ) are the fractional parts of the hybrid states when τ first reaches 110, 100, 101, 001, 011 respectively.

Based on the above equations, the Poincaré map can be calculated to describe the relation between (x 1 , 0, z 1 ) and (x 7 , 0, z 7 ) as follows. One dimension is missing in the matrix below; indeed, this dimension is useless in the computation of the stability, for more details see Chapter 3 Section 3.3.2.

  x 7 z 7   =   b 1 c 1 b 2 c 2     x 1 z 1   +   a 1 a 2   (5.13)
In the above equation, a 1 , a 2 , b 1 , b 2 , c 1 , c 2 are nonlinear combinations of the celerity values. Their expressions are given in the follows. We can easily derive that b 1 and c 2 are strictly positive, while b 2 and c 1 are strictly negative.

b 1 = B 1 D 0 b 2 = B 2 D 0 c 1 = C 1 D 0 c 2 = C 2 D 0 (5.14) B 1 = C ac0a1 C ac1a0 C ba0b1 C ba1b0 C cb0c0 C cb1c1 + C ac0a1 C ac1a0 C ba0b1 C ba1b0 C cb1c0 C cb1c1 + C ac0a1 C ac1a0 C ba1b0 C ba1b1 C cb1c0 C cb1c1 + C ac1a0 C ac1a1 C ba0b1 C ba1b0 C cb1c0 C cb1c1 + C ac1a0 C ac1a1 C ba1b0 C ba1b1 C cb1c0 C cb1c1 (5.15) B 2 = -C ac0a1 C ba0b0 C ba0b1 C cb0c0 C cb0c1 C cb1c1 -C ac0a1 C ba0b0 C ba0b1 C cb0c1 C cb1c0 C cb1c1 -C ac0a1 C ba0b0 C ba1b1 C cb0c1 C cb1c0 C cb1c1 -C ac0a1 C ba0b1 C ba1b0 C cb0c0 C cb0c1 C cb1c1 -C ac0a1 C ba0b1 C ba1b0 C cb0c1 C cb1c0 C cb1c1 -C ac0a1 C ba1b0 C ba1b1 C cb0c1 C cb1c0 C cb1c1 -C ac1a1 C ba0b1 C ba1b0 C cb0c1 C cb1c0 C cb1c1 -C ac1a1 C ba1b0 C ba1b1 C cb0c1 C cb1c0 C cb1c1
(5. 16)

C 1 = -C ac0a0 C ac0a1 C ac1a0 C ba0b1 C ba1b0 C cb0c0 -C ac0a0 C ac0a1 C ac1a0 C ba0b1 C ba1b0 C cb1c0 -C ac0a0 C ac1a0 C ac1a1 C ba0b1 C ba1b0 C cb1c0 -C ac0a1 C 2 ac1a0 C ba0b1 C ba1b0 C cb0c0 -C ac0a1 C 2 ac1a0 C ba0b1 C ba1b0 C cb1c0 -C ac0a1 C 2 ac1a0 C ba1b0 C ba1b1 C cb1c0 -C 2 ac1a0 C ac1a1 C ba0b1 C ba1b0 C cb1c0 -C 2 ac1a0 C ac1a1 C ba1b0 C ba1b1 C cb1c0
(5.17)

C 2 = C ac0a0 C ac0a1 C ba0b0 C ba0b1 C cb0c0 C cb0c1 + C ac0a0 C ac0a1 C ba0b0 C ba0b1 C cb0c1 C cb1c0 + C ac0a0 C ac0a1 C ba0b1 C ba1b0 C cb0c0 C cb0c1 + C ac0a0 C ac0a1 C ba0b1 C ba1b0 C cb0c1 C cb1c0 + C ac0a0 C ac1a1 C ba0b1 C ba1b0 C cb0c1 C cb1c0 + C ac0a1 C ac1a0 C ba0b0 C ba0b1 C cb0c0 C cb0c1 + C ac0a1 C ac1a0 C ba0b0 C ba0b1 C cb0c1 C cb1c0 + C ac0a1 C ac1a0 C ba0b0 C ba1b1 C cb0c1 C cb1c0 + C ac0a1 C ac1a0 C ba0b1 C ba1b0 C cb0c0 C cb0c1 + C ac0a1 C ac1a0 C ba0b1 C ba1b0 C cb0c1 C cb1c0 + C ac0a1 C ac1a0 C ba1b0 C ba1b1 C cb0c1 C cb1c0 + C ac1a0 C ac1a1 C ba0b1 C ba1b0 C cb0c1 C cb1c0 + C ac1a0 C ac1a1 C ba1b0 C ba1b1 C cb0c1 C cb1c0
(5.18)

D 0 = C ac0a0 C ac1a1 C ba0b0 C ba1b1 C cb0c0 C cb1c1 (5.19)

Symbolic Eigenanalysis

The stability of the characteristic state depends on the two eigenvalues of

  b 1 c 1 b 2 c 2   ,
which are:

λ 1 = b 1 + c 2 + (b 1 -c 2 ) 2 + 4c 1 b 2 2 (5.20) λ 2 = b 1 + c 2 -(b 1 -c 2 ) 2 + 4c 1 b 2 2 (5.21)
Property 6 These two eigenvalues are real and strictly positive.

Proof: From the expressions of b 1 , b 2 , c 1 , c 2 , we know that b 1 > 0, b 2 < 0, c 1 < 0, c 2 > 0, so we have (b 1 -c 2 ) 2 + 4c 1 b 2 > 0,
therefore these two eigenvalues are real. Moreover, the expression of the product of these two eigenvalues is:

λ 1 × λ 2 =b 1 c 2 -c 1 b 2 = D λ 1 ×λ 2 D 0 (5.22)
where

D λ 1 ×λ 2 = C ac0a1 C ac1a0 C ba0b1 C ba1b0 C cb0c1 C cb1c0 (5.23) D 0 = C ac0a0 C ac1a1 C ba0b0 C ba1b1 C cb0c0 C cb1c1 (5.24)
Therefore, λ 1 × λ 2 > 0. Since λ 1 is strictly positive, so is λ 2 . □ Suppose that two eigenvectors which are related to λ 1 and λ 2 respectively are

v 1 = (v 1 1 , v 2 1 ) and v 2 = (v 1 2 , v 2 2
). We have the following property on v 1 and v 2 .

Property 7 v 1 1 × v 2 1 < 0 and v 1 2 × v 2 2 > 0.
Proof: Since v 1 is a eigenvector related to λ 1 , we have:

λ 1   v 1 1 v 2 1   =   b 1 c 1 b 2 c 2     v 1 1 v 2 1   (5.25)
From which, we have:

v 2 1 = λ 1 -b 1 c 1 v 1 1 (5.26)
By developing the expression λ 1 -b 1 c 1 , we have:

λ 1 -b 1 c 1 = c 2 -b 1 + (c 2 -b 1 ) 2 + 4c 1 b 2 2c 1 (5.27)
We can see that whether

c 2 -b 1 is positive or negative, λ 1 -b 1 is always strictly positive because |c 2 -b 1 | < (c 2 -b 1 ) 2 + 4c 1 b 2 . Since c 1 is strictly negative, we have v 1 1 × v 2 1 < 0. Similarly for v 2 , we can have: v 2 2 = λ 2 -b 1 c 1 v 1 2 (5.28) and λ 2 -b 1 c 1 = c 2 -b 1 -(c 2 -b 1 ) 2 + 4c 1 b 2 2c 1 (5.29)
So, similarly, whether c 2 -b 1 is positive or negative, λ 2 -b 1 is always strictly negative. And since c 1 is strictly negative, we have v 1 2 × v 2 2 > 0. □ Based on these properties, we develop the following theorem to verify the stability of the characteristic state.

Theorem 7

The characteristic state is unstable if and only if λ 1 ≥ 1.

Proof: For any hybrid trajectory τ sufficiently close to the characteristic state but different from the characteristic state, which reaches the chosen Poincaré section at a hybrid state s = ((x 1 , 0, z 1 ), (0, 1, 1)), and then reaches this Poincaré section again for the first time at another hybrid state s ′ = ((x 2 , 0, z 2 ), (0, 1, 1)), such that there is no sliding mode between s and s ′ , since the hybrid state ((1, 0, 0), (0, 1, 1)) (which is related to the characteristic state) is a fixed point of the Poincaré map, we can have the following relation:

  x ′ 2 z ′ 2   =   b 1 c 1 b 2 c 2     x ′ 1 z ′ 1   (5.30) where (x ′ 1 , z ′ 1 ) and (x ′ 2 , z ′ 2 )
are the new coordinates of (x 1 , z 1 ) and (x 2 , z 2 ) by taking (1, 0) as the new origin, which means that Vectors (x ′ 1 , z ′ 1 ) and (x ′ 2 , z ′ 2 ) can be decomposed by:

x ′ i = x i -1, z ′ i = z i -0,
  x ′ 1 z ′ 1   = α 1 v 1 + β 1 v 2 (5.31)   x ′ 2 z ′ 2   = α 2 v 1 + β 2 v 2 (5.32)
where

α 1 , β 1 , α 2 , β 2 ∈ R and α 1 is not null as x ′ 1 < 0, z ′ 1 > 0, v 1 1 × v 2 1 < 0 and v 1 2 × v 2 2 > 0.
By multiplying Equation 5.31 on the left by

  b 1 c 1 b 2 c 2 
 and from Equation (5.30), we can derive that:

α 2 = λ 1 α 1 (5.33)
β 2 = λ 2 β 1 (5.34)
Consider the case where λ 1 ≥ 1. Since |α 1 | is strictly positive, then we have |α 2 | ≥ |α 1 |, which means that s ′ is not closer to the characteristic state than s in the direction of v 1 .

Therefore, any sequence of the intersection points of τ with this Poincaré section will never converge to the characteristic state, which indicates that the characteristic state is unstable. This proves that "the characteristic state is unstable" is a necessary condition for "λ 1 ≥ 1". Now consider the case where the characteristic state is unstable, meaning that there is no damped oscillation that converges to it. We suppose that λ 1 < 1. Since

λ 1 > λ 2 > 0, we have |α 2 | < |α 1 | and |β 2 | < |β 1 |,
which means that if τ is sufficiently close to the characteristic state, then it converges to the characteristic state, which contradicts the hypothesis stating that the characteristic state is unstable. This proves that "the characteristic state is unstable" is a sufficient condition for "λ 1 ≥ 1". □ Based on Hypothesis 1 and Theorem 7, the condition λ 1 ≥ 1 is a sufficient and necessary condition for the existence of sustained oscillations in this HGRN of canonical repressilator.

Since our final objective is to provide practical information for the construction of synthetic networks, conditions like λ 1 ≥ 1 might not be a good result, because the set of models under this constraint is not easy to figure. In the next section, a method to compute separable constraints based on the condition λ 1 ≥ 1 is presented.

Computation of Sufficient Separable Constraints on Parameters

In this section, we propose a method to compute separable constraints on parameters based on the condition λ 1 ≥ 1 which is developed in the previous section. In this chapter, "separable constraints" mean constraints with separable form: Each parameter is constrained by an interval, for example C ac0a0 ∈ [C ac0a0 , C ac0a0 ]. The reason why we choose constraints of separable form is that they can be easily interpreted and used. In other words, these separable constraints represent a n-dimensional bounding box in the space of parameters (celerities). What we want to do is to find such a n-dimensional box so that any model in this bounding box satisfies the condition λ 1 ≥ 1, which means that any model in this box has sustained oscillations.

Firstly, we present a simplification of the condition λ 1 ≥ 1. Secondly, we introduce a method to verify if all models in a given bounding box satisfy this simplified condition. At last, using the method in the second part, we propose a search algorithm to find some separable constraints.

Before introducing the Bernstein coefficients, we firstly introduce the notion of multiindice. A multi-indice is a vector of non-negative integers. Given two multi-indice i = (i 1 , i 2 , ..., i n ) and j = (j 1 , j 2 , ..., j n ), we write i ≤ j if ∀k ∈ {1, 2, ..., n}, i k ≤ j k . We also write i j for ( i 1 j 1 , i 2 j 2 , ..., in jn ) and

  j i   for   j 1 i 1     j 2 i 2   ...   j n i n   which is the multiplication of all binomial coefficient   j k i k   , k ∈ {1, 2, ..., n}.
Using the multi-indice, a polynomial f : R n → R can be represented as follows:

f (x) = i∈I d a i x i (5.43)
where a i ∈ R, i and d are multi-indices, I d is the set of all multi-indices i such that i ≤ d, and

x i = x i 1 1 x i 2 2 ...x in n is the product of all x i j j
, where x j is the j th variable of polynomial f . f can also be expressed by Bernstein expansion as follows:

f (x) = i∈I d b i B d,i (x) (5.44)
where One fundamental property of Bernstein coefficients for our approach is the range enclosure property, which can be derived from the convex hull property. The convex hull of a set S, noted Conv(S), is the smallest convex set that contains S. Since this HGRN of the canonical repressilator has 12 parameters (see Table 5.1), and if we assume that the number of possible smallest intervals for each parameter are the same, noted m, then there are at most m 12 smallest bounding boxes to check. Verifying all these boxes can be time consuming. In our implementation, we assume that the intervals of these three genes are identical, which means that we search for bounding boxes such that C aciaj = C baibj = C cbicj and C aciaj = C baibj = C cbicj for any i, j ∈ {0, 1}, where, for example, [C ac0a0 , C ac0a0 ] is the interval of the parameter C ac0a0 , so that we only need to consider 4 independent intervals when searching for bounding boxes. This assumption is only applied here to decrease the number of possible bounding boxes. A similar assumption about the symmetry between these three genes was also made in works based on differential equations, see for example [START_REF] Buşe | Dynamical properties of the repressilator model[END_REF]. We also assume that the minimal length of interval is greater or equal to 0.5. A value smaller than 0.5 could naturally be chosen, but this might exponentially increase the number of possible bounding boxes, which could also exponentially increase the execution time. With these assumptions, we obtain 5 bounding boxes which satisfy the condition2. In the results below, (y, x) represents (c, a), (a, b) and (b, c), where y inhibits x, for example C xy0x0 presents C ac0a0 , C ba0b0 and C cb0c0 .

B d,i (x) = β d 1 ,i 1 (x 1 )...β dn,in (x n ) (5.45)
β d k ,i k (x k ) =   d k i k   x i k k (1 -x k ) d k -i k (5.46) b i = j≤i   i j     d j   a j ( 5

Lemma 1 (Convex hull property)

Conv {(x, f (x)) | x ∈ B} ⊆ Conv {(i/d, b i ) | i ∈ I d }.

Lemma 2 (Range enclosure property) min {b

i | i ∈ I d } ≤ f (x) ≤ max {b i | i ∈ I d },
Bounding box 1:

C xy0x0 ∈ [0, 0.5], C xy0x1 ∈ [0, 0.5], C xy1x0 ∈ [0.5, 1], C xy1x1 ∈ [0, 1] Bounding box 2: C xy0x0 ∈ [0, 0.5], C xy0x1 ∈ [0.5, 1], C xy1x0 ∈ [0, 0.5], C xy1x1 ∈ [0, 0.5] Bounding box 3: C xy0x0 ∈ [0, 0.5], C xy0x1 ∈ [0.5, 1], C xy1x0 ∈ [0.5, 1], C xy1x1 ∈ [0, 1] Bounding box 4: C xy0x0 ∈ [0.5, 1], C xy0x1 ∈ [0.5, 1], C xy1x0 ∈ [0, 0.5], C xy1x1 ∈ [0, 0.5] Bounding box 5: C xy0x0 ∈ [0.5, 1], C xy0x1 ∈ [0.5, 1], C xy1x0 ∈ [0.5, 1], C xy1x1 ∈ [0, 1]
We can see that these constraints are easy to interpret and some intuitive results can be derived from them: For instance, from Bounding box 2, we can get that sustained oscillations exist if the values of C xy0x1 are close to each other and are sufficiently larger than any C xy0x0 , C xy1x0 and C xy1x1 .

Summary

In this chapter, we constructed a HGRN of the canonical repressilator. By computing and analyzing analytically a Poincaré map and based on a hypothesis, we develop a sufficient and necessary condition for the existence of sustained oscillations. Then we adapt the range enclosure property of Bernstein coefficients to find some bounding boxes in parameters space which satisfy this sufficient and necessary condition. These bounding boxes (intervals of parameters) can provide useful information for the design of synthetic circuits. Moreover, an intermediate result implies some new control strategies for sustained oscillations: Controlling the absolute values of the derivatives of one gene under certain regulations such that these values are sufficiently small, while keeping other parameters of the system unchanged. The work of this chapter can be considered as an extension of the limit cycle analysis method in Chapter 3 on a specific gene regulatory network because both works rely on the method of Poincaré map. The difference between the work of this chapter and the work of Chapter 3 is that, in the work of this chapter we analyze symbolically a HGRN while in the work of Chapter 3 the exact values of parameters are assumed to be known. This symbolic analysis approach can be further applied to other networks. This chapter focuses on the analysis of a HGRN of the canonical repressilator, while another important problem is the identification of parameters based on data which is investigated in the next chapter.

Chapter 6

PARAMETER IDENTIFICATION BASED ON TIME SERIES DATA

In the previous chapter, we studied the dynamical properties of a specific HGRN under certain constraints on parameters. An another challenge is how to find these constraints on parameters from data. In this chapter, we study the problem of parameter identification of HGRNs. In the literature, there are works focusing on identifying parameters of HGRNs based on biological knowledge [START_REF] Behaegel | Constraint identification using modified Hoare logic on hybrid models of gene networks in 24th International Symposium on Temporal Representation and Reasoning[END_REF][START_REF] Michelucci | Evolutionary continuous optimization of hybrid Gene Regulatory Networks[END_REF]. This thesis studies parameter identification methods based on time series data. We consider the parameter identification problem as an optimization problem and we mainly focus on meta-heuristic optimization methods like the genetic algorithm. As opposed to many parameter identification problems, which aim to find a specific choice of parameters, in this chapter we want to find a set of parameters that can provide more information about the possible qualitative properties of the system. We use the HGRN model of the repressilator studied in the previous chapter as a case study, although the method can be generalized on other HGRNs.

Definition of the Parameter Identification Problem

This section defines the parameter identification problem concerned in this chapter. Firstly, the parameters and the form of the parameter set to identify are presented in Section 6.1.1. Then the type of data is presented in Section 6.1.2.

Representation of the Parameter Set to Identify

In this chapter, we consider the same HGRN studied in Chapter 5, but the method is generic and can be extended for other HGRNs. Chapter 5 focused on the analysis of qualitative behaviors so the absolute values of the thresholds that separate discrete levels are not considered, while in this chapter, since the objective is to identify parameters from continuous expression data, we not only need to consider the parameters describing the time derivatives in each discrete state but also need to consider these thresholds.

Therefore, there are two classes of parameters considered in the work of this chapter: The parameters describing the time derivatives of each discrete state are denoted by

V ac0a0 , V ac0a1 , V ac1a0 , V ac1a1 , V ba0b0 , V ba0b1 , V ba1b0 , V ba1b1 , V cb0c0 , V cb0c1 , V cb1c0 , V cb1c1
; and the parameters describing the thresholds are denoted by θ ca , θ ab , θ bc . By assuming that all these parameters are positive, the dynamics of gene X is described by the following equations if X does not reach its maximum or minimum value.

dX dt =                V XY 0X0 if x Y < θ Y X ∧ x X < θ XZ V XY 0X1 if x Y < θ Y X ∧ x X ≥ θ XZ -V XY 1X0 if x Y ≥ θ Y X ∧ x X < θ XZ -V XY 1X1 if x Y ≥ θ Y X ∧ x X ≥ θ XZ (6.1)
where (Y, X, Z) ∈ {(c, a, b), (a, b, c), (b, c, a)}; x a /x b /x c is the continuous expression of gene a/b/c over time. If X reaches its maximum or minimum value, then its time derivative is 0 which corresponds to the existence of a sliding mode. In this chapter, we impose that the maximum value of each gene is 1 and the minimum value of each gene is 0 (the reason for this choice will be presented in the next subsection), so θ ca , θ ab , θ bc are real numbers between 0 and 1. In fact, the HGRN studied in Chapter 5 can be considered as a normalization over two levels (between 0 and 2) of the model considered in this chapter. Many works about parameter identification aim to find a specific choice of parameters, in other words, assign a specific value to each parameter. In the work of this chapter, instead of finding a specific choice of parameters, we want to find a parameter set which can be described by constraints on parameters under the following form:

∀(Y, X) ∈ {(c, a), (a, b), (b, c)}, ∀i, j ∈ {0, 1}, V XY iXj ∈ I XY iXj = [I XY iXj , I XY iXj ] and ∀(Y, X) ∈ {(c, a), (a, b), (b, c)}, θ Y X ∈ I Y X = [I Y X , I Y X ] where ∀(Y, X) ∈ {(c, a), (a, b), (b, c)}, ∀i, j ∈ {0, 1}, I XY iXj , I XY iXj , I Y X , I Y X ∈ R, 0 ≤ I XY iXj < I XY iXj , 0 ≤ I Y X < I Y X ≤ 1.
In other words, each parameter is constrained individually by an interval. The form of these constraints is the same as the form of separable constraints in Chapter 5. The difference is that we also consider here the constraints on thresholds. In this chapter, this kind of parameter set constrained by intervals is called regular parameter set. A regular parameter set is denoted by P s .

In order to simplify the parameter identification problem, we choose to fix the length of these intervals to identify; more precisely, we want to find a regular parameter set under the following form:

∀(Y, X) ∈ {(c, a), (a, b), (b, c)}, ∀i, j ∈ {0, 1}, V XY iXj ∈ [max(0, VXY iXj -ϵ XY iXj ), VXY iXj + ϵ XY iXj ] and ∀(Y, X) ∈ {(c, a), (a, b), (b, c)}, θ Y X ∈ [max(0, θY X -ϵ Y X ), min(1, θY X + ϵ Y X )]
where ∀(Y, X) ∈ {(c, a), (a, b), (b, c)}, ∀i, j ∈ {0, 1}, ϵ XY iXj and ϵ Y X , which are positive real numbers, are chosen (fixed), so we only need to identify the values of VXY iXj and θY X .

The reason why we choose to identify a parameter set is that a parameter set can cover more possible qualitative properties of the system compared to a single choice of parameters, and the reason why we choose to use a regular parameter set is that it is easier to sample models from such a parameter set.

Data Representation

In this chapter, we investigate parameter identification methods based on time series data. The inference of gene regulatory networks from time series data is in fact an active research area [START_REF] Byron | A comparative review of recent bioinformatics tools for inferring gene regulatory networks using time-series expression data[END_REF][START_REF] Ribeiro | Learning any memory-less discrete semantics for dynamical systems represented by logic programs[END_REF]. The time series data used in this chapter is formally represented as a sequence: (x t 0 X , x t 1 X , ..., x tm X ) where X ∈ {a, b, c} represents one gene of the system; the sequence (t 0 , t 1 , ..., t m ) is a sequence of time representing the moments when we know the continuous expressions of genes; x t i X represents the expression of gene X at time t i . For now, we do not have biological time series data that correspond to the canonical repressilator, thus simulated data is used in this chapter to verify the parameter identification methods. More precisely, we randomly choose parameters of this model and make a simulation from a randomly chosen initial state to obtain time series data. In this chapter, we simply choose that, for each gene, the maximum value is 1 and the minimum value is 0. Indeed, a priori, the maximum and minimum values are unknown; in that case, we can also consider these values as parameters to identify. In order to make the data more realistic, we also add noise independently on each data point following a normal distribution.

The data are presented in Table 6.1 and their illustration is given in Fig 6 .1. In Table 6.1, there are some values which are less than 0 or bigger than 1 because of the added noise. The parameters of the HGRN used to obtain these data are given as follows. Although these values are chosen randomly without any biological consideration, they can be used to verify the parameter identification methods at the first stage. We consider the parameter identification problem as an optimization problem. More precisely, we use an objective function to describe how well a regular parameter set fits the data and we want to find a regular parameter set that minimizes this objective function. This objective function is given in Eq (6.2).

V ac0a0 = 1 V ac0a1 = 1.9 V ac1a0 = 1.3 V ac1a1 = 0.4 V ba0b0 = 3.8 V ba0b1 = 2.5 V ba1b0 = 2.7 V ba1b1 = 3.3 V cb0c0 = 1.5 V cb0c1 = 0.8 V cb1c0 = 1.9 V cb1c1 =
lim ns→∞ 1 n s M i ∈Sample(Ps) f (M i ) (6.2)
where n s models are randomly sampled from a given regular parameter set P s (here n s models mean n s different choices of parameters) following the uniform distribution. These n s sampled models are denoted by Sample(P s ) = {M 1 , M 2 , ..., M ns }. The expression of function f is given in Eq (6.3).

f (M i ) = 1 3(m + 1) X∈{a,b,c} j∈{0,1,2,...,m} w(x t j X )(x i,t j X -x t j X ) 2 (6.3)
where, for each sampled model M i , we compute the simulated time series data from the same initial states (x t 0 a , x t 0 b , x t 0 c ), and these simulated data are denoted by (x i,t 0 X , x i,t 1 X , ..., x i,tm X ), X ∈ {a, b, c}. The mean squared error is used to measure the difference between the input data and these simulated data. We also add weights to the mean squared error because some data points are more important than others for parameter identification. The weight of data point x t j X is denoted by w(x t j X ), which is proportional to the estimated time derivative of x t j X . More precisely, w(x

t j X ) = max(| x t j X -x t j-1 X t j -t j-1 |, | x t j+1 X -x t j X t j+1 -t j |) if 0 < j < m, w(x t j X ) = | x t j+1 X -x t j X t j+1 -t j | if j = 0 and w(x t j X ) = | x t j X -x t j-1 X t j -t j-1 | if j = m.
In fact, this objective function describes the average of function f on a regular parameter set.

In practice, it is impossible to let n s converge to infinity, so in the implementation we choose a finite value of n s to compute the estimated value of the objective function.

Estimation of Potential Parameter Set

The final objective of this work is to find a regular parameter set that optimizes globally this objective function (a global minimum). However, finding the global minimum is non-trivial. So what we can do is to find some local minimums.

In the parameter space, there exist some local minimums which are not reasonable, see for example Fig 6 .2. This local minimum is not reasonable because the frequencies of the models that belong to this local minimum are too high. Such unreasonable local minimum exists because there is a finite number of data points in the real sequence. The number of such unreasonable local minimums decreases when the number of data points increases. In order to avoid finding these unreasonable local minimums, before applying the optimization algorithms, we propose a method to estimate a potential parameter set that does not contain these unreasonable local minimums. This idea is shown in Fig 6 .3. The grey area represents the whole parameter space, the red squares represent some regular parameter sets which are local minimums. In this illustration, for example, the left red square with the illustration of simulated data represents an unreasonable local minimum and the right one with the illustration of simulated data represents a reasonable local minimum. What we want to do is to find a potential parameter set (the blue area) that does not contain the unreasonable local minimums. Then we could search for local minimums starting from this potential parameter set, which could improve the efficiency of the searching process.

Two types of constraints are considered to compute this potential parameter set. The first type consists in constraints on the maximum absolute value of the time derivatives of genes, which can be computed by the estimated time derivatives. More precisely, the maximum absolute value of the estimated time derivative of one gene is increased by an order of magnitude and this increased value is used as the constraint on the maximum absolute value of the time derivatives of this gene. The second type is based on the number of boundaries reached between t 0 and t m . We assume that this sequence corresponds to an oscillation, and that it shows more than one period and less than two periods. So, based on the analysis of this HGRN in Chapter 5, we know that the number of boundaries reached between t 0 and t m must be between 6 and 36, where 6 is the minimum number of boundaries reached by one period and the number of boundaries reached by two periods must be less than 36. Details about how this potential parameter set is used are presented in the following section.

Parameter Identification based on Genetic Algorithm

Generally, the methods to optimize an objective function can be classified in two categories based on whether the information of gradient is used or not. For now, we have not yet found a method to compute the gradient of this objective function, so we chose to use methods that are not based on gradient.

The most used class of optimization methods that is not based on gradient is metaheuristics [START_REF] Burke | Search methodologies: introductory tutorials in optimization and decision support techniques[END_REF][START_REF] Liang | Comprehensive learning particle swarm optimizer for global optimization of multimodal functions[END_REF][START_REF] Liang | Classified perturbation mutation based particle swarm optimization algorithm for parameters extraction of photovoltaic models[END_REF][START_REF] Eberhart | A new optimizer using particle swarm theory in MHS'95[END_REF][START_REF] Kennedy | Particle swarm optimization in Proceedings of ICNN'95-international conference on neural networks[END_REF][START_REF] Hansen | The CMA evolution strategy: a comparing review, Towards a new evolutionary computation: Advances in the estimation of distribution algorithms[END_REF][START_REF] Storn | Differential evolution-a simple and efficient heuristic for global optimization over continuous spaces[END_REF]. Among different meta-heuristics, in this chapter, we choose to use the genetic algorithm [START_REF] Reeves | Genetic algorithms[END_REF][START_REF] Katoch | A review on genetic algorithm: past, present, and future[END_REF][START_REF] Such | Deep neuroevolution: Genetic algorithms are a competitive alternative for training deep neural networks for reinforcement learning[END_REF] which is one of the most studied meta-heuristics. In fact, genetic algorithms have also been used to identify parameters of HGRNs in [START_REF] Michelucci | Evolutionary continuous optimization of hybrid Gene Regulatory Networks[END_REF], which, contrary to this work, uses biological knowledge to identify the parameters. The idea of genetic algorithm comes from the process of natural selection. The main pipeline of genetic algorithm is illustrated in Fig 6 .4. Details about each step are presented as follows. In the step of initial population, some regular parameter sets that intersect the potential parameter set are randomly generated. In practice, to verify if a regular parameter set intersects the potential parameter set, we simply verify if the center of this regular parameter set (for example, if V XY iXj ∈ [max(0, VXY iXj -ϵ XY iXj ), VXY iXj + ϵ XY iXj ], then the center of V XY iXj is VXY iXj ) belongs to the potential parameter set. These randomly generated regular parameter sets are called the initial population (the first generation).

In the step of selection, we compute the estimated values (as stated before, the exact values of the objective function are hard or maybe impossible to compute, so we estimate the values by choosing a finite value of n s ) of the objective function of each regular parameter set in the current population and choose the top n top regular parameter sets, which are the n top regular parameter sets with the minimum n top estimated values of the objective function.

In the step of generation of new population (new generation), we compute the new population based on the top n top regular parameter sets from the previous population. The generation of new population is based on two operators: Crossover and mutation. In the literature, many sorts of crossover and mutation have been proposed [START_REF] Lim | Crossover and mutation operators of genetic algorithms[END_REF]. In this chapter, we adapt the classic crossover and mutation which are presented as follows.

Consider two regular parameter sets P p s and P q s . The intervals of P p s and P q s are denoted by

V XY iXj ∈ [max(0, V p XY iXj -ϵ XY iXj ), V p XY iXj + ϵ XY iXj ] θ Y X ∈ [max(0, θp Y X -ϵ Y X ), min(1, θp Y X + ϵ Y X )]
and

V XY iXj ∈ [max(0, V q XY iXj -ϵ XY iXj ), V q XY iXj + ϵ XY iXj ] θ Y X ∈ [max(0, θq Y X -ϵ Y X ), min(1, θq Y X + ϵ Y X )]
respectively, where (Y, X) ∈ {(c, a), (a, b), (b, c)}, i, j ∈ {0, 1}. The operator crossover can generate a new regular parameter set, denoted by P pq s , from P p s and P q s . The intervals of P pq s are

V XY iXj ∈ [max(0, V pq XY iXj -ϵ XY iXj ), V pq XY iXj + ϵ XY iXj ] θ Y X ∈ [max(0, θpq Y X -ϵ Y X ), min(1, θpq Y X + ϵ Y X )]
where

V pq XY iXj = α V p XY iXj + (1 -α) V q XY iXj θpq Y X = α θp Y X + (1 -α) θq Y X
where α is a randomly generated real number from [0, 1] following a uniform distribution. The mutation operator can generate a new regular parameter set P p′ s from P p s . The intervals of P p′ s are

V XY iXj ∈ [max(0, V p′ XY iXj -ϵ XY iXj ), V p′ XY iXj + ϵ XY iXj ] θ Y X ∈ [max(0, θp′ Y X -ϵ Y X ), min(1, θp′ Y X + ϵ Y X )]
where

V p′ XY iXj = max(0, β XY iXj V p XY iXj ) θp′ Y X = min(1, max(0, β Y X θp Y X ))
where β XY iXj , β Y X are randomly generated real numbers from [1 -δ, 1 + δ] where δ is a small positive real number. In other words, the mutation operator adds some small noises on the original intervals.

In order to compute the new population, firstly we apply the crossover on pairs of top n top regular parameter sets from the previous populations to get some new regular parameter sets, and then we apply the mutation on these new regular parameter sets to get finally the new population.

Different criteria can be used to determine when the algorithm stops; in this work, we choose to fix the number of generations, in other words the loop in Fig 6 .4 is repeated for a fixed number of times.

Application

The minimum value of the objective function of each generation (the estimated value of the objective function of the regular parameter set that minimizes the objective function in the current generation) is shown in Fig 6 .5. The execution time of the algorithm is 17 minutes. Computations were performed on a standard laptop computer, with an Intel Core I7-8550U 1.80GHz processor and 16.0GB RAM. Globally, we can see that new generations are generally better with regard to the estimated values of the objective function and it seems that this curve converges. So we can assume that a local minimum is found which means that the best regular parameter set in the last population is very close to a local minimum.

Simulated data of this best regular parameter set in the last population are shown in Fig 6 .6. We can see that the simulated data are close to the input data, which indicates that this parameter set could potentially cover the real qualitative properties of the system. This regular parameter set is given as follows. We can see that this regular parameter set does not include the real parameters. In fact, because of the facts that we search for a set of parameters instead of a specific choice of parameters and noises are added to data, data, does there exist another local minimum which generates data that visually fit better the input data? Or can we verify if the found local minimum is the global minimum? We note that the global minimum does not necessarily generate data that visually fit best the input data, but we can assume that it is true if the objective function is chosen properly.

Summary

In this chapter, we propose a method based on a genetic algorithm to identify parameter sets from time series data. The results show that the obtained parameter sets could potentially cover the qualitative properties of the real system. The major limit of this work is that only simulated data are used. To further evaluate this method, real data from biological experiments are required, because the distribution of real data could be different from simulated data which could potentially influence the effectiveness of the method. There are also open questions related to the identification of the global minimum that can be further investigated. Generally, this work is only a beginning of the study of parameter identification of HGRN based on time series.

The study of this chapter and the previous one were based on the 3-dimensional repressilator. In the next chapter, we will investigate the dynamical properties of the 4-dimensional repressilator with a discrete formalism.

Chapter 7

CONDITION FOR A DISCRETE PERIODIC ATTRACTOR IN 4-DIMENSIONAL REPRESSILATORS

In this chapter, we study the 4-dimensional repressilators, which is a continuation of the work in Chapter 5. Initially, the idea was to develop, like in Chapter 5, constraints on celerities of 4-dimensional repressilators for the existence of sustained oscillations. However, it is a non-trivial task, because this extra dimension makes the use of Poincaré map to analyze the stability of characteristic states more complicated. Thus, in this chapter, we turn to discrete models and aim to find the condition for the existence of a discrete periodic attractor based on the asynchronous semantics, because the existence of a discrete periodic attractor based on the asychronous semantics in a discrete model is a sufficient condition for the existence of a limit cycle in the associated HGRN. We aim to find such a condition described by the topological features of influence graphs, which follows some theoretical works [START_REF] Paulevé | Static analysis of Boolean networks based on interaction graphs: a survey[END_REF][START_REF] Richard | Necessary conditions for multistationarity in discrete dynamical systems[END_REF][START_REF] Richard | Negative circuits and sustained oscillations in asynchronous automata networks[END_REF][START_REF] Remy | Graphic requirements for multistability and attractive cycles in a Boolean dynamical framework[END_REF][START_REF] Richard | Attractor separation and signed cycles in asynchronous Boolean networks[END_REF]. Our major contributions include: 1) discovering that, with one exception, the relations between gene regulation thresholds do not impact the existence of discrete periodic attractors in any of the influence graphs considered in this study; 2) identifying a sufficient and necessary condition of simple form for the existence of a discrete periodic attractor when the same exception is ignored; 3) identifying new topological features of influence graphs that are necessary for predicting the existence of discrete periodic attractor in 4-dimensional repressilators. The work of this chapter has been presented at CMSB 2023 as regular paper [START_REF] Sun | Condition for Periodic Attractor in 4-Dimensional Repressilators in Computational Methods in Systems Biology[END_REF].

4-dimensional Repressilators

In the literature, there are works studying extensions of the 3-dimensional canonical repressilator by adding more genes into the network [START_REF] Page | Oscillations in well-mixed, deterministic feedback systems: Beyond ring oscillators[END_REF][START_REF] Perez-Carrasco | Combining a toggle switch and a repressilator within the AC-DC circuit generates distinct dynamical behaviors[END_REF][START_REF] Goh | Sustained oscillations in extended genetic oscillatory systems[END_REF][START_REF] Tomazou | Computational redesign of synthetic genetic oscillators for independent amplitude and frequency modulation[END_REF][START_REF] Zhang | Independent control of amplitude and period in a synthetic oscillator circuit with modified repressilator[END_REF]. Most of these works focus on 123 Consider the influence graph of Fig 7 .1, G 2 inhibits G 0 and G 3 , so G 2 has two distinct thresholds because of Assumption 7, meaning that it has three discrete levels: 0, 1, 2. Since each of the other genes only influences one other gene, they only have two discrete levels: 0, 1. A similar assumption can be found in [START_REF] Abou-Jaoudé | From structure to dynamics: frequency tuning in the p53-mdm2 network: I. logical approach[END_REF] for example.

G 0 G 1 G 2 G 3 - - - - -
The relations between the two thresholds of G 2 can have impact on the dynamical properties of the system. To show the relations between thresholds on an influence graph, we introduce the notion of influence graph with thresholds which is defined as IGS = (V, A, s) where V and A are the sets of genes (nodes) and regulations (edges) between genes respectively, as in the definition of an influence graph (see Chapter 2), and the function s assigns an integer to each regulation that represents the minimum discrete level of the source gene necessary to inhibit the target gene. Thus, the function s also characterizes the relationship between thresholds. G 3 ) = 2 (which is the number on the arc) means that G 3 is inhibited by G 2 if the discrete level of G 2 is bigger or equal to 2. For the regulation G 2 → G 0 in the same graph, equal to 1. We can see that, in the left influence graph with thresholds, the threshold of G 2 triggering the inhibition of G 0 is smaller than the threshold triggering the inhibition of G 3 , while in the right influence graph with thresholds, the situation is reversed. We can also see that the function s gives all possible discrete levels of the system.

From the influence graph of

s(G 2 → G 0 ) = 1 means that G 0 is inhibited by G 2 if the discrete level of G 2 is bigger or G 0 G 1 G 2 G 3 -1 -1 -1 -1 -2 G 0 G 1 G 2 G 3 -2 -1 -1 -1 -1
A priori, different discrete models (logic programs) can be associated to the same influence graph with thresholds, particularly when one gene is inhibited by several genes. Moreover, different choices of logic programs lead to different transition graphs, in other words different dynamical properties. In this chapter, we make an assumption about the dynamics when one gene is inhibited by several genes.

Assumption 8

In an influence graph with threshold IGS = (V, A, s), for any gene G, its discrete level can decrease by 1 if there exists a regulation from G ′ to G and the current discrete level of G ′ is bigger or equal to s(G ′ → G), otherwise its discrete level can increase by 1.

In fact, Assumption 8 is equivalent to assume that the inhibitions are disjunctive, meaning that only one inhibitor is enough to decrease the target gene. Similar assumptions about the disjunction or conjunction of gene regulation can be found in [START_REF] Melkman | Determining a singleton attractor of an AND/OR Boolean network in O (1.587 n) time[END_REF][START_REF] Akutsu | Finding a periodic attractor of a Boolean network[END_REF]. Consider the influence graph with thresholds on the left of Fig 7 .2; Assumption 8 leads to the following logic program (for reminder, the discrete level of gene G i is noted by a i ):

a 0 = 0 ← (a 2 ≥ 1) a 0 = 1 ← (a 2 < 1) a 1 = 0 ← (a 0 ≥ 1) ∨ (a 3 ≥ 1) a 1 = 1 ← (a 0 < 1) ∧ (a 3 < 1) a 2 = 0 ← (a 1 ≥ 1) ∧ (a 2 ≤ 1) a 2 = 1 ← (a 1 ≥ 1) ∧ (a 2 = 2) a 2 = 1 ← (a 1 < 1) ∧ (a 2 = 0) a 2 = 2 ← (a 1 < 1) ∧ (a 2 ≥ 1) a 3 = 0 ← (a 2 ≥ 2) a 3 = 1 ← (a 2 < 2) (7.1)
Note that, for instance, there needs to be two rules in order to make G 2 increase to the expression level 2: One to update it from level 0 to level 1 (line 7) and one to make it increase from 1 to 2 (line 8); this is because we did not constraint the dynamics to be unitary and we thus need to encode this property inside the rules. Using Assumption 8, we get a unique discrete model from any influence graph with thresholds, which simplifies the analysis.

In this work, we use the asynchronous semantics, which is formally defined as follows. We consider a system with N genes noted as G 0 , G 1 , G 2 , ..., G N . For any discrete state d s at time t, if there exists a logic rule a i 0 = k ← ϕ where i 0 ∈ {0, 1, 2, ..., N }, such that d s satisfies ϕ and d i 0 s ̸ = k, then at time t + 1, the system can reach the new discrete state d ′ s such that d ′i s = d i s for i ∈ {0, 1, 2, ..., N } \ {i 0 } and d ′i 0 s = k. The transition graph of discrete states based on the asynchronous semantics corresponding to the discrete model of Eq (7.1) is illustrated in Fig 7 .3.

The general logic rules for an arbitrary IGS are given in Eq (7.2).

a i = k + 1 ← (a i = k) ∧ (k < M ax G i ) ∧ (∀G ∈ reg(G i ), G < s(G → G i )) a i = k -1 ← (a i = k) ∧ (k > 0) ∧ (∃G ∈ reg(G i ), G ≥ s(G → G i )) (7.2)
where M ax G i is the maximum discrete level of G i and reg(G i ) is the set of all genes that inhibit G i . Obviously, the rules of Eq (7.1) can be derived from the ones of Eq (7.2) by simplification. Some simplifications also involve the knowledge of the asynchronous

Feature Selection and Search of Candidate Condition based on Decision Tree

The definition of a discrete periodic attractor is given as follows. One example of discrete periodic attractor is given by the red discrete states in Fig 7 .3.

In order to find a condition for the existence of a discrete periodic attractor, we firstly construct a decision tree model following the five steps below. The reason why we want to construct a decision tree model is that if we could obtain a decision tree with a classification accuracy of 1.0 to predict the existence of a discrete periodic attractor, then this decision tree is equivalent to a sufficient and necessary condition for the existence of a discrete periodic attractor in 4-dimensional repressilators. This work indeed only considers a finite number of discrete models, and a decision tree can be intuitively explained.

1. Generate all influence graphs with thresholds of 4-dimensional repressilators.

2. For each influence graph with thresholds, check the existence of a discrete periodic attractor using an attractor identification algorithm. Here, an influence graph having a discrete periodic attractor means that the associated discrete model (which is unique because of Assumption 8) has a discrete periodic attractor.

3. Compute manually some features which could be potentially used to predict the existence of a discrete periodic attractor.

4. Construct a decision tree model which uses the features in the previous step to predict the existence of a discrete periodic attractor.

5. Manually drop some features which do not influence the prediction result.

In step 1, there are 50625 influence graphs with thresholds in total without removing the graphs that are equivalent.

In step 2, we use the function attracting_components of the Python library NetworkX to verify the existence of a discrete periodic attractor. We find that any influence graph with thresholds which has a discrete periodic attractor has only one discrete periodic attractor.

In step 3, we compute two classes of features on the influence graph to describe the topology of the influence graph: The number of cycles of length n and the total out-degree of cycles of length n.

For the first class, since the system has 4 genes, there are only cycles of length 2, 3 and 4. We use C2, C3 and C4 to represent the numbers of cycles of length 2, 3 and 4, respectively. For example, for the influence graph in Fig 7 .1, there is no cycle of length 2 or 4 (C2 = 0, C4 = 0) and two cycles of length 3 (C3 = 2). It is logical to use these features to predict the existence of a discrete periodic attractor because, in this class of repressilators, the length of a cycle determines whether it is a negative feedback loop or a positive feedback loop as there are only inhibition regulations and the presence of loops is related to the existence of attractor(s). For example, it is already known that the presence of a negative feedback loop is a necessary condition for sustained oscillations [START_REF] Richard | Negative circuits and sustained oscillations in asynchronous automata networks[END_REF] and the presence of positive feedback loop is a necessary condition for multistability [START_REF] Richard | Necessary conditions for multistationarity in discrete dynamical systems[END_REF].

The second class of features is a new class of features introduced in this work which is defined formally as follows.

Definition 16 (Total out-degree of cycles of length n) The total out-degree of cycles of length n is the total number of arcs which go from a vertex which belongs to a cycle of length n to a vertex which does not belong to this cycle.

For example, for the influence graph in Fig 7 .1, the arc G 2 → G 3 goes from the cycle of length 3: G 0 → G 1 → G 2 , to G 3 , which does not belong to this cycle. There are two arcs like this in this influence graph:

G 2 → G 3 (for the cycle G 0 → G 1 → G 2 ) and G 2 → G 0 (for the cycle G 1 → G 2 → G 3 )
. So the total out-degree of cycles of length 3 is 2. Since, in this influence graph, there is no cycle of length 2, the total out-degree of cycles of length 2 is 0. Since any graph considered in this chapter has only 4 genes, the total out-degree of cycles of length 4 is always 0. We use OD2 and OD3 to denote the total out-degree of cycles of length 2 and 3, respectively.

To explain the motivation about these features describing the total out-degree of cycles, let's consider the two influence graphs in Fig 7 .4. The dynamical properties of these two influence graphs are different: Any influence graph with thresholds associated to the left influence graph has a discrete periodic attractor while any influence graph with threshold associated to the right influence graph does not have a discrete periodic attractor. However, the topologies of these two influence graphs are similar: The numbers of cycles of length 2, 3 and 4 of these two graphs are identical and they both have 6 arcs.

G 0 G 1 G 2 G 3 - - - - - - G 0 G 1 G 2 G 3 - - - - - - Figure
In order to find a condition for discrete periodic attractor, we need to find a way to exhibit the topological difference between these two graphs, and these new features are effective: For the left graph, OD2 = 1, OD3 = 3; for the right one, OD2 = 2, OD3 = 2. Note that these features do not depend on the relations between thresholds.

In step 4, we construct a decision tree to predict the existence of a discrete periodic attractor based on the 5 features C2, C3, C4, OD2 and OD3 using all influence graphs with thresholds considered in this chapter. This decision tree is constructed automatically using the decision tree model of the Python library Scikit-learn.

The accuracy of prediction of this decision tree is nearly 0.9990. Initially, we wished this accuracy to be 1 because in that case, the decision tree would provide a sufficient and necessary condition for the existence of a discrete periodic attractor. This small lack of accuracy is actually caused by a few influence graphs with thresholds all related to the same influence graph. By analyzing this influence graph, a very interesting result arises:

There exists one particular influence graph such that, for any influence graph with thresholds that is not associated to this influence graph (or any isomorphism), the existence of a discrete periodic attractor does not depend on the relations between thresholds and can be predicted by this decision tree with an accuracy of 1.

This particular influence graph is shown in Fig 7 .5. This figure also presents all influence graphs with thresholds, associated to this influence graph, having a discrete periodic attractor. In fact, with the exception of the relation between the thresholds presented in this figure (the only arcs assigned with numbers), the relations between the thresholds of G 2 do not influence the existence of a discrete periodic attractor, meaning that for any

G 0 G 1 G 2 G 3 -2 -2 -2 -1 -1 -1
Figure 7.5 -The particular influence graph whose the relations between thresholds influence the existence of a discrete periodic attractor. Amongst all influence graphs with thresholds associated to this influence graph, only a subset has a discrete periodic attractor; this subset is characterized by the thresholds depicted in the figure .   order of thresholds of G 2 added to this figure, it always has a discrete periodic attractor.

In step 5, we manually drop features that do not influence the accuracy of the decision tree model. To do so, we re-train the decision tree each time after dropping one feature and observe if the accuracy decreases. Finally, we find that only keeping the features OD2 and OD3 ensures the same accuracy. We have also verified that this is the only couple of features that can maintain this accuracy. The final decision tree is shown in Fig 7 .6. In this tree, the blue leaves predict the existence of a discrete periodic attractor, and the other leaves predict the non-existence of a discrete periodic attractor. We can see that except the second leaf from the right, models in all other leaves are classified correctly. In fact, this second leaf from the right contains all models associated to the particular influence graph of Fig 7.5. This means that apart from this particular influence graph, this tree describes a sufficient and necessary condition for the existence of a discrete periodic attractor in 4-dimensional repressilators.

Condition Simplification

In this subsection, we compute a simplified sufficient and necessary condition for the existence of a discrete periodic attractor based on the decision tree of Fig 7.6. The four paths which end at blue leaves, which are the leaves related to the existence of a discrete periodic attractor, are equivalent to the following logic rules: Moreover, for all influence graphs, OD2 and OD3 are not independent and they are linked by the following constraints. Since there is a finite number of models, these constraints can be easily obtained by enumerating all models and comparing the values of OD2 and By combining Eq (7.4) and Eq (7.5), we get the following result, which is a sufficient and necessary condition for the existence of a discrete periodic attractor in case that the influence graph is not equivalent to the one in This result is of simple form and we can also find some patterns in it: The values of OD2 are powers of 2 (2 0 , 2 1 , 2 2 ) except 0, and OD3 increases as OD2 increases. These patterns might lead to some general theoretical results for N-dimensional repressilators.

(OD2 ≤ 4 ∧ OD3 ≤ 4 ∧ OD2 ≤ 1 ∧ OD3 ≤ 2 ∧ OD3 > 1 ∧ OD2 ≤ 0) ∨ (OD2 ≤ 4 ∧ OD3 ≤ 4 ∧ OD2 ≤ 1 ∧ OD3 > 2) ∨ (OD2 ≤ 4 ∧ OD3 ≤ 4 ∧ OD2 > 1 ∧ OD2 ≤ 2 ∧ OD3 > 3) ∨ (OD2 ≤ 4 ∧ OD3 > 4) (7.

Number of Oscillatory Dimensions in a Discrete

Periodic Attractor

In this subsection, we also investigate the number of oscillatory dimensions in the discrete periodic attractors. For a discrete periodic attractor, oscillating in 3 dimensions means that there exists one dimension i 0 and an integer a such that for any discrete state d s in this discrete periodic attractor, d i 0 s = a, and for any dimension i which differs from i 0 , we can find two discrete states d s 1 , d s 2 in this discrete periodic attractor, such that d i s 1 ̸ = d i s 2 . An example of a discrete periodic attractor that oscillates in 3 dimensions is given in Fig 7 .7 where there is no oscillation in the first dimension. Oscillating in 4 dimensions means that for any dimension i, we can find two discrete states d s 1 , d s 2 in this discrete periodic attractor, such that d i s 1 ̸ = d i s 2 . For example, the discrete periodic attractor in Fig 7 .3 oscillates in 4 dimensions.

By automatically verifying the isomorphisms of all influence graphs with discrete periodic attractors except the influence graph of Fig 7 .5, we find that there are, in total, only 8 different (non-isomorphic) influence graphs which always have discrete periodic attractors. Any influence graph with thresholds corresponding to these 8 influence graphs has only one discrete periodic attractor. Among these 8 influence graphs, 2 of them (Fig 7

Summary

In this chapter, we study the condition for the existence of a discrete periodic attractor in 4-dimensional repressilators under some dynamical assumptions. With the guide of decision tree models, we find a special influence graph for which the relations between thresholds influence the existence of a discrete periodic attractor. For all other influence graphs, we show that the existence of a discrete periodic attractor does not depend on the relations between thresholds and we find a sufficient and necessary condition with a simple form, describing the topology of the influence graph, for the existence of a discrete periodic attractor.

G 0 ) = s(G 2 → G 0 ) = s(G 3 → G 0 ) = 1. G 0 G 1 G 2 G 3 G 0 G 1 G 2 G 3
G 0 G 1 G 2 G 3 G 0 G 1 G 2 G 3 G 0 G 1 G 2 G 3 G 0 G 1 G 2 G 3 G 0 G 1 G 2 G 3 G 0 G 1 G 2 G 3
We use an exhaustive and computational approach to find this condition and we find some patterns in this condition. In our next step, we would like to prove this result in a more mathematical way, and try to extend this result for repressilators in N dimensions.

The topological feature used in this work, that is, the total out-degree of cycles of length n, could be potentially simplified based on other more common features. If it is possible, the simplified result could be more easily used for the design of new synthetic circuits.

Meanwhile, only the topology of influence graphs is considered in this work. For future works, we will also investigate how the relations between thresholds influence some complex dynamical properties.

Chapter 8

CONCLUSION AND PERSPECTIVES

Summary of Contributions and Limits

With a lower dynamical complexity compared to ODEs and other hybrid systems, HGRN is a promising formalism to understand the global dynamical properties of gene regulatory networks. Methods related to HGRNs have not been exhaustively explored in the literature. In this thesis, we mainly focus on developing new analysis methods of HGRNs. There are two major research directions: The first one is to develop general analysis methods and the second one is to analyze some specific networks modeled as HGRNs.

For the general analysis methods, we propose a limit cycle analysis method which can identify all limit cycles inside simple cycles of discrete states (limit cycles that do not cross the same discrete state for more than once in one period) and analyze their stability; we also propose a reachability analysis method which can verify if the hybrid trajectory from a given singular hybrid state can reach certain region (a set of hybrid states). Among the state-of-the-art HGRNs of gene regulatory networks, these methods can fulfill successfully their tasks. While for general HGRNs (including other HGRNs which can model systems that are not gene regulatory networks), these methods might still have some limits.

For both methods, we ignore the existence of non-deterministic hybrid states. Indeed, for the current existing HGRNs of gene regulatory networks, ignoring these states does not influence the analysis results, mainly because such cases are very rare based on current observations. But in the future, we might have to consider non-deterministic behaviors. In that case, these algorithms should become more complex (for the identification of limit cycles, there will be more connections in the graph of discrete domains; for the reachability analysis method, a singular hybrid state could lead to several distinct hybrid trajectories), so some optimizations of these algorithms are required; and the method to analyze the stability of limit cycles should be adapted, because we might need to consider several Poincaré maps.

For the limit cycle analysis method, for now we only consider limit cycles inside simple cycles of discrete states, while, a priori, there could exist limit cycles that cross the same discrete state or the same discrete domain for more than once in one period. In that case, the enumeration of all possible cycles that could potentially contain limit cycles, becomes more difficult. In order to solve this problem, more mathematical tools are likely to be required.

For the reachability analysis method, for now we have only solved the reachability problem of hybrid trajectories with regular behaviors. Whether chaotic hybrid trajectories are decidable or not is still an open question. Therefore, exploring the boundary between decidable hybrid trajectories and undecidable ones should be an interesting theoretical study. There are also open questions related to the discrimination of chaotic hybrid trajectories: Whether we can find a sufficient and necessary condition for a hybrid trajectory to be chaotic.

In our works of analyzing some specific networks, we focus on a class of networks called repressilators. We first study the 3-dimensional repressilator using HGRNs and then study 4-dimensional repressilators using discrete models. There are also some improvements that can be done for these works.

The analysis of the 3-dimensional repressilator is based on the analysis of the stability of the characteristic state, which, in this case, has only one cycle of discrete states around it. For a HGRN of a 4-dimensional repressilator, a characteristic state can have several cycles of discrete states around it, so the application of Poincaré map becomes more complex. For now, we have not yet found a method to analyze such a characteristic state. In order to find such a method, we might need to consider other methods from control theory.

In the analysis of periodic attractors in discrete 4-dimensional repressilators, we find new topological features to describe a sufficient and necessary condition for the existence of a periodic attractor. We use an exhaustive strategy to prove that it is a sufficient and necessary condition because there is a finite and not too large number of influence graphs. If we could further prove it in a more mathematical way, then this result could be potentially extended for n-dimensional repressilators.

The final objective of these works on repressilators is to implement a synthetic repressilator with a controllable period into a living cell. Thanks to this thesis, some conditions about oscillation have been found. The next step is to compare these theoretical results with practice. So a continuation of this work is to collaborate with a team of biologists to verify whether these conditions are compatible with the real systems.

We have also investigated the parameter identification methods based on time series data. An open question of this work is that whether we can ensure that the global minimum of the objective function is found, which could be an interesting continuation of this work.

Perspectives

Following the works of this thesis, there are other interesting directions that can be explored.

Predict Possible Long-term Continuous Behaviors

There is a strong link between HGRNs and discrete models: The transition graph of discrete states of a HGRN corresponds to the transition graph of discrete states of the associated discrete model following the asynchronous semantics. The major difference between these two frameworks is that the celerities of HGRNs give extra information about the long-term continuous behaviors. Indeed, in discrete models, for example, we can identify discrete periodic attractors, but we do not know whether the continuous oscillations in these discrete periodic attractors are sustained oscillations, damped oscillations or chaotic attractors. By sampling many HGRNs from a discrete model under certain hypothesis of the distribution of celerities and by using our methods to classify trajectories (introduced in Chapter 4), we could predict possible long-term continuous behaviors of a discrete model, for example, we can know if a system can have sustained oscillations or a stable fixed point, depending on the parameters. Furthermore, we can also predict the possible co-existence of attractors, for example, a co-existence of a stable limit cycle and a stable fixed point. Theoretically this approach can be applied on HGRNs in any dimension. But for models in higher dimensions, a large amount of simulations is required, which can be time consuming. To deal with this problem, it could be useful to explore new methods for the simulation of HGRNs.

Improve the Algorithm to Identify Discrete Attractors

In our work of 4-dimensional repressilators, we find a condition, which is described by topological features of the influence graphs, for the existence of a discrete periodic attractor. If we can extend these results for n-dimensional repressilators (the form of the condition can be different in n-dimensional cases) and understand better how these topological features lead to periodic attractors, then these new results could be potentially used to accelerate the identification of attractors in more general discrete models. For instance, for certain discrete models, they could potentially give extra information about the existence or the location of certain periodic attractors.

Networks Inference based on HGRNs

In the parameter identification part of this thesis, we assume that we already know the influence graph. In some cases, such information is however unknown. It could be interesting to investigate inference methods of influence graphs based on HGRNs, which can lead to automatic construction of HGRNs with little biological knowledge (by combining the parameter identification approaches proposed in this thesis). In the literature, there are methods that can infer influence graphs based on discrete models [START_REF] Ribeiro | Learning any memory-less discrete semantics for dynamical systems represented by logic programs[END_REF][START_REF] Folschette | GULA: Learning (From Any) Semantics of a Biological Regulatory Network in BIOSS-IA 2020 workshop[END_REF]. These methods can be potentially extended for HGRNs by not only considering the discrete levels, but also considering the signs of the time derivatives estimated from data.

Applications on Other Networks

Overall, some methods of this thesis have not yet been applied on specific biological problems. If some suitable applications can be found, it could be interesting to further investigate the merits of these proposed methods. Such applications can be, for example, to find all periodic trajectories of a given system, or to estimate the continuous basin of attraction of certain attractors. Also, even though this framework is proposed to model gene regulatory networks, it has potential to solve problems of other classes of networks, for example, ecological networks [100, 101], because the dynamics of HGRNs can be considered as a simplification of certain types of ordinary differential equations. Abstract: Using dynamical models to reveal dynamical properties of gene regulatory networks can help us better understand the nature of these biological systems and develop new medical treatments. In this thesis, we focus on a class of hybrid dynamical systems called Hybrid Gene Regulatory Network (HGRN) and aim to analyze long-term dynamical properties. We propose methods to find limit cycles and analyze their stability, and to analyze the reachability in HGRNs. This is followed by a deeper study of some networks of interest for Systems Biology: The repressila-tors, and we find conditions for the existence of sustained oscillations in the 3-dimensional canonical repressilator, and conditions, which are described by topological features of the networks, for the existence of a periodic attractor in discrete 4-dimensional repressilators. In summary, this thesis proposes new methods to analyze some properties of HGRNs that were not investigated before, for instance, the stability of N-dimensional limit cycles, the reachability, etc. The results can be further developed in the future to study other large complex networks.

3 . 1 r

 31 (h): The reduction vector of hybrid state h, introduced in Section 3.3.2 S r : A reduction compatible zone, introduced in Section 3.3.2 Chapter 2

Figure 2 . 1 -

 21 Figure 2.1 -Left: A 3-dimensional influence graph. Middle (resp. Right): A transition graph of discrete states with asynchronous (resp. synchronous) semantics, based on the influence graph on the left and the logic program of Eq. 2.5.

Figure 2 . 2 -Definition 2 (

 222 Figure 2.2 -Example of a HGRN in 2 dimensions. Left: Influence graph (negative feedback loop with 2 genes, A = G 1 , B = G 2 ). Middle: Example of corresponding parameters (celerities). Right: Corresponding example of dynamics; abscissa represents gene A and ordinate represents gene B.

Fig 2 . 2

 22 only reaches one new boundary at a time, in general a hybrid trajectory can reach several new boundaries at the same time. When a hybrid trajectory reaches several output boundaries at the same time, it can cross any of them but can only cross one boundary at a time, which causes non-deterministic behaviors, see for example Fig 2.3.The simulation of HGRNs is presented more formally as follows.Consider a hybrid state h = (π, d s ) and a hybrid trajectory τ which reaches h at t > 0.

Figure 2 . 3 -

 23 Figure 2.3 -Illustration of a non-deterministic behavior.

Figure 3 . 2 -

 32 Figure 3.2 -A limit cycle (the blue periodic hybrid trajectory) with other hybrid trajectories (for instance, the red hybrid trajectory) reaching it.

Figure 3 .

 3 Figure 3.3 -A limit cycle with other hybrid trajectories converging to it in finite time.

Figure 3 . 4 -

 34 Figure 3.4 -Periodic hybrid trajectories that are not limit cycles.

Figure 3 . 5 -

 35 Figure 3.5 -Illustration of discrete domains and sequence of discrete domains of the HGRN in Fig 2.2.

Fig 3 . 6 and

 36 Fig 3.5. -Only considering the sign of celerities, it is possible that there is a hybrid trajectory which begins from D i and reaches D j without going through another boundary; see for example 00 + and 0 + 0 in Fig 3.6 and Fig 3.5: Since the celerity of 00 is positive in the first dimension and negative in the second, it is possible that there

Figure 3 . 6 -

 36 Figure 3.6 -Graph of discrete domains of the HGRN in Fig 2.2.

Theorem 1 A

 1 hybrid state h = (π, d s 0 ) belongs to the compatible zone S of T = (D 0 , D 1 , D 2 , ..., D m ) if and only if (π, d s 0 ) ∈ D 0 , (s -1 (M (D 0 ,D 1 ) s(π)), d s 1 ) ∈ D 1 , (s -1 (M (D 0 ,D 1 ,D 2 ) s(π)), d s 2 ) ∈ D 2 , ...., (s -1 (M (D 0 ,D 1 ,...,D m-1 ) s(π)), d s m-1 ) ∈ D m-1 and (s -1 (M (D 0 ,D 1 ,...,Dm) s(π)), d sm ) ∈ D m , where M (D 0 ,D 1 ,...,D i ) is the transition matrix of (D 0 , D 1 , ..., D i ) and D i is inside discrete state d s i (i ∈ {0, 1, ..

Figure 3 . 7 -

 37 Figure 3.7 -Illustration of the depth first algorithm on discrete domain (0, 0 + ).

Figure 3 . 8 -

 38 Figure 3.8 -Examples of hybrid trajectories inside different discrete trajectories of the HGRN in Fig 2.2.

Figure 3 . 9 -

 39 Figure 3.9 -Illustration of the Poincaré map of a periodic hybrid trajectory of the HGRN in Fig 2.2, where G() = s -1 (M T s()) is the Poincaré map.

Figure 3 . 10 -

 310 Figure 3.10 -Illustration of Definition 9 (Stability of limit cycles in HGRNs).

Figure 3 . 11 -

 311 Figure 3.11 -The sequence of turning states of a periodic hybrid trajectory τ : T s (τ ) = (h 1 , h 2 , ..., h 9 ).

Figure 3 . 13 -

 313 Figure 3.13 -Illustration of Proposition 2. h i and h j are two turning states on a periodic hybrid trajectory. Black boxes represent neighborhoods in the same discrete domain of h i /h j that correspond to Proposition 2.

Theorem 3 . 3

 33 Theorem Consider a periodic hybrid trajectory τ inside the closed discrete trajectory T = (D 1 , D 2 , ...D m , D 1 ), and λ 1 , λ 2 , ..., λ p the eigenvalues of the reduction matrix A of T .

Figure 3 .

 3 Figure 3.15 -Example of a limit cycle inside a sequence of discrete domains where each discrete domain contains only a singular hybrid state.

Figure 3 .

 3 Figure 3.16 -Left: Influence graph of a negative feedback loop with 3 genes, used to build three models given in Table 3.1. Right: Influence graph of a cell cycle model with 5 genes from [10]; the multiplex (m) expresses constrains on the joint activation of En and Ep on B.

Figure 3 .

 3 Figure 3.17 -Illustration of stable limit cycles and stable fixed point in HGRNs in 3 dimensions. A: Stable limit cycle in the first HGRN. B: Stable limit cycle in the second HGRN. C: Stable fixed point in the third HGRN.

  Fig 3.17 C.

Figure 3 .

 3 Figure 3.18 -Abstracted representations of the chosen Poincaré sections in the HGRNs in 3 dimensions, illustrating all possible discrete trajectories which start from and return to this Poincaré section. The blue dots represent the discrete domains and each arrow depicts one or several different discrete trajectories (each following a unique sequence of discrete domains). A: First HGRN. B: Second HGRN. C: Third HGRN.

Figure 3 .

 3 Figure 3.19 -Illustration of a hybrid trajectory inside the closed discrete trajectory 10 in Fig 3.18 A on the Poincaré section. Red triangle represents the compatible zone (not including boundaries). v 1 and v 2 represent the two eigenvectors. Hybrid trajectory from hybrid state h 1 on the compatible zone returns to the Poincaré section at hybrid state h 2 .

Figure 4 . 1 -

 41 Figure 4.1 -Illustration of Problem 1 and hybrid trajectory halting in finite time. Red arrows represent the hybrid trajectory from h ini . The blue rectangle represents R target .

  Does the hybrid trajectory τ from h ini enter the region R target ? In other words, does there exist t 0 such that τ (t 0 ) ∈ R target ? Problem 1 is illustrated in the examples of Fig 4.1, where the initial hybrid state of the hybrid trajectory (the hybrid trajectory is illustrated by red arrows) is h ini and the blue rectangle represents R target .

Figure 4 . 2 -

 42 Figure 4.2 -Illustration of hybrid trajectories attracted by cycles of discrete domains and predecessor in the same discrete state. Blues rectangles represent R target of Problem 1 and blue boxes represent their predecessors in the same discrete state.

Fig 4 . 2 .

 42 For better illustration, the examples in Fig 4.2 has only two dimensions, but this idea can be generalized for n-dimensional models.

Figure 4 . 3 -

 43 Figure 4.3 -Illustration of the idea of Theorem 6. The surrounding rectangle represents a discrete domain which has two free dimensions and the zone surrounded by dashed lines represents the compatible zone of a certain cycle of discrete domains. v 1 and v 2 represent the two eigenvectors of the reduction matrix of this cycle of discrete domains. h 1 , h 2 , ... represent the intersection points between a hybrid trajectory and this discrete domain.

  k are the eigenvectors of the reduction matrix of this cycle (D i , D i+1 , ..., D p , D 0 , ..., D i ). The idea of the Non-intersection condition is illustrated in Fig 4.4 by a 2dimensional example. Initially, we want to verify ∀j ∈ {1, 2, ..

2

 2 and which are reached before these hybrid trajectories reaching new discrete states. This set of hybrid states is illustrated by the red zone in Fig 4.4 middle (note that Fig 4.4 is just a illustration, the shape of the real red zone is slightly different from the one in this figure). If this red zone does not intersect the blue rectangle, then Condition1 is satisfied. However, this intersection is still difficult to compute. So secondly, we move the "thickness" of the red zone to the blue rectangle in Fig 4.4 right. Now to verify if Condition1 is satisfied, we only need to verify if the red trajectory in Fig 4.4 right reaches the blue and red zone, which can be done automatically and is called Non-intersection condition.

Figure 4 . 4 -

 44 Figure 4.4 -Illustration of the Non-intersection condition in Algorithm 2.

Figure 4 . 5 -

 45 Figure 4.5 -Illustration of the idea of computing the sequence l j .

Figure 4 . 6 -

 46 Figure 4.6 -Estimation results of P LC and P F P and illustration of the transition graph of discrete states

Figure 5 . 1 -

 51 Figure 5.1 -The influence graph of the canonical repressilator.

Figure 5 .

 5 Figure 5.2 gives two simulations with two different choices of parameters. The simulation on the left represents a sustained oscillation while the simulation on the right represents a damped oscillation. In the simulation on the left, gene C continues to increase from t = 0 until it reaches the maximal value, which also means that the hybrid

Figure 5 . 2 -

 52 Figure 5.2 -Simulations of the HGRN of the canonical repressilator with two different choices of parameters (Abscissa represents time and ordinate represents the sum of the fractional part and the discrete state of each gene). Parameters of the model on the left: C ac0a0 = 1, C ac0a1 = 1.9, C ac1a0 = 1.3, C ac1a1 = 0.4, C ba0b0 = 3.8, C ba0b1 = 2.5, C ba1b0 = 2.7, C ba1b1 = 3.3, C cb0c0 = 1.5, C cb0c1 = 0.8, C cb1c0 = 1.9, C cb1c1 = 1.5. Parameters of the model on the right: C ac0a0 = 1.5, C ac0a1 = 0.7, C ac1a0 = 0.6, C ac1a1 = 1.6, C ba0b0 = 2.1, C ba0b1 = 0.4, C ba1b0 = 0.3, C ba1b1 = 3.3, C cb0c0 = 1.25, C cb0c1 = 0.25, C cb1c0 = 0.23, C cb1c1 = 1.23.

Figure 5 . 3 -

 53 Figure 5.3 -Transition graph of discrete states of the HGRN of the canonical repressilator.

3 111Figure 5 . 4 -

 354 Figure 5.4 -Illustration of the closed hybrid trajectory with only instant transitions (green arrows), two special hybrid trajectories which can reach directly the characteristic state (black arrows), and a hybrid trajectory without sliding mode (red arrows) in this HGRN of the canonical repressilator.

Figure 5 . 5 -

 55 Figure 5.5 -Illustration of different qualitative behaviors in a HGRN of negative feedback loop in 2 dimensions. The three subfigures represent three different choices of parameters.

a 1 a 2 

 2 where i ∈ {1, 2}. In fact, this change of coordinates allows to remove the affine vector   in Equation 5.13.

. 47 )

 47 where d and i are multi-indices and k ∈ {1, 2, ..., n}. The values b i , for i ∈ I d are called Bernstein coefficients.

  non-linear constraints[START_REF] Ziat | Combination of Boxes and Polyhedra Abstractions for Constraint Solving in International Symposium on Formal Methods[END_REF][START_REF] Pelleau | A constraint solver based on abstract domains in International Workshop on Verification[END_REF].

Figure 5 . 6 -

 56 Figure 5.6 -Illustration of the algorithm to search for some bounding boxes.

Figure 6 . 1 -

 61 Figure 6.1 -Illustration of the time series data used to verify the parameter identification methods.

1 . 5 θ 2 6

 152 ca = 0.8 θ ab = 0.4 θ bc = 0.

Figure 6 . 2 -

 62 Figure 6.2 -Illustration of the real sequence of gene A and some simulated sequences of gene A by models belonging to a local minimum which is not reasonable.

Figure 6 . 3 -

 63 Figure 6.3 -Illustration of the idea to estimate a potential parameter set that does not contain the unreasonable local minimums. The grey area represents the whole parameter space, the red squares represent some regular parameter set which are local minimums and the blue area represents this potential parameter set.

Figure 6 . 4 -

 64 Figure 6.4 -Illustration of the major steps of a genetic algorithm.

Figure 7 . 1 -Assumption 7

 717 Figure 7.1 -An influence graph of a 4-dimensional repressilator.

Fig 7 . 1 ,

 71 by considering all different relations between thresholds, we can get two different influence graphs with thresholds as illustrated in Fig 7.2. For the regulation G 2 → G 3 in the left influence graph with thresholds, s(G 2 →

Figure 7 . 2 -

 72 Figure 7.2 -Different influence graphs with thresholds, corresponding to the influence graph of Fig 7.1.

Figure 7 . 3 -

 73 Figure 7.3 -Transition graph of discrete states of a model of 4-dimensional repressilator (corresponding to Eq (7.1)) based on asynchronous semantics. Red discrete states represent a discrete periodic attractor.

7 . 4 -

 74 Figure 7.4 -Left: Influence graph always having a discrete periodic attractor. Right: Influence graph never having a discrete periodic attractor.

3 )Since

 3 OD2 and OD3 are integers, these logic rules can be simplified as follows.(OD2 = 0 ∧ OD3 = 2) ∨ (OD2 ∈ {0, 1} ∧ OD3 ∈ {3, 4}) ∨ (OD2 = 2 ∧ OD3 = 4) ∨ (OD2 ∈ {0, 1, 2, 3, 4} ∧ OD3 ∈ {5,6, 7, ...})(7.4)

If

  OD2 = 0 then OD3 ∈ {0, 1, 2, 3} If OD2 = 1 then OD3 ∈ {0, 1, 2, 3, 4} If OD2 = 2 then OD3 ∈ {0, 1, 2, 3, 4, 5} If OD2 = 3 then OD3 ∈ {0, 1, 2, 3}If OD2 = 4 then OD3 ∈ {0, 1, 2, 3, 4, 6, 8}(7.5)

Fig 7. 5 :(

 5 OD2 = 0 ∧ OD3 ∈ {2, 3}) ∨ (OD2 = 1 ∧ OD3 ∈ {3, 4}) ∨ (OD2 = 2 ∧ OD3 ∈ {4, 5}) ∨ (OD2 = 4 ∧ OD3 ∈ {6, 8}) (7.6)

  8) can have discrete periodic attractors which oscillate in both 3 and 4 dimensions depending on different relations between thresholds, and the other 6 (Fig 7.9) only have discrete periodic attractors which oscillate in 4 dimensions.

Figure 7 . 7 -

 77 Figure 7.7 -Example of a discrete periodic attractor (red discrete states) that oscillates in 3 dimensions. The discrete model corresponds to the left influence graph in Fig 7.8, with s(G 1 → G 0 ) = s(G 2 → G 0 ) = s(G 3 → G 0 ) = 1.

  Figure 7.7 -Example of a discrete periodic attractor (red discrete states) that oscillates in 3 dimensions. The discrete model corresponds to the left influence graph in Fig 7.8, with s(G 1 → G 0 ) = s(G 2 → G 0 ) = s(G 3 → G 0 ) = 1.

Figure 7 . 8 -

 78 Figure 7.8 -The two influence graphs which have discrete periodic attractors oscillating in 3 or 4 dimensions. All arcs represent inhibitions between genes.

Figure 7 . 9 -

 79 Figure 7.9 -The six influence graphs which have discrete periodic attractors oscillating only in 4 dimensions. All arcs represent inhibitions between genes.

Titre:

  Identifier et analyser les comportements dynamiques à long terme des réseaux de régulation génétique à l'aide de modélisation hybride Mot clés : Modélisation hybride, Cycle limite, Accessibilité, Répressilateur, Attracteur, Réseaux de régulation des gènes Résumé : Utiliser des modèles dynamiques pour révéler les propriétés dynamiques des réseaux de régulation des gènes peut nous aider à mieux comprendre la nature de ces systèmes biologiques et à développer nouveaux traitements médicaux. Dans cette thèse, nous nous concentrons sur une classe de systèmes dynamiques hybrides appelés réseaux de régulation des gènes hybrides (HGRN) et visons à analyser les propriétés dynamiques à long terme. Nous proposons des méthodes pour trouver des cycles limites et analyser leur stabilité, et pour analyser l'accessibilité dans HGRNs. Ceci est suivi d'une étude plus approfondie de certains réseaux d'intérêt pour la biologie des systèmes : Les répressilateurs, et nous trouvons des conditions pour l'existence d'oscillations soutenues dans le répressilateur canonique en dimension 3, et des conditions, décrites par les caractéristiques topologiques des réseaux, pour l'existence d'un attracteur périodique dans les répressilateurs discrets en dimension 4. En résumé, cette thèse propose de nouvelles méthodes pour analyser certaines propriétés des HGRNs qui n'ont pas été étudiées auparavant, par exemple la stabilité des cycles limites à N dimensions, l'accessibilité, etc. Les résultats pourront être développés à l'avenir pour étudier d'autres grands réseaux complexes. Title: Identifying and Analyzing Long-term Dynamical Behaviors of Gene Regulatory Networks with Hybrid Modeling Keywords: Hybrid modeling, Limit cycle, Reachability, Repressilator, Attractor, Gene regulatory networks

  For instance, the hybrid state M in Fig 2.2 belongs to an attractive boundary. ■ If c i s > 0, d i s is not the maximal discrete level of the i th gene, and the i th component of c s ′ is positive, then the hybrid trajectory from the current hybrid state will cross instantly the boundary e and enter the discrete state d s ′ . e is called an output boundary of d s . dτ (t + )

	dt

dτ (t)

dt i = 0, and τ (t) = (π, d s ).

  4 s ) and we also suppose that c 1′ s and c 2′ s are positive. Consider the neighborhood in the same discrete domain of h 4 : N D (h 4 , r 4 ). Any hybrid state in the neighborhood in the same discrete state of h 1 can be represented by h ′ 1

  and the restriction of τ on [0, t 0 ] is inside the discrete state d s (that is, ∀t ∈ [0, t 0 ], τ (t) ∈ d s ). In this case, we also say that τ reaches τ (t 0 ) ∈ R inside d s .-For any hybrid state h that belongs to an input boundary of d s but does not belong to any Z i , the hybrid trajectory from h cannot reach any hybrid state of R insided s .Examples of predecessors in the same discrete state are illustrated in Fig 4.2 where blues rectangles represent R target and blue boxes present their predecessors in the same discrete state. Consider Problem 1, if the hybrid trajectory τ from h ini has already crossed at least one discrete state (we say τ has already crossed a discrete state at t 0 if there exists t < t 0 such that τ (t 0 ) and τ (t) do not belong to the same discrete state) without reaching the region R target , then Problem 1 is equivalent to "Does τ reach P re ds target (R target )?".

	Proposition 5

describes the last part of the simulation inside Cycle

  the current simulation is attracted by a cycle of discrete domains or if the current simulation is probably a chaotic hybrid trajectory, we use Theorem 6 or Property 5 respectively.The main idea of the function Stop_condition is the following. Its second argument, Stop_condition is, knowing that this hybrid trajectory is attracted by a cycle of discrete domains, to determine if the hybrid trajectory can reach R target after an infinite number of transitions (see Fig 4.2 right). If it is the case, the function returns "Reached". Otherwise, if from the current hybrid state, there is no more chance to reach R target (see Fig 4.2 left), then the function returns "Yes". For both cases, this function can give the right answer in finite time, and the result stops the algorithm. However, if both cases don't apply, the function returns "No" and the algorithm continues. The idea of the function Stop_condition is based on the proof of Theorem 5. The details of the function Stop_condition are presented in Algorithm 2.

	objective of the function
	else else if Current simulation is probably a chaotic hybrid trajectory then end if return "unknown result" end if end if 18: end while 7: 13: 14: 15: 16: 17: Cycle 6: 19: return "R target is not reached"

D , represents the cycle of discrete domains attracting the current hybrid trajectory. Its first argument, Cycle h , represents the last part of the simulation inside Cycle D . The Algorithm 1 Reachability analysis algorithm Input 1: A hybrid state

h ini = (π ini , d s ini ) Input 2: A region R target = (π, d starget ) | π i ∈ [a i , b i ], i ∈ {1, 2, ..., N } Output: "R target is reached", "

R target is not reached" or "unknown result" 1: Current hybrid state h := h ini 2: while h is not a fixed point do 3:

h ′ := next hybrid state so that h → h ′ is a transition 4: if Transition h → h ′ reaches R target then 5:

return "R target is reached"

h := h ′ 8:

if Current simulation is attracted by a cycle of discrete domains then 9:

if Stop_condition(Cycle h , Cycle D , R target ) returns "Yes" then 10:

return "R target is not reached" 11: else if Stop_condition(Cycle h , Cycle D , R target ) returns "Reached" then 12: return "R target is reached" Algorithm 2 Stop condition Input 1: A list of hybrid states Cycle h = (h 0 , h 1 , h 2 , ...., h p , h ′ 0 ) which D Input 2: A list of discrete domains Cycle D = (D 0 , D 1 , D 2 , ..., D p , D 0 ) such that ∀i ∈ {0, 1, 2, ..., p} , h i ∈ D i , h ′ 0 ∈ D 0 and the hybrid trajectory from h 0 is attracted by

  h ∞ p such that the intersection hybrid states between this hybrid trajectory and D 0 , D 1 , D 2 , ..., D p converge to or reach these hybrid states after an infinite number of transitions.

	10:	var_stop ← var_stop ∧ F alse
	11:	end if
	12:	end if
	13: end for
	14: if var_stop = T rue then
	15:	return "Yes"
	16: else if var_stop = F alse then
	17:

2: var_stop ← T rue 3: for i ∈ {0, 1, 2, ..., p} do 4: if D i is on an input boundary of d starget then 5: if Hybrid trajectory from h ∞ i reaches R target inside d starget then 6: return "Reached" 7: else if Non-intersection condition is satisfied then 8: var_stop ← var_stop ∧ T rue 9:

else if Non-intersection condition is not satisfied then

Table 6 .

 6 1 -Time series data used to verify the parameter identification methods.

	t i	0	0.24 0.49 0.73 0.98 1.22 1.47 1.72 1.96 2.21 2.45 2.71 2.95 3.19 3.44 3.68 3.93 4.18 4.42 4.67 4.91
	x ti a x ti b	0.51 0.96 0.95 0.84 0.71 0.63 0.54 0.40 0.23 0.35 0.74 1.01 0.97 0.95 0.83 0.78 0.61 0.56 0.44 0.21 0.33 0.66 0.00 -0.02 0.00 -0.02 0.00 -0.02 -0.01 0.39 0.95 0.35 0.03 -0.02 0.00 -0.01 0.03 0.01 0.00 -0.01 0.36 1.03
	x ti c	0.82 0.62 0.92 1.00 0.99 0.98 0.98 1.01 0.90 0.44 0.00 0.26 0.67 0.94 0.99 1.01 1.01 0.97 1.00 0.88 0.43

  Definition 15 A discrete periodic attractor is a set of discrete states E a such that: -E a contains at least two discrete states. -For any two discrete states d s , d ′ s ∈ E a , there exists a path in the transition graph of discrete states from d s to d ′ s . -For any discrete state d s ∈ E a and for any discrete state d ′ s / ∈ E a , there is no path in the transition graph of discrete states from d s to d ′ s .

Computations were performed on a standard laptop computer, with an Intel Core I7-8550U 1.80GHz processor and 16.0GB RAM.

I gratefully acknowledge the funding provided by China Scholarship Council (CSC) that made my PhD work possible.

Condition Simplification

The condition:

can be reformulated as:

which is equivalent to:

or:

This last condition is equivalent to (P 1 ≥ 0) ∨ (P 2 ≥ 0) where P 1 and P 2 are polynomials on parameters. The expressions of P 1 and P 2 are given as follows.

(5.39)

(5.40)

Condition (P 1 ≥ 0) ∨ (P 2 ≥ 0) seems preferable to λ 1 ≥ 1 because it only contains polynomials. In fact, one can easily prove that solutions for (P 1 ≥ 0) ∨ (P 2 ≥ 0) exist. For example, by only considering Cac0a0 and Cac0a1, which are two parameters describing the derivative of gene A when gene C does not inhibit gene A, P 1 and P 2 can be expressed by:

(5.41)

where p ij (with i ∈ {1, 2}, j ∈ {1, 2, 3, 4}) are expressions of parameters which do not include C ac0a0 and C ac0a1 . We can see that if C ac0a0 and C ac0a1 converge to 0 while other parameters remain unchanged, P 1 and P 2 converge to p 14 and p 24 respectively, which are both positive. This indicates that the solutions of (P 1 ≥ 0) ∨ (P 2 ≥ 0) exist and also implies a new control strategy for the existence of sustained oscillations: Controlling the derivatives of gene A when A is not inhibited by C such that these derivatives are sufficiently small, while keeping other parameters unchanged.

Satisfiability under Separable Constraints

In this subsection, we adapt the range enclosure property of Bernstein coefficients to verify if all models in a given bounding box satisfy the condition (P 1 ≥ 0) ∨ (P 2 ≥ 0). The Bernstein coefficients have been used in the literature to, for example, compute images for polynomial dynamical system [START_REF] Dang | Reachability Analysis for Polynomial Dynamical Systems Using the Bernstein Expansion[END_REF][START_REF] Dang | Image computation for polynomial dynamical systems using the Bernstein expansion[END_REF], or compute affine lower bound functions for polynomials [START_REF] Garloff | A comparison of methods for the computation of affine lower bound functions for polynomials in International Workshop on Global Optimization and Constraint Satisfaction[END_REF], etc. ∀x ∈ B, where B = [0, 1] n is the unit box.

The range enclosure property over-approximates the range of the image of f on B and can be used to verify if all models in a given bounding box satisfy the condition (P 1 ≥ 0) ∨ (P 2 ≥ 0). To do so, we need to firstly make a change of variables of polynomials P 1 and P 2 , such that all variables are included in [0, 1]. For example, the variable

. By doing so, we get two new polynomials P ′ 1 and P ′ 2 . So now, to verify if (P 1 ≥ 0) ∨ (P 2 ≥ 0) is always true in a given bounding box, we only need to verify if (P ′ 1 ≥ 0) ∨ (P ′ 2 ≥ 0) is always true in the unit box. To do so, we compute the Bernstein coefficients {b 1,i } and {b 2,i } (where i ∈ I d ) of P ′ 1 and P ′ 2 respectively. A sufficient condition for the condition "(P ′ 1 ≥ 0) condition2), according to the range enclosure property. In fact, since the minimum value of the image of P ′ 1 on the unit box is always larger or equal to the minimum value of {b 1,i }, "{b 1,i } are not negative" (∀i ∈ I d , b 1,i ≥ 0) indicates that "P ′ 1 is not negative on the unit box", and the same holds for P ′ 2 . Since there is a finite number of Bernstein coefficients, condition2 can be verified. Therefore, in this chapter, condition2 is used to verify if all models in a given bounding box always have sustained oscillations.

Search of Separable Constraints

Based on the method introduced in the previous subsection, we propose a depth first algorithm to find some bounding boxes which satisfy condition2. This algorithm is illustrated in Figure 5.6. Initially, each parameter is included in an interval. In this implementation, without loss of generality, we assume that each parameter is included initially in [0, 1]. Then, we verify if condition2 is satisfied for this bounding box using the method proposed in the previous subsection. If it is satisfied, then it is a bounding box such that all models in it have sustained oscillations. If condition2 is not satisfied, then there might be some models in this bounding box which do not have sustained oscillations, in this case the bounding box is split into two smaller bounding boxes (by splitting the largest interval into two) which have the same volume and the process is repeated on each of these two new bounding boxes. Each path in this algorithm will stop, either when a bounding box which satisfies condition2 is found or when the length of the largest interval is smaller than a certain threshold. In fact, similar ideas are widely used to find solution sets under even the global minimum does not necessarily include the real parameters.

The above result is only the result of one execution of the algorithm. By executing the algorithm for several times, different regular parameter sets (local minimums) can be found. Some regular parameter sets can generate data that are visually close to the input data, while others cannot. This leads to an open question of this work: For a local minimum found by this method which cannot generate data that visually fit the input Blue leaves represent the models classified as having a discrete periodic attractor and orange leaves represent the models classified as not having a discrete periodic attractor. "gini" describes the purity of models in a node regarding the two classes considered here: Models having a discrete periodic attractor and models not having a discrete periodic attractor; if all models in a node belong to the same class then gini = 0, otherwise gini > 0 (gini = 1 -( number class1 number total ) 2 -( number class2 number total ) 2 ). "sample" represents the number of models in a node. The first value of "value" represents the number of models not having a discrete periodic attractor and the second value of "value" represents the number of models having a discrete periodic attractor.