Professeur

Dina Irofti

Delay Effects: a Journey from Multi-agent Systems to Genetic Networks

Keywords: Physical Constants xi Time-delay systems, Consensus problem, Stability analysis, Synthetic biology, Multiagent systems, Spectral analysis

This thesis discusses diverse types of interconnected systems through networks. We address networks of agents with cooperative tasks and propose a new consensus protocol with delays and anticipatory agents. We study the consensus reaching conditions for networks organized under the proposed model. Moreover, we derive some theoretical results, which can apply to a more general class of systems, concerning stability issues when the considered system has multiple imaginary roots. In terms of networks, this situation can correspond to the case of switching topology networks, when the network can even be disconnected at some point. We separately discuss the case of zero characteristic roots, and roots laying on the imaginary axis, except the origin. Finally, we propose a gene network model with a functionality similar to a multiplexer circuit. Thus, we control two outputs with three input signals, and we carry out a stability analysis.

We prove the uniqueness and the stability of the network steady states, and validate the continuous and deterministic model with a stochastic model.

and to be very open to my suggestions and initiatives at the same time. He supported and encouraged me to apply for and obtain two important research grants during my thesis, a DAAD Research Grant in 2015, and a Mitacs Globalink Research Award in 2016. This funding had a substantial contribution to my thesis and to my growth as a research scientist.

I would like to express my special appreciation and thanks to Marc Roussel, who not only introduced me to the wonderful new world of synthetic biology, but also was a tremendous mentor for me. He showed me how to conduct research in an efficient and organized way: properly writing a proposal for a project, learning new techniques, finding solutions when encountering unexpected problems etc. I am also indebted to him for the corrections he suggested that helped me to improve the quality of this thesis, and for all his guidance throughout the last couple of years. I would also like to thank to Fatihcan Atay for the collaboration we had during my thesis. He stimulated my interest in consensus problem and showed me how important it is to state even the most challenging problems in simple terms in order to find efficient solutions.

I am extremely grateful to Keqin Gu for all the discussions we had during the last three years. He stimulated my interest in multiple roots problems and his constructive remarks significantly improved my writing. He also showed me how a distinguished research professor can find the time to have very close professional and personal discussions with the students, which I find extremely admirable.

I would like to thank the doctoral committee members: I deeply appreciate their careful reading of my thesis and their constructive remarks. I would like to express my appreciation to many other academics and colleagues that inspired my research. My time spent at the L2S laboratory would have been by far less enjoyable without my colleagues, and it cheers me up to remember the scientific and personal discussions we had, especially at lunch time or during coffee breaks.

I would like to express my appreciation to my family for all their love and support. Last but not least, I am extremely grateful to Nicolas Delerue, whose assistance and care kept me going, even in the trickiest situations.

Contents

Abstract i

List of Figures

Physical Constants

The molar mass of TetR = 23 328 g/mol [START_REF] Postle | Nucleotide sequence of the repressor gene of the tn10 tetracycline resistance determinant[END_REF] Cell volume (E. coli) = 1.7 un motif très typique. Une partie de ces résultats a été publiée et présentée lors de deux conférences internationales [START_REF] Atay | A delayed consensus algorithm in networks of anticipatory agents[END_REF][START_REF] Irofti | On the delay margin for consensus in directed networks of anticipatory agents[END_REF] et un article de journal et en cours de préparation.

Si la première partie de la thèse se concentre sur les réseaux avec une topologie fixée et le Laplacien normalisé, nous considérons par la suite une étude théorique qui peut correspondre par exemple aux réseaux ayant une topologie changeante, avec un Laplacien normalisé. Un tel type de réseau peut même être déconnecté. Le nombre de sous-graphes connexes correspond au nombre de valeurs propres égales à zéro du Laplacien. Cela montre l'intérêt d'étudier la valeur maximale des racines multiples à l'origine, comme présenté dans le chapitre 4. Cette étude sort du cadre des réseaux ; elle est plutôt générale pour une classe de systèmes à retards. La contribution principale de ce chapitre consiste en une formule pour trouver la limite supérieure du nombre de racines multiples à l'origine. Quand cette limite maximale est atteinte, nous démontrons que les paramètres du système satisfont quelques conditions algébriques. Nous comparons notre résultat avec des résultats similaires dans la littérature. Cette comparaison indique premièrement que notre limite est plus précise et deuxièmement que les méthodes utilisées pour dériver cette limite sont différentes par rapport à un autre résultat trouvé dans la littérature : nous utilisons des techniques de l'algèbre linéaire et non pas de

Résumé étendu en français xiv l'analyse complexe. Les résultats de cette étude ont été publiés dans un chapitre d'ouvrage [START_REF] Irofti | On the Codimension of the Singularity at the Origin for Networked Delay Systems[END_REF] et publiés et présentés lors d'une conférence internationale [START_REF] Boussaada | Computing the codimension of the singularity at the origin for delay systems in the regular case: A vandermondebased approach[END_REF].

Si le Laplacien n'est pas normalisé, nous pouvons avoir des racines imaginaires, qui ne sont pas à l'origine. Nous attaquons dans le chapitre 5 un problème encore plus général, mais qui peut concerner aussi les réseaux : le cas de racines caractéristiques doubles sur

Context and motivation

In general, a multi-agent system (MAS) consists of a finite number of agents that are interconnected under a given topology and aim to achieve a certain cooperative behaviour.

Thus, we say that the agents are interconnected on a network (see Figure 1.1). Examples for such a network in engineering are smart grids, multi-agent learning, wireless communication, sensor networks and robotics. It is very common that the dynamics of such networks are affected by delays. The origin of these delays has different sources: transmission delays, information processing delays, measurement delays, scheduling delays, queuing delays. The communication aspects play a central role in such networks, aiming to fulfill a cooperative control task. There are many different cooperative control applications for large groups of individuals that can be found in nature. Common examples are schools of fish, flocks of birds and collective food-gathering in ant colonies [START_REF] Strogatz | SYNC: The Emerging Science of Spontaneous Order[END_REF].

A fundamental property of such cooperation tasks is that the group behaviour is the result of the interactions between all the individuals in the group, and it is generally not dictated by only one of the individuals. A well-known illustration of this behaviour is a school of fish, where every fish only knows where the neighbouring fish are heading, but have no idea towards which direction the whole school is heading [START_REF] Couzin | Effective leadership and decision-making in animal groups on the move[END_REF]. This situation is depicted in Figure 1.2.

Nevertheless, practical applications show other properties of multi-agent systems [START_REF] Olfati-Saber | Consensus and cooperation in networked multi-agent systems[END_REF].

In power networks, the power generators have to provide a constant voltage no matter how many consumers are connected to the network [START_REF] Pavella | Transient stability of power systems : theory and practice[END_REF]. Transport systems are made up of airplanes, trains and buses, and their common task for instance is to bring people and food from one place to another, irrespective on the number of available transportation vehicles [START_REF] Fax | Information flow and cooperative control of vehicle formations[END_REF]. Communication on the Internet consists on many routers aiming to Example of how a network can be depicted. This network can illustrate a transport system made up of airplanes. Their common task for instance is to bring people and food from one place to another within European Union countries.

transmit the information from millions of sources to equally many receivers, no matter how the transmission capacity changes. All these examples show that it is quite common for a network to change the topology in time. Moreover, this illustrates that in a large network the number of agents can be time-varying, as individuals join or leave the network. We call such networks complex systems.

Complex systems are often decomposed into subsystems. If the subsystems are identical or similar, they can be modeled as nodes of a network. The interactions between these subsystems represents the links connecting the nodes. One example describing this situation can be the Internet: the servers and routers represent the nodes of the network, and the connection between them is described by the links. Here the delays could model the communication time between the nodes, or the time requested by a node to process the information before sending it to its neighbouring nodes. For this example, the analysis of the network dynamics is crucial for understanding the process of information transfer. This analysis usually relies on the study of network behaviour with respect to some of its parameters (the communication delays, the information processing delays, or other parameters). The study of such complex systems allows not only to investigate the conditions under which the information transfer is stable or robust, but also to draw some conclusions on the relation between the dynamics and the topology of the network.

As an example for the latter, we may consider how computer viruses are spread on the Internet. We suppose that the routers are connected under a certain topology, and that the probability that an infected computer infects its neighbours needs to exceed a threshold in order to allow the virus to persist [START_REF] Brauer | Mathematical Epidemiology[END_REF]. Then, the network's topology is directly linked to the probability that the virus persists. This example also shows the relevance of network science in the fields of biology.

Network models are used in ecology to study the population dynamics in complex ecosystems. The interaction between populations of all kinds are studied, starting from the relation between a potential infected population and a virus, between predator and prey to more mutual interactions, like the interactions between plants and pollinators. A natural question arising in ecology concerns the stability and robustness of such ecosystems. When studying such population dynamics, delays arise naturally by modeling growth, decay, maturation and regeneration in ecological network models. Knowledge about interactions in these networks may be used to maintain biodiversity and to prevent diseases spreading [START_REF] Montoya | Ecological networks and their fragility[END_REF]. Another example illustrating why the study of time-delay dynamics on networks in biology is important concerns the activation and inhibition of genes in gene regulatory networks. Indeed, synthetic genetic networks have been recently created to implement processing units such as toggle-switches [START_REF] Gardner | Construction of a genetic toggle switch in Escherichia coli[END_REF] and basic logical operations [START_REF] Wang | Engineering modular and orthogonal genetic logic gates for robust digital-like synthetic biology[END_REF][START_REF] Tamsir | Robust multicellular computing using genetically encoded nor gates and chemical 'wires[END_REF][START_REF] Siuti | Engineering genetic circuits that compute and remember[END_REF]. These networks have by far a less complex topology, when compared to large networks such as the Internet, but their study is certainly challenging, due to the biochemical processes involved and the tiny volume a cell commonly used for this study has. Because it is a relatively recent area of research, synthetic biology focuses on the design and fabrication of biological systems that generally do not exist in the natural world. New circuits and biological pathways have already been constructed by synthetic biologists. Given this background, an immediate extension of this relatively recent field will arise naturally, by putting together some (sub)systems already implemented to create larger, complex systems, with new functionalities. Actually, a few steps have already been made in this direction [START_REF] Seok Moon | Genetic programs constructed from layered logic gates in single cells accessed terms of use detailed terms[END_REF].

Related work

As stated in [START_REF] Bell | The horologium oscillatorium of christian huygens[END_REF], three great works laid the foundation of modern mechanics: Galilei's "Discourses on Two New Sciences"(1638), Huygens's "Horo-logium Oscillatorium"(1673), and Newton's "Philosophiae Naturalis Principia Mathematica"(1687). Of these, cooperative behaviour has been first studied in Christian Huygens' work [START_REF] Huygens | Horologium oscillatorium. Apud F. Muguet[END_REF], when he discovered synchronization of pendulum clocks. Even if this study took place more than 300 years ago, the mechanisms behind this type of synchronization, but also other types, have been analyzed only recently, in the last decades. For instance, pioneering definitions of cooperative systems can be found in [START_REF] Enslow | What is a "distributed" data processing system?[END_REF]. Enslow introduced the cooperative nomenclature in distributed data processing. However, the term of cooperative system focal fish influential neighbours Example of cooperative behaviour in nature: in a school of fish, every fish (red fish denoting the focal fish) only knows where the neighbouring fish (depicted in blue) are heading, but have no idea towards which direction the whole school (all fish, including the ones represented in gray) is heading.

was redefined in 1981 [START_REF] Victor | Functionally accurate, cooperative distributed systems[END_REF] in the context of systems with common goals, like problem solving tasks. An overview of other early works on cooperative systems is given in [START_REF] Vamos | Cooperative systems based on non-cooperative people[END_REF],

and one of the first PhD theses on this topic is [START_REF] Tsitsiklis | Problems in decentralized decision making and computation[END_REF].

Pioneering work on the consensus problem dates a half-century back in the context of horse race bettors [START_REF] Eisenberg | Consensus of subjective probabilities: The pari-mutuel method[END_REF]. A formal study of consensus investigating linear multi-agent systems in the discrete-time domain dates back to the 1960s (see [START_REF] Degroot | Reaching a consensus[END_REF] and references therein). In recent years, there has been a tremendous amount of renewed interest in cooperative control tasks (see [START_REF] Münz | Delay robustness in cooperative control[END_REF] and references therein), as consensus (micro-controllers with local clocks agreeing on a common time), rendezvous (robots meeting at a point or achieving a certain formation), flocking (cars moving with the same velocity and with a certain distance between each other), and synchronization (phase synchronization of generators in power networks).

However, it is beyond the scope of this thesis to review all important publications in this area. The key papers and the work we focus on are cited throughout the chapters of this thesis. Here we cite only a few studies on the consensus problem [START_REF] Olfati-Saber | Consensus and cooperation in networked multi-agent systems[END_REF][START_REF] Olfati | Consensus problems in networks of agents with switching topology and time-delays[END_REF][START_REF] Olfati-Saber | Ultrafast consensus in small-world networks[END_REF][START_REF] Olfati-Saber | Consensus filters for sensor networks and distributed sensor fusion[END_REF][START_REF] Lin | Local control strategies for groups of mobile autonomous agents[END_REF][START_REF] Vincent | Convergence in multi-agent coordination consensus and flocking[END_REF][START_REF] Bauso | Non-linear protocols for optimal distributed consensus in networks of dynamic agents[END_REF][START_REF] Cortes | Achieving coordination tasks in finite time via nonsmooth gradient flows[END_REF][START_REF] Mehyar | Distributed averaging on asynchronous communication networks[END_REF][START_REF] Bliman | Average consensus problems in networks of agents with delayed communications[END_REF][START_REF] Sandhu | Relative sensing networks: Observability, estimation, and the control structure[END_REF] related to multi-agent systems and cooperative control: cooperative behaviour of flocks and swarms [START_REF] Olfati-Saber | Flocking for multi-agent dynamic systems: Algorithms and theory[END_REF][START_REF] Olfati-Saber | Swarms on sphere: A programmable swarm with synchronous behaviors like oscillator networks[END_REF][START_REF] Savkin | Coordinated collective motion of groups of autonomous mobile robots: analysis of vicsek's model[END_REF][START_REF] Moshtagh | Distributed geodesic control laws for flocking of nonholonomic agents[END_REF][START_REF] Xi | A stochastic algorithm for self-organization of autonomous swarms[END_REF][START_REF] Freeman | Distributed estimation and control of swarm formation statistics[END_REF], synchronization of coupled oscillators [START_REF] Preciado | Synchronization in generalized erdos-renye networks of nonlinear oscillators[END_REF][START_REF] Sepulchre | Collective motion and oscillator synchronization[END_REF][START_REF] Jadbabaie | On the stability of the kuramoto model of coupled nonlinear oscillators[END_REF][START_REF] Papachristodoulou | Synchronization in oscillator networks: Switching topologies and non-homogeneous delays[END_REF][START_REF] Chopra | On synchronization of kuramoto oscillators[END_REF], algebraic connectivity of complex networks [START_REF] Olfati-Saber | Ultrafast consensus in small-world networks[END_REF][START_REF] Cao | A lower bound on convergence of a distributed network consensus algorithm[END_REF][START_REF] Kim | On maximizing the second smallest eigenvalue of a state-dependent graph laplacian[END_REF],

Time-delay systems are widely used to model systems in various disciplines, including cooperative control. Numerous classic examples of such systems can be found, for instance, in [START_REF] Kolmanovskii | Introduction to the Theory and Applications of Functional Differential Equations[END_REF]. More recent applications in networked system abound [START_REF] Atay | The consensus problem in networks with transmission delays[END_REF][START_REF] Seuret | Consensus under communication delays[END_REF]. Other studies on the analysis and control of time-delay systems are [START_REF] Pilbauer | Delayed resonator with distributed delay in acceleration feedback: Design and experimental verification[END_REF][START_REF] Hengster-Movric | Cooperative synchronization control for agents with control delays: A synchronizing region approach[END_REF][START_REF] Michiels | Control design for timedelay systems based on quasi-direct pole placement[END_REF]. An accurate estimate of delays is often rather difficult, which made it especially important to analyze how the system stability may change as the delays deviate from the nominal values.

One such study for systems with two delays can be found in [START_REF] Gu | On stability crossing curves for general systems with two delays[END_REF]. Another study on complex time-delay systems is [START_REF] Fatihcan | Complex Time-Delay Systems[END_REF], where various problems, such as coupled oscillators, the suppression of oscillations by time-delays, the stabilization of periodic solutions, the stability of neural networks, and traffic flow when time-delays are present are discussed.

Such networks require a wide range of methods and techniques to address major challenges of modelling, analyzing, and understanding the link between the self-organization and the evolution of the network. A collection of recent work on delays and networked control systems, discussing different approaches and cooperative control techniques can be found in [START_REF] Seuret | Delays and Networked Control Systems[END_REF].

Focus and contributions of the thesis

We address the consensus problem in both directed and undirected networks with fixed topology. We propose a new consensus algorithm with delays, where all agents in the network use past information to estimate future opinions (or states) of their neighbours.

One conceptual contribution of the proposed model consists in the different source of delays: in the consensus problem literature, delays are commonly modelling either the time needed by the network agents to process the information (for instance, the reaction time), or the transmission delay (i.e. the necessary time for the information to be transmitted among the agents). In our model, the delays are conceptually different by their anticipatory nature. Moreover, the proposed algorithm is simple: it does not involve sophisticated techniques of predictive or adaptive control, when compared to other results in the literature. Even better, it turns out that our model is also efficient when a good choice of parameters is made. More precisely, we use two parameters describing how far in the past the agents should look to estimate future states. One other contribution of Chapter 3 is that we derive (necessary and sufficient) conditions on these two parameters, and we study the delay margin for consensus reaching. The results we find have multi-faceted interpretations. The first important interpretation is that the anticipation can improve performance in networks of interacting agents, but when the agents try to anticipate too far in the future the system can diverge. This interpretation is something that we would have expected: using very old information can destabilize the system. However, a second important interpretation is that, in the case of directed networks, there is a non-trivial link between how far back in the past the agents should look and the network topology. This is another difference when compared to other consensus protocols in the literature: the condition for consensus reaching in the case of undirected networks is necessary and sufficient, whereas for directed networks it is not sufficient. We give simulation results from networks with various topologies to show that the consensus convergence speed can be considerably improved when using our algorithm (without changing the network topology, as in other results from the literature). This improvement is especially notable in networks with poor connectivity. One important contribution is that we find a formula for parameters choice that guarantees ultrafast convergence speed for any type of networks, when using the proposed algorithm. Moreover, we show that this choice is a local optimum in terms of network performances, as it can guarantee a convergence speed corresponding to an algebraic connectivity equal to 1. Additionally, we also discuss some interesting observations on directed circulant networks' eigenvalues: we consider some circulant topologies and show that their Laplacian eigenvalues follow a specific pattern, depending on the topology type. A part of these results has been published and presented at two international conferences [START_REF] Atay | A delayed consensus algorithm in networks of anticipatory agents[END_REF][START_REF] Irofti | On the delay margin for consensus in directed networks of anticipatory agents[END_REF], and another part is in preparation for a further submission.

If we focus in the first part of this thesis on connected networks with fixed topology and with normalized Laplacian, we also consider a theoretical study that can correspond, for instance, to switching topology networks, with normalized Laplacian. When considering such a network, we can end up with a disconnected graph. The number of connected subgraphs in a disconnected graph corresponds to the number of Laplacian zero eigenvalues.

This motivates the theoretical study on the bound of maximal multiplicity at the origin in systems with time-delays, presented in Chapter 4. This study exceeds the networks context, and the results can apply to a rather general class of systems with more than one time-delay (but still, with a finite number of time-delays). The main contribution relies on a derived formula for the theoretical upper bound of algebraic multiplicity of the zero characteristic root. Moreover, when this bound is reached, algebraic conditions on a system's parameters are given. We compare our results to other results in the literature. This comparison has two implications. First, it points out that the bound we find is more precise and it cannot be exceeded. Second, the methods we use to derive the bound formula are conceptually different (we basically use linear algebra techniques, and not complex analysis) when compared to the other result. This work has been published in a book chapter [START_REF] Irofti | On the Codimension of the Singularity at the Origin for Networked Delay Systems[END_REF] and published and presented at an international conference [START_REF] Boussaada | Computing the codimension of the singularity at the origin for delay systems in the regular case: A vandermondebased approach[END_REF].

If the Laplacian is not normalized, it can even have imaginary non-zero eigenvalues. We address in Chapter 5 a more general problem, which can also concern networks, of double imaginary characteristic roots in time-delay systems with two constant delays. This study partially solves some issues related to time-delay systems considered to be an open problem since 2005. We focus on studying the multiple imaginary characteristic roots (of multiplicity two, three and four) behaviour under a small perturbation of parameters.

The parameters can be the two delays, but also a delay and another parameter, or two other system parameters. One important contribution of this work is that the results can be extended and applied to a quite general class of systems, with or without delays, even with distributed time-delays. The techniques involved rely on what is called a geometric approach in the literature, whereas other results addressing this problem in the time-delay systems community made use of Puiseux series. Another important contribution is that we derive an algebraic criterion to decide whether the pair of complex conjugated double roots move towards stability or instability under a small deviation of parameters. In other words, we propose a criterion involving a limited number of relatively simple computations (only partial derivatives of the quasi-polynomial with respect to the parameters, up to the third order). We illustrate the proposed method on various examples, and we even show how to extend this study to the case of imaginary characteristic roots of multiplicity three and four. Moreover, for the double non semisimple roots case, we also compare the geometric approach we use to so-called algebraic approach. We verify that the two approaches give coherent results, and discuss their advantages and limitations, but also the open questions they raise. These results have been published and presented at three international conferences [START_REF] Gu | Migration of double imaginary characteristic roots under small deviation of two delay parameters[END_REF][START_REF] Irofti | Migration of imaginary roots of multiplicity three and four under small deviation of two delays in time-delay systems[END_REF][START_REF] Irofti | Geometric vs. algebraic approach: A study of double imaginary characteristic roots in time-delay systems[END_REF], and part of the

Thesis outline

We give a set of definitions and notations used throughout all other chapters of this Chapter 2

Definitions, notations and prerequisites

Throughout this chapter we introduce a set of definitions, concepts and notations relevant for other chapters of this thesis. More precisely, we give a non-exhaustive, but necessary background for graph theory, time-delay systems and synthetic biology, and discuss preliminary concepts used throughout the following chapters of the thesis.

Graph theory

We briefly review some relevant concepts and notations from graph theory. For details, the reader is referred to standard texts such as [START_REF] Fan | Spectral graph theory[END_REF] or [START_REF] Godsil | Algebraic Graph Theory[END_REF].

A graph G can be defined by using a vertex set V and an edge set E. We view the edges as a set of connections between the vertices. Graphically, a graph can be represented by drawing a point for each vertex in the V set, and a curve joining the endpoints for each edge in the E set. The corresponding notation is G = (V, E). A finite graph is a graph where the V and E sets are finite. In this thesis, we consider finite graphs. A path is a simple graph whose vertices can be ordered such that two vertices are adjacent if and only if they are consecutive in the ordering. A cycle is similar to a path, but where a cyclic order on the vertices can be imposed such that the first vertex is also the last vertex. A subgraph of a graph is also a graph defined by a subset of the larger graph vertex set and a subset of the larger graph edge set. Note that the paths and the cycles can be subgraphs of a given graph. We say that a graph edges, and define a cycle as a non-trivial circuit in which the only repeated vertex is the first/last one. Figure 2.4 depicts a summary of these definitions: a walk is a sequence of alternating vertices and edges; the trails set is a subset of the walks set, where no edge is repeated; paths are a subset of the trails set (no vertex is repeated); the circuits represent a subset of the trails (they are closed trails) and don't intersect the path set because at least one vertex is repeated; finally, cycles are a subset of the circuits, and don't intersect the paths subset, as they have only one repeated vertex. If a graph has any two of the following properties, it has all three.

G = (V, E) is connected if for every u, v ∈ V there exist a path from u to v in the graph. Otherwise, G is called disconnected. A walk is a list v 0 , e 1 , v
1. the graph is connected.

2. the graph has no cycles.

|E|= |V |-1.

Therefore, any graph with any two of these properties is called a tree. Examples of trees are depicted in Depending on whether the edges are directed or not, we can have undirected graphs, oriented graphs and directed graphs (also called digraph), as illustrated in Figure 2.7.

Note that in an undirected graph an edge is an unordered pair of vertices. An ordered pair of vertices is called a directed edge. Edges in an oriented graph are directed edges that can be traversed on only one direction. For directed graphs, the directed edges can be traversed in one or two directions.

A graph G is called a bipartite graph (or a bigraph) if its vertices set can be decomposed into two disjoint sets such that no two vertices within the same set are adjacent. Some examples of bipartite graphs are given in Figure 2.8. A circulant graph is a graph of n vertices in which the i th vertex is adjacent to the (i + j) th and (i -j) th vertices for some fixed values of j. Note that if j takes values in {1, 2, . . . , n} set, then the graph is a complete graph.

Time-delay systems

The delay concept models a physical transfer phenomenon, i.e. the transport or propagation of energy, information or even material that take place when an interaction between two physical systems occurs. The time-delay systems are also called systems with aftereffect or dead-time of hereditary systems [START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF]. As many dynamic processes contain this very aftereffect phenomenon, the worldwide engineers present a great interest for investigating time-delay systems, provided that they need their models to fit as much as possible with the real processes.

This type of systems' modelling applies to many scientific areas like physics, mechanics, biology, chemistry, economics or population dynamics, but there is a particular interest for time-delay systems in the control engineering field.

Modeling

Time-delay systems of retarded type belong to the class of functional differential equations (FDE) which are infinite dimensional, opposing to ordinary differential equations (ODE), which are finite. In physical processes modelling, a classical hypothesis assumed is that we can predict the future behavior of a deterministic system by taking into account the present state of the system. But to some extent we can suppose that the past states also affect the system future responses.

Thus, instead of an n-dimensional vector moving through Euclidean space R n , in the case of time-delay systems of retarded type the state vector is a function x t corresponding to a past time-interval [t -τ, t] and not to a single discrete value of time t (as in the case of ODEs).

A more rigorous mathematical approach for describing time-delay systems is by aid of delay differential equations (DDE). This is a type of differential equation where the derivative of the unknown function depends not only on present time and the unknown variable at the present time, but also on the values of the function at previous times.

Generally, a delay differential equation is written as

d dt x(t) = f (t, x(t), x t ), (2.1)
where x(t) ∈ R n and x t = {x(τ )|τ ≤ t} is the trajectory of the solution in the past.

Considering this general form, f is a functional operator mapping R × R n × C into R.

We denote by C the Banach space of continuous functions [-τ, 0] → C n with norm

φ = sup -τ ≤θ≤0 |φ(θ)|, where |•| is any norm in C n .
The representation (2.1) corresponds to the case of time-delay systems of retarded type.

Even if it is beyond the scope of this thesis, we briefly mention the time-delay systems of neutral type that involve a more complex mathematical approach, where the right-hand side also contains the same highest derivative order for some components of x(t) at both present time t or past times.

d dt x(t) = f (t, x t , ẋt , u t ), where, if x : [-τ, α] → C n , α > 0, is a continuous function, then x t , u t ∈ C, 0 ≤ t ≤ α are defined by x t (θ) = x(t + θ) u t (θ) = u(t + θ) , for -τ ≤ θ ≤ 0.
In the following paragraphs, we consider DDEs in which x(t) and x(t-τ i ) appear linearly, giving time-delay systems of retarded type written as

ẋ(t) = A 0 x(t) + m i=1 A i x(t -τ i ), (2.2)
where the state variable at time t, x(t), has the dimension n and A i ∈ R n×n , i = 0, m are real matrices. Generally, the time-delays τ i are constants, but there is also a great interest in systems with time-dependent delays τ i (t) or state-dependent delays τ i (t, x(t)).

Another case is the case where the delays are distributed over some delay interval.

Throughout this thesis, we simply suppose that the problem is non-singular, meaning that the delays are constant and bounded by a maximal positive constant. Thus, in the equation (2.2) we can consider that the time-delays τ i , i = 0, m are distinct positive numbers, sorted in ascending order (0

< τ 1 < τ 2 < • • • < τ m ).
This case is studied into more detail in Chapter 4.

A general form for modelling linear, time-invariant (LTI) systems containing delays, as written in [START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF], is

ẋ(t) = q l=1 D l ẋ(t -τ D l ) + m i=0 (A i x(t -τ i ) + B i u(t -τ i )) + r j=1 t t-τ G j (G j (θ)x(θ) + H j (θ)u(θ))dθ,
where the matrices D i corresponds to the neutral part, τ 0 = 0, A 0 represents the instantaneous feedback gain, the matrices A i , i = 1, m represent discrete-delay phenomena and the sum of integrals corresponds to distributed-delay effects, weighted by the matrices G j , j = 1, r over the time intervals [t -τ G j , t]; B i and H j are input matrices, and u(t) is the control signal.

Next we discuss the initial values problem and the spectral proprieties w.r.t. the equation (2.2), i.e. we consider the retarded case (D l = 0), where the time-delays are point-wise or discrete, meaning that the memory effect is "selective" (G j = 0). Note that there are no input matrices (B i = 0 and H j = 0).

The initial values problem

Unlike the ODEs, in the case of DDEs the initial values x(t) is not enough so determinate the solution for t ≥ 0. It is necessary to specify a history function, also called a segment function, i.e. it is necessary to know x(t) for the interval τ ≤ t ≤ 0 for the DDE to be defined for the interval 0 ≤ t ≤ τ . Thus, the initial condition for a time-delay system of retarded type described by (2.2) is a history (or segment) function φ belonging to a Banach space of continuous functions mapping [-τ m , 0] into R n , equipped with a supremum norm defined as [START_REF] Michiels | Stability and Stabilization of Time-Delay Systems (Advances in Design and Control)[END_REF].

f (x) s = sup{|f (x)| : x ∈ C([τ m , 0], R n )}.
We consider the general form of a DDE ẋ(t) = f (t, x t ).

(2.3)

A function x is said to be a solution on the interval [σ -τ, σ + a] of the (2.3), where f is defined on a subset Ω of R × C with values in the n-dimensional state space,

f : Ω ⊂ R × C → R n , if there exist σ ∈ R and a > 0 such that conditions • x ∈ C([σ -τ, σ + a], R n ), • (t, x t ) ∈ Ω, • x(t) satisfies (2.3) for t ∈ [σ -τ, σ + a]
hold. A solution x(σ, φ, f ) of (2.3) with initial values φ at σ is called a solution through (σ, φ) if there is a positive a such that the following conditions are satisfied:

• x(σ, φ, t) is a solution of (2.3) on t ∈ [σ -τ, σ + a], a > 0, • x σ (σ, φ, f ) = φ. Supposing that Ω is an open subset of R × C and f is a continuous function, f ∈ C(Ω, R), we call x ∈ C([σ -τ, σ + α], R n ) with α > 0 a backward continuation of the solution through (σ, φ) if • x σ = φ and • (σ 1 , x σ 1 ) ∈ Ω and x is a solution of (2.3) on [σ 1 -τ, σ] through (σ 1 , x σ 1 ) for any σ 1 [σ -α, σ], α > 0 [67].
Concerning the techniques involved for solving a system of DDEs, usually it is possible to reduce the DDE to a sequence of ODEs and to adapt the methods used for solving the initial values problem of ODEs. The method of steps [START_REF] Michiels | Stability and Stabilization of Time-Delay Systems (Advances in Design and Control)[END_REF], for example, is used to build the solution step by step beginning with a given initial condition. The solution on [0, τ ] (σ = -τ and a = 2τ ) is computed using the function φ defined on [-τ, 0], by solving the ODE

ξ1 (t) = A 0 ξ 1 (t) + m i=1 A i φ(t -τ i ),
with the initial condition ξ 1 (0) = x(0) = φ(0). The next iterations are based on the same idea: at the k th step the evolution of (2.2) is determined by solving the ODE

ξk (t) = A 0 ξ k (t) + m i=1 A i ξ k-1 (t -τ i ), with the initial condition ξ k ((k -1)τ ) = x((k -1)τ ) = ξ k-1 ((k -1)τ ).
Finally a unique, globally defined forward solution of (2. This method guarantees by construction the existence and uniqueness of solutions. However, the existence of solutions is commonly proven using fixed point theorems [START_REF] Michiels | Stability and Stabilization of Timedelay Systems. An Eigenvalue Based Approach[END_REF].

Similar to ODEs, we can say that there is a solution as defined in (2.

Spectrum properties

The solutions of (2.2) are based on the characteristic equation

det M (s) = 0, (2.4)
where M (s) is the characteristic matrix, defined by

M (s) := sI -A 0 - m i=1 A i e -sτ i . (2.5)
The spectrum σ(A) corresponds to the transcendental equation (2.4) that has an infinite number of roots (called as in the finite-dimensional case characteristic roots), meaning that #σ(A) = ∞.

Similarly to the finite-dimensional case, a complex number s as an eigenvalue of A has an algebraic multiplicity (the multiplicity as a root of the characteristic matrix M (s))

and a geometric multiplicity (the dimension of the nlsull space of M (s)). Moreover, an eigenvalue is called simple if the corresponding algebraic multiplicity is one and multiple otherwise. More precisely, a multiple eigenvalue is called semi-simple if its algebraic multiplicity is (greater than 1 and) equal to its geometric multiplicity, and non semisimple otherwise.

Like in the finite-dimensional case, an eigenvalue of A, s ∈ σ(A), has a corresponding eigenvector that can be either a right eigenvector u ∈ C\{0} corresponding to a right eigenfunction ue sθ , with θ ∈ [-τ, 0], such that M (s)u = 0, or a left eigenvector v = 0 corresponding to a left eigenfunction, such that

v * M (s) = 0. (2.6)
In other words, the maximal number of linearly independent eigenvectors corresponding to a multiple eigenvalue is called geometric multiplicity of the eigenvalue.

If there is a single eigenvector corresponding to a multiple eigenvalue, then the eigenvalue is called nonderogatory. For instance, the nonderogatory case of zero eigenvalue, s = 0, with algebraic multiplicity two and geometric multiplicity one is known as Bogdanov-Takens singularity, and it has been extensively studied in the literature, including in the models with time delays (see [START_REF] Bogdanov | Versal deformations of a singular point of a vector field on the plane in the case of zero eigenvalues[END_REF], [START_REF] Bogdanov | Bifurcations of a limit cycle for a family of vector fields on the plane[END_REF], [START_REF] Boussaada | Computing the codimension of the singularity at the origin for delay systems in the regular case: A vandermondebased approach[END_REF], [START_REF] Campbell | Zero singularities of codimension two and three in delay differential equations[END_REF], [START_REF] Faria | On the study of singularities for a planar system with two delays. Dynamics of Continuous[END_REF], [START_REF] Takens | Forced oscillations and bifurcations[END_REF]).

We denote by quasi-polynomial the left-hand side of equation (2.4), i.e. a function of the form

q(s) = n i=1 p i (s)e τ i s , (2.7)
where for 1 ≤ i ≤ n, p i (s) are polynomials of s with complex coefficients and τ i are complex numbers, also called the delays of the quasi-polynomial. Supposing that the delays are distinct and the complex polynomials are not identically null, we call the degree of the quasi-polynomial q(z) (see also [START_REF] Wielonsky | A rolle's theorem for real exponential polynomials in the complex domain[END_REF]) the number:

deg q(s) := n -1 + n i=1 d i ,
where d i , for 1 ≤ i ≤ n represents the degree of p i (s) polynomial:

d i = deg p i (s).
However, this thesis discusses real quasi-polynomials, that is the delays τ i are real numbers such that a ik s k , with 1 ≤ i ≤ n and ∀ k = 0, . . . , deg(p i ) the condition a ik = 0 holds. Several other definitions can be encountered in the literature describing sparsity. Among others, we mention the lacunary polynomials [START_REF] Marden | Geometry of Polynomials[END_REF].

0 ≤ τ 1 < τ 2 < • • • < τ n ,
When referring to spaces and subspaces, (but also manifolds and sub-manifolds, or algebraic varieties and subsets) we define the codimension of a subspace Y in a space X as the difference between its dimensions:

codim(Y ) = dim(X) -dim(Y ).
The following definition is important, as it plays a crucial role in Chapter 5.

Definition 2.1. We say that the characteristic equation (2.4) (thus, the quasi-polynomial (2.7)) has a multiple root s = s 0 of multiplicity m if

∂ k p ∂s k s=s 0 = 0, for k = 0 . . . m -1 ∂ m p ∂s m s=s 0 = 0.
In the sequel, we provide some examples illustrating the definitions and prerequisites above. More precisely, we consider three finite dimensional cases where s 1 = 1 is a simple, non semi-simple, respectively, semi-simple root.

Example 2.1 (s 1 = 1 simple eigenvalue). We consider the matrix

A = 0 -1 2 3
.

Its characteristic polynomial (the quasi-polynomial) is

q(s) = det(A -sI 2 ) = -s -1 2 3 -s = (s -1)(s -2).
Note that the root s 1 = 1 has the algebraic multiplicity 1. Next, we find the eigenvectors corresponding to s 1 = 1, by writing

-s 1 -1 2 3 -s 1 u 1 u 2 = 0 0 .
Therefore, the root s 1 = 1 has the algebraic multiplicity 1 and the geometric multiplicity also 1, since there is only one linearly independent (right) eigenvector, which is

u = 1 -1 .
Example 2.2 (s 1 = 1 non semi-simple eigenvalue). Given the matrix

A = 1 1 0 1
, we solve the characteristic equations to find the characteristic eigenvalues and eigenvectors. The characteristic polynomial can be written det(A -sI 2 ) = (1 -s) 2 ; thus, s 1 = 1 has the algebraic multiplicity 2. There is only one linearly independent (right) eigenvector corresponding to this eigenvalue, which is u = 1 0

. The geometric multiplicity is 1, different from the algebraic multiplicity, so in this case s 1 = 1 is a non semi-simple root.

Example 2.3 (s = 1 semi-simple eigenvalue). Consider that the matrix A is the identity matrix

A = 1 0 0 1 .
The characteristic polynomial is the same as in Example 2.2, but here we find two linearly independent (right) eigenvectors corresponding to s 1 = 1. This means that the geometric multiplicity of this characteristic root is 2, equal to its algebraic multiplicity. That being the case, we call s 1 = 1 a semi-simple root.

Geometric approach

Control systems often depend on parameters, such that their characteristic equation may be written as q(s, p) = 0, (

where s is the Laplace variable and p ∈ R n is a vector of n parameters. We may have parameters due to internal dynamics. For instance, modelling in physical, biological or social sciences sometimes requires taking into account the time delays inherent in the phenomena. Depending on the model complexity, but also on how much information is known, we may choose a model with continuous constant delays, or a model with distributed delays [START_REF] Cushing | Time delays in single species growth models[END_REF][START_REF] Macdonald | Biological delay systems[END_REF]. For instance, in the case of a time-delay system with two constant delays, the characteristic equation is written of the form

q 1 (s, τ 1 , τ 2 ) = r 0 (s) + r 1 (s)e -τ 1 s + r 2 (s)e -τ 2 s , (2.9)
where r k (s), k = 0, 1, 2 are polynomials of s with real coefficients, and the delays τ 1 , τ 2 are the two parameters.

Also common is the case when p contains controller parameters. Classical example are PI, PD and PID controllers, where the continuous time domain controller is expressed in the Laplace domain as

q 2 (s) = K P 1 + 1 T i s , (2.10)
q 3 (s) = K P (1 + T d s) , (2.11)
q 4 (s) = K P 1 + 1 T i s + T d s , (2.12)
respectively, where K P is the proportional gain, T i the integral time, and T d the derivative time. Furthermore, in many practical applications, a time delay of the process model τ m may be involved [START_REF] O'dwyer | Pi and pid controller tuning rules: an overview and personal perspective[END_REF][START_REF] Morȃrescu | Stability crossing boundaries and fragility characterization of pid controllers for siso systems with i/o delays[END_REF]. These include proportional plus delay (2.13), integrator plus delay model (2.14), first order lag plus delay (2.15), first order lag plus integral plus delay (2.16) expressed below:

q 5 (s) = K m (1 + e -sτm ) (2.13) q 6 (s) = K m e -sτ s (2.14) q 7 (s) = K m e -sτ 1 + sT m (2.15) q 8 (s) = K m e -sτm s(1 + sT m ) (2.16)
If in equation (2.13) there are two different proportional gains, then we obtain the model of a proportional retarded controller:

q 9 (s) = K p + K r e -sτm .
(2.17)

Furthermore, [START_REF] Villafuerte | Tuning of proportional retarded controllers: Theory and experiments[END_REF] showed that proportional retarded controller outperforms a PD controller on an experimental DC-servomotor setup. Obviously, any control among (2.10) to (2.17) results in a characteristic equation that depends on the control parameters.

Many studies have been conducted on the stability of systems that depend on parameters. For example, for systems with a single delay as the parameter, methods of identifying all the stable delay intervals are given in [START_REF] Lee | On the tau-decomposition method of stability analysis for retarded dynamical systems[END_REF] and [START_REF] Walton | Direct method for tds stability analysis[END_REF]. For system with two parameters, a rich collection of stability charts (the parameter regions where the system is stable) for time delay systems are presented in [START_REF] Stépán | Retarded dynamical systems: stability and characteristic functions[END_REF]. For systems with two delays as the parameters, a geometric approach is introduced in [START_REF] Gu | On stability crossing curves for general systems with two delays[END_REF]. The analysis is based on the continuity of the characteristic roots as functions of parameters (which needs to be carefully evaluated in the case of time delay systems of neutral type [START_REF] Michiels | Stability and Stabilization of Time-Delay Systems (Advances in Design and Control)[END_REF][START_REF] Gu | A review of some subtleties of practical relevance for time-delay systems of neutral type[END_REF]), and consists of identifying the parameters that correspond to imaginary characteristic roots and judging the direction of crossing of these roots as parameters change. Such an analysis is often known as D-subdivision method [START_REF] Gryazina | D-decomposition technique stateof-the-art[END_REF].

Challenges due to non-differentiability arise when the imaginary roots concerned are multiple roots. Such problems have traditionally been solved using Puiseux series [START_REF] Kato | Perturbation Theory for Linear Operators[END_REF][START_REF] Knopp | Theory of functions : parts I and II[END_REF], see, for example, [START_REF] Chen | An eigenvalue perturbation approach to stability analysis, part i: Eigenvalue series of matrix operators[END_REF][START_REF] Chen | An eigenvalue perturbation approach to stability analysis, part ii: When will zeros of time-delay systems cross imaginary axis[END_REF][START_REF] Li | On computing puiseux series for multiple imaginary characteristic roots of lti systems with commensurate delays[END_REF] for systems with one parameter.

In Chapter 5, we study systems with two delays, and the extension to two parameterdependent systems can be found in Chapter 7. More precisely, we present a method to analyze the migration of roots in a neighborhood of the parameter corresponding to a multiple imaginary characteristic root. The method of analysis uses traditional complex analysis, and does not require Puiseux series. A preliminary version of this method, which is restricted to the case of two delays as the parameters, was presented in [START_REF] Gu | Migration of double imaginary characteristic roots under small deviation of two delay parameters[END_REF].

Synthetic biology

Synthetic biology is a relatively recent interdisciplinary branch of biology and engineering. Even if its name has been first used at the beginning of twentieth century [START_REF] Le | Théorie physico-chimique de la vie et générations spontanées[END_REF],

the first description close to our understanding of synthetic biology today was given by Polish geneticist Wac law Szybalski in 1974 [START_REF] Kohn | Control of Gene Expression[END_REF]. However, defining synthetic biologic is a debatable topic even nowadays. Still, it is commonly accepted that the central aim of synthetic biology is to design new biological components and systems that do not already exist in the natural world, but also to re-design existing biological systems.

As mentioned in chapter 2 of [START_REF] Freemont | Synthetic Biology -A Primer[END_REF], an important aspect of synthetic biology is the systematic design based on engineering principles such as modularity, characterization and standardization. The modularity means the reduction of a biological device or system in its component parts, also called bioparts. Each part need to be characterized in detail.

This means for instance that a biopart behaviour has to be rigorously studied under different conditions. The standardization means that these parts should be build so that they can be easily combined into devices with well defined functions (e.g. logic gates).

Thus, new devices can be created from existing standard parts. Another concept widely used in synthetic biology, as in control theory is referring to systems. One example of a system realization in synthetic biology field is a counter [START_REF] Smolke | It's the dna that counts[END_REF], but researchers aim to be able to build intracellular control functions in the future.

In order to engineer new biological systems, it is important not only to understand how information flows in biological systems, but also how the information flow is controlled.

Biology has highly complex aspects, rather different when compared to engineering complexity issues. A biological system is nevertheless different from an electrical circuit: it can be subject to modifications based on natural selection, which is not the case of and T (thymine). In order to understand how the information flows in biological systems, we need to know how the two strands interact. Geneticists know that each of the four bases can only pair with one of the others. This gives a complementary feature that allows the cell to always have two copies of the information stored by the DNA. Thus, in case one strand is damaged, no information is lost.

The information flow in biological system is usually represented by the central dogma, as depicted in Figure 2 A smart consensus algorithm for multi-agent systems

Chapter overview

In this chapter we introduce and analyze a delayed consensus algorithm as a model for interacting agents using anticipation of their neighbors' states to improve convergence.

The contribution of the proposed model mainly relies on the different nature of the delays used in the model, namely the anticipatory nature: every agent in the network uses past information to estimate future states of its neighbours. Other models in the literature include delays coming form other sources, as information transmission or information processing delays. We analyze the proposed model for directed and undirected networks of all type of fixed topologies, and derive conditions for the systems to reach consensus.

We explicitly calculate the dominant characteristic root of the consensus problem as a measure of the speed of convergence.

Another important contribution of this chapter is that we illustrate how choosing the model parameters can be critical for consensus reaching. We give simulation results from networks with various topologies and show that the speed of consensus can be dramatically improved, especially in networks with poor connectivity. Moreover, we show how the two parameters of the proposed model can be chosen to guarantee an ultrafast convergence speed for any type of network, and prove that this choice is even a local optimum in terms of network performances.

Part of the research presented in this chapter was carried out during a visit to Max Planck Institute for Mathematics, Leipzig, Germany, during the summer of 2015, under financial support from DAAD Research Grant. Part of the work (the case α = 1) has been published and presented at two international conferences [START_REF] Atay | A delayed consensus algorithm in networks of anticipatory agents[END_REF][START_REF] Irofti | On the delay margin for consensus in directed networks of anticipatory agents[END_REF], and another part of the work (the general case, where α = 1, together with Sections 3.5 and 3.6) are in preparation for a further submission.

Problem statement

Consensus problems arise in a wide range of applications in distributed computing [START_REF] Lynch | Distributed Algorithms[END_REF],

management science [START_REF] Degroot | Reaching a consensus[END_REF], flocking and swarming theory [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF], distributed control [START_REF] Fax | Information flow and cooperative control of vehicle formations[END_REF],

and sensor networks [START_REF] Olfati-Saber | Consensus filters for sensor networks and distributed sensor fusion[END_REF], among others. In these applications, multi-agent systems interact to agree on a common value for a certain cooperative behaviour. The interaction rule that specifies the information exchange between the agents is called the consensus protocol or consensus algorithm. The consensus problem on networks in continuous time can be formulated as

ẋi = u i (t), i = 1, . . . , n, (3.1)
where n is the number of agents in the network, x i ∈ R is the state of the agent i at time t, which changes under the interaction with other agents as described by the consensus protocol u i (t).

Definition 3.1. The system (3.1) is said to reach consensus if for any set of initial

conditions {x i (0)} there exists c ∈ R such that lim t→∞ x i (t) = c
for all i, in which case the number c is called the consensus value.

In the classical linear case, the consensus protocol typically has the form

u i (t) ∼ n j=1 a ij (x j (t) -x i (t)). (3.2)
where a ij are nonnegative numbers describing the strength of the influence of agent j on agent i. It can be shown, under certain conditions regarding the network connectivity, that the system (3.1) under the protocol (3.2) reaches consensus from arbitrary initial conditions. There are several consensus protocols studied in the literature, some of them including delays from different sources (transmission delays, information processing delays).

In this chapter, we study another delayed consensus protocol where the delay comes from a quite different source, namely from the anticipatory nature of the agents. In other words, we consider a network of intelligent agents who try to anticipate the future states of their neighbors in their interaction, which is a common situation, for instance, in economic systems.

We therefore consider an anticipatory algorithm of the form

u i (t) ∼ n j=1 a ij (x j (t + δ) -x i (t)), (3.3)
where xj (t + δ) is the anticipated state of the neighbor x j at some time (t + δ) in the future. In other words, every node in the network will adjust its state by anticipating its neighbors' states, based on the present and past states. Note that a ij = a ji in the case of undirected networks, i.e. the communication is bidirectional between every two interconnected nodes in the network in this case.

Obviously, the future state of a node is an unknown, and we need a mechanism to estimate it. One possible solution for the agents is to use a first order estimation derived from past observations (we say that they employ memory), as follows. Agent i, knowing the current state x j (t) of a neighbor j and remembering its past state x j (t -τ ) as well, uses a first-order estimation to get the future state xj (t + δ) by linear extrapolation:

xj (t + δ) = x j (t) + x j (t) -x j (t -τ ) τ δ = (1 + α)x j (t) -αx j (t -τ ), (3.4)
where

α = δ τ . (3.5)
The idea is graphically summarized in Figure 3.1. By comparison, the classical consensus algorithm (3.2) can be viewed as a zero-order estimation where an agent's expectation of the short-term future of its neighbors is represented simply by the present states,

xj (t + δ) = x j (t).
Substituting (3.4) into (3.3), we obtain the consensus algorithm

u i (t) ∼ n j=1 a ij [(1 + α)x j (t) -αx j (t -τ ) -x i (t)], (3.6)
on which we focus in the sequel. More precisely, in this chapter we study the convergence of the following consensus problem

ẋi (t) = 1 d i n j=1 a ij [(1 + α)x j (t) -αx j (t -τ ) -x i (t)], (3.7)
time

x j t -τ t t + δ
first-order estimate zero-order estimate where a ij ≥ 0 ∀i, j and d i = n j=1 a ij is the (generalized) degree of node i. Dividing the summation by the degrees d i gives rise to a normalized Laplace operator, which is a natural choice in several applications and has some advantageous properties, as will be briefly reviewed in Section 3.3. In particular, the normalization bounds the spectrum of Laplacian regardless of the network size, thus allowing comparison of networks of very different sizes. Equation (3.7) is a delay differential equation; however, the source of the delay and the form of the resulting equation is quite different here from the existing literature. For example, in the well-known delayed consensus protocol considered in [START_REF] Olfati | Consensus problems in networks of agents with switching topology and time-delays[END_REF], namely

u i (t) ∼ n j=1 a ij (x j (t -τ ) -x i (t -τ )) , (3.8)
the delay can be viewed as modeling delayed information processing, since the protocol feeds back the same information (x j -x i ) to the system, as in (3.2), but only after a delay τ ≥ 0. In this case, it is known that there exists an upper limit τ max such that the system (3.1) under the protocol (3.8) reaches consensus from arbitrary initial conditions if and only if τ < τ max [START_REF] Olfati | Consensus problems in networks of agents with switching topology and time-delays[END_REF]. Another type of interaction, which models delayed information transmission, is given by the consensus protocol [START_REF] Atay | The consensus problem in networks with transmission delays[END_REF][START_REF] Seuret | Consensus under communication delays[END_REF][START_REF] Moreau | Stability of continuous-time distributed consensus algorithms[END_REF][START_REF] Atay | Consensus in networks under transmission delays and the normalized Laplacian[END_REF]]

u i (t) ∼ n j=1 a ij (x j (t -τ ) -x i (t)) . (3.9)
Here the interpretation is that the information coming from a neighbouring node j takes some time τ to reach site i. It was shown that the system (3.1) under the protocol (3.9) reaches consensus from arbitrary initial conditions regardless of the value of the delay τ as long as the network contains a spanning tree [START_REF] Atay | The consensus problem in networks with transmission delays[END_REF][START_REF] Atay | Consensus in networks under transmission delays and the normalized Laplacian[END_REF].

As mentioned above, in both (3.8) and (3.9), the delay τ is imposed by system constraints. By contrast, the delay in (3.7) comes from a quite different source, namely from the anticipatory nature of the agents. As such, τ in (3.7), together with α, is to be seen as a design parameter. Note that the anticipative consensus protocol (3.6) reduces to the classical undelayed case (3.2) by setting either α = 0 or τ = 0, whereas putting α = -1 yields the consensus protocol under signal transmission delay (3.9).

In this chapter, we consider (3.7) for general α, τ ∈ R + . We investigate if and under what conditions the system (3.7) can reach consensus from arbitrary initial states, and furthermore, whether the introduction of a positive τ , namely anticipation, really brings any advantages into the consensus dynamics. Nevertheless, we note that system (3.7)

is still described by delay-differential equations and its analysis is subject to the same difficulties one faces in studying infinite-dimensional systems on Banach spaces.

Preliminaries

We briefly present some concepts and notations from graph theory, relevant for this chapter. For more details, the reader is referred to standard texts such as [START_REF] Fan | Spectral graph theory[END_REF] or [START_REF] Godsil | Algebraic Graph Theory[END_REF].

Let G = (V, E) be a directed graph with a finite set V of vertices and a set of directed edges E ⊂ V × V consisting of ordered pairs of vertices. We denote by A = [a ij ] the (weighted) adjacency matrix of the graph, where a ij > 0 if there is a link from j to i, and a ij = 0 otherwise. For a network with n agents, A is a square matrix of dimension n and a ij describes the strength of the influence of agent j on agent i. Note that a ij = a ji ∀i, j in the case of undirected networks, which means that A is a symmetric matrix.

The in-degree d i of node i is defined as d i = n j=1 a ij , i.e., the sum of the elements of the i th row of A, and D = diag(d 1 , . . . , d n ) denotes the diagonal degree matrix. We consider simple non-trivial graphs without self-loops; in particular, d i is assumed to be nonzero for all i.

The normalized Laplacian matrix is defined as

L = I -D -1 A, (3.10)
where I denotes the identity matrix. The normalized Laplacian naturally arises in a class of important problems, in particular in random walks on networks, as D -1 A is the transition matrix for probability distributions arising from such walks [START_REF] Fan | Spectral graph theory[END_REF]. The eigenvalues {λ i } of L lie in the unit-diameter disc centered at 1 on the complex plane, that is,

|1 -λ k |≤ 1, k = 1, 2, . . . , n, (3.11)
which can be seen by an application of Gershgorin's circle theorem (see e.g. [START_REF] Horn | Matrix Analysis[END_REF]).

The sum of each row of L equals zero by the definition of in-degree d i ; hence, L always has a zero eigenvalue, which we denote by λ 1 , corresponding to the eigenvector

v 1 = (1, 1, . . . , 1) .
The eigenvectors of L form a complete set that span R n , and the eigenvalues λ i are real numbers and can be ordered as

0 = λ 1 ≤ λ 2 ≤ . . . ≤ λ n ≤ 2.
The second eigenvalue λ 2 , which is also called the spectral gap, is positive if and only if the graph is connected. In fact, the multiplicity of the zero eigenvalue equals the number of connected components of the graph. In a connected graph, λ 2 gives an indication of how difficult it is to disconnect the graph into two large pieces by removing a small number of edges, and is thus directly related to graph connectivity.

Remark 3.2. Among undirected networks, only complete graphs have all their eigenvalues (except λ 1 ) greater than 1, i.e. λ 1 = 0 and

λ 2 = • • • = λ n = n n-1 > 1
where n is the number of nodes in the network [START_REF] Fan | Spectral graph theory[END_REF].

In matrix form, (3.7) becomes

ẋ(t) = D -1 A[(1 + α)x(t) -αx(t -τ )] -x(t), (3.12)
with x = (x 1 , x 2 , . . . , x n ) . For simplicity, we make the standing assumption that L has a complete set of eigenvectors {v k }, a condition which generically holds for matrices in R n×n . Then we can write x(t) = n k=1 u k (t)v k for some scalar coefficients u k , which transforms (3.12) into a system of n decoupled scalar equations

uk (t) = [(1 -λ k )(1 + α) -1]u k (t) -α(1 -λ k )u k (t -τ ), k = 1, . . . , n. (3.13)
The characteristic equation corresponding to the eigenmode (3.13) is

ψ k (s) := s -(1 + α)(1 -λ k ) + 1 + α(1 -λ k )e -sτ = 0. (3.14)
Consequently, the characteristic equation for the whole system (3.12) can be written as

Ψ(s) := n k=1 ψ k (s) = 0. (3.15)
The linear equation (3.13) has solutions of the form e st corresponding to characteristic roots s of (3.15). Note that s = 0 is always a characteristic root for the first factor ψ 1 (s) = s -α + αe -sτ corresponding to the first eigenmode, λ 1 = 0. Thus, points on the synchronization subspace spanned by v 1 = (1, 1, . . . , 1) can be at best neutrally stable.

If that is the case, and in addition lim t→∞ u k (t) = 0 for all k ≥ 2, then the system converges to a point on v 1 , i.e. it reaches consensus, from arbitrary initial conditions.

This clearly happens if and only if zero is a simple root of Ψ and all other roots of Ψ have negative real parts. Moreover, in this case the speed of convergence from general initial conditions depends on the slowest of these modes u k , k ≥ 2. Hence, we factor (3.15) into directions along and transverse to the synchronization subspace span(v 1 ) as

Ψ(s) = ψ 1 (s) Ψ(s)
, where

Ψ(s) := n k=2 ψ k (s), (3.16)
and use the transverse part to quantify the speed of convergence, which motivates the following definition. Note that zero is simple root of ψ 1 and all its other roots have negative real parts.

Although s * may not be unique, its real part is unique and essentially determines the speed of convergence to consensus in case it is negative.

When τ = 0, system (3. (H) Zero is a simple eigenvalue of the Laplacian matrix L.

Condition (H) is a requirement on the connectivity of the network. In fact, for undirected networks, i.e., when a ij = a ji ∀i, j, the eigenvalues λ k of L are all real and can be ordered

as 0 = λ 1 ≤ λ 2 ≤ . . . ≤ λ n ≤ 2.
Thus, condition (H) is equivalent to λ 2 > 0, which is a well-known necessary and sufficient condition for the network to be connected. In fact, the second eigenvalue λ 2 , which is called the algebraic connectivity of the graph, is also a quantitative measure of how well the network is connected [START_REF] Fiedler | Algebraic connectivity of graphs[END_REF]. Since the slowest mode e -λ 2 t determines the speed of convergence to consensus, we have a direct relation between consensus speed and graph connectivity in the undelayed/non-predictive case (3.2). Note that λ 2 can be arbitrarily small even if the network is connected, so consensus can be a very slow process depending on the connection structure.

In this chapter we study the conditions under which the predictive consensus algorithm (3.7) can yield much faster convergence independent of the network structure, and we show how to choose appropriate values for the design parameters α and τ to achieve this.

Convergence of the consensus algorithm

We begin this section by some observations on the roots of a certain complex function.

We introduce here some results on the roots of the complex function

ψ(s) = s -β 1 -e -s , β ∈ C, (3.18)
that are used in the proofs of results in subsections 3.4.1 and 3.4.2.

When β is real, ψ(s) becomes a special case of the function

ψ(s) := s -a 1 -a 2 e -s , a 1 , a 2 ∈ R, (3.19)
which has been studied in the classical paper of Hayes [START_REF] Hayes | Roots of the transcendental equation associated with a certain difference-differential equation[END_REF]. The properties for (3.19) are therefore well-known; here we recall the stability region depicted in Figure 3.

β < 1 (3.20) if β ∈ R, or Re(β) < Im(β) cot(Im(β)) (3.21) if Im(β) = 0.
Proof. Clearly ψ(0) = 0; so zero is always a root of ψ. Moreover, explicit curve on the complex plane described by

ψ (0) = 1 -β is nonzero if and only if β = 1,
Re(β) = Im(β) sin(2 Im(β)) 1 -cos(2 Im(β)) = Im(β) cot(Im(β)), (3.22)
which is depicted in Figure 3.3. It can be seen (for instance by taking the limit Im(β) → 0) that this curve encloses the ray R defined above, on which the statement of the lemma was shown to hold. Hence for the values of β described by the condition (3.21), all roots of ψ have negative real parts except for a simple root at zero.

Undirected networks

We can now state a convergence result for the consensus problem (3.7) in undirected networks.

Theorem 3.5. The system (3.7) defined on a connected undirected graph reaches consensus from arbitrary initial conditions if and only if ). We note that the roots s of ψ k satisfy

ατ < 1. ( 3
s = a k + b k e -sτ , (3.24)
with

a k = (1 + α)(1 -λ k ) -1 b k = -α(1 -λ k ). (3.25) (3.26)
For τ = 0, (3.24) reduces to

s = a k + b k = -λ k < 0 for k ≥ 2, ( 3
1 ≤ |b k τ |= |ατ (1 -λ k )|≤ ατ,
where we have substituted for b k from (3.26) and used (3.11). This shows that, as long as condition (3.23) holds, no roots can cross the imaginary axis, and so all roots of ψ k for k = 2, . . . , n have negative real parts. This completes the proof of the theorem.

The performance of the consensus algorithm depends on the transverse dominant root of the problem, as defined in Definition 3.3. We use Lambert W function to solve for the characteristic roots.

Recall that the Lambert W function is defined as the inverse function of the mapping

z → ze z , f (z) = ze z for z ∈ C [103]; in other words, W (z) satisfies W (z)e W (z) = z. (3.29)
Since f is not one-to-one, W (z) is multi-valued. We let W 0 denote the principal branch of the Lambert function.

Proposition 3.6. The root of the characteristic equation (3.14) having the largest real part is given by which can be re-written as

s = 1 τ W 0 (τ b k e -a k τ ) + a k , (3.30
s -a k = b k e -(s-a k )τ e -a k τ . A change of variables s → (s -a k )τ gives s = τ b k e -a k τ e -s ,
from which, by the definition (3.29) of Lambert function, we immediately have s = W (τ b k e -a k τ ). Going back to the original variable s shows that the roots of ψ k are given by

s = 1 τ W (τ b k e -a k τ ) + a k , (3.32)
It follows by a recent result from [START_REF] Shinozaki | Robust stability analysis of linear time-delay systems by Lambert W function: Some extreme point results[END_REF] that the root with largest real part is given by the principal branch of the Lambert function, which establishes (3.30). Then, Definition 3. Note that this result can be depicted for every possible generic value of eigenvalues, λ, as we can solve the characteristic equation of the system in view of (3.32), taking into account different values of the delay τ . In the color map of Figure 3.4, lighter colors correspond to more negative real parts for the roots. Hence we note, from the change in color as we move vertically upwards from the horizontal axis, that a positive value of τ can indeed improve convergence speed for a given eigenmode.

In Section 3.5 we confirm this observation through actual simulation of several networks.

Directed networks

Unlike the case of undirected networks, the condition (3.23) is not sufficient for consensus in directed networks. The distinction clearly lies in the fact that the Laplacian eigenvalues are real for undirected networks but not necessarily for directed ones. It is therefore of interest to have some understanding of the network structures and the corresponding Laplacian eigenvalues that give rise to reduced delay margins 1 . We take up this task in the sequel: we state a necessary condition for consensus, which serves as an upper bound for the delay margin. This necessary condition is given by the following theorem.

Theorem 3.7. The inequality

ατ < 1 (3.33)
is a necessary condition for (3.7) to reach consensus from arbitrary initial conditions.

In other words, the delay margin for the stability of consensus in (3.7) is at most 1 α .

Proof. Consider the first factor

ψ 1 (s) = s -(1 -λ 1 )(1 -α) + 1 + α(1 -λ 1 )e -sτ in (3.15),
where λ 1 = 0. Since ψ 1 (0) = 0 and ψ 1 (0) = 1 -ατ , zero is always a root of ψ 1 and is a simple root unless ατ = 1. For ατ > 1, on the other hand, ψ 1 is unstable; in fact, it has a real and positive root, which can easily be seen by plotting the real functions α -s and αe -sτ and observing that they must intersect at a positive value of s if ατ > 1.

Therefore, (3.33) is a necessary condition that ψ 1 , and hence the characteristic equation (3.15), has a simple zero root and all the remaining roots have strictly negative real parts.

As we have seen in Subsection 3.4.1, the root having the largest real part is given by the principal branch W 0 of the Lambert function. Hence, for the stability of (3.14) it suffices to calculate

s = 1 τ W 0 (τ be -aτ ) + a,
and check whether its real part is negative.

As a first step in numerical investigations, we randomly generate complex numbers λ satisfying

|1 -λ|≤ 1 (3.34)
and check the stability of (3.14) using (3.30) for the range of delay values τ ∈ (0, 1) with α = 1. The results of the experiment are shown in Figure 3.5 for one million randomly generated λ. It can be seen that only a small fraction of λ actually yield an unstable ψ.

To gain further insight into the nature of instability, we plot in Figure 3.6 the location of the stable and unstable Laplacian eigenvalues λ in the complex plane. As expected, for sufficiently small delays consensus is stable, and as the delay increases, two regions of unstable eigenvalues (shown with red color) grow from the circle boundary, eventually meeting at λ = 0 as τ ↑ 1 (

Ultrafast consensus with predictive algorithm

This section puts some light on one of the main advantages of the proposed model. We

show that by well choosing τ and α, we not only avoid instability, but also guarantee a very fast convergence speed for any network. More precisely, the convergence speed corresponds to λ 2 = 1, for τ chosen equal to

τ * := ln 1 + α α . ( 3
     λ k > 1 - 1 (1+α) ln( 1+α α ) if Im(λ k ) = 0 Re(λ k ) > 1 -Im(λ k ) cot Im(λ k )(1 + α) ln 1+α α if Im(λ k ) = 0 . ( 3
Φ k (s) = s -(1 -λ k )(1 + α) + 1 + α(1 -λ k ) 1 + α α -s
.

Evidently, Φ k (-1) = 0, meaning that s = -1 is always a real root of the characteristic equation (3.15). In order to prove the second part of the theorem we consider s = s + 1.

Clearly, Re(s) For λ k with Im(λ k ) = 0, we can rewrite condition (3.20) of the form of the first equation of (3.36). This applies to real eigenvalues of directed networks and to all eigenvalues of undirected networks. As for the directed networks where ∃λ k such that Im(λ k ) = 0, the condition (3.21) of Lemma 3.4 can be written of the form

≤ -1 ⇔ Re(s ) ≤ 0. We write Φ k (s -1) = 0 as s -(1-λ k )(1+α)+α(1- λ k ) 1+α α 1-s = 0, which means s + (1 -λ k )(1 + α) -1 + 1+α α -s = 0. We replace
(1 + α)(1 -Re(λ k ))τ * < -(1 + α) Im(λ)τ * cot (-(1 + α) Im(λ k )τ * ) .
As τ * > 0, we can cancel it without changing the sign. Then, taking into account the odd parity of function cot, we obtain exactly the second relation of (3.36).

λ 2 > 1 - 1 (1 + α) ln( 1+α α ) , (3.37)
then all characteristic roots satisfy Re(s) ≤ -1 (we say that s = -1 is the "dominant root" and it corresponds to τ = τ * ). Re(s) ≤ -1 (we say that s = -1 is the "dominant root" and it corresponds to τ = τ * ).

Proof. Using l'Hôpital rule we remark that lim α→∞

(1 + α) ln 1 + α α = 1. (3.38)
In the case of undirected networks, this means that α > 0 can be as large as we need, since condition λ k > 0 always holds ∀k ≥ 2. We can conclude thus that for any given For the case of directed networks, given equation (3.38), we can see that α > 0 can be chosen as large as we need if all Laplacian eigenvalues λ k satisfy

Re(λ k ) > 1 -Im(λ k ) cot(Im(λ k )). (3.39)
We do a change of variables λ = λ -1 = σ + iω. Then equation (3.39) can be written as σ > -ω cot(ω). Assume, without loss of generality, that ω > 0. This means that ω cot(ω) > 0. If σ < 0, condition ω cot(ω) > σ is immediately satisfied. On the other hand, if σ > 0, instead of condition (3.39) we can check condition ω 2 cot 2 (ω) > σ 2 . The fact that the eigenvalues are in the unit circle is translated into the inequality σ 2 +ω 2 ≤ 1.

So we can write σ 2 ≤ 1 -ω 2 . But because we assumed ω > 0, we have |ω| |sin(ω)| > 1, which leads to 1 -ω 2 < ω 2 cot 2 (ω). This means that condition (3.39) is true for all eigenvalues of the normalized Laplacian. Thus, we conclude by saying that if we set τ = τ * , for any given directed network we can always choose α > 0 satisfying condition (3.36) such that s = -1 the dominant root of the characteristic equation (3.15). Remark 3.11. For any given undirected and connected network we can find α such that equation (3.37) is satisfied. Indeed, the right-hand side of (3.37) is a continuous function going down to zero, as Figure 3.8 shows. This is to say that for any given undirected and connected network we can guarantee dominant root s = -1.

This result is remarkable because it means for instance that we can reach consensus 10 times faster in the case of a sparse network (with the algebraic connectivity 0.1), without rewiring the network (see [START_REF] Olfati-Saber | Ultrafast consensus in small-world networks[END_REF], where the proposed algorithm leads to a significant increase of the speed of convergence -even up to 100 times faster -but by using a slightly more complex procedure and by actually changing the network structure).

In the sequel, we present a couple of examples to compare the convergence speed of the proposed model with the classical consensus problem. We shall see that the performances' improvement is even more dramatic in the case of sparse networks.

Example 3.1. Consider a Erdös-Rényi random oriented graph with 100 nodes. We generate it using a probability p = 0.05. We choose α = 1 and τ = τ * for the proposed model (3.12). We compare the performances of our model with the case of classical consensus problem. The comparison between the standard deviation for the two cases is shown in Figure 3.9. We can see that using the proposed algorithm (3.12) and carefully selecting the parameters τ and α as above, the value of consensus is reached faster than in the case of classical consensus problem: using the proposed algorithm (3.12) we obtain dominant root s = -1, whereas in the case of classical consensus problem the dominant root is 0.5706. This means that for this network, algorithm (3.12) is more than 20 times faster. Now we know how to choose τ and α in order to have convergence for the case of directed networks. Next section is dedicated to the optimality question.

Optimality of the delay parameter

This section gives further evidence to show that choosing τ = τ * is even a local optimum choice in terms of network performances. Consider a directed network of the form (3.12)

and let α ∈ R + . Consider only a factor of the characteristic equation (3.15):

s -(1 -λ)(1 + α) + 1 + α(1 -λ)e -sτ = 0.
In the equation above, s can be seen as a function of τ . We would like to know what happens to the real part of s when τ = τ * changes. This is why we first differentiate with respect to τ the equation above, in order to find an expression for the derivative of s with respect to τ evaluated in τ = τ * :

∂s ∂τ -α(1 -λ)e -sτ ∂s ∂τ τ + s τ =τ * = 0.
As we know that s = -1 when τ = τ * , we can replace e -sτ = 1+α α in the last equation, and obtain ∂s ∂τ

[1 -τ * (1 -λ)(1 + α)] + (1 -λ)(1 + α) = 0, which leads to ∂s ∂τ = (1 -λ)(1 + α) (1 -λ)(1 + α)τ * -1 .
Note that in the case of undirected networks there is no difference between the real part of the right-hand part of the equation above and the right-hand side itself, as all involved variables are real. So, in order to see what happens with the sign of the real part of the left-hand side we shall study two cases.

Case 1. If λ > 1, then we immediately obtain ∂s ∂τ > 0.

Case 2. If λ < 1, then the numerator (1 -λ)(1 + α) is positive, but the denumerator is negative, as we had already chosen (1 -λ)(1 + α)τ * < 1 (see Lemma 3.4, where

β = (1 -λ)(1 + α)τ *
). Thus, we obtain ∂s ∂τ < 0.

For λ = 1 we can easily see that the sign does not change.

This is to say that a perturbation on τ = τ * will make the characteristic roots to move one way (towards left or right), corresponding to eigenvalues of the Laplacian L, λ > 1, and to move in the opposite direction for λ < 1. This situation is depicted in Figure 3.12.

For the case of directed networks, we have a similar situation, only that instead of λ we have Re(λ):

Re ∂s ∂τ τ =τ * = Re (1 -λ)(1 + α) (1 -λ)(1 + α)τ * -1 .
If we write λ = Re(λ) + i Im(λ), we multiply the fraction with the complex conjugate, and separate real and imaginary part, then we obtain

= 0 0.5 1 1.5 2 Re(s) -2 -1 0 1 = = 1 = = ln 2 s = -1 Figure 3
.12: Evolution of real part of dominant characteristic roots as a function of τ for α = 1. We note that for τ = ln 2, the dominant root is -1. Moreover, we note that the network can reach consensus up to τ = (1/α) = 1.

Re

∂s ∂τ τ =τ * = (1 -Re(λ)) 2 (1 + α) 2 τ * -(1 -Re(λ))(1 + α) + Im 2 (λ)(1 + α) 2 τ * [(1 -Re(λ))(1 + α)τ * -1] 2 + [Im(λ)(1 + α)τ * ] 2 .
It is straightforward that if Re(λ) > 1 in the equation above, then Re ∂s ∂τ τ =τ * > 0, ∀α > 0. Further, we claim that for any given network there always exists an eigenvalue λ k of the normalized Laplacian L such that Re(λ k ) > 1. If not, suppose there is an eigenvalue λ k with Re(λ k ) ≤ 1. Then we have n k=2 Re(λ k ) ≤ n-1. But this contradicts the property of the normalized Laplacian n k=2 Re(λ k ) = n, and thus proves the claim. Given that every network has an eigenvalue of L with real part greater that 1, for which Re ∂s ∂τ τ =τ * > 0, we have now the sufficient condition to prove that Re(s) = -1 is a local optimum for τ = τ * . This argument shows that, in networks where Re(λ k ) > 1 for all k ≥ 2, the derivative Re(∂s/∂)τ > 0; so, decreasing τ from τ * will improve convergence speed. An example of such a network is given in Example 3.3, and the accompanying Figure 3.13 demonstrates the situation. This difficulty does not occur in undirected networks, since Re(∂s/∂τ ) takes both negative and positive values over the set of Laplacian eigenvalues (except for complete graphs, in case of binary a ij ).

As stated above, τ * is only locally optimum. There are cases where the classical consensus protocol (3.2) performs better than using τ * in the proposed algorithm (3.7), as depicted in Example 3.3.

Example 3.3. We consider a directed network made up of 9 nodes. The adjacency matrix is

A =                     0 1 1 0 0 0 0 1 0 0 0 1 1 1 1 0 1 0 1 0 0 0 1 0 0 1 1 0 1 1 0 0 0 0 0 0 1 0 1 0 0 1 0 1 1 0 1 1 1 1 0 0 0 0 1 1 1 1 0 1 0 1 0 0 0 0 1 0 1 0 0 1 0 1 0 1 1 0 1 0 0                     .
The eigenvalues of the corresponding Laplacian L are: 0, 1.0740 ± 0.3691i, 1.3590, 1.1899 ± 0.1679i, 1.0576, and 1.0278 ± 0.2463i.

It is easy to see here that we obtain s = -1 dominant root, but the consensus value is reached slower than in the case of classical consensus problem. Figure 3.13 depicts the value of the real part of the dominant root corresponding to each (pair of ) eigenvalue(s), as a function of τ , for α = 1. The figure shows that for τ = ln 2 we obtain Re(s) = 1 for all eigenvalues, but this value is a bit "slower" compared to the case τ = 0. Example 3.3 shows that the local optimum given by τ * may not always perform better than the undelayed algorithm. Nevertheless, even in such cases we can still find (under some conditions on the network) a delay τ for protocol (3.7) in order to obtain better performances that in the case of classical protocol (3.2).

Concluding remarks

We have introduced, in Sections 3.2 and 3. Anticipation can improve performance in networks of interacting agents if the parameter τ is chosen judiciously, otherwise the system can diverge as agents try to anticipate too aggressively into the future. We have given simulation results from various networks and shown that a positive τ can indeed improve speed of consensus, in some cases dramatically, especially in networks with poor connectivity.

For both undirected and directed networks, we show that by well choosing τ and α the proposed model guarantees a very fast convergence speed for any network, and moreover, this choice is a local optimum in terms of network performance.

Chapter 4

A theoretical bound for maximal multiplicity at the origin

Chapter overview

This chapter focuses on a rather general class of time-delay systems, with a finite number of delays. We study the case where such systems have multiple characteristic roots at the origin, and show the link between the system parameters and the maximal multiplicity at the origin. Part of the work discussed in this chapter has been published in a book chapter [START_REF] Irofti | On the Codimension of the Singularity at the Origin for Networked Delay Systems[END_REF] and published and presented at an international conference [START_REF] Boussaada | Computing the codimension of the singularity at the origin for delay systems in the regular case: A vandermondebased approach[END_REF].

The chapter is structured as follows. The considered class of systems is described in Section 4.2. An exact formula for the theoretical bound of the maximal algebraic multiplicity of the zero characteristic root is given in Section 4.3. Moreover, algebraic conditions on a system's parameters are given for this maximal bound to be reached. These results are briefly compared to other results in the literature in Section 4.4, where our contribution is highlighted: we use two simple examples to illustrate that our bound is more precise.

Two mechanical examples are discussed in Section 4.5 in order to illustrate the results.

Problem formulation

In the case of an n-dimensional linear system of ordinary differential equations ẋ = Ax, where x ∈ R n and A ∈ M n (R), the n eigenvalues of the matrix A are at the same time the spectral values of the system. Since the corresponding characteristic equation is a polynomial of degree n in the Laplace variable s, it has at most n complex roots -the eigenvalues of the system. So, we can say that the (algebraic) multiplicity of a given spectral value can at most be equal to the dimension of the state-space.

On the other hand, the case of time-delay systems is slightly different as, in this case, the characteristic equation has an infinite number of roots. Consider a system described by the following delay-differential equation:

ẋ = N i=0 A i x(t -τ i ), (4.1)
under appropriate initial conditions belonging to the Banach space of continuous func-

tions C([-τ N , 0], R n ).
Here, x = (x 1 , . . . , x n ) denotes the state-vector, matrices A i ∈ M n (R) for i = 0 . . . N (see [START_REF] Ahlfors | An introduction to the theory of analytic functions of one complex variable[END_REF]106]). The N constant delays

τ i , i = 1 . . . N satisfy 0 = τ 0 < τ 1 < τ 2 < • • • < τ N .
In the corresponding characteristic equation (a transcendental equation in the Laplace variable) some exponential terms appear due to the delays. In this chapter, we consider matrices A i , i = 1 . . . N , such that the system (4.1)

has the characteristic function q : C × R N + → C defined as:

q(s, τ ) = r 0 (s) + N i=1 r i (s) e -τ i s , (4.2)
where τ is the vector containing all delays, and r i (s) are polynomials of s with real coefficients.

Without any loss of generality, assume that the polynomial r 0 is a monic polynomial of degree n in s and the polynomials r i are such that deg(

r i ) ≤ (n -1), ∀1 ≤ i ≤ N.
The study of zeros of the quasi-polynomial (4.2) plays a crucial role in the analysis of asymptotic stability of zero solution of system (4.1). Indeed, the zero solution is asymptotically stable if all the zeros of (4.2) are in the open left half-complex plane [START_REF] Michiels | Stability and Stabilization of Time-Delay Systems (Advances in Design and Control)[END_REF].

According to this definition, the parameter space which is spanned by the coefficients of the polynomials r i , can be split into stability and instability domains (the so-called D-decomposition, see [START_REF] Michiels | Stability and Stabilization of Time-Delay Systems (Advances in Design and Control)[END_REF][START_REF] Gryazina | D-decomposition technique stateof-the-art[END_REF] and references therein). These two domains are separated by a boundary corresponding to a spectrum consisting of roots with zero real parts.

Moreover, under appropriate algebraic restrictions, a given root corresponding to the boundary can have high multiplicity.

A typical example of a non-simple (multiple) zero spectral value is the Bogdanov-Takens singularity which has an algebraic multiplicity two and a geometric multiplicity one.

Cases with higher order multiplicities of the zero spectral value are known as generalized Bogdanov-Takens singularities and could be involved in concrete applications [START_REF] Campbell | Zero singularities of codimension two and three in delay differential equations[END_REF]. For instance, the Bogdanov-Takens singularity is identified in [START_REF] Hale | Period doubling in singularly perturbed delay equations[END_REF] where the case of two coupled scalar delay equations modeling a physiological control problem is studied. In [START_REF] Saldivar Marquez | Analysis and Control of Oilwell Drilling Vibrations[END_REF], this type of singularity is also encountered in the study of coupled axial-torsional vibrations of an oilwell rotary drilling system.

Commonly, the time-delay induces desynchronizing and/or destabilizing effect on the dynamics. However, new theoretical developments in control of finite-dimensional dynamical systems suggest the use of delays in the control laws as controller parameters for stabilization purposes. For instance, the papers [START_REF] Sieber | Bifurcation analysis of an inverted pendulum with delayed feedback control near a triple-zero eigenvalue singularity[END_REF][START_REF] Boussaada | Inverted pendulum stabilization: Characterization of codimension-three triple zero bifurcation via multiple delayed proportional gains[END_REF] are concerned with the stabilization of the inverted pendulum by delayed control laws. In [START_REF] Sieber | Bifurcation analysis of an inverted pendulum with delayed feedback control near a triple-zero eigenvalue singularity[END_REF], the authors prove that some appropriate delayed proportional-derivative (PD) controller stabilizes the inverted pendulum by identifying a singularity of multiplicty three for a system of two coupled delayed equations. In [START_REF] Boussaada | Inverted pendulum stabilization: Characterization of codimension-three triple zero bifurcation via multiple delayed proportional gains[END_REF], the same singularity is studied by using a particular delay block configuration.

In this chapter we investigate this type of singularity (multiple zero singularity) and give an answer to the question concerning its maximal multiplicity (as a function of the coefficients of the polynomials r i (s), i = 1 . . . N ). This work is motivated by the importance of the multiplicity of singularities in nonlinear stability analysis [START_REF] Carr | Application of Center Manifold Theory[END_REF][START_REF] Hale | Introduction to functional differential equations[END_REF][START_REF] Guckenheimer | Nonlinear oscillations, dynamical systems, and bifurcation of vector fields[END_REF][START_REF] Kuznetsov | Elements of applied bifurcation theory[END_REF].

Maximal algebraic multiplicity of zero singularity

We consider N non-zero distinct delays 0 < τ 1 < τ 2 < . . . < τ N in equation (4.2), and the polynomials r i such that r 0 is a monic polynomial with deg(r 0 ) = n and deg(r i ) ≤ (n-1)

for 1 ≤ i ≤ N . Let n 0 = max 1≤i≤N deg(r i ).
We denote by a i,k the coefficient of s k for the polynomial r i , with 0 ≤ i ≤ N . Note that a 0,n = 1. Using these notations, we write the following proposition. Moreover, such a bound is reached if and only if for 0 ≤ k ≤ 0 m -1, the parameters of (4.2) simultaneously satisfy : We need to introduce some notations. Denote by q (k) (s) the k-th derivative of q(s)

a 0,k = - N i=1 a i,k + k-1 l=0 (-1) l+k+1 a i,l τ i k-l (k -l) ! . ( 4
with respect to the variable s. Recall that zero is an eigenvalue of algebraic multiplicity m ≥ 1 for (4.1) if q(0) = q (k) (0) = 0 for all k = 1, . . . , (m -1) and q (m) (0) = 0.

We first prove condition (4.3), which follows directly from the Lemma 4.2.

Lemma 4.2. Zero is a root of q (k) (s) for k ≥ 0 if and only if all coefficients of r i , with

0 ≤ i ≤ N , satisfy condition (4.3).
Proof. We define the family ∇ k for all k ≥ 0 by

∇ k (s) = N i=1 d k ds k r i (s) + k-1 j=0 (-1) j+k k j N i=1 τ k-j i d j ds j r i (s) , (4.4)
where d 0 ds 0 f (s) := f (s) and k j denote the binomial coefficient equal to k! j!(k-j)! . We prove by induction that zero is a root of q (k) (s) for k ≥ 0 if and only if zero is a root of ∇ k (s). Note that the defined family (4.4) is polynomial since r i and their derivatives are polynomials. As for the induction, we prove that differentiating k times q(s) we obtain a recursive formula:

q (k) (s) = N i=1 d k ds k r i (s)e τ i s + k-1 j=0 (-1) j+k k j N i=1 τ k-j i d j ds j r i (s)e τ i s ,
where we can set e -τ i s = 1, since we study the multiplicity at s = 0. Thus, we obtain the relation (4.4) that finally leads to the condition (4.3).

Example 4.1. To illustrate Lemma 4.2, consider the scalar delay differential equation:

ẋ(t) + a 0,0 x(t) + a 1,0 x(t -τ 1 ) + a 2,0 x(t -τ 2 ) = 0,
where a 0,0 , a 1,0 , a 2,0 ∈ R. The corresponding characteristic quasi-polynomial function is

given by: q(s, τ ) = s + a 0,0 + a 1,0 e -τ 1 s + a 2,0 e -τ 2 s .

For k = 0, equality (4.3) gives the first sufficient condition guaranteeing a multiplicity at least one for s = 0. Indeed, q(0, τ ) = 0 = a 0,0 +a 1,0 +a 2,0 which gives a 0,0 = -(a 1,0 +a 2,0 ).

A sufficient condition for an algebraic multiplicity at least two for s = 0 is obtained by computing the first partial derivative of q with respect to s. This leads to condition a 0,1 = 1 = τ 1 a 1,0 + τ 2 a 2,0 , which satisfies equality (4.3) with k = 1.

In a similar way, one can easily show that if equality (4.3) is satisfied for all k = 0, 1, 2, then the zero singularity is of multiplicity three (which is the maximal multiplicity) if and only if

a 0,0 = - τ 1 + τ 2 τ 1 τ 2 , a 1,0 = - τ 2 (τ 1 -τ 2 ) τ 1 , a 2,0 = τ 1 τ 2 (τ 1 -τ 2 )
.

In the sequel, we consider the variety corresponding to polynomials ∇ k vanishing, i.e.

∇ 0 (0) = . . . = ∇ m-1 (0) = 0 and ∇ m (0) = 0. Recall that we aim to find the maximal algebraic multiplicity m of zero singularity. The first elements from the family ∇ k for s = 0 are

∇ 0 (0) = 0 ⇔ N i=0 a i,0 = 0, ∇ 1 (0) = 0 ⇔ N i=0 a i,1 - N i=1 a i,0 τ i = 0, ∇ 2 (0) = 0 ⇔ 2! N i=0 a i,1 -2! N i=1 a i,0 τ i + N i=1 a i,0 τ 2 i = 0.
If we consider a i,k and τ i as variables, the obtained algebraic system is nonlinear and solving it without giving values for n and N becomes a very difficult task. Even by using Gröbner basis methods [START_REF] Cox | Ideals, varieties, and algorithms. An introduction to computational algebraic geometry and commutative algebra[END_REF], this task is still complicated because the set of variables depends on N , n and n 0 . This is why we additionally assume that all the polynomials r i satisfy the condition deg(r i ) = n 0 for all 1 ≤ i ≤ N . We consider a i,k as variables and τ i as parameters and we use the following notation: a 0 = (a 0,0 , a 0,1 , . . . , a 0,n-1 ) T and

a i = (a i,0 , a i,1 , . . . , a i,n 0 ) for 1 ≤ i ≤ N , τ = (τ 1 , τ 2 , .
. . , τ N ), and a = (a 1 , a 2 , . . . , a N ) T .

Consider the ideal I 1 generated by the n polynomials < ∇ 0 (0), ∇ 1 (0), . . . , ∇ n-1 (0) >.

As it can be seen from Lemma 4.2, the variety V 1 associated with the ideal I 1 has the linear representation a 0 = M 1 (τ ) a, where the matrix

M 1 ∈ M n,N (n 0 +1) (R[τ ]).
In other words, we write a matrix equation for equality (4.3), where a 0 contains all a 0,k coefficients with 0 ≤ k < n corresponding to r 0 polynomial, and a contains all coefficients corresponding to r i polynomials, with 1 ≤ i ≤ N .

Note that in this variety there is no restriction on the components of a when a 0 is left free. Since a 0,k = 0 for all k > n, the remaining equations consist in an algebraic system only in a and parametrized by τ .

Consider now the ideal I 2 generated by the N (n 0 + 1) polynomials:

I 2 =< ∇ n+1 (0), ∇ n+2 (0), . . . , ∇ n+N (n 0 +1) (0) > .
We remark that the variety V 2 corresponding to ideal I 2 can be written as M 2 (τ )a = 0, which is nothing else that a homogeneous linear system with the matrix

M 2 ∈ M N (n 0 +1) (R[τ ]
). More precisely, M 2 is a functional Vandermonde matrix:

M 2 (τ ) =(V (τ 1 ), d dτ 1 V (τ 1 ), . . . , d n 0 dτ n 0 1 V (τ 1 ), . . . , V (τ N ), . . . , d n 0 dτ n 0 N V (τ N ))
, where

V (x) = ((-x) n+1 , (-x) n+2 , . . . , (-x) n+N (n 0 +1) ) T . (4.5)
We note that every subset of vectors

F k = (V (τ k ), . . . , d n 0 dτ n 0 k V (τ k )) is a family of vectors in R N (n-1) ([τ k ]). Denote by V k the k-th component
(n 0 + 1 -k)! N ˙ 1≤i<l≤N (τ i -τ l ) (n 0 +1) 2 1≤h≤N τ (n 0 +1)(n+1) h .
Recall that the distinct delays are strictly positive. Thus, the determinant of matrix M 2 cannot vanish. This means that the only solution for the subsystem M 2 (τ )a = 0 is the zero solution, that is a = 0.

Next, we can write

∇ n (0) = 0 ⇔ 1 = - N i=1 n-1 l=0 (-1) l+n+1 a i,l τ i n-l (n -l) ! .
Substituting the unique solution of V 2 into the last equality leads to a contradictory result. This shows that the maximal algebraic multiplicity of the zero singularity is less or equal to N (n 0 + 1) + n.

Discussion on results

We begin this section with two examples.

Example 4.2 (A vector disease model). We consider the scalar delay-differential equation

ẋ(t) = -bx(t) + cx(t -τ ) -cx(t)x(t -τ ), (4.6)
presented by Cooke in [START_REF] Cooke | Stability analysis for a vector disease model[END_REF] to describe the spread of a communicable disease carried by a vector. This very simple deterministic biological model describes the evolution of infectious persons. More precisely, x designates the infected host population governed by the law (4.6), where b is the recovery rate (assuming that the host recovery is an exponential process) and c is the contact rate between infected and uninfected populations.

We linearize equation (4.6) near the point x = 0, to obtain

ẋ(t) = -bx(t) + cx(t -τ ), (4.7)
where b, c and τ are positive constants. The corresponding characteristic (transcendental) function is

q(s) = s -b + ce -sτ . (4.8)
Zero is a spectral value for (4.7) if and only if (4.8) vanishes at zero, i.e. condition -b + c = 0 holds. We compute the first derivative of (4.7) with respect to s, and obtain

q (s) = 1 -τ ce -sτ .
Thus, it is easy to see that s = 0 has an algebraic multiplicity two (i.e. q(0) = 0, q (0) = 0, and q (0) = 0) under conditions b = c and τ = 1/c. Note that for satisfying the vanishing of the second derivative of (4.7),

q (s) = τ 2 ce -sτ ,
the domain in which the positive parameters b, c and τ can take values is empty. We can thus conclude that the maximal multiplicity of the zero spectra value, s = 0, is at most two. This upper bound is less than the number of free parameters (b, c, and τ ) of the system. Using the QPmR software package (see appendix B), we choose b = 0.5 and represent the spectrum of system (4.7) in Figure 4.1.

Example 4.3 (Inverted pendulum on a cart).

There are a few ways of modeling an inverted pendulum on a cart (see Figure 4.2). This system is an example commonly found in the control theory literature due to its various applications. The motion equations usually refer to the motion of both the pendulum itself and the cart, under some assumptions (the absence or not of the friction or other resistance to movement). Here we consider the simple case of an inverted pendulum on a cart governed by a second-order differential equation (as in [START_REF] Boussaada | Computing the codimension of the singularity at the origin for delay systems: The missing link with birkhoff incidence matrices[END_REF]):

1 - 3ε 4 cos 2 θ θ + 3ε 8 θ2 sin(2θ) -sin θ + F cos θ = 0,
where ε = m m+M is the relative mass, M is the mass of the cart, m is the mass of the pendulum, and F is the driving force, representing the control law. If we choose F of the form

F = k 1 θ(t -τ 1 ) + k 2 θ(t -τ 2 )
and set ε = 3 4 , then the characteristic transcendental function is

q(s) = s 2 - 16 7 + 16k 1 7 e -sτ 1 + 16k 2 7 e -sτ 2 . (4.9)
Some simple computations show us that for

k 1 = - 7 8 τ 1 2 -7 k 2 = 8 τ 1 2 8 τ 1 2 -7 τ 2 = 7 8 τ 1 (4.10) (4.11) (4.12)
we have q(s = 0) = 0, q (s = 0) = 0, q (s = 0) = 0 and q (3) (s = 0) = 0, i.e. for this system the zero root has multiplicity 3. The spectrum corresponding to τ 1 = 1 is depicted in

d 1 ≥ 1, d 2 ≥ 1 . . . d N ≥ 1, d 1 + d 2 + . . . + d N = D + N.
Let f i,j (s) stands for the function f i,j (s) = s j-1 e τ i s , for 1 ≤ j ≤ d i and 1 ≤ i ≤ N .

Let be the number of zeros of the function

f (s) = 1≤i≤N,1≤j≤d i c i,j f i,j (s),
that are contained in the horizontal strip α ≤ Im(s) ≤ β.

Assuming that

1≤k≤d 1 |c 1,k |> 0, . . . , 1≤k≤d N |c N,k |> 0, then (τ N -τ 1 ) (β -α) 2 π -D -N + 1 ≤ ≤ (τ N -τ 1 ) (β -α) 2 π + D + N -1.
Idea of the proof. We apply the argument principle (see Appendix A) on a rectangle {s : α ≤ Im(s) ≤ β and -a ≤ Re(s) ≤ a} and then let a tend to ∞.

We apply Proposition 4.3 on the inverted pendulum example, and obtain that the maximal multiplicity of the zero characteristic root s = 0 is 6. One of the main contributions of the result given by Proposition 4.1 is that it provides a sharper upper bound. Indeed, as proven in [START_REF] Boussaada | Inverted pendulum stabilization: Characterization of codimension-three triple zero bifurcation via multiple delayed proportional gains[END_REF], the maximal multiplicity of the zero root s = 0 in the case of the inverted pendulum is three. Obviously, this bound is more precise than the result given by Pólya-Szegö: 3 < 6. This is because the approach presented in this chapter is a constructive one: it is obvious that the maximal multiplicity of zero root decreases below the limit given by the Pólya-Szegö result if a component of the vector parameters is not left free (because we'll have more restrictions on the parameters).

Moreover, according to Proposition 4.1, this bound is reached when the parameters of characteristic equation (4.9) satisfy condition (4.3). We note that formula (4.3) can be actually translated into equations (4.10)-(4.12).

We can thus conclude that the upper bound given by Proposition 4.1 is the same as the maximal number of roots in the horizontal strip found by Pólya-Szegö, but only when all the parameters of the quasi-polynomial are left free. In all other cases this bound is sharper than Pólya-Szegö's one.

Illustrative examples

In the sequel, we have chosen two mechanical systems, benchmarks in the field of control theory and robotics, to illustrate previous results: an inverted pendulum on a slope and a double inverted pendulum on a cart.

Inverted pendulum on a slope

We consider an inverted pendulum on a cart, moving on a slope with angle α, as in Under the assumptions above, we have the following model for the system depicted in Figure 4.4 (see [START_REF] Kokame | State difference feedback for stabilizing uncertain steady states of non-linear systems[END_REF]):

   θ -g L sin(θ -α) + 1 L d cos α = 0 M d + b ḋ = F -M g sin α , (4.13)
where L = J+mL 2 mL .

We choose a control law with 2 delays

F = k 1 θ(t -τ 1 ) + k 2 θ(t -τ 2 ), (4.14)
where τ 1 and τ 2 are two positive distinct delays and the gains k 1 and k 2 real numbers.

Using the set of equations (4.13), we obtain a state-space representation of our system: we choose the state vector x = d L θ ḋ L θ T and we linearize around the equilibrium point (0, L α, 0, 0).

Thus, we can write the state-space representation

ẋ(t) = A 0 x(t) + A 1 x(t -τ 1 ) + A 2 x(t -τ 2 ), (4.15)
where

A 0 =        0 0 1 0 0 0 0 1 0 0 -b M 0 0 g L b cos α M 0        ; A 1 =        0 0 0 0 0 0 0 0 0 k 1 L M 0 0 0 -k 1 cos α M 0 0        ; A 2 =        0 0 0 0 0 0 0 0 0 k 2 L M 0 0 0 -k 2 cos α M 0 0       
.

We assume the measured physical parameters are: m = 0.05kg, M = 0.5kg, L = 2.5m,

J = 1 3 mL 2 [kg m 2 ], α = π 6 , b = 1Ns/m, g = 9.8m/s 2 .
The characteristic equation of (4.15) is det q(s, τ ) = 0, where the characteristic transcendental function q(s, τ ) = sI 6 -A 0 -A 1 e -sτ 1 -A 2 e -sτ 2 can be written as:

q(s, τ ) = P 0 (s) + 2 i=1 P i (s)e -sτ i
, where

P 0 (s) = s 4 + 2 s 3 - 147 50 s 2 - 147 25 s, P 1 (s) = √ 3 k 1 s 2 + 7 √ 3 k 1 5 s, P 2 (s) = √ 3 k 2 s 2 + 7 √ 3 k 2 5 s.
Proposition 4.4. Under the effect of the delayed feedback given by (4.14) the multiplicity of the zero spectral value s = 0 for system (4.15) is at most 4.

Proof. It is easy to see that q(0) = 0, since s is a common factor for all P i , with i ∈ {0, 1, 2}. The first derivative of q(s) is

q (s) = 6 s 2 - 147 s 25 + 4 s 3 + 7 √ 3 k 1 e -s τ 1 5 + 7 √ 3 k 2 e -s τ 2 5 + 2 √ 3 k 1 s e -s τ 1 + 2 √ 3 k 2 s e -s τ 2 - √ 3 k 1 s 2 τ 1 e -s τ 1 - √ 3 k 2 s 2 τ 2 e -s τ 2 - 7 √ 3 k 1 s τ 1 e -s τ 1 5 - 7 √ 3 k 2 s τ 2 e -s τ 2 5 - 147 25 .
If k 1 is of the form

k 1 = 7 √ 3 5 -k 2 (4.16)
with the parameter k 2 left free, then q (0) = 0 and the maximal multiplicity of the zero spectral value is at least 2.

Next, we compute the second derivative of q(s) and find that under conditions

   k 1 = 7 √ 3 5 - 42 τ 1 -9 10 √ 3 τ 1 -10 √ 3 τ 2 k 2 = 42 τ 1 -9 10 √ 3 τ 1 -10 √ 3 τ 2 , (4.17)
and assuming that the delays are positive and distinct, the multiplicity at the origin is

at least 3. Furthermore, if                  k 1 = 1729 √ 3 20 (-147 τ 1 2 +63 τ 1 +55) k 2 = -7 √ 3 (14 τ 1 -3) (42 τ 1 -9) 20 (-147 τ 1 2 +63 τ 1 +55) τ 2 = 63 τ 1 +110 294 τ 1 -63 τ 1 > 0.215 (4.18)
holds, then the third derivative of q(s) vanishes at zero, i.e. the multiplicity at origin is 4. This is also the maximal multiplicity of the zero spectral values, since there is no set of real values (k 1 , k 2 , τ 1 , τ 2 ) with τ 1 and τ 2 positive and distinct, that satisfies q (4) (0) = 0. Remark 4.5. Note that in conditions set (4.18) τ 1 is a parameter left free, such that τ 1 > 0.215. If the choose for instance τ 1 = 1, we obtain k 1 = -5.16, k 2 = 7.59 and τ 2 = 0.75. We replace those values in the polynomials P 0 , P 1 and P 2 ; the coefficients obtained verify Remark 4.6. Using a constructive approach we found an upper bound of the multiplicity at the origin. Since a component of vector of parameters a (containing the coefficients of P 0 ) is not left free, we have more constraints on the gains and delays, and consequently the maximal multiplicity of the zero spectral values decreases under the bound given by Pólya-Szegö result (4 < 6).

a 0,k = - N i=1 a i,k + k-1 l=0 (-1) l+k+1 a i,l τ i k-l (k -l) ! . for 1 ≤ k ≤ 4 and 1 ≤ i ≤ 2. Choosing τ 1 = 1, k 1 = -5.
Remark 4.7. The upper bound of the multiplicity at the origin is linked not only to the vanishing coefficients of a, but also to system's physical parameters. For instance, for the same choice of parameters except for the length L of the inverted pendulum (the distance from the pivot to the center of gravity of the pendulum) we can obtain a sharper bound. More precisely, if we consider that L decreases under 2.5m, then the maximal multiplicity of the zero spectral value is 3.

Double inverted pendulum

We consider a double inverted pendulum described in Figure 4.6, where θ 1 and θ 2 are the angular position of the lower and upper pendulum with respect to the vertical axis, and r is the position of the cart with respect to a reference point.

The motion of the system can be described by the Lagrange equation

d dt ∂L ∂ θ - ∂L ∂θ = Q (4.19)
where L = T -V is the Lagrangian, T is the kinetic energy, V is the potential energy and Q is a vector of generalized forces. The kinetic and potential energies, T and V, of the system can be written as the sum of energies of each individual component, the cart and the two inverted pendulums. The kinetic and potential energy are written in the form

T = T 0 + T 1 + T 2 V = V 0 + V 1 + V 2 ,
where

T 0 = 1 2 m 0 ṙ2 T 1 = 1 2 m 1 ṙ + l 1 θ1 cos θ 1 2 + l 1 θ1 sin θ 1 2 + 1 2 I 1 θ2 1 = 1 2 m 1 ṙ2 + 1 2 m 1 I 2 1 + I 1 θ2 1 + m 1 l 1 ṙ θ1 cos θ 1 T 2 = 1 2 m 2 ṙ + L 1 θ1 cos θ 1 + l 2 θ2 cos θ 2 2 + L 1 θ1 sin θ 1 + l 2 θ2 sin θ 2 2 + 1 2 I 2 θ2 2 = 1 2 m 2 ṙ2 + 1 2 m 2 L 2 1 θ2 1 + 1 2 m 2 l 2 2 + I 2 θ2 2 + m 2 L 1 ṙ θ1 cos θ 1 + m 2 l 2 ṙ θ2 cos θ 2 + m 2 L 1 l 2 θ1 θ2 cos (θ 1 -θ 2 ) V 0 = 0 V 1 = m 1 gl 1 cos θ 1 V 2 = m 2 g (L 1 cos θ 1 + l 2 cos θ 2 ) .
Then, the Lagrangian of the system can be written

L = 1 2 (m 0 + m 1 + m 2 ) ṙ2 + 1 2 m 1 l 2 1 + m 2 L 2 1 + I 1 θ2 1 + 1 2 m 2 l 2 2 + I 2 θ2 2 + (m 1 l 1 + m 2 L 1 ) cos (θ 1 ) ṙ θ1 + m 2 l 2 cos (θ 2 ) ṙ θ2 + m 2 L 1 l 2 cos (θ 1 -θ 2 ) θ1 θ2 -(m 1 l 1 + m 2 L 1 ) g cos θ 1 -m 2 l 2 g cos θ 2
We differentiate the Lagrangian with respect to θ and θ and rewrite the equation

m i r + (m 1 l 1 + m 2 L 1 ) cos (θ 1 ) θ1 + m 2 l 2 cos (θ 2 ) θ2 -(m 1 l 1 + m 2 L 1 ) sin (θ 1 ) θ2 1 -m 2 l 2 sin (θ 2 ) θ2 2 0 = (m 1 l 1 + m 2 L 1 ) cos (θ 1 ) r + m 1 l 2 1 + m 2 L 2 1 + I 1 θ1 + m 2 L 1 l 2 cos (θ 1 -θ 2 ) θ2 + m 2 L 1 l 2 sin (θ 1 -θ 2 ) θ2 -(m 1 l 1 + m 2 L 1 ) g sin θ 1 0 = m 2 l 2 cos (θ 2 ) r + m 2 L 1 l 2 cos (θ 1 -θ 2 ) θ1 + m 2 l 2 2 + I 2 θ2 -m 2 L 1 l 2 sin (θ 1 -θ 2 ) θ1 2 -m 2 l 2 g sin θ 2 .
Thus, the system model described by three nonlinear differential equations of second order can be translated into a state-space representation, consisting of six first order ODEs (see for more details [START_REF] Furuta | Computer control of a double inverted pendulum[END_REF] and [START_REF] Bogdanov | Optimal control of a double inverted pendulum on a cart[END_REF]).

We chose θ = r θ 1 θ 2

T

, where r is the position of the cart with respect to the origin; m 1 = 0.5kg and m 2 = 0.5kg are the masses of the lower, respectively upper pendulum and L 1 = 0.5m, L 2 = 0.5m are the corresponding lengths. We consider that only the position of the cart is affected by the command and we linearize around θ 1 = θ 2 = 0 and r = 0. Finally we obtain a state-space representation

ẋ(t) = A 0 x(t) + A 1 x(t -τ 1 ) + A 2 x(t -τ 2 ), (4.20)
where

x = θ θ T ∈ R 6 , A 0 =        0 3 I 3 0 -441 46
u = γ 11 θ 1 (t -τ 1 ) + γ 12 θ 2 (t -τ 1 ) + γ 21 θ 1 (t -τ 2 ) + γ 22 θ 2 (t -τ 2 ), (4.21)
and τ 1 and τ 2 , satisfying 0 < τ 1 < τ 2 are the delays. We denote by 0 3 and I 3 the zero matrix and the identity matrix of dimension three, respectively.

The characteristic equation of (4.20) is det q(s, τ ) = 0, where the characteristic transcendental function q(s, τ ) = sI 6 -A 0 -A 1 e -sτ 1 -A 2 e -sτ 2 can be written as:

q(s, τ ) = P 0 (s) + Proof. It is easy to see that q(0) = q (0) = 0, since s 2 is a common factor for all P i , with then the multiplicity at the origin is at least 3, since q (2) (0) = 0. Next, we compute the third derivative of q(0), under the condition (4.22), in order to check if there is a real set of (γ 11 , γ 12 , γ 21 , γ 22 , τ 1 , τ 2 ) with 0 < τ 1 < τ 2 that satisfies q (3) (0) = 0. Thus, we find the condition

i ∈ {0,
γ 21 = - 5 6 2 γ 22 τ 1 -2 τ 2 γ 22 -49 τ 1 τ 1 -τ 2 . (4.23)
Now we know that if the system's parameters satisfy conditions (4.22) and (4.23), γ 12 , γ 22 , τ 1 and τ 2 are left free and 0 < τ 1 < τ 2 , then the multiplicity at the origin is at least 4 (since q(0) and its first three derivatives at origin are zero). Following the same reasoning, we are looking for constraints on γ 12 , γ 22 , τ 1 and τ 2 such that the next derivative of q at the origin vanishes, and we obtain hold, the multiplicity at the origin is at least 6,

γ 12 = 2401 4 τ 1 τ 2 -γ 22 - 539
γ 22 = 49 60
τ 1 490 τ 1 τ 2 -245 τ 2 2 -44 τ 1 -τ 2 , (4.25)
and in addition, for the sixth derivative of q at the origin to vanish, the condition (4.26) must be also satisfied, Next, we can easily check that q (7) = 0, i.e. there is no set of real values for (γ 11 , γ 12 , γ 21 , γ 22 , τ 1 , τ 2 ) with 0 < τ 1 < τ 2 , that satisfies q (7) (0) = 0, so s = 0 multiplicity at the origin is 7.

- 2593080 τ 3 1 τ 2 23 + 7779240 τ 1 2 τ 2
a 0,k = - N i=1 a i,k + k-1 l=0 (-1) l+k+1 a i,l τ i k-l (k -l) ! .

Concluding remarks

The work in this chapter focuses on the study of spectral values at the origin, in the general case of time-delay systems. We have found a bound for maximal multiplicity of zero spectral value.

In the case of a regular quasi-polynomial (i.e. all coefficients are non-zero) we have actually recovered the result of Pólya-Szegö. The difference lies in the essential elements of the proof. Our approach is a rather constructive one: we have made use of linear algebra, Vandermonde matrices, and mathematical induction, whereas Pólya-Szegö's proof is based on some elements of complex analysis like the argument principle (see Appendix A). Moreover, the real bound is sharper than the one given by the Pólya and Szegö result. For instance, for the case of a vector disease model, the quasi-polynomial is a regular one, so the bound given by the Pólya and Szegö result is equal to the one we obtain (which is 2). On the other hand, in the inverted pendulum example, where the corresponding quasi-polynomial is sparse (i.e. there are some null coefficients), the Pólya and Szegö result claims that the maximal multiplicity at the origin is 4, but in fact it cannot be greater than 3. The results has been illustrated on other two systems, the inverted pendulum on a slope and double inverted pendulum, and the relation between the quasi-polynomial coefficients has been verified.

Chapter 5

Insights into geometric approach

Chapter overview

This chapter discusses the behaviour of multiple imaginary characteristic roots under a small deviation of parameters. More precisely, we study the cases of multiplicity two, three and four, by considering a class of control systems depending on two parameters.

In this chapter, we discuss only the case where the two parameters are the delays; a more general case, with a larger parameter domain, is presented in Chapter 7.

One of the main contributions of the geometric approach discussed in this chapter is that it does not use Puiseux series, as compared to other results in the literature. Moreover, we show how an algebraic criterion can be derived in order to judge what happens to a pair of complex conjugated double roots, when the parameters move into a small sector. The work in this chapter has been published and presented at three international conferences (see [START_REF] Gu | Migration of double imaginary characteristic roots under small deviation of two delay parameters[END_REF][START_REF] Irofti | Migration of imaginary roots of multiplicity three and four under small deviation of two delays in time-delay systems[END_REF][START_REF] Irofti | Geometric vs. algebraic approach: A study of double imaginary characteristic roots in time-delay systems[END_REF]) and a paper has been accepted for publication in Automatica journal.

The chapter is structured as follows. The considered class of systems is discussed in more details in Section 5.2. The double root case is presented in Section 5.3. A concrete example of a time-delay system is given, illustrating how to apply the theory. The same example is also used in Subsection 5.3.5, to illustrate the comparison between the geometric approach and an algebraic approach to the same problem. This comparison sheds some light on the differences between the two methods, their limitations and the open questions to which they give rise. An extension of the proposed method is made in Section 5.4 to study characteristic roots with higher multiplicity: three and four. The concluding remarks can be found in Section 5.5.

Problem statement

Consider a system with the characteristic equation of the form q(s, p) = 0, (

where s is the Laplace variable and p ∈ R n is the vector parameter of dimension n. In the first part of this chapter, we discuss the case n = 2. Note that the characteristic equation (5.1) can be alternatively written as q(s, p 1 , p 2 ) = 0. For p 1 = p 10 and p 2 = p 20 , we assume q(s, p 1 , p 2 ) has an imaginary root s 0 = iω 0 of m th order. In other words,

∂ k q ∂s k s=s 0 p 1 =p 10 p 2 =p 20 = 0, for k = 0 . . . m -1. (5.2)
We further assume that s 0 is not a root of order (m + 1), i.e.

∂ m q ∂s m s=s 0 p 1 =p 10 p 2 =p 20 = 0. (

We make the following additional non-degeneracy assumption: where Re(•) and Im(•) denote the real and imaginary part of a complex number, respectively. Note that D in (5.4) may also be written as

D = det
D = Im ∂q * ∂p 1 • ∂q ∂p 2 s=s 0 p 1 =p 10 p 2 =p 20
, where (•) * denotes the complex conjugate of a complex number.

Equations (5.2)-(5.4) will be the standing assumptions in the remaining part of this chapter. Assumption (5.4) contains the first-order partial derivatives of q with respect to the two parameters p 1 and p 2 . If it had had zero-order instead of first-order partial derivatives, then the condition D = 0 would have been a degenerate case. We can refer to a condition evoking all second-order derivatives as a "more degenerate" case for instance. The following definition sets this terminology.

Definition 5.1. We call a system of the form (5.1) least degenerate if assumptions (5.2)-(5.4) are satisfied.

Obviously, if (5.1) satisfies (5.2), then (5.3) and (5.4) represents the least degenerate case. Therefore, we will refer to (5.2) and (5.4) as the least degeneracy assumptions, and such a point (s 0 , p 0 ) as the least degenerate root of order m th .

In this chapter, we will study the behaviour of multiple imaginary characteristic roots as (p 1 , p 2 ) varies in a small neighborhood of (p 10 , p 20 ) under the least degeneracy assumptions. In the sequel, we present double roots case (m=2) in Section 5.3 and we discuss imaginary characteristic roots of multiplicity three (m = 3) and four (m = 4) in Section 5.4.

Double root case

We assume that q(s, p) is analytic with respect to s, and continuously differentiable with respect to (s, p) up to a third order. In view of the implicit function theorem, a consequence of inequality (5.4), which is part of the non-degeneracy assumption, is that the characteristic equation (5.1) defines the pairs (p 1 p 2 ) in a small neighbourhood of the critical point p 0 = (p 10 p 20 ) as a function of s, in a sufficiently small neighbourhood of s 0 .

Introduce the notation:

N ε (x 0 ) = {x | |x -x 0 | < ε} , N • ε (x 0 ) = {x | 0 < |x -x 0 | < ε} ,
to describe the concept of neighbourhood of a point x 0 in an n-dimensional space. Then the above can be more precisely stated as follows.

Proposition 5.2. There exists an ε > 0 and a sufficiently small δ > 0 such that for all s ∈ N δ (s 0 ), we may define p 1 (s) and p 2 (s) as the unique solution of (5.1) with (p 1 (s), p 2 (s)) ∈ N ε (p 10 , p 20 ). The functions so defined are differentiable up to the third order.

It should be pointed out that in general, for s ∈ N δ (s 0 ), characteristic equation (5.1) may have other solutions outside of N ε (p 10 , p 20 ).

Recall that stability crossing curves are defined in [START_REF] Gu | On stability crossing curves for general systems with two delays[END_REF] as the set of all points (p 1 , p 2 ) ∈ R 2 + such that q(s) has at least one zero on the imaginary axis. Therefore, the set

T (ω 0 ,p 10 ,p 20 ) = {(p 1 (iω), p 2 (iω)) ∈ N ε (p 10 , p 20 ) | iω ∈ N δ (iω 0 )} ,
which is a curve in the p 1 -p 2 space that passes through the point (p 10 , p 20 ), is the restriction of stability crossing curves to a neighborhood of (p 10 , p 20 ). Therefore, T (s 0 ,p 10 ,p 20 )

will be known as the local stability crossing curve. Roughly speaking, it is a separation curve that divides the parameter space into regions, such that the number of characteristic roots on the right half complex plane remain constant as the parameters vary within each such region.

We will also denote two curves as follows

T + (ω 0 ,p 10 ,p 20 ) = {(p 1 (iω), p 2 (iω)) ∈ N ε (p 10 , p 20 ) | iω ∈ N δ (iω 0 ), ω > ω 0 } , and
T - (ω 0 ,p 10 ,p 20 ) = {(p 1 (iω), p 2 (iω)) ∈ N ε (p 10 , p 20 ) | iω ∈ N δ (iω 0 ), ω < ω 0 } .
The curves T + (ω 0 ,p 10 ,p 20 ) and T - (ω 0 ,p 10 ,p 20 ) will be known as the positive and negative local stability crossing curves, respectively.

Geometric insights: cusp and local bijection

Parameterize a neighborhood of s 0 by a radial variable u and an angle θ. Then, a point in the neighbourhood of s 0 can be written as s = s 0 + ue iθ .

(5.5)

Moreover, p 1 and p 2 can be considered as functions of u and θ. Note that ∂s ∂u = e iθ . For the sake of convenience, we write γ = e iθ = ∂s ∂u .

( = 0, (5.8) in view of (5.4). Equation (5.8) has two important implications.

First, if we set γ = i, the equation (5.8) indicates that the local stability crossing curve T (ω 0 ,p 10 ,p 20 ) may have a cusp at (p 10 , p 20 ) in the parameter space (see [START_REF] Guggenheimer | Differential Geometry[END_REF]). Indeed, as will be confirmed by considering the second-order derivative in the next subsection, T (ω 0 ,p 10 ,p 20 ) partitions a sufficiently small neighborhood of (p 10 , p 20 ) into a great sector (or G-sector) and a small sector1 (or S-sector) as shown in Figure 5.1. We investigate how the double roots at iω 0 migrate as (p 1 , p 2 ) moves from (p 10 , p 20 ) to the G-sector or the S-sector.

To obtain the second implication, we first show the following.

Lemma

∂p 2 ∂u u=0 u + o(u) = 0 + ∂ 2 q ∂s 2 s=s 0 p 1 =p 10 p 2 =p 20 γu + 0 + 0 + o(u),
from which we may conclude (5.9) in view of (5.3) for m = 2.

The implicit function theorem allows us to conclude the following from Lemma 5.3.

Proposition 5.4. Let s a , p 1a and p 2a be defined as in Lemma 5.3. Then there exists a sufficiently small neighborhood of (p 1a , p 2a ) such that the equation (5.1) defines a unique function s(p 1 , p 2 ) with the function value restricted in a small neighborhood of s a . Thus, the second implication of the equation (5.8) may be stated as the following corollary, which is a consequence of Propositions 5.2 and 5.4.

Corollary 5.5. Let s a , p 1a and p 2a be defined as in Lemma 5.3. Then equation (5.1)

defines a bijection between s in a small neighborhood of s a and (p 1 , p 2 ) in a small neighborhood of (p 1a , p 2a ).

Obviously, the small neighborhoods referred in Proposition 5.4 and Corollary 5.5 above should not include s 0 and (p 10 , p 20 ) in view of condition (5.2) for m = 1. Moreover, in view of continuity of solutions of (5.1) with respect to the parameters (p 1 , p 2 ), Corollary 5.5 may be equivalently stated as follows.

Corollary 5.6. For all (p 1 , p 2 ) ∈ N • ε (p 10 , p 20 ) with ε > 0 sufficiently small, the characteristic equation (5.1) has exactly two simple roots in a small neighborhood of s 0 .

Geometric insights: mapping in a neighborhood of a double root

In this section, it will be shown that we can very clearly describe the mapping between s and (p 1 , p 2 ) in the neighborhood of s 0 based on the second order derivative when s -s 0 is restricted to one quadrant. From this description, we may obtain the information on how the double root migrates as (p 1 , p 2 ) moves from (p 10 , p 20 ) to the G-sector or the S-sector in Figure 5.1 according to the sign of D, and whether the negative local stability crossing curve T - (ω 0 ,p 10 ,p 20 ) is on the clockwise side or on the counterclockwise side of T + (ω 0 ,p 10 ,p 20 ) in the S-sector.

Taking a derivative of (5.7) with respect to the radial variable u, we obtain

∂ 2 q ∂p 2 1 ∂p 1 ∂u 2 + 2 ∂ 2 q ∂p 1 ∂p 2 ∂p 1 ∂u ∂p 2 ∂u + 2 ∂ 2 q ∂p 1 ∂s ∂p 1 ∂u γ + ∂q ∂p 1 ∂ 2 p 1 ∂u 2 + ∂ 2 q ∂p 2 2 ∂p 2 ∂u 2 +2 ∂ 2 q ∂p 2 ∂s ∂p 2 ∂u γ + ∂q ∂p 2 ∂ 2 p 2 ∂u 2 + ∂ 2 q ∂s 2 γ 2 = 0. (5.10)
We set u = 0, and apply (5.8)

  -1   Re ∂ 2 q ∂s 2 γ 2 Im ∂ 2 q ∂s 2 γ 2   s=s 0 p 1 =p 10 p 2 =p 20 , (5.11)
which may also be written in a complex form

∂ 2 p 1 ∂u 2 ∂ 2 p 2 ∂u 2 s=s 0 p 1 =p 10 p 2 =p 20 = 1 D   Im ∂q * ∂p 2 ∂ 2 q ∂s 2 γ 2 -Im ∂q * ∂p 1 ∂ 2 q ∂s 2 γ 2   s=s 0 p 1 =p 10 p 2 =p 20 .
(5.12)

In view of (5.8), the tangent of the curve describing (p 1 , p 2 ) as a function of u at (p 10 , p 20 ) is determined by the second order derivative given in (5.11) or (5.12).

Before proceeding further, it is helpful to recall the following well known fact. It can be found in various elementary books that deal with geometry, see for example [START_REF] Gonzalez | A First Course in Continuum Mechanics[END_REF].

Lemma 5.7. Let x (0) ∈ R 2 and M ∈ R 2×2 be fixed. For any x ∈ R 2 , let θ be the angle to rotate x (0) to the direction of x in the counterclockwise direction. Let φ(θ) be the angle to rotate M x (0) to the direction of M x in the counterclockwise direction if det(M ) > 0, and in the clockwise direction if det(M ) < 0. Then the function φ(θ) satisfies the following: i) φ(θ) is a continuous and increasing function of θ ii) 0 < φ(θ) < π if and only if 0 < θ < π.

We now make the following two observations about the second order derivative expression (5.11).

First, set γ = i and γ = -i, the expression determines the tangent of T (ω 0 ,p 10 ,p 20 ) as ω → ω 0 from each side. As Second, as γ rotates through an angle of π/2 radians in a counterclockwise direction, ∂ 2 q ∂s 2 γ 2 rotates through a 180 • angle in the same direction; and

∂ 2 p 1 ∂u 2 ∂ 2 p 2 ∂u 2
T given in (5.11) also rotates through a 180 • angle in a direction determined by the sign of D, which is the determinant of the matrix inverted: the rotation is counterclockwise if D > 0, and it is clockwise if D < 0 (according to Lemma 5.7). With the above observations, and the fact that

p 1 (s) p 2 (s) = p 10 p 20 + u 2 2 ∂ 2 p 1 ∂u 2 ∂ 2 p 2 ∂u 2 s=s 0 p 1 =p 10 p 2 =p 20 + o(u 2 )
we may describe the local mapping (p 1 (s), p 2 (s)) in a very informative manner when s-s 0 is restricted to one quadrant. The situation for s-s 0 in the first quadrant

Q 1 = s = s 0 + ue iθ | 0 < u < δ, 0 ≤ θ ≤ π/2
with D > 0 is illustrated in Figure 5. The complete mapping (p 1 (s), p 2 (s)) with s-s 0 in all four quadrants may be divided into four possible cases depending on the sign of D and whether T - (ω 0 ,p 10 ,p 20 ) is on the counterclockwise or on the clockwise side of T + (ω 0 ,p 10 ,p 20 ) in the S-sector. The migration of the double roots in all cases is summarized in the following theorem. III to III when restricted to III, or from IV to IV when restricted to IV . As the S-sector (in a sufficiently small neighborhood) is contained in II ∩ III , we may conclude that for any (p 1 , p 2 ) in the S-sector, one of the two characteristic roots in the neighborhood of s 0 must be in region II, the other must be in region III, and obviously both in the left half-plane. Similarly, the G-sector (in a sufficiently small neighborhood) is contained in (I ∪ IV ) ∩ (II ∪ IV ). Therefore, for any (p 1 , p 2 ) in the G-sector, one of the two characteristic roots in the neighborhood of s 0 must be in I ∪ IV (in the right half-plane), and the other must be in II ∪ III (in the left half-plane). Case ii is illustrated in Figure 5.4. In this case, the S-sector is contained in I ∩ IV , and therefore, the two characteristic roots in the neighborhood of s 0 must be in regions I and IV , both in the right half-plane. The G-sector can still be expressed as (I ∪ IV ) ∩ (II ∪ IV ).

Case iii is illustrated in Figure 5.5, and Case iv is illustrated in Figure 5.6, and the conclusions can be drawn in a similar manner.

Geometric insights: algebraic S-sector condition

Theorem 5.8 indicates that the migration pattern of the two roots in the G-sector is always the same under the least degeneracy assumptions, which is the only case discussed in this chapter. However, judging the migration pattern of the two roots in the S-sector requires knowing the sign of D and which side of T + (ω 0 ,τ 10 ,τ 20 ) the curve T - (ω 0 ,τ 10 ,τ 20 ) is in the S-sector. Fortunately, by considering the third order derivatives, an explicit algebraic condition is possible. where Proof. Differentiate equation (5.10) with respect to u, we obtain

κ = Re ∂ 2 q ∂s 2 - ∂ 3 q ∂s 3 + 3 ∂ 2 q ∂p 1 ∂s ∂ 2 p 1 ∂u 2 + 3 ∂ 2 q ∂p 2 ∂s ∂ 2 p 2 ∂u 2 s=s 0 p 1 =p 10 p 2 =p
∂ 3 q ∂p 3 1 ∂p 1 ∂u 3 + 3 ∂ 3 q ∂p 2 1 p 2 ∂p 1 ∂u 2 ∂p 2 ∂u + 3 ∂ 3 q ∂p 2 1 ∂s ∂p 1 ∂u 2 γ + 3 ∂ 2 q ∂p 2 1 ∂p 1 ∂u ∂ 2 p 1 ∂u 2 +3 ∂ 3 q ∂p 1 ∂p 2 2 ∂p 1 ∂u ∂p 2 ∂u 2 + 6 ∂ 3 q ∂p 1 ∂p 2 ∂s ∂p 1 ∂u ∂p 2 ∂u γ + 3 ∂ 2 q ∂p 1 ∂p 2 ∂ 2 p 1 ∂u 2 ∂p 2 ∂u +3 ∂ 2 q ∂p 1 ∂p 2 ∂p 1 ∂u ∂ 2 p 2 ∂u 2 + 3 ∂ 3 q ∂p 1 ∂s 2 ∂p 1 ∂u γ 2 + 3 ∂ 2 q ∂p 1 ∂s ∂ 2 p 1 ∂u 2 γ + ∂q ∂p 1 ∂ 3 p 1 ∂u 3 + ∂ 3 q ∂p 3 2 ∂p 2 ∂u 3 + 3 ∂ 2 q ∂p 2 ∂s ∂ 2 p 2 ∂u 2 γ + ∂q ∂p 2 ∂ 3 p 2 ∂u 3 + ∂ 3 q ∂s 3 γ 3 = 0. (5.15)
We set u = 0 and use (5.8) where

B = ∂ 3 q ∂s 3 γ 3 + 3 ∂ 2 q ∂p 1 ∂s ∂ 2 p 1 ∂u 2 γ +3 ∂ 2 q ∂p 2 ∂s ∂ 2 p 2 ∂u 2 γ s=s 0 p 1 =p 10 p 2 =p 20 . Let ∂ k p 1 ∂u k ∂ k p 2 ∂u k ± = ∂ k p 1 ∂u k ∂ k p 2 ∂u k s=s 0 p 1 =p 10 p 2 =p 20 γ=±i , k = 1, 2, 3,
and

p 1 p 2 ± = p 1 (s 0 ± δi) p 2 (s 0 ± δi)
.

Then the Taylor series gives

p 1 p 2 ± = p 10 p 20 + δ ∂p 1 ∂u ∂p 2 ∂u ± + δ 2 2 ∂ 2 p 1 ∂u 2 ∂ 2 p 2 ∂u 2 ± + δ 3 6 ∂ 3 p 1 ∂u 3 ∂ 3 p 2 ∂u 3 ± + o(δ 3 ).
But according to (5.8) and (5.11), we have

∂p 1 ∂u ∂p 2 ∂u ± = 0, ∂ 2 p 1 ∂u 2 ∂ 2 p 2 ∂u 2 + = ∂ 2 p 1 ∂u 2 ∂ 2 p 2 ∂u 2 - . Therefore, ∆p 1 ∆p 2 ∆ = p 1 p 2 + - p 1 p 2 - = δ 3 6   ∂ 3 p 1 ∂u 3 ∂ 3 p 2 ∂u 3 + - ∂ 3 p 1 ∂u 3 ∂ 3 p 2 ∂u 3 -   + o(δ 3 ) = - δ 3 6   Re ∂q ∂p 1 Re ∂q ∂p 2 Im ∂q ∂p 1 Im ∂q ∂p 2   -1 s=s 0 p 1 =p 10 p 2 =p 20
Re(∆B)

Im (∆B) + o(δ 3 ), (5.17)
where Let -

∆B = B| γ=i -B| γ=-i = 2i - ∂ 3 q ∂s 3 + 3 ∂ 2 q ∂p 1 ∂s ∂ 2 p 1 ∂u 2 + + 3 ∂ 2 q
∂ 2 q ∂s 2 0 = -∂ 2 q ∂s 2 s=s 0 p 1 =p 10 p 2 =p 20
. Comparing the expressions (5.17) and (5.11) and using Lemma 5.7, we can see that the above can be achieved if we can reach the direction of -∂ 2 q ∂s 2 0 by rotating ∆B counterclockwise through an angle of θ ∈ (0, π) if D > 0 (which is Case i in Theorem 5.8). The rotation from ∆B to -∂ 2 q ∂s 2 0 needs to be clockwise if D < 0 (which is Case iii). The counterclockwise rotation from ∆B to -∂ 2 q ∂s 2 0 may be expressed as

Re(∆B) Im - ∂ 2 q ∂s 2 0 -Im(∆B) Re - ∂ 2 q ∂s 2 0 > 0,
which is equivalent to (5.13), and the conclusion is valid in this case in view of Case i in Theorem 5.8. It can be similarly shown that if we can rotate ∆B to the direction of -∂ 2 q ∂s 2 0

clockwise through an angle of θ ∈ (0, π), then (5.14) is satisfied, and the conclusion is valid in this case also in view of Case iii in Theorem 5.8. Similarly, we can show that κ > 0 and D > 0, or κ < 0 and D < 0 can guarantee that we can reach the direction of If κ = 0, higher order derivatives may be used to evaluate conditions in Theorem 5.8.

∂ 2 p 1 ∂u 2 , ∂ 2 p 2
It should be noticed that the roots of the characteristic equation discussed in Theorem 5.8 and Corollary 5.9 are restricted to a sufficiently small neighborhood of s 0 = jω 0 . Because characteristic roots are distributed symmetrically with respect to the real axis, there is also a double root at s * 0 = -jω 0 when p 1 = p 10 , p 2 = p 20 . When (p 1 , p 2 ) deviates from (p 10 , p 20 ), the migration of the two roots in the neighborhood of s * 0 follows the same pattern as those in the neighborhood of s 0 .

There may also be roots on the imaginary axis outside of the neighborhoods of s 0 and s * 0 . The migration of these imaginary roots needs to be analyzed separately.

Finally, the roots on the right half-plane remain on the right half-plane as long as (p 1 , p 2 ) stay within a sufficiently small neighborhood of (p 10 , p 20 ). Similarly, the roots on the left half-plane remain on the left half-plane when the deviation of (p 1 , p 2 ) is sufficiently small.

Illustrative example

In this subsection, we consider an example of time-delay system with two delay parameters and we show how to apply the results discussed in this chapter on such a system. (5.18)

Note that the function given in (5.18) is of the form (2.9), with polynomials r 0 , r 1 and r 2 of order 2, 1, and 0, respectively. For (τ 1 , τ 2 ) = (1, 2), system (5.18) has double imaginary roots at s = s 0 = ±iω 0 with ω 0 = 1. It can be computed that

D 1.74159 > 0 κ 30.7082 > 0 .
From the sign of D and κ, it can be concluded that this system belongs to Case ii of Theorem 5.8, i.e., T - (ω 0 ,τ 10 ,τ 20 ) is on the clockwise side of T + (ω 0 ,τ 10 ,τ 20 ) in the S-sector.

The stability crossing curve T (which contains both positive and negative local stability crossing curves) is plotted in Figure 5.8. It can be seen that T divides the region into three regions: region A is connected to the origin, region B is the small region on the upper side, and region C is the small region on the lower side. For τ 1 = 0, τ 2 = 0, the characteristic quasipolynomial is reduced to a polynomial, and it can be easily calculated that both roots are on the right half-plane. Therefore, the region connected to the origin has two right half-plane roots.

According to Corollary 5.9 or Theorem 5.8, both imaginary roots move to the right halfplane as (τ 1 , τ 2 ) moves to the S-sector (which is connected to region B). According to Theorem 5.8, as (τ 1 , τ 2 ) moves to the G-sector (which is connected to region A), one of the two imaginary roots moves to the right half-plane, the other moves to the left halfplane. In other words, as (τ 1 , τ 2 ) moves from region B to region A through (1, 2), one root moves from the right half-plane to the left half-plane passing through the point i on the imaginary axis, another root on the right half-plane moves to touch the imaginary axis at i then returns to the right half-plane.

Due to symmetry, another left half-plane root moves to the right half-plane through the point -i. Therefore, there are two more right half-plane roots when (τ 1 , τ 2 ) is in region B as compared to the case when (τ 1 , τ 2 ) is in region A. Thus, we conclude that there are four roots on the right half-plane when (τ 1 , τ 2 ) is in region B.

For region C, it can be easily calculated using the method described in [START_REF] Gu | On stability crossing curves for general systems with two delays[END_REF] that the two right half-plane roots cross the imaginary axis to the left half-plane as (τ 1 , τ 2 ) moves from region A. Therefore, there is no right half-plane root in this region, and the system is stable for (τ 1 , τ 2 ) in region C.

Link with the algebraic approach

We have presented, throughout this chapter, a geometric approach allowing us to analyze the behaviour of double characteristic roots for time-delay systems with two delays, when a small perturbation occurs on the parameters. However, we can find other studies in the literature aiming to identify the direction in which double roots cross the imaginary axis when the delay parameters change. We cite here, for instance, Chapter 5 of [START_REF] Michiels | Stability and Stabilization of Time-Delay Systems (Advances in Design and Control)[END_REF] and [START_REF] Morȃrescu | Qualitative Analysis of Distributed Delay Systems: Methodology and Algorithms[END_REF], where an algebraic approach is derived. In this subsection, we compare the two approaches and discuss how they are conceptually different. We check in Example 5.1 that both methods provide consistent results, and highlight the corresponding limitations and open questions.

Consider a linear time-delay system

ẋ(t) = A 0 x(t) + A 1 x(t -τ 1 ) + A 2 x(t -τ 2 ), (5.19)
where x ∈ R n is the state vector, the constant delays τ 1 and τ 2 are real and positive, with A 0 , A 1 , A 2 ∈ R n×n . The characteristic matrix of the system (5. [START_REF] Montoya | Ecological networks and their fragility[END_REF]) is given by (2.5).

In particular, the characteristic matrix for the considered time-delay system with two delays is

M (s, τ ) = sI -A 0 -A 1 e -sτ 1 -A 2 e -sτ 2 , (5.20)
where τ = (τ 1 τ 2 ) T , I is the n × n identity matrix, and s is the Laplace variable.

Here we suppose that real matrices A 0 , A 1 , A 2 are such that the characteristic function p(s, τ ) := det M (s, τ ) (also called quasi-polynomial) of system (5.19) has the form q(s, τ ) = p 0 (s) + p 1 (s)e -sτ 1 + p 2 (s)e -sτ 2 , (5.21)

where p k (s), k ∈ {0, 1, 2} are polynomials of s with real coefficients. is satisfied, the characteristic function of system (5.19) can be written of the form (5.21).

Note that conditions (C1)-(C4) are satisfied when matrices A 1 and A 2 are of rank 1.

Recall that, similar to the finite-dimensional case, the characteristic equation of system (5.19) is given by (2.4). The roots of the characteristic equation are called characteristic roots of system (5.19). Thus, the finite-dimensional nonlinear eigenvalue problem of system (5.19) can be written as M (s, τ )u = 0, (5.22) where the vector u ∈ C n \ {0} is called a right eigenvector corresponding to the characteristic root s. Note that, in a similar manner, we can construct a left eigenvector By considering (as, for instance, in [START_REF] Hale | Introduction to functional differential equations[END_REF]) a linear operator A = A 0 x(t) + A 1 x(t -τ 1 ) +

v T ∈ C n \ {0} corresponding to the eigenvalue s, satisfying v T M (s, τ ) = 0.
A 2 x(t -τ 2 ), with the domain

D(A) = φ φ, dφ dθ ∈ C ([-max{τ 1 , τ 2 }, 0] , R n ) , φ(0) = A 0 φ(0) + A 1 φ(-τ 1 ) + A 2 φ(-τ 2 )} , such that Aφ = dφ dθ , where C ([-max{τ 1 , τ 2 }, 0] , R n ) is the space of continuous func- tions from [-max{τ 1 , τ 2 }, 0] to R n ,
we can rewrite system (5.19) as an abstract ordinary differential equation ẋt = Ax t , where

x t (φ) ∈ C ([-max{τ 1 , τ 2 }, 0] , R n ) is the func- tion segment defined by x t (φ)(θ) = x(φ)(t + θ), θ ∈ [-max{τ 1 , τ 2 }, 0] . Note that
the characteristic roots are the eigenvalues of A. This allows us to write an infinitedimensional and linear eigenvalue problem equivalent to (5.22), (sI -A)u = 0, with

u ∈ C ([-max{τ 1 , τ 2 }, 0] , C n ).
In addition, the corresponding eigenfunctions of eigenvalue s take the form ue sθ , with θ ∈ [-max{τ 1 , τ 2 }, 0]. However, the algebraic approach described in the sequel is based on the finite-dimensional and nonlinear form of eigenvalue problem (5.22).

In the sequel, we consider the characteristic matrix (5.23) in Example 5.1. A right eigenvector u 0 satisfying (5.24) is of the form

u 0 = α 1 + iα 2 -(cos(2) -i(-1 + sin(2))) (α 1 + iα 2 )
, where α 1 , α 2 ∈ R. We compute the generalized right eigenvector u 1 = (β 1 + iβ 2 u 12 ) T , verifying (5.25), with β 1 , β 2 ∈ R and

u 12 = (cos(2) -i sin(2))α 1 + (i cos(2) + sin(2))α 2 - -(i + cos(2) -i sin(2)) (β 1 + iβ 2 ) .
In the same manner we write the left eigenvector v T 0 satisfying equation (5.26)

v 0 = (1 -i) (γ 1 + iγ 2 ) γ 1 + iγ 2 ,
and the left generalized eigenvector v T 1 such that condition (5.27) holds:

v 1 = -γ 1 -iγ 2 + (1 -i)δ 1 + (1 + i)δ 2 δ 1 + iδ 2 , with γ 1 , γ 2 , δ 1 , δ 2 ∈ R.
The normalization condition (5.42) leads to the following constraint on α 1 , α 2 , γ 1 , and γ 2 :

α 1 = 0, α 2 = - 1 γ 2 (cos 2 + sin 2 tan 2) , γ 1 = -γ 2 tan 2, γ 2 = 0.
Therefore, the normalized eigenvectors satisfy

u 0 = 1 γ 2 -i cos 2 cos 2(-1 + i cos 2 + sin 2) , u 12 = cos 2(-i cos 2 -sin 2) γ 2 -(cos 2 -i(-1 + sin 2)) (β 1 + iβ 2 ) , v 0 = γ 2 (1 + i)(1 + i tan 2) -tan 2 + i , v 1 = (-i + tan 2)γ 2 + (1 -i) (δ 1 + iδ 2 ) δ 1 + iδ 2 .
Now we write M 0 2 and M 1 1 of the form

M 0 2 = ∂ 2 ∂s 2 M (s, τ ) s=s 0 τ =τ 0 , M 1 1 = τ 11 ∂ ∂τ 1 M (s, τ ) s=s 0 τ =τ 0 + τ 21 ∂ ∂τ 2 M (s, τ ) s=s 0 τ =τ 0 ,
and recover the explicit formulae (5.40) and (5.41), respectively. More precisely, we obtain

M 0 2 = 2e -i cos 1 0 -i -(4 -i)e -2i 0 , M 1 1 = -2ie -i cos 1τ 11 0 e -2i iτ 21 --1 + e 2i τ 11 0 .
We compute s 1 as a function of τ 11 and τ 21 , using equation (5.45):

s 1 = ± (2i + (4 + 2i)e 2i ) τ 11 -2iτ 21 -1 + (1 -2i)e 2i .
Using formula (5.46), we write s 2 as a function of τ 11 and τ 21 We notice that if we fix τ 1 and increase τ 2 , then the double root splits up into two simple roots moving towards instability (with positive real part), as depicted in Figure 5.9, left.

On the other hand, if we fix τ 1 and decrease τ 2 , then the two simple roots will move one towards instability, and the other one towards stability, as illustrated in Figure 5.9, right. Figure 5.9 was obtained by using the QPmR algorithm developed by [START_REF] Vyhlídal | Mapping based algorithm for large-scale computation of quasi-polynomial zeros[END_REF]. Characteristic roots s(ε) are represented by coloured points. Yellow points s = ±i correspond to ε = 0, τ 1 = 1 and τ 2 = 2. We remark that these points become red as τ 2 increases, i.e. as perturbation ε increases. Fix τ 1 = 1 (i.e. τ 11 = 0). For case (a), the characteristic roots of the quasi-polynomial are computed by using the QPmR algorithm, for each value of τ 2 , from 2 to 3 using a step of 10 -3 . We can see that double root s = ±i splits up into two simple roots that will move towards the right half-plane. This is coherent to our computation of s(ε) of the form (5.34), using equations (5.45) and (5.46). For case (b), we decrease τ 2 from 2 (yellow points) to 1.9 (red points) using a step of 10 -3 , and we notice in Figure 5.9 right that the double root splits up into two simple roots going in opposite directions, one to the left half-plane, and the other one to the right half-plane. This is also consistent with our prediction based on computing s(ε) as in (5.34).

Re(s) case (a)

We notice that despite the fact that the two methods (algebraic and geometric approaches) are rather different, they both point out that characteristic roots follow the same pattern when the parameters are perturbed. The former method is based on eigenvector computation. Even if in this section we have only discussed the non semisimple double roots case, a similar study can be made for double semi-simple characteristic roots. The latter method involves the computation of partial derivatives of the quasi-polynomial with respect to the delays, up to the third order. Using the geometric approach implies that condition (5.4) must be satisfied. Nonetheless, the geometric approach applies to both semi-simple and non semi-simple characteristic roots. In other words, the main limitation when using the algebraic approach is that we have to separately treat semi-simple and non semi-simple double characteristic roots, and the geometric approach's main limitation is that assumption (5.4) has to be satisfied.

Both approaches show that there are typically two types of situations when a deviation on the delays occurs in the parameter space: either the two imaginary roots move to the same half-plane of the complex plane, or one of the imaginary roots moves to the left half-plane and the other one moves to the right half-plane. This is rather surprising in the case of double non semi-simple characteristic roots: the trajectory of such a root when a perturbation ε occurs is described by the equation (5.34), where s 0 lies on the imaginary axis, and s 1 is given by equation (5.45). Note that the double root splits up into two simple roots, the positive real part s 1 corresponding to a root, and the s 1 with negative real part corresponding to the other simple root. Some publications in the literature suggest that this is why non semi-simple double roots always split up into two simple roots following only one type of tendency, which is to move towards different half-planes of the complex plane. The main message of this section is that we have to pay attention to the quantification of the small neighbourhood, which remains an open problem at least for the algebraic approach. We can briefly illustrate this in Example 5.1, as follows.

We consider a point in S-sector, (τ 1 , τ 2 ) = (1.1, 2.6), as depicted in Figure 5.8. We have seen in this chapter that in this case both imaginary roots move to the right halfplane. We can verify that we obtain the same result by using the algebraic approach, by choosing ε = 0.1, τ 11 = 1 and τ 21 = 6 in equations (5.32)-(5.33). Thus, we compute s(ε)

written of the form (5.34), yielding the value of the two roots s = 0.975521 -0.254494i and s = 0.826199 -1.58361i. As both roots have positive real part, they are both in the right half-plane. Of course, we can always find an extremely small value for ε such that we end up with two simple roots lying on different sides of the imaginary axis, but this choice gives rise to a natural question regarding the relevance of considering such a value for the control community. We conclude by saying that when choosing values for ε, τ 11 and τ 21 we have to decide upon what a small neighbourhood means, which remains an open question.

Extension: multiplicity three and four

The geometric approach presented throughout this chapter can be extended for characteristic roots with multiplicity greater than two. In the sequel, we briefly discuss the case of multiplicity three and four for characteristic roots lying on the imaginary axis.

Multiplicity three

In this subsection, we study the migration of triple roots. and on the counterclockwise side of T - (ω 0 ,p 10 ,p 20 ) . In this case, two characteristic roots of (2.8) move to the right-half complex plane, and the third root moves to the left-half plane.

Case ii. D > 0 and (p 1 p 2 ) moves in the region on the clockwise side of T - (ω 0 ,p 10 ,p 20 ) and on the counterclockwise side of T + (ω 0 ,p 10 ,p 20 ) . In this case, two characteristic roots of (2.8) move to the left-half complex plane, and the third root moves to the right-half plane.

Case iii. D < 0 and (p 1 p 2 ) moves in the region on the clockwise side of T - (ω 0 ,p 10 ,p 20 ) and on the counterclockwise side of T + (ω 0 ,p 10 ,p 20 ) . In this case, two characteristic roots of (2.8) move to the right-half complex plane, and the third root moves to the left-half plane.

Case iv. D < 0 and (p 1 p 2 ) moves in the region on the clockwise side of T + (ω 0 ,p 10 ,p 20 ) and on the counterclockwise side of T - (ω 0 ,p 10 ,p 20 ) . In this case, two characteristic roots of (2.8) move to the left-half complex plane, and the third root moves to the right-half plane.

Proof. In the complex plane consider a point s in the neighbourhood of s 0 , as described by equation (5.5). Similar to the way we obtained (5.8) from (5.7), we may conclude from (5.10) using (5.2) for k = 2 and equation (5.8) that

  ∂ 2 p 1 ∂u 2 ∂ 2 p 2 ∂u 2   u=0 = 0.
(5.47) Differentiating (5.10) again with respect to u yields

∂ 3 q ∂p 3 1 ∂p 1 ∂u 3 + 3 ∂ 2 q ∂p 2 1 ∂ 2 p 1 ∂u 2 ∂p 1 ∂u + 3 ∂ 3 q ∂p 2 1 ∂p 2 ∂p 1 ∂u 2 ∂p 2 ∂u + 3 ∂ 3 p ∂p 2 1 ∂s ∂p 1 ∂u 2 γ+ +3 ∂ 2 p ∂p 1 ∂p 2 ∂ 2 p 1 ∂u 2 ∂p 2 ∂u + 3 ∂ 3 p ∂p 1 ∂p 2 2 ∂p 1 ∂u ∂p 2 ∂u 2 + 3 ∂ 2 q ∂p 1 ∂p 2 ∂p 1 ∂u ∂ 2 p 2 ∂u 2 + +6 ∂ 3 q ∂p 1 ∂p 2 ∂s ∂p 1 ∂u ∂p 2 ∂u γ + 3 ∂ 3 q ∂p 1 ∂s 2 ∂p 1 ∂u γ 2 + ∂q ∂p 1 ∂ 3 p 1 ∂u 3 + 3 ∂ 2 q ∂p 1 ∂s ∂ 2 p 1 ∂u 2 γ+ + ∂ 3 p ∂p 3 2 ∂p 2 ∂u 3 + 3 ∂ 2 q ∂p 2 2 ∂p 2 ∂u ∂ 2 p 2 ∂u 2 + 3 ∂ 3 p ∂p 2 2 ∂s ∂p 2 ∂u 2 γ + 3 ∂ 2 p ∂p 2 ∂s ∂ 2 p 2 ∂u 2 γ+ +3 ∂ 3 q ∂p 2 ∂s 2 ∂p 2 ∂u γ 2 + ∂q ∂p 2 ∂ 3 p 2 ∂u 3 + ∂ 3 q ∂s 3 γ 3 = 0. (5.48)
If we set u = 0 and use (5.8) and (5.47) in equation ( 5.15), we obtain

∂q ∂p 1 ∂ 3 p 1 ∂u 3 + ∂q ∂p 2 ∂ 3 p 2 ∂u 3 + ∂ 3 q ∂s 3 γ 3 s=s 0 p 1 =p 10 p 2 =p 20 = 0 or ∂q ∂p 1 ∂ 3 p 1 ∂u 3 + ∂q ∂p 2 ∂ 3 p 2 ∂u 3 s=s 0 p 1 =p 10 p 2 =p 20 = - ∂ 3 q ∂s 3 γ 3 s=s 0 p 1 =p 10 p 2 =p 20
.

We separate real and imaginary part to obtain

   Re ∂q ∂p 1 Re ∂q ∂p 2 Im ∂q ∂p 1 Im ∂q ∂p 2    s=s 0 p 1 =p 10 p 2 =p 20   ∂ 3 p 1 ∂u 3 ∂ 3 p 2 ∂u 3   s=s 0 p 1 =p 10 p 2 =p 20 = -    Re ∂ 3 q ∂s 3 γ 3 Im ∂ 3 q ∂s 3 γ 3    s=s 0 p 1 =p 10 p 2 =p 20 . Thus   ∂ 3 p 1 ∂u 3 ∂ 3 p 2 ∂u 3   s=s 0 p 1 =p 10 p 2 =p 20 = -    Re ∂q ∂p 1 Re ∂q ∂p 2 Im ∂q ∂p 1 Im ∂q ∂p 2    -1 s=s 0 p 1 =p 10 p 2 =p 20    Re ∂ 3 p ∂s 3 γ 3 Im ∂ 3 p ∂s 3 γ 3    s=s 0 p 1 =p 10 p 2 =p 20 . (5.49)
Using Lemma 5.7 and in view of (5.49), we know that a π/2 counterclockwise rotation of γ in the complex plane will generate a 3π/2 rotation in τ 1 -τ 2 parameter space, in the counterclockwise direction if D > 0, and in the clockwise direction if D < 0.

Accounting for higher order terms, the situation is illustrated in Figure 5.10 for Cases i and ii (D > 0), and in Figure 5.11 for Cases iii and iv (D < 0). In both Figures 5.10 (ω 0 ,p 10 ,p 20 ) ) in the neighbourhood of C (or (p 10 , p 20 )) may be expressed as I ∩ (II ∪ III ) ∩ IV . Therefore, for any (p 1 , p 2 ) in this region, there must be one root in region I, one root in either region II or region III, and one root in region IV . In other words, there must be two roots on the right half plane, and one root on the left half plane. This proves Case i. Case ii can be shown by noticing that the region on the clockwise side of Im(-) (or T - (ω 0 ,p 10 ,p 20 ) ) and on the counterclockwise side of Im(+) (or T + (ω 0 ,p 10 ,p 20 ) ) in the neighborhood of (p 10 , p 20 ) may be expressed as (I ∪ IV ) ∩ II ∩ III . Cases iii and iv may be shown similarly.

Remark 5.19. Note that

∂ 3 p 1 ∂u 3 γ=-i = - ∂ 3 p 1 ∂u 3
γ=i in view of (5.49). This means, in view of (5.8) and (5.47), that T + (ω 0 ,p 10 ,p 20 ) has the same tangent as T - (ω 0 ,p 10 ,p 20 ) at (p 10 p 20 ). Thus, T (ω 0 ,p 10 ,p 20 ) is a smooth curve. In other words, unlike the double root case discussed in Section 5.3, the stability crossing curve is smooth without a cusp at (p 10 p 20 ).

Multiplicity four

In this section we study the migration of quadruple roots. For system (5.1), s 0 is a quadruple root if conditions (5.2) and (5.3) hold for m = 4.

Parameterize s by u and θ (or γ) as in (5.5). From (5.8), (5.47) and (5.49), we immedi-

ately conclude   ∂ k p 1 ∂u k ∂ k p 2 ∂u k   u=0 = 0 for k = 1, 2, 3. (5.50)
The above is true for k = 3 due to (5.49) and equation ( 5 .

This can be solved to obtain

  ∂ 4 p 1 ∂u 4 ∂ 4 p 2 ∂u 4   s=s 0 p 1 =p 10 p 2 =p 20 = -    Re ∂q ∂p 1 Re ∂q ∂p 2 Im ∂q ∂p 1 Im ∂q ∂p 2    -1 s=s 0 p 1 =p 10 p 2 =p 20    Re ∂ 4 p ∂s 4 γ 4 Im ∂ 4 p ∂s 4 γ 4    s=s 0 p 1 =p 10 p 2 =p 20 . (5.51)
Similar to the triple root case, the last equation above shows that

∂ 4 p 1 ∂u 4 ∂ 4 p 2 ∂u 4
T rotates four times as fast as γ does. To understand this case, we shall divide the circle in s domain in 45 • pieces in the complex plane, in order to work with singly connected regions (see Figures 5.12 to 5.15, left).

Considering (5.50) and (5.51) for γ = i and γ = -i, we see that the local stability crossing curve T (ω 0 ,p 10 ,p 20 ) have a cusp at (p 10 , p 20 ) [START_REF] Guggenheimer | Differential Geometry[END_REF]. Indeed, T (ω 0 ,p 10 ,p 20 ) partitions a sufficiently small neighborhood of (p 10 , p 20 ) into a great sector (or G-sector) and a and V III, which are all in the left-half plane. So we may conclude that when (p 1 p 2 ) is in the S-sector, the characteristic equation (5.1) has a root in the right-half plane, and three others in the left-half plane. As for the G-sector, Figure 5.12 shows that it can be expressed as

(I ∪ II ) ∩ (III ∪ IV ) ∩ (V ∪ V I ) ∩ (V II ∪ V III ).
Thus, the characteristic equation (5.1) has two unstable roots in G-sector, within the regions (I ∪ II) and

(III ∪ IV ), and two stable roots, within the regions (V ∪ V I) and (V II ∪ V III).

Case ii: The S-sector can be expressed as

I ∩ (II ∪ III ) ∩ IV ∩ (V I ∪ V II ), as
shown in Figure 5.13. Therefore, for any (p 1 p 2 ) in S-sector, one characteristic root must be in (V I ∪ V II) (in the left half-plane), and the remaining three roots in right half-plane (one in I, one in II ∪ III, and one in IV ). Next, G-sector can be expressed

as (I ∪ II ) ∩ (III ∪ IV ) ∩ (V ∪ V I ) ∩ (V II ∪ V III
). Therefore, we can conclude that there are two roots on the left-half plane and two roots on the right-half plane.

For case iii and case iv, the conclusions can be drawn in a similar manner. Case iii is illustrated in Figure 5.14. S-sector can be expressed as I ∩(II ∪III )∩IV ∩(V I ∪V II ), and G-sector as

(I ∪ II ) ∩ (III ∪ IV ) ∩ (V ∪ V I ) ∩ (V II ∪ V III ).
Case iv is depicted in Figure 5.15, S-sector can be expressed as

(II ∪ III ) ∩ V ∩ (V I ∪ V II ) ∩ V III , and G-sector as (I ∪ II ) ∩ (III ∪ IV ) ∩ (V ∪ V I ) ∩ (V II ∪ V III ).

Illustrative example

In this subsection we consider an example of a systems with two delays as parameters, τ 1 and τ 2 . The corresponding characteristic equation has the form (2.9). In other words, consider the quasi-polynomial q(s, τ 1 , τ 2 ) = s 4 + a 03 s 3 + a 02 s 2 + a 01 s + a 00 + + a 12 s 2 + a 11 s + a 10 e -sτ 1 + + (a 21 s + a 20 ) e -sτ 2 .

The system has a triple imaginary root at s = s 0 = iω 0 , with ω 0 = 1, for (τ

1 τ 2 ) = (3, 5), a 03 = 1, a 12 = 1, a 21 = 1
, and the values of other coefficients are given in As depicted in Figure 5.16, the local stability crossing curves divide the neighbourhood of [START_REF] Postle | Nucleotide sequence of the repressor gene of the tn10 tetracycline resistance determinant[END_REF][START_REF] Kennell | Transcription and translation initiation frequencies of the Escherichia coli lac operon[END_REF] in the τ 1 -τ 2 plane into two regions. Next, it can be calculated that D > 0.

Therefore, for (τ 1 τ 2 ) taking values in the region below the curve, which is on the clockwise side of T - (ω 0 ,τ 10 ,τ 20 ) and on the counterclockwise side of T + (ω 0 ,τ 10 ,τ 20 ) , two roots

Coeff. Exact value will move in the left-half plane, and one root in the right-half plane. Similarly, for (τ

1 τ 2 )
taking values above the curve, two roots will move on the right-half plane, and one root on the left-half plane.

Concluding remarks

The migration of double roots of characteristic equations that depend on two parameters (which are two delays) is studied in Section 5.3 under the least degeneracy assumption.

It is shown that in the parameter space, the local stability crossing curve has a cusp and divides the neighbourhood of the critical point into two regions: an S-sector and a G-sector. As the parameter pair moves to the G-sector, one root moves to the left half-plane and the other moves to the right half-plane. If the parameter pair moves to the S-sector, a simple algebraic criterion may be used to judge whether both roots move to the right half-plane or the left half-plane. The proposed approach is also compared to

Motivation and functionality

The engineering of synthetic gene circuits with predicted or desired functionality in living cells remains an outstanding problem of the synthetic biology field. In the last two decades, an increasing number of studies have focused on designing biological circuits implementing new behaviours, such as controlled oscillations [START_REF] Michael | A synthetic oscillatory network of transcriptional regulators[END_REF][START_REF] Danino | A synchronized quorum of genetic clocks[END_REF] and counters [START_REF] Smolke | It's the dna that counts[END_REF][START_REF] Friedland | Synthetic gene networks that count[END_REF]. Recent studies abound in designing logic gates: NOT gates [START_REF] Wang | Engineering modular and orthogonal genetic logic gates for robust digital-like synthetic biology[END_REF][START_REF] Tamsir | Robust multicellular computing using genetically encoded nor gates and chemical 'wires[END_REF], AND gates [START_REF] Seok Moon | Genetic programs constructed from layered logic gates in single cells accessed terms of use detailed terms[END_REF][START_REF] Wang | Engineering modular and orthogonal genetic logic gates for robust digital-like synthetic biology[END_REF][START_REF] Tamsir | Robust multicellular computing using genetically encoded nor gates and chemical 'wires[END_REF][START_REF] Siuti | Engineering genetic circuits that compute and remember[END_REF], NAND gates [START_REF] Wang | Engineering modular and orthogonal genetic logic gates for robust digital-like synthetic biology[END_REF][START_REF] Tamsir | Robust multicellular computing using genetically encoded nor gates and chemical 'wires[END_REF], OR, NOR, and XOR gates [START_REF] Tamsir | Robust multicellular computing using genetically encoded nor gates and chemical 'wires[END_REF]. Moreover, it has been shown that gates can be layered in order to implement more complex circuits; for instance an AND gate with four inputs was obtained using three AND gates with two inputs each in [START_REF] Seok Moon | Genetic programs constructed from layered logic gates in single cells accessed terms of use detailed terms[END_REF], and circuits with five inputs were designed in [START_REF] Xie | Multiinput rnai-based logic circuit for identification of specific cancer cells[END_REF]. Designing new, more complex behaviours is important in synthetic biology: in [START_REF] Fernandez-Rodriguez | Post-translational control of genetic circuits using potyvirus proteases[END_REF] there is an attempt to control two signals with one single input. In this work we design a new circuit: we control two signals using three inputs. One of the inputs is the one that selects the expected output.

This new device can have the functionality of a multiplexer (or the selector part of a multiplexer, see for instance 74LVC1G157 data sheet, page 3). There are a few multiplexing and demultiplexing logic circuits in the literature, but usually they have only one single controlled output. For instance, in [START_REF] Pasotti | Multiplexing and demultiplexing logic functions for computing signal processing tasks in synthetic biology[END_REF] both inputs generate the same product.

The multiplexer from [START_REF] Pasotti | Multiplexing and demultiplexing logic functions for computing signal processing tasks in synthetic biology[END_REF] is also considered in a (theoretical) modular design in [START_REF] Huynh | Sbrome: A scalable optimization and module matching framework for automated biosystems design[END_REF].

Other multiplexer implementations can be found in [START_REF] Seok Moon | Construction of a genetic multiplexer to toggle between chemosensory pathways in Escherichia coli[END_REF] and [START_REF] Sergi Regot | Distributed biological computation with multicellular engineered networks[END_REF]. The contribution of our device is that it has two different outputs, designed to select only one of two incompatible biochemical pathways that could be carried out in a single cell. The idea is that at most one output can be selected at a time.

In other words, we would like to design a gene network similar to a multiplexer circuit.

In Figure 6.1, the gene network is represented as a black box. There are three input signals (signal A, signal B and signal S), and two outputs (signal a and signal b). The expected behaviour is as follows: output a is activated when signal A is activated, and more than this, when signal S is off; output b is activated when input B and input S are both on. Otherwise, no output is selected. The complete truth table is:

A B S output 0 0 0 0 0 0 1 0 0 1 0 0 0 1 1 b 1 0 0 a 1 0 1 0 1 1 0 a 1 1 1 b
In this truth table, 1 means that a signal is activated (or on), and 0 means that the signal is off. Note that signal S is used to select an output (signal a when S is off, and signal b when S is on). A short version of the true table would then be case A B S output

(i) 1 0 a (ii) 1 1 b
where can be either 1, or 0. Any other combination of input signals than case (i) or case (ii) would produce 0.

In terms of utility, the device can be used to select one of two pathways, but does not allow both to operate at the same time. This would be useful to allow two incompatible pathways to coexist in a cell. In the natural world, this happens for example in some cyanobacteria where it is necessary to alternate between photosynthesis and nitrogen fixation because one of the key enzymes that carries out nitrogen fixation is damaged by oxygen [START_REF] Adams | Heterocyst formation in cyanobacteria[END_REF]. Thus, the circuit can be use to segregate two incompatible metabolic activities, like nitrogen fixation and photosynthesis, or other pairs of biochemically incompatible activities that could be carried out in a single cell.

Another possible application would be the block copolymerization (the formation of block copolymers, consisting of two or more blocks of different polymers clustered together [START_REF] Mcnaught | Compendium of Chemical Terminology[END_REF]). The circuit could be used to choose which of two polymerization processes is active, and therefore build up a block copolymer, essentially by switching from one polymerization process to the other. However, the ultimate goal of this study would be to put this gene network into living cells. because of their properties: the fast growth rate (one generation per twenty minutes on average), the genetic simplicity (E. coli cells have about 4400 genes), and their reduced cost.

Proposed model

The model is depicted in Figure 6.2. The three inputs are Ara (corresponding to signal B), aTc (signal S) and IPTG (signal A). The outputs are RFP (signal b) and GFP (signal a). Genes araC, tetR, LuxI, and LacI are always on. Gene RFP will produce RFP when both P tet and P BAD promoters are activated. As tetR is always produced, it inhibits the P tet promoter. But the presence of aTc can release this repression. This means that the presence of aTc (signal S) will activate promoter P tet . On the other hand, AraC is also always produced and it inhibits in turn the P BAD promoter of the RFP gene. Still, the presence of Ara (signal B) releases the repression and thus the P BAD promoter will be activated. In other words, the presence of both Ara (signal B) and aTc (signal S) allows the expression of RFP (signal b).

At the same time, tetR represses the P tet promoter of the λcI gene. Having aTc on releases this repression. Thus, the production of λcI will inhibit the expression of the LuxR gene. So, the absence of aTc will allow the expression of LuxR. As LuxI is always on and produces the AHL enzyme, LuxR, if activated by the absence of aTC, will bind AHL and activate the P lux promoter of the GFP gene. LacI represses the P T ac promoter of the GFP gene, but IPTG can release this repression, and activate the P T ac promoter. This means that the presence of IPTG (signal A) and the absence of aTc (signal S) activate both promoters, and let gene GFP to be expressed (i.e. activate signal a). Note that we don't consider uptake or decay kinetics for input signals.

Kinetics

The activity of promoters P BAD , P tet , P lux , and P T ac can be expressed as (see Supplementary Information of [START_REF] Seok Moon | Genetic programs constructed from layered logic gates in single cells accessed terms of use detailed terms[END_REF]; note that here we did not take into account the RNA polymerase reactions):

P BAD = K bad 1 f bad T L 1 + K bad 1 f bad T L + K bad 2 f bad T P tet = 1 1 + 2K tet 1 f tet T + (K tet 1 ) 2 (f tet T ) 2 P lux = K lux 1 f lux T L 1 + K lux 1 f lux T L P T ac = 1 1 + K tac 1 f tac T (6.1) (6.2) (6.3) (6.4)
where

f k T L = (L k ) n k (K k D ) n k + (L k ) n k , (6.5)
and

f k T = 1 -f k T L , (6.6)
with k ∈ {bad, tet, lux, tac}. We have chosen the parameters for equations (6.1)- (6.6) as in Table 6.1.

The value of K tet D given in [START_REF] Seok Moon | Genetic programs constructed from layered logic gates in single cells accessed terms of use detailed terms[END_REF] is 3.6 ng/ml. The molar mass of TetR is 23 328 g/mol. Using this information, we write the value of K tet D in M.

Note that for each promoter, P BAD , P T ac , P lux , P tet and P λcI , we consider the states depicted in Figure 6.3. Thus, P λcI activity is given by

P λcI = 1 1 + K λ 1 D + K λ 2 D + αK λ 1 K λ 2 D 2 , (6.7
K λ 2 • 10 -8 M α 132 K λ 1 6.9 • 10 8 M -1 K λ 2 2.7 • 10 7 M -1 D [(λcI) 2 ]
Table 6.2: P λcI parameters [START_REF] Golovkin | Statistical fluctuations of the level of operator filling by repressor determine the level of noise of reporter gene expression[END_REF].

and it depends on the dimers concentration D. The reaction describing the formation of dimers (λcI) 2 is: (λcI) 2 K λ ----2λcI.

We can then find the dimers' concentration, D =[(λcI) 2 ], by solving the system:

     [λcI] total = [λcI] + 2 [(λcI) 2 ] K λ = [λcI] 2 [(λcI) 2 ] . (6.8)
Parameter values for equations (6.7)-(6.8) are given in Table 6.2.

Using the law of mass action from chemical kinetics, we can write the following set of differential equations describing the chemical species of the proposed model: where the dissociation constant k d can generally be written as

d[RF P ] dt = k rf p max P tet (t -τ rf p )P BAD (t -τ rf p ) -k rf p d [RF P ] d [λcI] total dt = k λ max P tet (t -τ λ ) -k λ d [λcI] total d[LuxR] dt = k lxr max P λcI (t -τ lxr ) -k lxr d [LuxR] -k cmp a [LuxR][AHL] + k cmp d [LuxR : AHL] d [LuxI] dt = k lux syn -k lux d [LuxI] d[AHL] dt = k ahl max [LuxI] -k ahl d [AHL] -k cmp a [LuxR][AHL] + k cmp d [LuxR : AHL] d[LuxR : AHL] dt = k cmp a [LuxR][AHL] -k cmp
k j d = ln 2 t j 1/2
, where j ∈ {rf p, λ, lxr, lux, ahl, cmp, gf p}.

Note that even if we assumed that LuxI is always on, we also consider the differential equation (6.12) describing the time evolution of the concentration of LuxI, in order to make sure that this concentration is reasonable, given the volume of an E. coli bacterium, which is about 1. Note also that both transcription and translation require a significant amount of time to form the completed product, hence the necessity to include a delay term τ j in the modeling. This delay term is simulated by considering the past ligand concentration L j corresponding to equations (6.9)-(6.12) and (6.15).

The half-life constant t j 1/2 or the degradation rate k j d , the maximum protein synthesis rate k j max , the synthesis rate k j syn , the association constant k j a , and the delay τ j for equations (6.9)-(6.15) are given in Table 6.3.

Choice of parameters

Estimating gene expression times mainly relies on the length of the genes and the average transcription/translation speeds. Using the values given in [START_REF] Kennell | Transcription and translation initiation frequencies of the Escherichia coli lac operon[END_REF], we compute the transcription and translation delays for λ repressor and LuxR. Moreover, we estimate the RFP and GFP transcriptional and translational delays at about 2 min.

The values chosen for maximum protein synthesis rates and for half-life time or degradation rates are justified by some similar values found in Supplementary Information of [START_REF] Pai | Optimal tuning of bacterial sensing potential[END_REF] and [START_REF] Lynch | The bioenergetic costs of a gene[END_REF]. For k cmp d we have chosen the value estimated in [START_REF] Weber | Dynamics of the quorum sensing switch: stochastic and non-stationary effects[END_REF], as for k cmp a we have chosen a value of the same order of values found in [START_REF] Takayama | Switch of spnr function from activating to inhibiting quorum sensing by its exogenous addition[END_REF][START_REF] Urbanowski | Reversible acyl-homoserine lactone binding to purified vibrio fischeri luxr protein[END_REF].

Parameter k lux syn has the same order as the value of LuxI protein synthesis rate given in [START_REF] Basu | A synthetic multicellular system for programmed pattern formation[END_REF]. We have initially designed a version of the model where the LuxI was induced by λcI, and LuxR was constitutive (always on). But because this version was two times slower than the case where LuxR is induced by λcI, and LuxI is constitutive, we have chosen the latter for our design.

The half-lives of GFP and RFP have an influence on the amplitude (the steady-state) of the concentrations of the two proteins. More precisely, increasing the half-life will also τ λ λcI transcription and translation delay 0.55 min [START_REF] Kennell | Transcription and translation initiation frequencies of the Escherichia coli lac operon[END_REF] τ lxr LuxR transcription and translation delay 0.59 min [START_REF] Kennell | Transcription and translation initiation frequencies of the Escherichia coli lac operon[END_REF] τ gf p GFP transcription and translation delay 2 min [START_REF] Kennell | Transcription and translation initiation frequencies of the Escherichia coli lac operon[END_REF] Table 6.3: Parameters values for equations (6.9)-(6.15).

increase the amplitude (and the fall and rise time 1 as depicted in Figures 6.7 and 6.8).

Additionally, we have chosen GFP and RFP half-lifes such that when the output signals are on, their amplitudes are similar. The optimal values would be 15.5 and 26 minutes, respectively, which is reasonable considering references [START_REF] Pai | Optimal tuning of bacterial sensing potential[END_REF][START_REF] Lynch | The bioenergetic costs of a gene[END_REF].

Another important parameter is k lxr d , because it plays a crucial role in GFP fall time. Even if in [START_REF] Goryachev | Systems analysis of a quorum sensing network: Design constraints imposed by the functional requirements, network topology and kinetic constants[END_REF] we can find the value 0.0231 min -1 , we have chosen k lxr d = 0.15 min -1 in order to decrease the GFP fall time and to improve the circuit design. This chosen value depends on the half-life of LuxR and can be engineered. A more detailed explanation of this choice can be found in Section 6.6.

1 By fall time we understand the time required for a signal to fall from 90% to 10% of its final value. By rise time we understand the time required for a signal to rise from 10% to 90% of its final value.

Steady-states study

We will generally denote steady-state values by an asterisk. For instance, [LuxI] * is the steady-state value of [LuxI]. In a similar way, P * lux denotes the promoter activity when the corresponding ligand, [LuxR : AHL] in this case, has its steady-state value.

Moreover, P * BAD corresponds to a constant [Ara] input concentration. Given these notations, the following proposition holds. Proposition 6.1. For any set of constant inputs ([IP T G], [Ara], [aT c]), the gene network described by equations (6.9)-( 6.15) has a unique steady-state given by We can also write, from equation (6.23), the equality (6.16), where P * tet and P * BAD can be expressed as a function of [aT c] and [Ara], respectively, in view of Table 6.1 and equations (6.1)-(6.2) and (6.5)-(6.6). In the same manner, the concentration [λcI] total can be written using equation (6.24). Thus, we obtain equation (6.17 We have implemented the delay differential equations with Hill functions (6.9)- (6.15) in Matlab/Simulink. For the simulation, we have considered constant inputs. Then, for different sets of constant input signals, we have compared the analytic steady-state values given by equations (6.16)-(6.22) with the steady-states given by Simulink, and we have obtained the same results. The steady-state values are given in Table 6. It is easy to see that the following conditions hold:

                                                               [RF P ] * = k rf p
d[λcI] total dt > 0 if [λcI] total < [λcI] * total , d[λcI] total dt = 0 if [λcI] total = [λcI] * total , and
d[λcI] total dt < 0 if [λcI] total > [λcI] * total ,
with the steady-state [λcI] * total defined as in (6.17). The phase-space flow on the [λcI] total axis therefore is depicted in Figure 6 We consider equation (6.9 is the equation of a parabola that passes through the origin, a strictly increasing function of [λcI]. The two must therefore intersect at a unique point in the positive quadrant, as depicted in Figure 6.5. Since [λcI] total → [λcI] * total asymptotically, and since system (6.8) has a unique solution for every value of [λcI] total , the solution to system (6.8) also approaches a unique asymptotic value. Specifically, let D * be the asymptotic concentration of dimers corresponding to [λcI] total = [λcI] * total . Thus, P λcI also approaches an asymptotic value P * λcI . We saw earlier that [LuxI] also tends to an asymptotic value [LuxI] * . Thus, equations (6.11), (6.13), and (6.14) asymptotically reduce to

d[LuxR] dt = k lxr max P * λcI -k lxr d [LuxR] -k cmp a [LuxR][AHL] + k cmp d [LuxR : AHL], d[AHL] dt = k ahl max [LuxI] * -k ahl d [AHL] -k cmp a [LuxR][AHL] + k cmp d [LuxR : AHL], d[LuxR : AHL] dt = k cmp a [LuxR][AHL] -k cmp d [LuxR : AHL].
We consider the state vector [LuxR AHL (LuxR : AHL)] T , and compute the Jacobian matrix

A =      -k lxr d -k cmp a [AHL] * -k cmp a [LuxR] * k cmp d -k cmp a [AHL] * -k ahl d -k cmp a [LuxR] * k cmp d k cmp a [AHL] * k cmp a [LuxR] * -k cmp d     
.

In view of steady-state equations (6.18) and (6.17), we note that the matrix A depends We note that these transition times are comparable, which gives an appropriate behaviour to our multiplexer design. For instance, we consider the time evolution of [GF P ] and [RF P ] concentrations when [IP T G] and [Ara] input signals are both on. The case when the select signal [aT c] changes from off to on is depicted in Figure 6.9, and the case when it switches from on to off is illustrated in Figure 6.10. We remark that, given that the rise and fall times of GF P and RF P are similar, the transition takes place so that there is only one output signal activated at a given time (either GF P , or RF P ).

Parameter scan

The half-lives are experimentally manipulable, which means that the degradation rates are manipulable [START_REF] Michael | A synthetic oscillatory network of transcriptional regulators[END_REF]. We can shorten the half-lives (i.e. increase the degradation rates)

to study how the proposed circuit behaviour changes. Moreover, we can also decrease the (maximum) protein synthesis rates and remark how the concentrations [GF P ] and [RF P ] vary.

In the sequel (see Figures 6.11 Thus, by saying that a molar concentration stays well below the threshold for physical relevance, we consider the order of magnitude of the concentration with respect to a typical E. coli cell volume. The interpretation of Figures The red dot corresponds to the value of the parameter considered in Table 6.3. The red dot corresponds to the value of the parameter considered in

Discussion on performance

We have mentioned in Section 6.4 that the k lxr d degradation rate plays a crucial role in the GFP fall time. Indeed, if we compare the performance of the proposed design with all the parameters chosen as in Table 6.3 to the same design where we choose k lxr d = 0.0231 min -1 (as in [START_REF] Goryachev | Systems analysis of a quorum sensing network: Design constraints imposed by the functional requirements, network topology and kinetic constants[END_REF]), t gf p 1/2 = 15.5 min and t rf p 1/2 = 26 min, we find that the latter design has a slower transient regime. We note that the half-lives of GFP and RFP affect their steady-state concentrations. Thus, the output steady-states are of same We remark that increasing the value of k lxr d decreases GFP fall time, but slightly increases the GFP rise time. Overall, we obtain better switching times.

We might even further improve the performances and decrease the GFP rise time by taking into account an eighth differential equation describing [LacI] evolution in time.

The differential equation corresponding to this concentration would depend on LacI degradation rate.

Concluding remarks

New logic circuits are important building blocks in synthetic biology. As very little analysis has been carried out for multiplexer genetic networks in the literature, we propose a multiplexer model where two outputs are controlled by three input signals.

Our design is described by a set of differential equations with delays. All the parameters involved in the model have been chosen according to experimental values found in the literature. We have seven variables in our design and have proven the uniqueness and stability of their steady states.

We show, as the simulations point out, that the temporal characteristics of the device are parameter-(in particular, delay-)dependent, even though the stability of the unique steady-state is unconditional. In effect, the device has three distinct behaviors for reasonable parameters values, even though they don't arise from a bifurcation. The delays considered in this model play a role only in the transient regime. Moreover, we have done a parameters scan and conclude that the circuit behaviour does not change when its parameters vary in a quite large range, which point out that the proposed design is relatively robust.

Chapter 7

Extensions and perspectives

Consensus problem in networks

The study of the proposed consensus algorithm (3.7) presented in Chapter 3 is only the tip of the iceberg. Further research should be carried out in connection with the optimal choice of τ and α parameters. For instance, an adaptive algorithm can be derived for a global optimum choice of parameters, based on the difference between the estimation and the real value of the neighbours' states. Another issue worthy to study is what happens if there is no normalization by d i in equation (3.7). In this case, an equivalent normalization for the network should be found.

For some applications, it would be interesting to derive the discrete case of this algorithm, and to compare it to the continuous case (3.7). Moreover, extensions of the proposed model can be made for other cooperative behaviours, like synchronization.

Another problem that has not been studied in this thesis is the consensus value. A natural question would be whether the consensus values given by the proposed model always coincide to the value obtained in the case of classical consensus problem, or other protocols in the literature. In this work, we have considered constant initial conditions. One question would be what happens in the case of sinusoidal or random initial conditions: will the consensus value change?

There are different issues that can be studied with respect to the robustness of the proposed algorithm. Are there networks where a small perturbation on τ = τ * parameter drastically changes the network performance or bring instability? How does our algorithm react when one agent in the network suddenly changes its opinion in a very unpredictable way? Furthermore, the case where the network has an uncertainty on its topology, or even the case on switching topology networks should be studied under the protocol (3.3).

Other perspective studies can focus on the analysis of the proposed protocol in some particular networks, like bipartite networks. In the sequel however, we present some observations on other types of network topology, namely circulant networks.

Some observations on directed circulant networks' eigenvalues as prototypes for worst delay tolerance

In this subsection, we first consider a network configuration that has all its eigenvalues on the circle, namely the directed cycle shown in Figure 7.1(a). We refer to this configuration as Network A. Indeed, it is easy to calculate that the characteristic polynomial for the Laplacian for Network A is

(1 -λ k ) n = 1;
thus the eigenvalues λ k are roots of unity shifted by one:

λ k = 1 -exp(2πi(k -1)/n), k = 1, . . . , n; see Figure 7.1(e).
It follows that for very small network sizes n, the eigenvalues are outside the red region of Figure 3.6 (although they are still all on the circle); so Network A will have a delay margin equal to 1. As the network size n becomes larger, the eigenvalues move into the red region and the delay margin decreases, as shown in Figure 7.2.

We next consider circulant networks with in-degree equal to two. Here we have two natural variations as extensions of Network A: either all connections follow the same direction, or connections of distance one and two have opposite directions. We refer to these networks as Network B and Network C, respectively; see For directed circulant networks with three connections per node, there are more possibilities for the choice of connection directions. We only show one possibility in Figure 7. 1(d), where all connections point in the same direction. We refer to this connection topology as Network D.

The Laplacian eigenvalues of network types B, C and D no longer all fall on the shifted unit circle (3.11), as those of Network A did. Nevertheless, the locus of eigenvalues still comes close to the circle near the origin, as shown in Figure 7.1(f)-(h). Hence, we have a similar conclusion: up to a certain network size, these networks have a delay margin Interestingly, Network C has a much higher delay tolerance than B, although all nodes in both networks have the same in-degrees.

We also test the accuracy of the curves in Figure 7.2 by direct simulation of the system (3.7). We pick τ = 0.9 and choose network sizes close to the curves in Figure 7.

Insights into geometric approach: extension to two parameter-dependent systems

The results presented in Chapter 5 apply not only to time-delay systems with two delays, but to a rather more general class of systems with two parameters, as discussed in Subsection 2.2.2. In the sequel, we give a few examples of various nature (linear, parameter-dependent polynomials, or non-linear, distributed delays systems) to illustrate the application of the theory synthesized in Theorem 5.8 and Corollary 5.9. These extensions also contain degenerate cases, that are identified and discussed.

Parameter-dependent polynomial

Polynomial characteristic equations are the most commonly seen in the undergraduate

textbooks. An algebraic study on the stability of parameter-dependent polynomial can be found in [START_REF] Mailybaev | On stability of polynomials depending on parameters[END_REF]. A representative study on robust stability of parameter dependent polynomials can be found in [START_REF] Barmish | New Tools for Robusteness of Linear Systems[END_REF]. A Puiseux series approach has been used to study the perturbation of the multiple roots under small parameter deviation in [START_REF] Chen | An eigenvalue perturbation approach to stability analysis, part i: Eigenvalue series of matrix operators[END_REF] and [START_REF] Chen | An eigenvalue perturbation approach to stability analysis, part ii: When will zeros of time-delay systems cross imaginary axis[END_REF].

In the following example, we will use the method we have arrived at to analyze such a system. According to Corollary 5.9, this means that both roots at i moves to the left half-plane as (p 1 , p 2 ) moves into S-sector. Furthermore, we may compute D = 3 > 0.

Distributed delays

Distributed delays also appear in many practical systems. An early example is given by Cushing [START_REF] Cushing | Time delays in single species growth models[END_REF] to model population dynamics as follows

where w(θ) is the kernel function. In chapter 2 of [START_REF] Insperger | Semi-discretization method for delayed systems[END_REF], the special case of w(θ) ≡ 1 was studied. In this case, the characteristic function becomes q(s) = s -a -b 1 -e -sσ s , s = 0.

The next example considers the case with two such distributed delays.

Example 7.2. Consider the "Cushing-like" system with the following characteristic quasi-polynomial: We note that in a neighbourhood of (p 10 , p 20 ) we can easily confirm that p 1 = p 2 for p 1 ≤ p 10 results in two imaginary roots in the neighbourhood of s 0 . This means that the positive and negative local stability crossing curves, T + (ω 0 ,p 10 ,p 20 ) and T - (ω 0 ,p 10 ,p 20 ) , coincide and the S-sector is empty. This situation is depicted in Figure 7.7. However, the conclusion about the G-sector still holds, i.e. there is a characteristic root in the right half-plane in the neighbourhood of s 0 , and another one in the neighbourhood of -s 0 , as (p 1 , p 2 ) moves to the G-sector. The aim of synthetic biology is to design and build new biological systems that respond to external signals by producing a desired chemical compound. More precisely, as specified in Section 2.3, the parts (bioparts) encoding biological functions have to be characterized, so that devices encoding human defined functions can be built. In our case, such a device implements a multiplexer. This circuit has to be implemented and tested in order to improve the design, as described in Figure 7.8. Of course, in the synthetic biology design cycle, a precise control of biological pathways is critical for a successful design. In our model (as presented in Chapter 6), we have basically considered a proportional (P) controller for the device: the input signals can be either on or off, and we neglect their uptake and decay kinetics (which is reasonable for a first implementation version).

q(s, τ 1 , τ 2 ) = s -a -b 1 -e -sτ 1 s -c 1 -e -sτ 2 s , (7.3
Our model specifications (basically the first step in the cycle) can be summarized under the two truth tables that can be found in Section 6.2. The design phase in Figure 7.8 can be depicted as in Figure 6.2 and the corresponding comments in Section 6.3. The device modeling can be described by the ordinary differential equations (6.9)-(6.15). To sum up, the work presented throughout Chapter 6 only reduces to three steps of the synthetic biology design cycle: specifications statement, design and modeling. Before the implementation step, it is important to rigorously study and analyze the model. We have completed the stability analysis, we have studied the transient regime, and done a parameter scan. However, the study of the stochastic model can also be carried out to complete the modelling phase. Such a study can also be motivated by the following.

The reaction-rate equations given in Section 6.4 are continuous and deterministic. However, the time evolution of the chemical network is not a continuous process, as the number of molecules can only change by integer amounts dictated by the stoichiometry of the elementary reactions. Moreover, chemical reactions are stochastic rather than deterministic processes [START_REF] Daniel | Exact stochastic simulation of coupled chemical reactions[END_REF]. Treating transcription and translation in detail for each gene would require a very large and unwieldy model, although it can be done [START_REF] Roussel | Stochastic kinetics description of a simple transcription model[END_REF][START_REF] Mäkelä | Stochastic sequence-level model of coupled transcription and translation in prokaryotes[END_REF].

As in differential equation models, the use of delays allows a considerable simplification of stochastic models of gene expression [START_REF] Bratsun | Delay-induced stochastic oscillations in gene regulation[END_REF][START_REF] Roussel | Validation of an algorithm for delay stochastic simulation of transcription and translation in prokaryotic gene expression[END_REF][START_REF] Ribeiro | Stochastic and delayed stochastic models of gene expression and regulation[END_REF]. Delay-stochastic simulation algorithms [START_REF] Bratsun | Delay-induced stochastic oscillations in gene regulation[END_REF][START_REF] Roussel | Validation of an algorithm for delay stochastic simulation of transcription and translation in prokaryotic gene expression[END_REF][START_REF] Gibson | Efficient exact stochastic simulation of chemical systems with many species and many channels[END_REF][START_REF] Vo | Efficient rejection-based simulation of biochemical reactions with stochastic noise and delays[END_REF] and software [START_REF] Lloyd-Price | Sgns2: A compartmental stochastic chemical kinetics simulator for dynamic cell populations[END_REF][START_REF] Moinat | Extending the stochastic simulation software package StochPy with stochastic delays, cell growth and cell division[END_REF] are available.

We recall here the continuous model, and we remove the differential equation corresponding to LuxI, and adjust the value of k ahl max accordingly. The stochastic simulation we discuss in the sequel are based on Markov-chain and Monte Carlo simulations. We note that if the DDE model (7.6)-(7.11) works with concentrations, the species in the case of stochastic simulations are expressed in numbers of Converting the protein synthesis terms from the deterministic model into chemical reactions, on the other hand, represents a complex process. Thus, we need to carefully make a set of assumptions that can simplify this process without compromising the kinetics of our model. The first assumption we make concerns the Hill coefficients: for the stochastic simulations we suppose that the Hill coefficients are integers. Indeed, a parameter scan we carried out on the deterministic model shows that the proposed circuit has the same qualitative behaviour when Hill coefficients are rounded to the closest integer. The second assumption we make is that the promoters only have two states: on and off. Thus, we consider that there is only one promoter controlling GFP expression (basically we take P lux and P T ac promoters as a single promoter), and this promoter can switch from on to off and the other way around, giving the process

G gfp off k gfp on ------ k gfp off G gfp on , (7.12)
where G gfp on/off represents the two states of the GFP gene promoter. When the gene is on, it produces proteins according to the delayed mass-action reaction [START_REF] Marc | The use of delay differential equations in chemical kinetics[END_REF] G gfp on (t) k gfp max ---→ G gfp on (t) + GFP(t + τ gfp ). (7.13) We note that considering the same value for the delay τ gfp in equation (7.13) as in the deterministic model represents another assumption we made, in order to keep consistency with equation (7.11). In other words, we consider the transcription and translation as a single process by assuming that the GFP gene is expressed with some mean delay τ gfp .

By using this assumption, we avoid the proliferation of reactions and parameters in the model, and expanding it too much. Reactions (7.12) and (7.13) together give the same kinetics for the synthesis of GFP as the first term on the right-hand side of equation (7.11) provided the rate constants conform to the equilibrium distribution between the on and off states of the gene, viz.

k
where k sw is a characteristic rate of switching, a new parameter added in the stochastic model when converting the first term of equation (7.11) to reactions (7.12) and (7.13).

In order to avoid an unnecessary increase in the number of parameters in the model, we also assume that all the genes have the same rate of switching k sw . When we carry out the stochastic simulations, we choose an arbitrary value for k sw such that its value is larger than the largest of the other stochastic rate constants. We remind the reader that in principle, each gene could have its own characteristic switching rate. However, to avoid a proliferation of parameters, we take the same k sw for all genes, keeping in mind that we only need this parameter to be sufficiently large in order to reproduce the kinetics of (7.8).

We also need consider decay of LuxR:

LuxR

k lxr d --→ ∅.
The remaining terms in (7.8) correspond to the binding equilibrium LuxR + AHL

k LA +c ------ k LA -c
LuxR • AHL.

The first term of equation (7.9) corresponds to AHL synthesis, which is assumed to proceed at constant rate. There is also a decay term:

∅ k ahl max ---→ AHL k ahl d --→ ∅.
For the case of λcI expression, the first term of equation ( 7 After writing all the equivalent chemical reactions, we also need to convert the rate constants into stochastic rate constants [START_REF] Daniel | Exact stochastic simulation of coupled chemical reactions[END_REF]. This can be done by using some formulas based on the conversion from concentration to number of molecules, which depends on the type of the rate constants. We can have zero-, first-, or second-order rate constants. Sketch of the proof. We consider z 0 ∈ Z(f ) a zero of f . Since f is non-null on Ω, there is m ≥ 0 such that f (z) = (z -z 0 ) m φ(z), where φ ∈ H(Ω) and φ(z 0 ) = 0. A trivial computation shows that f f (z) = m z-z 0 + φ φ (z), where φ φ is holomorphic on Ω.

( f g ) = f (z 0 ) g (z 0 ) .
The next proposition, known as the argument principle (see [START_REF] Polya | Problems and Theorems in Analysis[END_REF]), use the winding number of a closed, continuous, oriented curve in the plane around a given point which defines the total number of times that the curve turns in the given direction (and returns to its starting point) around the point. This graphical interpretation is equivalent with an analytical representation that can be written:

Z -P = 1 2iπ γ f f ,
where Z is the number of zeros and P is the number of poles of the function f inside the contour γ.

Proof. We denote by Z γ (f ) and P γ (f ) the sets of zeros and poles of f inside the contour γ. Proposition A.5 shows that for z ∈ Z γ (f ), the residue in z of f f is the order of the zero z. A complex number p belongs to P γ (f ) if and only if p is a zero of 1 f , thus p is a pole of order Res The formulation of the next theorem can also we found in [START_REF] Polya | Problems and Theorems in Analysis[END_REF]. Proof. We assume that f (z) and g(z) are regular on γ since f (z) and f (z) + g(z) are different from 0 and |f (z)|> |g(z)| sufficiently close to γ. Moreover, the last condition is equivalent to |g(z)| |f (z)| < 1, so we can say that the function 1 + g(z) f (z) has (on γ) a positive real part. That means that, as z moves on γ, the change of its argument is equal to 0. Furthermore f (z) + g(z) can be also written as f (z) 1 + g(z) f (z) . Given that the winding number of the image under the product w = f (z) 1 + g (z) f (z) of γ is equal to the sum of the winding numbers of the images of γ under w = f (z) and w = 1 + g (z) f (z) respectively, i.e. the image of γ under f (z)+g(z) has the same winding number as the one under f (z)

B.2 QPmR

QPmR (Quasi-Polynomial Mapping Based Rootfinder) is a Matlab file containing several functions for computation and analysis of the spectra of characteristic quasi-polynomials of both retarded and neutral time delay systems. It is easy to use (and it has quite a good degree of stability, even if it is based on a recursive algorithm) and is well-documented.

The main input is a quasi-polynomial represented either by two arrays, or by a function handle; the routine finds the zeros of the quasi-polynomial (in some specified region, with some specified precision), using the algorithms that can been found in the referenced bibliography (see http://www.cak.fs.cvut.cz/algorithms/qpmr).

B.3 TRACE-DDE

An alternative for plotting the spectral values is to use the TRACE-DDE package. It has a friendly graphical user interface, computes the characteristic roots, and provides stability charts (sets of asymptotically stable/unstable regions in the parameters plane) for the robust stability analysis. It is intuitive and its documentation can be found at http://sole.dimi.uniud.it/ ~dimitri.breda/research/software/.

B.4 SOSTOOLS

It is noteworthy to also mention this Matlab toolbox, even if its aim is rather different compared to DDE-Biftool, QPmR and TRACE-DDE. SOSTOOLS [START_REF] Papachristodoulou | SOSTOOLS: Sum of squares optimization toolbox for MATLAB[END_REF] is a free software designed for sum of squares optimization problems. It can determine whether a polynomial is the sum of squares (of other polynomials), it can be used to search Lyapunov functions, to find bounds for minimization problems, to decide upon the copositivity of a given matrix, to find an upper bound for structured singular values in robustness problems and many other problems.

3. 1

 1 Estimating the future state x j (t + δ) of an agent j from its present and past states. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 The stability region of the function ψ given in (3.19) in the parameter space. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 The locus of numbers β ∈ C satisfying condition (3.18) is the open set enclosed by the blue curve. . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Map of the dominant characteristic root of eigenmodes corresponding to a generic Laplacian eigenvalue λ and delay value τ for the case α = 1. . . 3.5 Fraction of complex numbers λ that make the characteristic equation (3.14) unstable, from among 10 6 randomly selected complex numbers inside the shifted unit circle (3.34). . . . . . . . . . . . . . . . . . . . . . . . 3.6 Values of random complex numbers λ inside the shifted unit circle, colored according to the stability of the characteristic equation (3.14). . . . . . . . 3.7 Standard deviation function of time in the case of classical consensus problem (dashed line) and algorithm (3.7) with α = 1 and τ = 0.8 (continuous line) for a directed network with 50 nodes. . . . . . . . . . . . . . . . . . . 3.8 Right-hand side of (3.37) as a function of α. . . . . . . . . . . . . . . . . . 3.9 Standard deviation function of time in the case of classical consensus problem (dashed line) and algorithm (3.12) with α = 1 and τ = τ * (continuous line) for Example 3.1. . . . . . . . . . . . . . . . . . . . . . . 3.10 Time evoluation of network's agents in the case of classical consensus problem (top) and algorithm (3.12) with α = 1 and τ = τ * (bottom) for Example 3.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Figures vii 3.11 Standard deviation function of time in the case of classical consensus problem (dashed line) and algorithm (3.12) with α = 1 and τ = τ * (continuous line) for Example 3.2. . . . . . . . . . . . . . . . . . . . . . . 3.12 Evolution of real part of dominant characteristic roots as a function of τ for α = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Real part of the dominant root function of τ , corresponding to eigenvalues 1.0740 ± 0.3691i, 1.3590, 1.1899 ± 0.1679i, 1.0576, and 1.0278 ± 0.2463i. . .

  l'axe des imaginaires pour les systèmes qui comportent deux retards. Cette étude résout partiellement quelques problèmes connexes aux systèmes à retard qui étaient considérés des questions ouverts depuis 2005. Nous étudions le comportement des racines caractéristiques imaginaires multiples (multiplicité deux, trois et quatre) lors d'une petite perturbation des paramètres. Ces paramètres peuvent être soit deux retards, soir un retard et un autre paramètre, ou bien deux autres paramètres du système. Une contribution importante de ce travail est que les résultats peuvent être généralisés et utilisés sur une classe beaucoup plus grande de systèmes, avec ou sans retards, ou même avec des retards distribués. Les techniques utilisées se basent sur l'approche géométrique, alors que les autres résultats dans la littérature utilisent les séries de Puiseux. Une contribution qu'il faut citer est le calcul d'un critère pour décider si la paire de racines doubles et complexe conjuguées se dirige vers la stabilité ou vers l'instabilité quand les paramètres sont soumis à une perturbation. Autrement dit, nous trouvons un critère qui implique un nombre limité de calculs qui sont relativement faciles à faire (des dérivés partiels d'un quasi-polynôme par rapport aux paramètres). Nous vérifions notre méthode sur des exemples très variés et nous montrons comment la méthode peut s'étendre aux case des racines de multiplicité trois et quatre. Pour le cas de racines doubles, nous comparons la méthode géométrique avec une méthode algébrique. Nous vérifions que les deux méthodes donnent des résultats cohérents et nous discutons leurs points forts, leurs limitations et les questions ouvertes qu'elles posent. Tous ces résultats ont été publiés et présentés lors de trois conférences internationales [10-12] et une partie des résultats a été acceptée pour publication dans le journal Automatica. La dernière partie de cette thèse présente une application fortement interdisciplinaire des réseaux génétiques. Ce travail se trouve à l'intersection de l'ingénierie et de la biochimie. Nous modélisons un circuit génétique qui fonctionne comme un multiplexeur. Plus précisément, il s'agit d'un modèle de réseau génétique qui réagit à trois signaux en entrée en répondant soit par un signal unique parmi les deux signaux possibles à la sortie, soit en ne répondant rien. Ce modèle est conçu pour être intégré dans une cellule Résumé étendu en français xv minuscule (d'un volume de quelques femtolitres). Le but de la biologie synthétique, un sujet relativement récent de recherche, est de finalement connecter plusieurs cellules de ce type pour construire des réseaux plus grands, avec des fonctionnalités plus complexes. Dans ce contexte, le travail présenté dans le chapitre 6 a pour but juste de poser une brique. Néanmoins, le modèle proposé dépasse en fonctionnalité et en complexité d'autres modèles de multiplexeurs dans la littérature, où un seul type de signal est produit à la sortie. Nous étudions et analysons la stabilité du modèle en prenant en compte les retards correspondant au temps nécessaire à la transcription et à la traduction des gènes. Une contribution notable, par rapport à d'autres circuits de la littérature, est que nous prouvons l'unicité et la stabilité des points d'équilibre. Pour simuler et tester notre modèle, nous utilisons des valeurs de paramètres expérimentales trouvées dans la littérature. Nous effectuons également un test de robustesse, pour vérifier si le comportement du circuit change lorsqu'une petite perturbation se produit sur l'un des paramètres. Le modèle est décrit en utilisant des diagrammes génétiques standards. Le modèle déterministe comporte des équations différentielles avec des retards. Néanmoins, la méthodologie utilisée dans la biochimie nécessite une validation du modèle déterministe par le modèle stochastique, ce qui implique 23 réactions chimiques avec 15 espèces en interaction. Cette validation ouvre les perspectives d'une mise en place d'une expérience. Ce travail est en préparation pour une soumission à un journal interdisciplinaire et les résultats préliminaires ont déjà été présentés lors d'une session de posters à l'Université de Lethbridge. Dedicated to my dear father.

Figure 1 . 1 :

 11 Figure 1.1:Example of how a network can be depicted. This network can illustrate a transport system made up of airplanes. Their common task for instance is to bring people and food from one place to another within European Union countries.

Figure 1 . 2 :

 12 Figure 1.2:Example of cooperative behaviour in nature: in a school of fish, every fish (red fish denoting the focal fish) only knows where the neighbouring fish (depicted in blue) are heading, but have no idea towards which direction the whole school (all fish, including the ones represented in gray) is heading.

  work has been accepted for publication in Automatica journal.The last part of this thesis focuses on a challenging interdisciplinary application of a gene network, lying at the intersection of engineering and biochemistry. We design a biochemical circuit made of genes that functions like a multiplexer. More precisely, we propose a model of a gene network that reacts to three input signals by producing either only one of two possible output signals, or none. This model is designed to be integrated in a tiny cell (of a volume of a few femtolitres). The aim of a relatively recent field called synthetic biology is to eventually connect several such cells into a larger network with a more complex functionality. In this sense, the work presented throughout Chapter 6 is a small piece in the puzzle. Nevertheless, the proposed model exceeds in functionality and complexity other multiplexer models in the literature, where only one type of output signal is produced. We study and analyze the stability of the model, and we take into account the delays corresponding to the required time of gene transcription and translation. One notable contribution, when compared to other such circuits in the literature, is that we prove the uniqueness and the stability of the steady states. To simulate and test the model, we use experimental parameter values found in the literature. Even so, given the biochemical nature of the model, we also do a parameter scan (similar to, but less rigorous than a robustness study) to check whether the circuit behaviour changes when a small perturbation occurs on any of the relevant parameters. The model is described by using standard gene diagrams. The deterministic model involves differential equations with time-delays. Nonetheless, the biochemistry methodology requires a deterministic model validation by the stochastic model, which implies 23 chemical reactions with 15 interacting species. This validation opens the perspectives of an experimental bench setup. This work is in preparation for a submission to an interdisciplinary journal, and preliminary results have been already presented during a poster session at the University of Lethbridge.

loopsFigure 2 . 1 :Figure 2 . 2 :

 2122 Figure 2.1: Example of loops in a graph.

Figure 2 . 3 :

 23 Figure 2.3: Example of a simple graph.

Figure 2 . 5 :

 25 Figure 2.5: Example of trees.

Figure 2 . 5 .

 25 A spanning tree of a graph is a subset of the graph that is a tree and that has the same vertex set as the graph. A random graph is a graph in which properties such as the number of graph vertices, graph edges, and connections between them are determined in some random way. A complete graph is a simple graph where the vertices are pairwise adjacent. A few examples of complete graphs are given in Figure2.6.

Figure 2 . 6 :Figure 2 . 7 :

 2627 Figure 2.6: Examples of complete graphs.

Figure 2 . 8 :

 28 Figure 2.8: Examples of bipartite graphs.

Figure 2 . 9 :

 29 Figure 2.9: The central dogma of molecular biology.

Figure 3 . 1 :

 31 Figure 3.1: Estimating the future state x j (t + δ) of an agent j from its present and past states.

Definition 3 . 3 .

 33 The number s * ∈ C is called the dominant transverse root of the consensus algorithm (or dominant root, for short) if Ψ(s * ) = 0 and all roots s of Ψ satisfy Re(s) ≤ Re(s * ).

  7) reduces to the classical consensus problem (3.1)-(3.2), and the characteristic equation (3.15) reduces to Ψ(s) = n k=1 (s + λ k );(3.17)whose roots are directly given by the negative Laplacian eigenvalues. In particular, the dominant root equals the negative of the eigenvalue having the smallest real part, among the set {λ k } k≥2 . Note that Re(λ k ) = 0 if and only if λ k = 0, by(3.11). Hence, reaching consensus under the undelayed protocol (3.2) is possible if and only if the following condition is satisfied.

Figure 3 . 2 :Lemma 3 . 4 .

 3234 Figure 3.2: The stability region of the function ψ given in (3.19) in the parameter space.

2 :Figure 3 . 3 :

 233 Figure 3.3: The locus of numbers β ∈ C satisfying condition (3.18) is the open set enclosed by the blue curve.

. 23 )

 23 Proof. To begin with, we show that the characteristic equation (3.15) has a simple root at zero and all the remaining roots have negative real parts if and only if condition(3.23) holds. Consider first the roots of the first factor ψ 1 (s) = s -α + αe -sτ in(3.15). By a change of variable s = sτ , we can equivalently consider the roots of ψ 1 (s ) = s -τ α + τ αe -s . Lemma 3.4 gives that ψ 1 , and therefore ψ 1 , has a simple root at zero and all the other roots have negative real parts if and only if(3.23) holds. It then suffices to show that all roots of the remaining factors ψ k , k = 2, . . . , n, in (3.15) have negative real parts under condition(3.23

  .27) and the roots are on the open left complex plane. We next look for roots that may cross the imaginary axis as τ is increased from zero. Letting s = iω, ω ∈ R, the imaginary part of (3.24) gives ω = -b k sin(ωτ ). Therefore, |ω|= |b k sin(ωτ )|≤ |b k ωτ |. (3.28) Notice from (3.24) and (3.27) that ω = 0. Hence, dividing (3.28) by |ω|,

) where W 0

 0 is the principal branch of the Lambert W function and a k and b k are defined in (3.25)-(3.26). Consequently, the dominant transverse root of the consensus problem equals max 2≤k≤n 1 τ W 0 (τ b k e -a k τ ) + a k . (3.31) Proof. As seen in the proof of Theorem 3.5, the roots of ψ k satisfy equation (3.24),

3 and

 3 the fact that the eigenvalues of the normalized Laplacian matrix L are real numbers imply (3.31). Using Proposition 3.6 and equation (3.30), we calculate the dominant root for any eigenmode (3.13) of system (3.7) corresponding to a generic Laplacian eigenvalue λ. This gives a universal map of the dominant roots of the consensus problem (3.7) for any graph topology and delay value, which is depicted in Figure 3.4 for the case α = 1.

Figure 3 . 4 :

 34 Figure 3.4: Map of the dominant characteristic root of eigenmodes corresponding to a generic Laplacian eigenvalue λ and delay value τ for the case α = 1. The colour represents the real part of the dominant root of a factor (3.14) of the characteristic equation. The region above the black line corresponds to the case when there is a root with positive real part. Note that, since zero is always an eigenvalue of the Laplacian in every network and α = 1, the condition τ < 1 actually defines the region for reaching consensus, in accordance with Theorem 3.5.

  Figure 3.6 (a) through (d)). The value τ = 1 is the upper value of allowable delay for any network, since at this value the characteristic equation (3.14) has a double zero root (see Theorem 3.7 and its proof).

Figure 3 . 5 :

 35 Figure 3.5: Fraction of complex numbers λ that make the characteristic equation (3.14) unstable, from among 10 6 randomly selected complex numbers inside the shifted unit circle (3.34).

Figure 3 . 6 :

 36 Figure 3.6: Values of random complex numbers λ inside the shifted unit circle, colored according to the stability of the characteristic equation (3.14). Red points represent unstable values and blue points represent stable values of λ, shown for four different values of the delay.

. 35 ) 3 . 8 .

 3538 Theorem Let α ∈ R + and choose τ = τ * , as defined in equation(3.35). Then the characteristic equation (3.15) has a root s = -1. Moreover, all characteristic roots satisfy Re(s) ≤ -1 (we say that s = -1 is the "dominant root" and it corresponds to τ = τ * ) if and only if the connected network has all eigenvalues satisfying

  Erdös-Rényi random oriented graph generated with the probability p = 0.06. Cyclic oriented network.

Figure 3 . 7 :

 37 Figure 3.7: Standard deviation function of time in the case of classical consensus problem (dashed line) and algorithm (3.7) with α = 1 and τ = 0.8 (continuous line) for a directed network with 50 nodes.

  1+α α -s by e -ln( 1+α α )s in the last equation, and we change the variables in order to obtain s + ln 1+α α (1 -λ k )(1 + α)(-1 + e -s ) = 0. This is nothing else but an equation of the form (3.18), with β = ln 1+α α (1 -λ k )(1 + α). So we can use Lemma 3.4 and say that Re(s ) ≤ 0 when condition (3.20) hold for undirected networks, and when condition (3.21) holds for directed networks.

Corollary 3 . 9 .

 39 Let α ∈ R + and choose τ = τ * , as defined in equation (3.35). Then the characteristic equation (3.15) has a root s = -1. Moreover, if and only if the undirected network has all eigenvalues satisfying

Proof.

  Since in the case of an undirected network the normalized Laplacian L has all the eigenvalues real numbers between 0 and 2, with 0 = λ 1 < λ 2 ≤ λ 3 ≤ . . . ≤ λ n , we can conclude that for any given connected undirected network by choosing τ = τ * we get s = -1 the dominant root of the characteristic equation (3.15) if and only if condition (3.37) holds. Theorem 3.10. Let τ = τ * , as defined in equation (3.35). If zero is a simple eigenvalue of the normalized Laplacian L, then we can always find α > 0 satisfying (3.36) ∀k ≥ 2, such that the characteristic equation (3.15) has all characteristic roots satisfying

Figure 3 . 8 :

 38 Figure 3.8: Right-hand side of (3.37) as a function of α.

Figure 3 . 9 :

 39 Figure 3.9: Example 3.1. Standard deviation function of time in the case of classical consensus problem (dashed line) and algorithm (3.12) with α = 1 and τ = τ * (continuous line).

Example 3 . 2 .

 32 Another type of topology commonly used in real-life networks is circular topology. We consider now an undirected circular network made up of 20 nodes. We choose α = 10 and τ = τ * . Figure3.10 shows the evolution in time of the states x i , i = 1, . . . , 20 for the two cases: classical consensus problem and consensus problem with prediction. Note that the initial conditions are the same. Obviously, the consensus in reached faster (less than 6 units of time) if we use our model. Furthermore, Figure3.11 depicts the comparison of the standard deviation in the two cases. We add that if for the classical consensus problem λ 2 = 0.048943, using the algorithm (3.12) we obtain λ 2 = 1.

Figure 3 . 10 :Figure 3 . 11 :

 310311 Figure 3.10: Example 3.2. Time evoluation of network's agents in the case of classical consensus problem (top) and algorithm (3.12) with α = 1 and τ = τ * (bottom).

Figure 3 . 13 :

 313 Figure 3.13: Real part of the dominant root function of τ , corresponding to eigenvalues 1.0740 ± 0.3691i, 1.3590, 1.1899 ± 0.1679i, 1.0576, and 1.0278 ± 0.2463i.

Proposition 4 . 1 .

 41 Consider the regular quasi-polynomial function (4.2). The maximal multiplicity of zero singularity cannot be larger than 0 m = N (n 0 + 1) + n, where N is the number of non-zero distinct delays, n 0 is the upper degree of the polynomial family r i , with 1 ≤ i ≤ N and n is the degree of r 0 .

. 3 ) 4 . 3 . 1

 3431 Proof of Proposition 4.1

  of the vector V , and by deg(V k ) the degree of the polynomial V k (x). Using equation (4.5), we can show that for any i = l, deg(V i ) = deg(V l ). This means that the vectors in family F k are linear independents. Moreover, no element from F l (the family of vectors in R N (n-1) ([τ l ])) can be written as a linear combination of elements of F k with l = k. The direct computation of det(M 2 ) gives det(M 2 ) = 1≤k≤n-2

Figure 4 . 1 :Figure 4 . 2 :

 4142 Figure 4.1: Spectrum distribution for a vector disease model. b = c = 0.5 and τ = 2.

Figure 4 Figure 4 . 3 :

 443 Figure 4.3: Spectrum distribution for an inverted pendulum. τ 1 = 1.

Figure 4 . 4 .Figure 4 . 4 :

 4444 Figure 4.4.The mass of the pendulum is m and is small enough comparing to the mass of cart M . We denote by θ the angle between the vertical direction of the slope and the inverted pendulum, and we assume that it can be measured or observed. The cart moves on the slope due to the external force F and we can measure its position d, as the distance between a fixed reference point and its center of gravity. We also denote by L the length of the inverted pendulum and by J the moment of inertia with respect to its center of gravity.

Figure 4 . 5 :

 45 Figure 4.5: Inverted pendulum on a slope: spectrum distribution for τ 1 = 1.

115 Figure 4 . 6 :

 11546 Figure 4.6: Double inverted pendulum

γ

  ih above are the gains of the command u of the form

2 23 - 2593080 τ 1 τ 2 3 23 -63504 τ 1 τ 2 Figure 4 . 7 :

 23247 Figure 4.7: Equation (4.26) solution in τ 1 -τ 2 domain

Remark 4 . 9 .

 49 If we chose a value for τ 1 contained in the restricted domain, τ 1 = 0.1 for instance, we obtain γ 11 = 709.783, γ 12 = -245.367, γ 21 = -620.637, γ 22 = 216.379 and τ 2 = 0.115. If we replace these values in polynomials P 0 , P 1 and P 2 , the obtained coefficients verify

for 1 ≤Figure 4 . 8 :

 148 Figure 4.8: Double inverted pendulum: spectrum distribution for τ 1 = 0.1 and τ 2 = 0.115

Figure 5 . 1 :

 51 Figure 5.1: G-sector and S-sector.

∂ 2 p 1 ∂u 2 ∂ 2 p 2 ∂u 2 T

 22 given in(5.11) for γ = i and -i have the same value, T - (ω 0 ,p 10 ,p 20 ) and T + (ω 0 ,p 10 ,p 20 ) (A C and C B in Figure5.1) are tangent to each other at the point (p 10 , p 20 ), thus forming a cusp.

Figure 5 . 2 :

 52 Figure 5.2: The mapping (p 1 (s), p 2 (s)) with s -s 0 in the first quadrant.

2 :Figure 5 . 2 )

 252 Figure 5.2) relative to the point C . The mapping with s-s 0 in the other three quadrants are similar.

Theorem 5 . 8 (

 58 Migration of Double Roots). If (p 1 , p 2 ) is in the G-sector in a sufficiently small neighborhood of (p 10 , p 20 ), then one root of (5.1) in the neighborhood of s 0 is in the right half-plane, the other is in the left half-plane.When (p 1 , p 2 ) is in the S-sector, then the two roots are either both in the left half-plane or both in the right half-plane. More specifically, Case i. If D > 0, and T - (ω 0 ,p 10 ,p 20 ) is on the counterclockwise side of T + (ω 0 ,p 10 ,p 20 ) in the S-sector, then both roots are on the left half-plane. Case ii. If D > 0, and T - (ω 0 ,p 10 ,p 20 ) is on the clockwise side of T + (ω 0 ,p 10 ,p 20 ) in the S-sector, then both roots are on the right half-plane. Case iii. If D < 0, and T - (ω 0 ,p 10 ,p 20 ) is on the counterclockwise side of T + (ω 0 ,p 10 ,p 20 ) in the S-sector, then both roots are on the right half-plane. Case iv. If D < 0, and T - (ω 0 ,p 10 ,p 20 ) is on the clockwise side of T + (ω 0 ,p 10 ,p 20 ) in the S-sector, then both roots are on the left half-plane. Proof. Consider Case i. The situation is illustrated in Figure 5.3. Let the region bounded by the arc EP B and line segments BC and CE be denoted as I, and the region bounded by the curves E P B , B C and C E be denoted as I . Similarly, region II is bounded by BQF , F C, CB, and region II is bounded by B Q F , F C , C B ; region III is bounded by F RA, AC, CF , and III is bounded by F R A , A C , C F ; region IV is bounded by ASE, EC, CA, and region IV is bounded by A S E , E C , C A . As discussed before the theorem, (p 1 (s), p 2 (s)) is a bijection from I to I when s is restricted to I. Similarly, (p 1 (s), p 2 (s)) is a bijection from II to II when restricted to II, or from

Figure 5 . 3 :Figure 5 . 4 :Figure 5 . 5 :

 535455 Figure 5.3: The mapping (p 1 (s), p 2 (s)) in a neighborhood of s 0 . Case i: D > 0, and T - (ω0,p10,p20) is on the counterclockwise side of T + (ω0,p10,p20) in the S-sector.

Figure 5 . 6 :

 56 Figure 5.6: The mapping (p 1 (s), p 2 (s)) in a neighborhood of s 0 . Case iv: D < 0, and T - (ω0,p10,p20) is on the clockwise side of T + (ω0,p10,p20) in the S-sector.

20 γ=i , and ∂ 2 p i ∂u 2

 2 may be evaluated by(5.11) or (5.12) with γ = i. If κ > 0 (5.14) instead, then both roots are in the right half-plane.

1 ,

 1 ∆p 2 ) counterclockwise through an angle θ ∈ (0, π) as is shown in Figure5.7.

Figure 5 . 7 : 2 ,

 572 Figure 5.7: T - (ω0,p10,p20) is on the counterclockwise side of T + (ω0,p10,p20) , and the angle φ needed to rotate A B to the direction of ∂ 2 p1 ∂u 2 , ∂ 2 p2 ∂u 2

∂u 2 T+

 2 by rotating (∆p 1 , ∆p 2 ) clockwise through an angle θ ∈ (0, π), and the conclusions are true in view of Case ii and Case iv in Theorem 5.8. We have exhausted all possibilities, and the corollary is proven.

Example 5 . 1 .

 51 Consider a time-delay system with the following characteristic quasipolynomial q(s, τ ) = s 2 -2s + 2 + [(2 cos 1) s -2 (cos 1 + sin 1)] e -τ 1 s + e -τ 2 s .

Figure 5 . 8 :

 58 Figure 5.8: τ 1 -τ 2 parameter space for Example 5.1. Characteristic equation (5.18) has two roots (region A), four roots (region B), and no roots (region C) with positive real part.

Remark 5 . 10 .= a 11 a 12 a 21 a 22 , A 1 = b 11 b 12 b 21 b 22 , A 2 = c 11 c 12 c 21 c 22 .

 51012221122221222 For n = 2, consider A 0 If at least one of the conditions (C1): {b 11 = 0, b 12 = 0, c 11 = 0, c 12 = 0} (C2): {b 11 = 0, b 21 = 0, c 11 = 0, c 21 = 0} (C3): {b 12 = 0, b 22 = 0, c 12 = 0, c 22 = 0} (C4): {b 21 = 0, b 22 = 0, c 21 = 0, c 22 = 0}

Remark 5 . 11 (

 511 An equivalent eigenvalue problem). The eigenvalue problem (5.22) corresponding to the functional differential equation (5.[START_REF] Montoya | Ecological networks and their fragility[END_REF]) is nonlinear and finite-dimensional.

s 2 = 3 80 + i 80 sec 2 ( 1 - 11 + 2 ( 4 -( 3 -

 32111243 cos 1 -i sin 1) (cos 1 -(1 -i) sin 1)((1 + i) cos 1 + sin 1) [(1 + i)((115 + 82i) cos 1 + (87 + 14i) cos 3 + (61 -2i) cos 5 + (9 + 18i) cos 7 -(73 -194i) sin (53 + 70i) sin 3 -(25 -70i) sin 5 -(45 + 18i) sin 7)τ cos 1 -i sin 1)((42 -40i) -(21 -4i) cos 2 + (22 -40i) cos 4i) cos 6 + (4 + 28i) sin 2 + (4 + 4i) sin 4 -(4 -24i) sin 6)τ 21 ]. Now that we have concrete expressions of s 1 and s 2 as functions of τ 11 and τ 21 for Example 5.1, we proceed with giving values to τ 11 , τ 21 , and ε in order to illustrate Proposition 5.13. Suppose γ 2 = β 1 = β 2 = γ 1 = γ 2 = 1. We consider two simple cases, as follows, and compute the value of s(ε) for each case:

Figure 5 . 9 :

 59 Figure 5.9: Double characteristic root behaviour for Example 5.1. Yellow points (lying on the imaginary axis) correspond to τ 10 = 1 and τ 20 = 2. Case (a): fix τ 1 = 1and increase τ 2 from 2 to 3. Case (b): fix τ 1 = 1 and decrease τ 2 from 2 to 1.9. For both cases, the double root s = ±i splits up into two simple roots. = 0 corresponds to yellow points, and these points become red as the value of τ 2 increases or decreases with a step of 10 -3 .

Theorem 5 . 18 .

 518 Suppose system (2.8) satisfies (5.2) and (5.3) for m = 3, and assumption (5.4) holds. Then, as (p 1 p 2 ) moves from (p 10 p 20 ) to one of the two regions of N ε (p 10 , p 20 ) divided up by T (ω 0 ,p 10 ,p 20 ) , at least one root moves to the right half-plane, and one other root moves to the left half-plane. The remaining root may move to either the left half-plane, or the right half-plane. Specifically: Case i. D > 0 and (p 1 p 2 ) moves in the region on the clockwise side of T + (ω 0 ,p 10 ,p 20 )

and 5 .Figure 5 . 10 :

 5510 Figure 5.10: The mapping (p 1 (s), p 2 (s)) in a neighborhood of s 0 . D > 0.

Figure 5 . 11 :

 511 Figure 5.11: The mapping (p 1 (s), p 2 (s)) in a neighborhood of s 0 . D < 0.

Figure 5 . 12 :

 512 Figure 5.12: The mapping (p 1 (s), p 2 (s)) in a neighborhood of s 0 . Case i: D > 0, and T - (ω0,p10,p20) is on the counterclockwise side of T + (ω0,p10,p20) in the S-sector.

Figure 5 . 13 :

 513 Figure 5.13: The mapping (p 1 (s), p 2 (s)) in a neighborhood of s 0 . Case ii: D > 0, and T - (ω0,p10,p20) is on the clockwise side of T + (ω0,p10,p20) in the S-sector.

Figure 5 . 14 :

 514 Figure 5.14: The mapping (p 1 (s), p 2 (s)) in a neighborhood of s 0 . Case iii: D < 0, and T - (ω0,p10,p20) is on the counterclockwise side of T + (ω0,p10,p20) in the S-sector.

Figure 5 . 15 :

 515 Figure 5.15: The mapping (p 1 (s), p 2 (s)) in a neighborhood of s 0 . Case iv: D < 0, and T - (ω0,p10,p20) is on the clockwise side of T + (ω0,p10,p20) in the S-sector.

a 20 630+193s 1 +219s 3 -583s 5 +108s 6 -s 9 +2s 11 -172c 1 +216c 3 -10c 5 +18c 6 +2c 11 223s 1 -337s 5 +2s 11 +105c 1 +435c 5 a 11 -58s 1 +1836s 2 +271s 3 -

 361131113 200s 5 +223s 7 +336s 8 -2s 11 +2s 13 +466c 1 -380c 2 -1515c 3 +250c 5 +15c 7 +20c 8 +4c 11 223s 1 -337s 5 +2s 11 +105c 1 +435c 5 a 02 -1439+606s 1 -418s 2 -5s 4 +1626s 5 -30s 6 +47s 8 +30s 10 +24s 11 -528c 1 +996c 2 -225c 4 +1620c 5 -100c 6 +84c 8 -36c 10 -12c 11 446s 1 -674s 5 +4s 11 +210c 1 +870c 5 a 10 -194s 1 +684s 2 +1119s 3 -679s 5 +167s 7 -96s 8 -s 11 -2s 13 +37c 1 -1742c 2 +486c 3 +145c 5 -188c 7 +302c 8 -2c 11 +2c 13 223s 1 -337s 5 +2s 11 +105c 1 +435c 5 a 00 -1849+1116s 1 -490s 2 -117s 4 +2112s 5 +66s 6 -73s 8 -2s 10 -12s 11 -850c 1 +1396c 2 -301c 4 +2430c 5 +76c 6 -38c 8 +44c 10 +4c 11 446s 1 -674s 5 +4s 11 +210c 1 +870c 5 a 01 808+261s 1 +290s 2 +3s 4 -1815s 5 -92s 6 +65s 8 -44s 10 -6s 11 +691c 1 -628c 2 +95c 4 +225c 5 +50c 6 -67c 8 -18c 10 -16c 11 223s 1 -337s 5 +2s 11 +105c 1 +435c 5

Figure 5 . 16 :

 516 Figure 5.16: The positive local stability crossing curve T + (ω0,τ10,τ20) is the curve on the left-hand side of the point (τ 10 , τ 20 ) = (3, 5). The negative local stability crossing curve T - (ω0,τ10,τ20) is the curve on the right-hand side of the point (τ 10 , τ 20 ) = (3, 5). These two curves divide τ 1 -τ 2 space into two regions.

Figure 6 . 1 :

 61 Figure 6.1: Multiplexer circuit.

Figure 6 . 2 :

 62 Figure 6.2: Proposed model.

  ), where P * tet depends on [aT c]. Note that equation (6.28) gives k cmp a [LuxR][AHL] = k cmp d [LuxR : AHL]; using this equality in (6.25) and (6.27) and replacing the expression of [LuxI] concentration (6.19), we obtain equalities (6.18) and (6.20) where P * λcI depends on ([λcI] total concentration, which depends on) [aT c] concentration. Knowing the concentration (6.18) of [LuxR] and (6.20) of [AHL], we can now write [LuxR : AHL] concentration from equation (6.28) as in (6.21). Finally, from equation (6.29) we write the steady state (6.22) where P * T ac can be expressed as a function of [IP T G], and P * lux depends on [LuxR : AHL], which depends on [λcI] total , which depends on [aT c].

  4. The concentration of input and output signals are given in M, as well as the corresponding boolean values. This simulation confirms the fact that the steady-state does not depend on time-delays.

. 4 .

 4 The arrows indicate the direction of motion along the axis, and the dot marks the stable steady state value of [λcI] total . Accordingly, [λcI] total → [λcI] * total asymptotically. The stability of the steady-state values of [RF P ] and [LuxI] can be similarly dealt with.

Figure 6 . 6 :

 66 Figure 6.6: Imposed time evolution of the inputs (above) vs. time evolution of the outputs (below).

Figure 6 . 7 :Figure 6 . 8 :Figure 6 . 9 :Figure 6 . 10 :

 676869610 Figure 6.7: GFP signal when [aT c] switches from off to on and [IP T G] switches from on to off.

  -6.24), we consider low concentrations of input signals (i.e. the concentration is zero) and high concentrations of these signals (i.e. [IP T G] = 2•10 -5 , [Ara] = 1 • 10 -5 and [aT c] = 5 • 10 -6 ), and we analyze how output concentrations [RFP] and [GFP] change when the parameters above mentioned vary. We recall that one perspective of this study would be to put this gene network into E. coli bacteria.

Figure 6 . 11 :

 611 Figure 6.11: Dependence of [RF P ] on t rf p 1/2 . Left: for low [IP T G] and high [Ara] and [aT c] concentrations, the increase of [RF P ] concetration is linear, as described by equation (6.9). Right: for high [Ara] and low [IP T G] and [aT c] concentrations, [RF P ] increases also in a linear way, but stays well below the threshold for physical relevance.The red dot corresponds to the value of the parameter considered in Table6.3.

Figure 6 . 12 :

 612 Figure 6.12: Dependence of [RF P ] on k rf p max . Left: for low [IP T G] and high [Ara] and [aT c] concentrations, the increase of [RF P ] concetration is linear, as described by equation (6.9). Right: for high [Ara] and low [IP T G] and [aT c] concentrations, [RF P ] increases also in a linear way, but stays well below the threshold for physical relevance.The red dot corresponds to the value of the parameter considered in Table6.3.

Figure 6 . 16 :

 616 Figure 6.16: Dependence of [GF P ] on k lxr d . Left: for high [IP T G] and low [Ara] and [aT c] concentrations, the [GF P ] concentration decreases as k lxr d increases. Right: for low [IP T G] and high [Ara] and [aT c] concentrations, the decrease in [GF P ] concentration saturates for sufficiently large k lxr d , and always stays well below the threshold for physical relevance. The red dot corresponds to the value of the parameter considered in Table 6.3.

Figure 6 . 17 :

 617 Figure 6.17: Dependence of [GF P ] on k cmp d . Left: for high [IP T G] and low [Ara] and [aT c] concentrations, the [GF P ] concentration decreases as k lxr d increases. Right: for low [IP T G] and high [Ara] and [aT c] concentrations, the decrease in [GF P ] concentration saturates for sufficiently large k cmp d, and always stays well below the threshold for physical relevance. The red dot corresponds to the value of the parameter considered in Table6.3.

Figure 6 . 18 :

 618 Figure 6.18: Dependence of [GF P ] on t gf p 1/2 . Left: for high [IP T G] and low [Ara] and [aT c] concentrations, [GF P ] concentration increases linearly as t gf p 1/2 increases. Right: for low [IP T G] and high [Ara] and [aTc] concentrations, the increase in [GF P ] concentration is also linear, as described by equation (6.15), but stays well below the threshold for physical relevance. The red dot corresponds to the value of the parameter considered in Table6.3.

Figure 6 . 19 :Figure 6 . 20 :

 619620 Figure 6.19: Dependence of [GF P ] on k λ max . Left: for high [IP T G] and low [Ara] and [aT c] concentrations, [GF P ] concentration does not vary that much. Right: for low [IP T G] and high [Ara] and [aT c] concentrations, the decrease in [GF P ] concentration saturates for sufficiently large values of k λmax , and always stays well below the threshold for physical relevance. The red dot corresponds to the value of the parameter considered in Table6.3.

Figure 6 . 21 :

 621 Figure 6.21: Dependence of [GF P ] on k ahl max . Left: for high [IP T G] and low [Ara] and [aT c] concentrations, the increase in [GF P ] concentration saturates for sufficiently large values of k ahl max . Right: for low [IP T G] and high [Ara] and [aT c] concentrations, [GF P ] increases when k ahl max increases, but stays well below the threshold for physical relevance. The red dot corresponds to the value of the parameter considered in Table 6.3.

Figure 6 . 22 :

 622 Figure 6.22: Dependence of [GF P ] on k lux syn . Left: for high [IP T G] and low [Ara] and [aT c] concentrations, the increase in [GF P ] concentration saturates for sufficiently large values of k lux syn . Right: for low [IP T G] and high [Ara] and [aT c] concentrations, [GF P ] increases when k lux syn increases, but stays well below the threshold for physical relevance. The red dot corresponds to the value of the parameter considered in Table 6.3.

Figure 6 . 23 :

 623 Figure 6.23: Dependence of [GF P ] on k cmp a . Left: for high [IP T G] and low [Ara] and [aT c] concentrations, the increase in [GF P ] concentration saturates for sufficiently large values of k cmp a . Right: for low [IP T G] and high [Ara] and [aT c] concentrations, [GF P ] increases when k cmp a increases, but stays well below the threshold for physical relevance.The red dot corresponds to the value of the parameter considered in Table6.3.

Figure 6 . 24 :

 624 Figure 6.24: Dependence of [GF P ] on k gf p max . Left: for high [IP T G] and low [Ara] and [aT c] concentrations, [GF P ] linearly increases when k gf p max increases, as described by equation (6.15). Right: for low [IP T G] and high [Ara] and [aT c] concentrations, [GF P ] also increases when k gf pmax increases, in a linear way, but stays well below the threshold for physical relevance. The red dot corresponds to the value of the parameter considered in Table6.3.

Figure 7 .

 7 1 (b)-(c).

  (a) Circulant network A: every node is connected to the next one counterclockwise. (b) Circulant network B: every node is connected to the next two nodes counterclockwise. (c) Circulant network C: similar to B, except the outer connection direction is counterclockwise and the inner direction is clockwise. (d) Circular network D: every node is connected to the next three nodes counterclockwise.

Figure 7 . 1 :

 71 Figure 7.1: Example of directed circulant networks (upper row) and their Laplacian eigenvalue patterns (lower row). The eigenvalues are plotted in the complex plane for the network sizes n = 20 (black dots), and n = 100 (blue crosses); the red circle depicts the shifted unit circle (3.11) in the complex plane.

Figure 7 . 2 :

 72 Figure 7.2: Delay margin as a function of the network size n for various types of circulant networks: type A (blue circles), type B (green points), type C (red crosses), and type D (purple triangles).

Figure 7 . 3 :Example 7 . 1 .Figure 7 . 4 ,

 737174 Figure 7.3: Time evolution of system (3.7) for directed circulant networks of type C of various sizes n. If we fix α = 1, at a chosen delay value τ = 0.9, the critical network size for consensus equals 24, as read off from Figure 7.2. For smaller sizes Network C reaches consensus, while for larger sizes it diverges. The trajectories of only 5 nodes are shown for clarity.

Figure 7 . 4 : 1 . 5 ≤ p 1 ≤ 4 .5 and 2 ≤ p 2 ≤ 4 .

 741514224 Figure 7.4: p 1 -p 2 parameter space for Example 7.1.

Figure 7 . 5 :

 75 Figure 7.5: τ 1 -τ 2 parameter space for Example 7.2. Characteristic equation (7.3) has no roots on the right half-plane when (τ , τ 2 ) is in region A and two roots with positive real part when (τ 1 , τ 2 ) is in region B.

Figure 7 . 6 :Example 7 . 4 . 5 ) 3 2 , -3 2 , system ( 7 . 5 ) 2 2.

 7674532752 Figure 7.6: p 1 -p 2 parameter space for Example 7.3. Point (p 10 , p 20 ) = (1, 2) corresponds to double root at ω 0 = 1. D = 0. The stability crossing curve T does not have a cusp.

Figure 7 . 7 : 2 2

 772 Figure 7.7: p 1 -p 2 parameter space for Example 7.4. Point (p 10 , p 20 ) = ( 3 2 , -3 2 ) corresponds to ω 0 = √ 2 2 . S-sector is empty. In G-sector, characteristic equation (7.5) has two unstable roots.

7. 3

 3 Designs of a biochemical multiplexer: the stochastic model Part of the research presented in this section was carried out during a visit to the University of Lethbridge, Alberta, Canada, in 2017, in the context of a further collaboration with Prof. Marc R. Roussel. The financial support was provided by Prof. Roussel from his Natural Sciences and Engineering Research Council (NSERC) grant, and by L2S laboratory, France.

Figure 7 . 8 :

 78 Figure 7.8: The synthetic biology design cycle.

  obtain the same decay kinetics as in equation (7.7), both λcI and (λcI) 2 must decay with the same rate constant:

  on, the rfp gene produces proteins according to the delayed mass-→ G rfp on (t) + RFP(t + τ rfp ).(7.27)The kinetics of the second term of this equation are given by the decay reaction RFP

Definition A. 3 .

 3 We consider an open subset Ω of the complex plane and a function f ∈ H(Ω\S), where S ⊂ Ω is a (finite or not) set containing the poles z i of f , at each of which the function must have a Laurent series. Then, we say that f is a meromorphic function on the open subset Ω.We notice that actually a meromorphic function on Ω is holomorphic on all Ω except a set of isolated points z i . Definition A.4. Let f be a meromorphic function on the open subset Ω of the complex plane, s ∈ Ω and f (z) = n∈Z c n (z -s) n the Laurent series of f in the neighborhood of s. Then, we call residue of f around the point s: Res s (f ) = c -1 . Proposition A.5. If Ω is a domain where the function f ∈ H(Ω) is non-null and z belongs to the set Z(f ) denoting the zeros of f , then z is a zero of order Res z f f , called logarithmic residue of f in z.

Proposition A. 6 ( 3 .

 63 Calculating residues).Let Ω be an open subset of the complex plane, f and g two meromorphic functions on Ω and z 0 ∈ Ω. Then:1. If z 0 is a simple pole of f , then Res z 0 (f ) = lim z→z 0 (z -z 0 )f (z); 2. If z 0 is a k th -order pole of f , then Res z 0 (f ) = lim z→z 0 d k-1 dz k-1 (z-z 0 ) k (k-1)! f (z) ; If z 0 is a simple pole of f g , then Res z 0

Proposition A. 7 (

 7 Argument Principle). Let γ denote a closed continuous curve without double points and Ω the closed interior of γ. The function f (z) is assumed to be regular in Ω, except possibly at finitely many poles, finite and non-zero on γ. As the point z moves along γ in the positive sense, the point w = f (z) describes a certain closed continuous curve the winding number of which is equal to the number of zeros inside γ minus the number of poles inside γ. The proposition remains true also when f (z) is only continuous and non-zero on γ.

Theorem A. 8 (

 8 Rouché's Theorem). We suppose that f (z) and g(z) are two functions that are regular in the interior of Ω, continuous in the closed domain Ω and that furthermore |f (z)|> |g(z)| on the boundary γ of Ω. Then the function f (z) + g(z) has exactly the same number of zeros inside Ω as f (z).

  1.1 Example of how a network can be depicted. . . . . . . . . . . . . . . . . . 1.2 Example of cooperative behaviour in nature. . . . . . . . . . . . . . . . . 2.1 Example of loops in a graph. . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Multiple edges in a graph. . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Example of a simple graph. . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 A summary of walk, trail, path, circuit and cycle definitions. . . . . . . .

2.5 Example of trees. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Examples of complete graphs. . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Example of an undirected, an oriented and a directed graph. . . . . . . . 2.8 Examples of bipartite graphs. . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 The central dogma of molecular biology. . . . . . . . . . . . . . . . . . . . 2.10 The basic structure of a gene. . . . . . . . . . . . . . . . . . . . . . . . . .

  • 10 -15 L [4] Dans cette thèse, j'aborde le problème de consensus pour les réseaux orientés et nonorientés ayant une topologie fixe. Nous proposons un nouvel algorithme de consensus résultats concerne les réseaux orientés, pour lesquels nous démontrons l'existence d'une liaison fortement non-triviale entre comment régler les deux paramètres et la topologie du réseau. Celle-ci représente une autre différence, si nous comparons notre algorithme avec d'autres protocoles du consensus dans la littérature : notre condition pour atteindre le consensus est une condition nécessaire et suffisante pour les réseaux non-orientés, mais pas suffisante pour les réseaux orientés. Nous testons notre algorithme sur des réseaux ayant des topologies différentes et nous présentons les résultats des simulations qui montrent que la vitesse de convergence peut être considérablement améliorée en utilisant notre protocole. Cette amélioration est notamment visible dans le cas des réseaux qui ne sont pas très bien connectés. Une contribution considérable de notre

	Transcription rate for E. coli	= 40 nucl/s [5]
	Translation rate for E. coli	= 15 amino acids/s
	λ repressor residue count	= 236 amino acids
	LuxR length	= 750 nucl.

xi Résumé étendu en français qui comporte des retards et où tous les agents du réseau utilisent de l'information passé pour estimer les opinions (ou les états) futurs des agents voisins. Une contribution conceptuelle du modèle proposé concerne la nature différente des retards : dans la littérature au sujet du problème du consensus, les retards modélisent généralement soit les temps morts nécessaires aux agents du réseau pour traiter l'information (le temps de réaction, par exemple), soit les temps nécessaires à la transmission effective de l'information (le temps nécessaire pour que l'information circule parmi les agents du réseau considéré). Dans notre modèle, les retards ont une nature différente, car ils peuvent anticiper. De plus, l'algorithme que nous proposons a l'avantage d'être facile à implémenter et il ne comporte pas de techniques sophistiqués de nature prédictive ou adaptative, habituellement retrouvées dans les autres résultats de la littérature spécialisée. Encore mieux, nous démontrons que le modèle proposé est extrêmement efficace quand un choix optimal des paramètres est fait. Plus exactement, nous manipulons deux paramètres décrivant dans quelle mesure les agents devraient chercher dans le passé pour estimer les états futurs. Une autre contribution du chapitre 3 concerne le calcul des conditions nécessaires et suffisantes qui dépendent de ces deux paramètres, pour que le réseau atteigne le consensus. Nous étudions également un indicateur appelé delay margin et qui décrit dans quelle mesure les retards peuvent augmenter dans le réseau. Les résultats obtenus ont plusieurs niveaux d'interprétation. Une première interprétation importante est que la nature anticipative de l'algorithme proposé peut considérablement améliorer les performances, mais si les agents essayent d'anticiper trop loin dans le futur, le système peut diverger. Cette conclusion n'est pas très surprenante : en utilisant de l'information qui n'est pas mise à jour, nous pouvons rendre le système instable. Cependant, une deuxième xii Résumé étendu en français xiii interprétation des algorithme consiste en un réglage des paramètres assurant une vitesse de convergence extrêmement rapide, qui est valable pour n'importe quelle topologie du réseau. En effet, ce choix des paramètres représente un choix optimal et assure une vitesse de convergence qui correspond à une connectivité algébrique égale à 1. Nous discutons également quelques propriétés et remarques intéressantes sur les réseaux orientés : nous considérons quelques topologies circulaires et montrons que les valeurs propres du Laplacien suivent

  3, the basic notation and review relevant background information from graph theory. Section 3.2 also introduces the proposed model and its contribution with respect to other models in the literature. All other sections of this chapter present our results: we have presented a graphical depiction of

the dominant roots for arbitrary eigenmodes of the system as a function of the delay value τ , thereby deriving a universal picture applicable for any undirected graph. By giving necessary and sufficient conditions, we have shown that (3.7) reaches consensus from arbitrary initial states if and only if ατ < 1. Another result focuses on directed networks, where we prove that the same condition is not sufficient.

  1, 2}. Moreover, if γ 11 has the form (4.22), with the parameters γ 12 , γ 21 and γ 22

	left free,
	γ 11 =	245 6	-	5 3	γ 22 -γ 21 -	5 3	γ 12 ,	(4.22)

  5.3. Consider s a ∈ N • δ (s 0 ), δ > 0 sufficiently small, and let p 1a = p 1 (s a ), p 2a = p 2 (s a ) as defined in Proposition 5.2. Then

			∂ ∂s	q(s, p 1a , p 2a )	s=sa	= 0.	(5.9)
	Proof. Let
			s a = s 0 + uγ, |γ| = 1,
	then,
	∂q ∂s s=sa p 1 =p 1a p 2 =p 2a	=	∂q ∂s s=s 0 p 1 =p 10 p 2 =p 20	+	∂ 2 q ∂s 2 s=s 0 p 1 =p 10 p 2 =p 20	γu +	∂ 2 q ∂s∂p 1	s=s 0 p 1 =p 10 p 2 =p 20	∂p 1 ∂u u=0	u
		+	∂ 2 q ∂s∂p 2 s=s 0 p 1 =p 10 p 2 =p 20

  As the tangent direction of the local stability crossing curve T (ω 0 ,p 10 ,p 20 ) at the cusp(p 10 , p 20 ) is ∂ 2 p 1 ∂u 2 , ∂ 2 p 2, it can be easily seen that the T - (ω 0 ,p 10 ,p 20 ) is in the counterclockwise side of T + (ω 0 ,p 10 ,p 20 ) if we may reach the direction of ∂ 2 p 1 ∂u 2 , ∂ 2 p 2

		T
	∂u 2	+
		∂p 2 ∂s	∂ 2 p 2 ∂u 2	+	p 2 =p 20 p 1 =p 10 s=s 0	.

  .2) for k = 3.

	Differentiate (5.15) again with respect to u, taking into account (5.50); we obtain
	∂q ∂p 1	∂ 4 p 1 ∂u 4 +	∂q ∂p 2	∂ 4 p 2 ∂u 4	u=0 s=s 0 p 1 =p 10 p 2 =p 20	= -	∂ 4 p ∂s 4 γ 4	s=s 0 p 1 =p 10 p 2 =p 20

Table 5 .

 5 1, where s k stands for sin k, and c k stands for cos k. The approximate values of the coefficients in Table 5.1 are: a 20 = -2.19272, a 11 = 6.27284, a 02 = -2.19272, a 10 = 2.03748, a 00 = -2.52733, and a 01 = 5.60094.

Table 5 . 1 :

 51 Coefficients' values for illustrative example subsection.

  The usual cells used in genetic engineering are Escherichia coli (E. coli) cells. Researchers usually prefer E. coli cells to other bacteria

	[IPTG]
	P tet	λcI gene	LuxR
		P	luxR gene
			LuxR AHL
	[aTc]		LuxI
			luxI gene		GFP	[GFP]
			P lux	P Tac	gfp gene
	TetR
	tetR gene
	AraC
	araC gene
					RFP	[RFP]
	[Ara]		P BAD	P tet	rfp gene

Table 6 . 1 :

 61 ) Parameters for promoters P BAD , P tet , P lux , and P T ac (see Supplementary Information of[START_REF] Seok Moon | Genetic programs constructed from layered logic gates in single cells accessed terms of use detailed terms[END_REF]).

		P BAD	P tet	P lux	P T ac
	L	[Ara]	[aTc]	[LuxR:AHL]	[IPTG]
	n	2.2	2	1.7	1.6
	K D 10 -4 M 1.5 • 10 -10 M	10 -8 M	3 • 10 -6 M
	K 1	12	51	8.3	1950
	K 2	4.4	-	-	-
					State
	P BAD	O2	O1	I1-I2
				AraC:Ara
				complex
		O2	O1	I1-I2
			O1	O2
				AraC
				I1-I2
	P Tac		Operator
			Lacl
			Operator
	P lux		Operator
			LuxR:AHL
			complex
			Operator
	P tet		O1	O2
				tetR
			O1	O2
			tetR
			O1	O2
			tetR	tetR
			O1	O2
	P		O1	O2
	λcI
				(λcI) 2
			O1	O2
			(λcI) 2
			O1	O2
			(λcI) 2	(λcI) 2
			O1	O2
		Figure 6.3: Promoters.

  7 • 10 -15 L. We have used this equation to set deterministic simulations, but when considering the stochastic model (as in Chapter 7) this equation can be dropped from the model, as there is no interest to consider the dynamics of a gene that is always on (as LuxI gene, or LacI gene) in a stochastic model. Nevertheless, we have considered the equation (6.12) in the preliminary model presented in this chapter, in order to verify that the steady-state value, and the transient regime value of LuxI make sense from a biophysical point of view.

  Remark 6.2. Note that the unique steady-state (6.[START_REF] Pavella | Transient stability of power systems : theory and practice[END_REF])-(6.22) does not depend on timedelays τ rf p , τ λ , τ lxr and τ gf p . Thus, P tet , P BAD , P λcI , P T ac , and P lux have constant values depending on the input signals. Remark that the concentrations [LuxI] and [AHL] do not depend on the input signals. This is due to the fact that the LuxI gene is always on.In the sequel, we solve the set of equations (6.23)-(6.29) and show that each concentration ([RF P ], [λcI] total , [LuxR], [LuxR : AHL], [LuxI], [AHL] and [GF P ]) can be written as a function of the parameters and depending on input concentrations ([aT c],[IP T G] and [Ara]). From equation (6.26), we can easily see that equation(6.19) holds.

	k rf p max P * tet P * BAD -k rf p d [RF P ] = 0	(6.23)
	k λ max P * tet -k λ d [λcI] total = 0	(6.24)
	k lxr max P * λcI -k lxr d [LuxR]
	-k cmp a [LuxR][AHL] + k cmp d [LuxR : AHL] = 0	(6.25)
	k lux syn -k lux d [LuxI] = 0	(6.26)
	k ahl max [LuxI] -k ahl d [AHL]
	-k cmp a [LuxR][AHL] + k cmp d [LuxR : AHL] = 0	(6.27)
	k cmp a [LuxR][AHL] -k cmp d [LuxR : AHL] = 0	(6.28)
	k gf p max P * T ac P * lux -k gf p d [GF P ] = 0	(6.16) (6.29)
		(6.17)
		(6.18)
		(6.19)
		(6.20)
		(6.21)
	* T ac P * lux	(6.22)
	Proof. At steady-state, the concentrations of all chemical species involved in the gene
	network are constant; this means that the left-hand side of equations (6.9)-(6.15) is
	zero. Moreover, we don't consider uptake kinetics or decay kinetics for input signals;
	this means that each promoter's activity is a constant at a given concentration of the
	input signals. Thus, we obtain the following set of nonlinear equations.

Table 6 . 4 :

 64 Steady-state response for constant input signals. For each signal we show the concentration order of magnitude in M (left sub-column), and the corresponding logical values in bold (right sub-column).Proposition 6.3. The unique steady-state (6.16)-(6.22) is stable.Proof. Parts of the stability analysis of this device reduce to one-dimensional analyses.We first consider equation (6.10). P tet depends on the select signal [aT c]. Since we don't consider uptake or decay kinetics for the input signals, P tet is in fact constant for a given

	6.5.1 Stability analysis
	set of inputs. The differential equation (6.10) reduces to
	d [λcI] total dt	= k λ max P tet -k λ d [λcI] total .

  ). P tet depends on the select signal [aT c], and P BAD on the input signal[Ara]. Since we don't consider uptake or decay kinetics for the input signals, Figure 6.4: Phase-space flow on the [λcI] total axis. ] > [RF P ] * , where the steady-state value [RF P ] * is defined as in equation (6.16). Thus, we can immediately show that [RF P ] → [RF P ] * asymptotically. The stability of the [LuxI] * steady-state can be proven in the same manner. We can where the steady-state value [LuxI] * is defined in equation (6.19). Following the same reasoning as in the case of the [λcI] total analysis, we can see that [LuxI] → [LuxI] * asymptotically. The stability analysis is different in the case of [LuxR], [AHL], and [LuxR : AHL], which are coupled to each other and to the value of [λcI] total . We can dispose of the latter dependence, as we already know that [λcI] total tends asymptotically to [λcI] *total , by showing that for any positive value of [λcI] total , the system of equations (6.8) has a unique positive solution. This can be proven geometrically. First, rewrite Geometry of the system (6.8). The parameters of the two equations affect the intercept of the line, which must however remain positive, and the steepness of the parabola. There is always exactly one intersection.([λcI], [(λcI) 2 ]]) plane. Note that equation (6.33a) is the equation of a line with a positive intercept ([λcI] total /2) and a negative slope (-1). At the same time, equation (6.33b)

		2 ]
		λ [( cI)
	P tet and P BAD are in fact constants for a given set of inputs. The differential equation
	(6.9) simplifies to
		d[RF P ] dt	= k rf p max P tet P BAD -k rf p d [RF P ].
	Similar to the case of [λcI] total stability analysis, we can write
	Figure 6.5:	d[RF P ] dt dt d[RF P ]	> 0 = 0	λcI] if [RF P ] < [RF P ] * , [ if [RF P ] = [RF P ] * ,
	and if [RF P analyze equation (6.12) by writing d[RF P ] dt < 0
		d[LuxI] dt	> 0	if [LuxI] < [LuxI] * ,
		d[LuxI] dt	= 0	if [LuxI] = [LuxI] * ,
	and if [LuxI] > [LuxI] system (6.8) to the form d[LuxI] dt < 0
		[(λcI) 2 ] =	1 2	([λcI] total -[λcI]) ,	(6.33a)
		[(λcI) 2 ] = [λcI] 2 /K λ .	(6.33b)
	We know that system (6.8) has a solution wherever these two curves intersect in the

total

[ λcI] * * ,

  6.11-6.24 is that, overall, the gene network's behaviour does not change when one half-life time or one synthesis rate varies: the proposed circuit still follows the truth table depicted in Section 6.2, which corresponds to the desired functionality.

Table 6 .

 6 3. Dependence of [GF P ] on t λ 1/2 . Left: for high [IP T G] and low [Ara] and [aT c] concentrations, [GF P ] concentration does not vary that much. Right: for low [IP T G] and high [Ara] and [aT c] concentrations, the decrease in [GF P ] saturates for sufficiently large t λ 1/2 , and always stays well below the threshold for physical relevance. The red dot corresponds to the value of the parameter considered in Table 6.3. Dependence of [GF P ] on t lux 1/2 . Left: for high [IP T G] and low [Ara] and [aT c] concentrations, the increase in [GF P ] concentration saturates for sufficiently large values of t lux 1/2 . Right: for low [IP T G] and high [Ara] and [aT c] concentrations, [GF P ] concentration increases, but stays well below the threshold for physical relevance. The red dot corresponds to the value of the parameter considered in Table 6.3. Figure 6.15: Dependence of [GF P ] on t ahl 1/2 . Left: for high [IP T G] and low [Ara] and [aT c] concentrations, the increase in [GF P ] concentration saturates for sufficiently large values of t ahl 1/2 . Right: for low [IP T G] and high [Ara] and [aT c] concentrations, [GF P ] concentration increases, but stays well below the threshold for physical relevance. The red dot corresponds to the value of the parameter considered in Table 6.3.

	1.0280	#10 -8			5	#10 -10
					4
	1.0279
	(M)			(M)	3
	[GFP]			[GFP]	2
	1.0279
					1
	1.0278	0	t 1/2 6 (min) 2	4	0	t 1/2 6 (min) 10 -3 10 -2 10 -1 10 0 4
	Figure 6.13:

Table 6 .

 6 3. order for the two designs. We can easily check this by using equations (6.22) and (6.16):we remark that increasing t gf p 1/2 and t rf p 1/2 , [GF P ] * and [RF P ] * will increase, respectively.If in the case of the former design the rise and fall times are as indicated in Subsection 6.5.2, switching times are different for the latter design. For GFP the rise and fall times are 80 minutes, and for RFP the rise and fall times are about 60 minutes, except

for the case when IP T G is on and aT c switches from off to on, when the GFP fall time is 140 minutes. In order to improve this switch time, the decay rate of LuxR should be engineered as in Table

6

.3.

  d[RFP] dt = k rfp max P tet (t -τ rfp )P BAD (t -τ rfp ) -k rfp d [RFP] d [λcI] total dt = k λ max P tet (t -τ λ ) -k λ d [λcI] total d[LuxR] dt = k lxr max P λcI (t -τ lxr ) -k lxr d [LuxR] Tac (t -τ gfp )P lux (t -τ gfp ) -k gfp

			(7.6)
			(7.7)
		-k LA +c [LuxR][AHL] + k LA -c [LuxR • AHL]	(7.8)
	d[AHL] dt	= k ahl max -k ahl d [AHL]
		-k LA +c [LuxR][AHL] + k LA -c [LuxR • AHL]	(7.9)
	d[LuxR • AHL] dt	= k LA +c [LuxR][AHL] -k LA -c [LuxR • AHL]	(7.10)
	d[GFP] dt	= k gfp max P d [GFP]	(7.11)

  The effective chemical reactions can be derived in a similar way for other genes in the model. The LuxR expression implies the process lxr on/off represents the two states of the λcI promoter. Also, the gene expression produces protein according to the delayed mass-action reaction Similarly to the case of GFP expression control, we can write k lxr on = k sw P lxr on (t), k lxr off = k sw P lxr off (t),

					(7.19a)
					(7.19b)
	where
		P lxr on (t) = P λcI (t),	(7.20a)
		P lxr off (t) = 1 -P λcI (t).	(7.20b)
		G lxr off	k lxr on ------off k lxr	G lxr on ,	(7.17)
	where G G lxr on (t)	k lxr max ---→ G lxr on (t) + LuxR(t + τ lxr ).	(7.18)

  .7) gives the same kinetics as P tet (t), P λ off (t) = 1 -P tet (t).

	and
		P λ on (t) =
	the process
		G λ off	k λ on -----off k λ	G λ on
	with
	G λ on (t)	k λ max ---→ G λ on (t) + λcI(t + τ λ ).
	The corresponding rate constants are selected as above:
		k λ on = k sw P λ on (t),	(7.21a)
		k λ off = k sw P λ off (t),	(7.21b)

Roughly speaking, we define here the delay margin as the maximum delay value in the proposed model (3.7), guaranteeing the consensus reaching.

We have used the word "small" in a sense analogous to "small solution": a small sector is contained by a sector with straight sides with arbitrarily small angle when the neighborhood is sufficiently small.

Acknowledgements

I would like to express my special appreciation and thanks to my supervisors, especially to Silviu Niculescu for being a good advisor. This thesis wouldn't have been possible without his support during these three years. He opened many doors to the research and academic world for me. What I appreciate most is his ability to give good advice Acknowledgements

Abbreviations

AHL

We consider the characteristic matrix M (s, τ ) = s + 2e -sτ 1 cos(1) -1 1

-1 -e -sτ 2 + 2e -sτ 1 sin(1) s -1 .

(5.23)

We compute the characteristic function corresponding to the characteristic matrix (5.23), and obtain the quasi-polynomial (5.18) of Example 5.1. We recall that for τ = τ 0 = (τ 10 , τ 20 ) T = (1, 2), the characteristic function (5.18) has a double imaginary root s = s 0 = ±iω 0 at ω 0 = 1. Thus, the eigenvalue s 0 = ±iω 0 = ±i in Example 5.1 is double and non semi-simple. This implies that there exist u 0 and v 0 , a right, respectively left, eigenvector, and u 1 and v 1 a right, respectively left, generalized eigenvector, such that conditions M 0 u 0 = 0 (5.24)

are simultaneously satisfied (see [START_REF] Campbell | Zero singularities of codimension two and three in delay differential equations[END_REF], Chapter 7 in [START_REF] Hale | Introduction to functional differential equations[END_REF]), where M 0 = M (s 0 , τ 0 ), (5.28)

∂M (s, τ ) ∂s s=s 0 τ =τ 0 .

(5.29) Proposition 5.12. Eigenvectors u 0 , u 1 , v 0 and v 1 satisfy v T 0 M 0 1 u 0 = 0 (5.30)

(5.31)

Proof. Multiply equation (5.25) by v T 0 on the left and obtain (5.30) in view of (5.26). Multiply equation (5.25) by v T 1 on the left and obtain v T 1 M 0 u 1 = -v T 1 M 0 1 u 0 . Multiply equation (5.27) by u 1 on the right and obtain v T 1 M 0 u 1 = -v T 0 M 0 1 u 1 . Thus, we have (5.31) because v T 1 M 0 1 = 0 and M 0 1 u 1 = 0.

We consider a simple case, where the delay parameters are under a small perturbation ε > 0:

τ 2 (ε) = τ 20 + ετ 21 .

(5.33)

As τ (ε) is smooth, we can write a Taylor series expansion of the characteristic matrix (5.20). Moreover, provided that the eigenvalue is non semi-simple, we can write a Puiseux series expansion of it and of the corresponding right eigenvector (see, for instance, [START_REF] Kato | Perturbation Theory for Linear Operators[END_REF]):

u 0 (ε) = u 0 + ε 1 2 w 1 + . . . (5.35) We can now replace equations (5.32), (5.33), (5.34), and (5.35) in the eigenvalue problem (5.22) and write

Note that we use series expansion of exponential functions, as ε → 0 (for instance, we write e -εs 0 τ 10 = 1 -εs 0 τ 11 +

). We collect the terms of equal powers of ε, more precisely we collect only the first three orders (ε 0 , ε

) and obtain

)

where the matrices M 0 , M 0 1 , M 0 2 , and M 1 1 are given by

M 0 1 = I + τ 10 A 1 e -s 0 τ 10 + τ 20 A 2 e -s 0 τ 20 , (5.39)

M 1 1 = τ 11 s 0 A 1 e -s 0 τ 10 + τ 21 s 0 A 2 e -s 0 τ 20 .

(5.41)

Notice that equations (5.38) and (5.39) are explicit expressions of (5.28) and (5.29), respectively. For the sake of uniqueness of u 0 and without any loss of generality, we can use the normalization

(5.42)

Divide equation (5.36) by s 1 = 0 and obtain 1 s 1 M 0 w 1 = M 0 u 1 in view of (5.25). Since M 0 = 0, we can write w 1 as

(5.43)

Replace (5.43) in equality (5.37) and multiply by v T 0 on the left. We obtain

Now we use the second equality from the normalization (5.42), and condition (5.30) from Proposition 5.12 in equation (5.44), and write the expression of s 1 :

Given the continuity properties of the spectrum (Chapter 1 in [START_REF] Michiels | Stability and Stabilization of Timedelay Systems. An Eigenvalue Based Approach[END_REF]) and the Puiseux series expansion (5.34), it is easy to see that under a small perturbation of delay parameters (5.32)-(5.33) the double characteristic root s 0 = iω 0 splits up into two simple characteristic roots, which will move to the right or left half-plane function of the sign of Re(s( )). These results can be summarized in the following proposition and remarks.

Proposition 5.13. Assume that a double and non semi-simple characteristic root of system (5.21) is located on the imaginary axis, but not at the origin. Then, under a small perturbation of delay parameters of the form (5.32)- (5.33) this double characteristic root splits up into two simple roots, and each one of them will move towards stability (instability) if Re(s(ε)) < 0 (Re(s(ε)) > 0).

Remark 5.14. Provided that Re(s 0 ) = 0, we might suppose, in general, that the sign of Re(s(ε)) is given by the sign of Re(s 1 ) as defined in equation (5.45), when ε is very small.

However, knowing the sign of Re(s 1 ) is not enough to conclude over a global tendency of the double root to move towards stability or instability, as we shall illustrate in the sequel. This is why we take a further step, multiply equation (5.37) with v T 1 on the left, and obtain the following expression of s 2 ,

in view of normalization (5.42). Other terms of s(ε) might be found in a similar way.

Remark 5.15. In order to approximate the value of s(ε), as a function of ε, τ 11 , and τ 21 , we use s(ε) defined as in equation (5.34), where s 1 is given by (5.45), s 2 is given by (5.46), with M 0 1 , M 0 2 , and M 1 1 given by (5.39)-(5.41), u 0 and v T 0 are normalized right, respectively left, eigenvectors, satisfying (5.42), and u 1 and v 1 are generalized right, respectively left, eigenvectors, satisfying equations (5.24)-(5.27).

Remark 5.16. Note that the double non semi-simple root will split up into two simple roots. The plus sign from equation (5.45) corresponds to one of these simple roots, and the minus sign corresponds to the other simple root.

Remark 5.17. Roughly speaking, the double root will split up into two simple roots that will follow one of the two tendencies: either one root moves towards a half-plane, and the other root towards the other half-plane, or both roots moves towards the same halfplane. This means that we have two types of qualitative behaviour when a perturbation ε arises. We illustrate both situations in Example 5.1. In the same manner, region II the sector ECF , region III the sector F CG, and so on. Thus, the neighbourhood of s 0 shown in left side of Figures 5.12 We shall only show the case where they are in the S-sector. Their locations do not affect the validity of the conclusion. When one or both points F and I are outside of the S-sector, the proof for the G-sector is slightly more involved, but still possible. Similar to the case discussed in Section 5.3 (see Corollary 5.5) we can show that (p 1 (s) p 2 (s)) is a bijection from R to R when s is restricted to R, with R a region from the set {I, II, . . . V III}, and R the corresponding region in the set {I , II , . . . , V III }.

Consider Case i. The S-sector (in a sufficiently small neighbourhood) can be expressed as (II ∪ III ) ∩ V ∩ (V I ∪ V II ) ∩ V III , as depicted in Figure 5.12 right. But the corresponding regions are (II ∩III), which is in the right-half plane, and V , (V I ∪ V II) a conceptually different, algebraic method. The differences between the two approaches are commented on, and it is shown that they give the same results.

The migration of imaginary characteristic roots of multiplicity three and four is discussed in Section 5.4. Note that both studies make use of a conventional approach, without using Puiseux series.

Under the least degeneracy assumption, neither the triple root case, nor the quadruple case may be stable under a small deviation of the parameters. More precisely, in the case of triple roots the local stability crossing curve divides the neighbourhood of (p 10 p 20 ) into two regions of roughly equal size. In one region, there are two roots on the right-half plane, and one root on the left-half plane. In the other region, there are two roots on the left-half plane, and one root on the right-half plane.

In the case of quadruple roots, the stability crossing curve has a cusp (similar to the double root case), and divides the neighbourhood of (p 10 p 20 ) into a S-sector and a G-sector. When the parameters move into the G-sector, there are two roots on the right-half plane and the other two on the left-half plane. When the parameters move into the S-sector, either there are three right-half plane roots and one left-half plane root, or there are three left-half plane roots and one right-half plane root. We can next write Remark 6.4. Note that the stability analysis is delay independent: the delays drop out of the stability calculations, so their values clearly affect neither the steady-state values, nor their stability. The stability is also independent of the functional forms of the promoter activities.

Remark 6.5. We have chosen a geometrical stability analysis, rather than computing a Jacobian matrix of dimension 7 due to the complexity of the latter task.

Transient regime

We vary the three input concentrations with time and observe the time dependence of in the S-sector.

The stability crossing curve T is depicted in Figure 7.5 for the range 2 ≤ τ 1 ≤ 4, 8.5 ≤ τ 2 ≤ 10.5. We remark that T divides this area into two regions: region A is the region containing S-sector and connected to the origin, and region B contains the G-sector. For τ 1 = τ 2 = 0, the characteristic quasi-polynomial reduces to a polynomial that has only one root equal to a. As a is negative, we conclude that for (τ 1 , τ 2 ) in region A the quasi-polynomial (7.3) has no root with positive real part, and the system is stable.

Next, according to Corollary 5.9 or Theorem 5.8, both imaginary roots at i move to the left half-plane as (τ 1 , τ 2 ) moves from the cusp to the S-sector (region A). Furthermore, according to Theorem 5.8, as (τ 1 , τ 2 ) moves to G-sector (region B), one of the imaginary roots moves to the right half-plane, and the other one to the left half-plane. In other words, as (τ 1 , τ 2 ) moves from region A to region B through (τ 10 , τ 20 ), one of the two imaginary roots moves from the left half-plane to the right half-plane passing through the point i of the imaginary axis, and the other root moves in the left half-plane to touch the imaginary axis at i and then returns to the left half-plane.

Due to symmetry reasons, another left half-plane root moves to the right half-plane through the point -i. Thus, in region B there are two more roots with positive real part, as compared to the region A. Thus, we conclude that there are two roots on the right half-plane when (τ 1 , τ 2 ) is in region B.

Degenerate cases

In this subsection, two examples will be presented to illustrate degenerate cases. The first example shows the local stability crossing curve may not have a cusp when one of the least degeneracy assumptions, D = 0, is violated. The second example shows that the S-sector may be empty.

Moreover, we can further improve the deterministic model, by taking into account the number of plasmids and the dimerisation of [LuxR • AHL], which adds a few other reactions and parameters to our model. The effective chemical reactions can be derived in a similar way, as described above. Thus, once we have all the equations needed for the stochastic model simulation, we can validate the continuous model.

Global conclusions and perspectives

We introduce multi-agent systems and consensus problem in Chapter 3, where a new consensus protocol for networks with anticipatory agents is proposed. The network stability under the proposed protocol is investigated and consensus conditions are given.

Chapters 4 and 5 discuss stability issues for a larger class of systems. Limit cases are considered, when multiple characteristic roots lay on the imaginary axis. More precisely, Chapter 4 presents a set of conditions on the system parameters such that the characteristic root at the origin reaches its maximal multiplicity. Chapter 5 focuses on systems with two parameters and multiple roots on the imaginary axis, except the origin. A criterion to decide if a small deviation on the parameters induces a tendency toward the stability or instability is given for the considered type of systems with two parameters. An extension to multiplicity three and four is also studied in Chapter 5. A gene network model is presented in Chapter 6. This design is similar to a multiplexer circuit: two outputs are controlled using three input signals. The differential equations describing the systems are given, as well as steady-state values obtained analytically and by using simulations. A brief robustness study is presented, showing how the circuit behavior changes when the parameters vary.

The analysis of time-delay dynamics on networks may give some insights into the understanding of complex systems, from engineering to life sciences and economics. Adapting to and addressing problems relevant to our changing world is of great importance for the control community. Thus, more research should be carried out to answer a challenging demand of understanding the self-organization and evolution of large networks with common tasks. As new opportunities to apply control principles and methods are exploding, it is high time for researchers from control theory to join researchers from different domains (e.g. computer science or biochemistry) and come up with new breakthroughs in engineering and technology that can considerably improve the world we live in.

Appendix A

Other mathematical tools relevant for Chapter 4

In this appendix, we briefly review some definitions from complex analysis, linear algebra and numerical analysis. For details, the reader is referred to standard texts such as [START_REF] Mathews | Complex Analysis for Mathematics and Engineering[END_REF] and [START_REF] Allaire | Numerical Linear Algebra[END_REF].

Complex analysis

Definition A.1. A holomorphic function is a complex-valued function of one or more complex variables that is complex differentiable in a neighborhood of every point in its domain.

Next we denote by H(D) the set of holomorphic functions that take values in the domain D. Moreover, we denote by C r,R the annulus centered at the origin, of internal radius r and external radius R, and by C ρ the circle centered at the origin with a radius of ρ > 0.

Theorem A.2 (Laurent). We consider z 0 ∈ C and two real numbers r and R, such that 0 < r < R. Every function f belonging to H(z 0 + C r,R ) can be written as a series, called a Laurent series, as:

with h ∈ C r,R and, for each ρ ∈ (r, R):

plus the one under 1 + g(z) f (z) (which is zero), we can conclude that f (z) and f (z) + g(z) has exactly the same number of zeros inside Ω.

Linear algebra and numerical analysis

A matrix with the terms of a geometric progression in each row (or column) is called a Vandermonde matrix. For instance, an m × n matrix

Sometimes we need a mathematical representation showing a relationship between two objects. We consider two objects belonging to different classes, A and B. An incidence matrix is a matrix having one row for each element of A and one column for each element of B. By convention, its elements are either 1 or 0, 1 if the corresponding row and column are related, and 0 if they are not.

Appendix B

Software Tools

We mention in this appendix a few open source Matlab packages very useful for spectra study in the case of time-delay systems. The results of this thesis have been illustrated not only by using code developed by the author, but also by making use of some Matlab packages named below.

B.1 DDE-Biftool

The second version of DDE-BIFTOOL is a free Matlab package for the stability and bifucation analysis of delay differential equations. It consists actually in a set of Matlab routines that allow to compute and analyze branches of solutions and bifurcations (see for more details the manual http://www.cs.kuleuven.ac.be/publicaties/rapporten/ tw/TW330.abs.html). Generally, the user has to write at least four files in order to define the system one would like to study: an initialization file containing the dimension and the name of the system, a second file describing the system of the form ẋ = f (coefficients, x(t), delayed x(t)),

the Jacobian file that can be replaced with df deriv.m, and another file called sys tau containing the position of the delays in the parameter vector. The manual also provides some step-by-step examples of how to define a time-delay system (with constant delays or state-dependent delays) and what function to use for stability analysis and branch manipulation.

Abstract: This thesis discusses diverse types of interconnected systems through networks. We address networks of agents with cooperative tasks and propose a new consensus protocol with delays and anticipatory agents. We study the consensus reaching conditions for networks organized under the proposed model. Moreover, we derive some theoretical results, which can apply to a more general class of systems, concerning stability issues when the considered system has multiple imaginary roots. In terms of networks, this situation can correspond to the case of switching topology networks, when the network can even be disconnected at some point.

We separately discuss the case of zero characteristic roots, and roots laying on the imaginary axis, except the origin. Finally, we propose a gene network model with a functionality similar to a multiplexer circuit. Thus, we control two outputs with three input signals, and we carry out a stability analysis. We prove the uniqueness and the stability of the network steady states, and validate the continuous and deterministic model with a stochastic model.