Camille Poignard

John Tyson

INDUCING CHAOS IN A GENE REGULATORY NETWORK BY COUPLING AN OSCILLATING DYNAMICS WITH A HYSTERESIS-TYPE ONE

In this paper, we investigate the chaotic behavior of a gene regulatory network modeled by four differential equations and seventeen parameters. This network, called V-system, has been designed to couple in a simple way an oscillating system with one having a bistable switch. After having studied it analytically, we exhibit (by a constructive proof) the mechanism responsible of chaos for a general differential system presenting such a coupling. Namely, given a generic one-parameter family of smooth vector fields on R n presenting a Hopf bifurcation, we prove that under an assumption on the Jacobian at the bifurcation point, we can create such a chaotic system by perturbing the parameter thanks to a hysteresis-type dynamics. Finally, we numerically show that the mechanism highlighted previously takes place in the V-system, for a particular set of values of its parameters.

Introduction

This article is motivated by the growing interest in understanding and controlling biological regulatory systems, two main topics in systems and synthetic biology, where identification and characterization of regulatory units with prescribed dynamical features are essential (see [START_REF] Tyson | Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell[END_REF] and [START_REF] Lu | Next-generation synthetic gene networks[END_REF]). Many metabolic and genetic intracellular regulations are homeostatic, that is they maintain constant some viability parameters (rate, concentration, level, etc...) by perpetually adapting the internal state of the cell to a changing environment. Temperature, pH, osmotic pressure are emblematic, but other examples are the control of internal concentration of metal ions which are essential as trace but become lethal at higher rate (as can be seen in [START_REF] Pécou | A mathematical model for copper homeostasis in enterococcus hirae[END_REF]). In general, breaking of homeostasis can lead important damages, like Wilson's disease in the context of copper homeostasis ( [START_REF] Pécou | A mathematical model for copper homeostasis in enterococcus hirae[END_REF]).

The question dealt with in this article is how to break homeostasis and induce a chaotic dynamics in a given model of homeostatic regulatory unit. Concrete motivation of such a general problem is the investigation of a way to destroy an organism having a stable dynamics by destabilizing its metabolism (see [START_REF] Pécou | Desynchronization of one-parameter families of stable vector fields[END_REF]). Other motivations are the investigations of the role played by chaos in the formation and the behavior of cells, which have begun to be supported by experiments since the last decade (see [START_REF] Furusawa | Chaotic expression dynamics implies pluripotency: when theory and experiment meet[END_REF]). In most cases, the homeostatic regulation is due to a negative feedback, that is to say a mechanism creating a response that inhibits the action starting this mechanism. More rarely, it can also be caused by a positive feedback (an example is the blood clotting process).

More precisely, we aim at inducing chaos in a gene regulatory network (G.R.N), issued from a continuous modeling: this means its dynamical behavior is described by differential equations, that translate in mathematical terms the set of chemical reactions defining this system. It has been invented by J.J Tyson 1 and E. Pécou, and was called V-system in recognition of its origins. Its particularity lies in the fact that it was designed so as to couple in a very simple way, an oscillating system with one having a hysteresis-type dynamics, the interest being that those two dynamics are the typical mechanisms causing respectively a negative and a positive loop (see [START_REF] Tyson | On the existence of oscillatory solutions in negative feedback cellular control processes[END_REF], [START_REF] Tyson | Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell[END_REF]). This guarantees the plausibility of our network and of the emergence of chaos in many other G.R.Ns. Usually, given one G.R.N represented by a differential system, we can reduce its set of chemical reactions (governed by the mass-action law) into a smaller system, thanks to the classical quasi-steady state approximation. This leads to the appearance of regulation functions expressing the nonlinearities effects present in this small system (see [START_REF] Pécou | Mathematical comments on basic topics in systems biology[END_REF], [START_REF] Murray | Mathematical biology[END_REF]). This is the case of the Goodwin model ( [START_REF] Bryan | Oscillatory behavior in enzymatic control processes[END_REF], [START_REF] Othmer | The qualitative dynamics of a class of biochemical control circuits[END_REF]) which describes the production of the mRNA m of a gene, coding for an enzyme that catalyses a chemical reaction of which final product z represses the production of m. The repression of the mRNA is described by the following equation:

ṁ = k 1 + z j 2 -γm,
in which k is a kinetics constant, j is a threshold parameter and γ designs a rate of degradation. The network we consider in this paper concerns four genes of which expressions are regulated at the transcriptional step. The four proteins coded by them are represented by the variables (A i ) i=1,••• , [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF] . The regulation functions are of the same kind as in the Goodwin model but they are more complex because there is here a competition among the proteins that block or activate the transcription of the genes associated to them. Here are the (adimensionalized) equations representing their evolution in concentration:

                                                           Ȧ1 = k 1 + k 11 A 1 j 11 2 + k 13 A 3 j 13 2 1 + A 1 j 11 2 + A 2 j 12 2 + A 3 j 13 2 -γ 1 A 1 Ȧ2 = k 21 A 1 j 21 2 1 + A 1 j 21 2 -γ 2 A 2 Ȧ3 = k 3 1 + A 4 j 34 2 -γ 3 A 3 Ȧ4 = k 4 + k 4 A 1 j 41 2 1 + A 1 j 41 2 + A 3 j 43 2 -γ 4 A 4 .
There are 17 parameters in this model: the terms j ik are concentrations assumed to be constant (they represent some thresholds), the quantities k i design kinetics constants and γ i the degradation terms associated to A i . Thus all our parameters are taken in R + . As explained above, the V-system is the matching of two sub-systems that we'll denote V A1,A2 and V A3,A4 which are respectively composed of the two first equations and of the last two ones. The first sub-system V A1,A2 (in which A 3 is this time considered as a parameter) admits a Hopf bifurcation, which creates an oscillating behavior near the associated bifurcation point. The second one V A3,A4 (in which A 1 is seen as a parameter) admits a hysteresis-type dynamics which makes it jump from a stable steady-state to the other one. Let us describe how the genes act on each other. In each equation the linear degradation of A i is counterbalanced by a nonlinear term. In the first equation, the production of A 1 is activated by itself and by A 3 while A 2 inhibites A 1 : A 2 and A 3 are competing with each other on the operator region of the gene encoding protein A 1 . The second equation tells us that A 1 activates the production of A 2 and bounds it from above by k 21 . The third one is the same as in the Goodwin's model, and shows the inhibition of A 3 by protein A 4 .

In the last equation, the same mechanism appears as in the first one: the more A 3 is produced the more A 4 tends to decrease, while A 1 tends to increase the production of A 4 .

Picture 1 shows the interaction graph that sums up the effects (positive or negative) of the proteins on each other (an edge that ends with an arrow shows a positive feedback, while a one that ends with a stroke stands for a negative feedback).

Our goal is to find a set of values of the 17 parameters that we'll denote P 1 for which the associated V-system V 1 exhibits chaos, with the only rule that the chosen values be strictly positive, and to explain the mechanism responsible of this behavior. Let us notice that our result is in coherence with Thomas' conjecture which states that a positive (the bistable switch) and a negative (the Hopf subsystem) loops are necessary conditions for chaos ( [START_REF] Thomas | Deterministic chaos seen in terms of feedback circuits: analysis, synthesis, "labyrinth chaos[END_REF]). Many other conjectures and results linking the properties of the graph interaction of a gene network, with its dynamical behavior have been stated for a few decades (see [START_REF] Kaufman | A new necessary condition on interaction graphs for multistationarity[END_REF], [START_REF] Thomas | On the relation between the logical structure of systems and their ability to generate multiple steady states or sustained oscillations[END_REF], [START_REF] Soulé | Graphic requirements for multistationarity[END_REF] and very recently [START_REF] Pécou | The interaction graph structure of mass-action reaction networks[END_REF] and [START_REF] Richard | Stable periodicity and negative circuits in differential systems[END_REF]).

A 1 A 2 A 3 A 4 1 Figure 1.
The interaction graph of the V-system.

In our setting, a chaotic dynamical system is a system having a strictly positive topological entropy, which guarantees a sensitive dependence on initial conditions. Usually, the method used to find chaos in a particular differential system is the computing of Lyapunov exponents in different directions: having a strictly positive Lyapunov exponent implies the chaotic dynamics of the system considered. Here such an approach fails because the chaos in the V-system is weak, as it is the case in a lot of models representing a physical phenomenon of the nature (see [START_REF] Dyson | Birds and frogs[END_REF], [START_REF] Kaneko | Chaos as a source of complexity and diversity in evolution[END_REF], [START_REF] Klages | Microscopic chaos, fractals and transport in nonequilibrium statistical mechanics[END_REF]): numerical experiments realized by Tyson's team did not get through in this direction. Moreover the main drawback of this numerical approach is that it does not explain where the chaos comes from: here we want to demonstrate the presence of chaos so that our result be applicable in other models of genetic regulation.

We thus aim at proving that, any general differential system obtained as the matching of the two dynamics described above has (under some generic assumptions on the way these are coupled) a chaotic dynamics. In a more mathematical setting, we want to prove that given a generic one-parameter family of smooth vector fields on R n admitting a local Hopf bifurcation, it is possible to impose a hysteresis dynamics on the parameter so as to obtain a (higher-dimensional) chaotic dynamical system. In the article [START_REF] Pécou | Desynchronization of one-parameter families of stable vector fields[END_REF], the following analogous result is proved: for a one-parameter family of smooth ordinary differential equations that all admit a globally stable asymptotic state, it is possible, under very mild conditions, to construct a feedback on the parameter to get a chaotic system. Unlike what is done in this paper, our result deals with the case of an exchange of stability, appearing at the critical bifurcation point. The other and main difference concerns the approaches chosen to prove the results. Indeed, the proof of E. Pécou's theorem relies on the construction of a homoclinic orbit which permits to meet or to reconstruct a well-known chaotic situation, namely those respectively described by the Shilnikov theorem (see [START_REF] Shilnikov | A case of the existence of a denumerable set of periodic motions[END_REF], and [START_REF] Guckenheimer | Nonlinear oscillations, dynamical systems, and bifurcations of vector fields[END_REF]), and the geometric Lorenz attractors ( [START_REF] Guckenheimer | Structural stability of lorenz attractors[END_REF], [START_REF] Guckenheimer | Nonlinear oscillations, dynamical systems, and bifurcations of vector fields[END_REF]). Similarly another paper in which this approach using the Shilnikov theorem is also adopted is the one of L. Chua (see [START_REF] Chua | The double scroll family[END_REF]). It concerns a system with three equations and four parameters modeling the behavior of a class of electronic circuits. The proof of the existence of a homoclinic orbit is very numerical, and specific to the particular class of piecewise linear differential equations composing this model. Here we could also try to apply the Shilnikov theorem thanks to the spiraling motion associated to the Hopf bifurcation, although we cannot directly apply it (see proposition (4.4) of subsection 4.1.c). But this way to do would not be relevant for the concrete models of genetic regulation, because it would require the finding of such an orbit (satisfying in addition the right assumptions concerning the associated stable and unstable manifolds) in those systems. This is usually impossible due to their too high number of equations and parameters. In particular, trying to find such an orbit in our V-system which admits seventeen ones would not be realistic. To sum up, as in [START_REF] Pécou | Desynchronization of one-parameter families of stable vector fields[END_REF] we also want to consider the theoretical aspect of the problem but with a more constructive approach than the one adopted in it and in [START_REF] Chua | The double scroll family[END_REF], that would fit better with the concrete differential models of genetic regulation, especially the V-system.

Our strategy consists in the construction of a Horseshoe map using the singular perturbations theory: indeed we impose a fast hysteresis dynamics on the parameter so that our dynamical system be a slow-fast one. The use of the slow-fast dynamics has the advantage of exhibiting a hierarchy among the parameters ( [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF]): concerning the V-system, it permits us to see which parameters have to be chosen great compared to the other ones. As usual, the Horseshoe-type dynamics is created from a Poincaré section of which return map permits the transition from a time-continuous dynamical system to a discrete one. We refer the reader to [START_REF] Guckenheimer | Nonlinear oscillations, dynamical systems, and bifurcations of vector fields[END_REF], [START_REF] Iooss | Comportement chaotique des systèmes déterministes[END_REF], [START_REF] Wiggins | Global bifurcations and chaos[END_REF] for more details on this classical way to proceed.

The paper is divided in four parts. In the first one (section 2), we recall the Hopf bifurcation and hysteresis phenomena and define precisely the notion of hysteresis curve. The second part (section 3) is the analysis of the V-system in which we exhibit the important parameters responsible of the presence of the two desired dynamics in this system. We also make a first hierarchy among the seventeen ones so as to reduce the dimension to three, according to the singular perturbations theory, in view of the geometrical study that will follow. In section 4 we state and prove the main result on chaos created by perturbing the parameter of a vector field on R n admitting a Hopf bifurcation (Theorem (4.2)). We begin by the two-dimensional case (see subsection 4.1). The mechanism creating chaos relies on a linear model (see 4.1.c) with two spirals having close but distinct centers, for which we construct a one-dimensional Poincaré section which is covered twice. Then we prove that chaos is kept up when the situation is no more linear by constructing a Horseshoe from this one-dimensional section (see 4.1.c). Let us notice that the idea of fitting two spirals well was already pointed out in [START_REF] Lozi | Etude des motifs dans les oscillations de la réaction de Bélousov-Zhabotinsky(I): Existence d'un confineur principal pour un modèle contraint particulier[END_REF] by R. Lozi, even though it was not prove in this article this can lead to chaos. We finish subsection 4.1 by proving a corollary of Theorem (4.2). Next, the general case n ≥ 2 is proved at the end of section 4, under the hypothesis that all the non purely imaginary eigenvalues of the Jacobian at the bifurcation point have a strictly negative real part. The exactly same reasoning as in dimension two can be applied, because this assumption forces the dynamics to locally restrict itself to a manifold of dimension two. Lastly, in the fourth part (section 5) we use three scales of time to be in the situation of section 4. This permits us to see the parameters of our system, from which the chaotic dynamics mostly depends on. We summarize in Theorem (5.1) the conditions on these key parameters, resulting from all our study. Finally, we exhibit a set of parameters P 1 defining a V-system V 1 , for which we numerically prove that the same mechanism as the one highlighted previously takes place.

2. Short recalls on local bifurcations theory and hysteresis-type dynamics. Notations and definitions.

2.1. Hopf bifurcation. Given a differential system representing a physical phenomenon, some qualitative changes may appear in the solutions of this system as we vary the parameters defining it: for instance the emergence of two fixed points or an exchange of stability among them. Bifurcations theory is the study of those sudden changes of phase portrait (see [START_REF] Guckenheimer | Nonlinear oscillations, dynamical systems, and bifurcations of vector fields[END_REF], [START_REF] Wiggins | Global bifurcations and chaos[END_REF]). The number of parameters defining the system is called the codimension of the bifurcation.

All this text deals with a typical local bifurcation of codimension one, that is which appear near a critical point (called the bifurcation point) of a one-parameter differential system for a particular value of the parameter (the bifurcation value): the Hopf bifurcation. It consists in the birth of a periodic orbit near a critical point for which the associated Jacobian has exactly two eigenvalues with zero real part. As the parameter vary near the bifurcation value, this critical point loses its stability and a periodic orbit appears. The precise mathematical assumptions required to get such a bifurcation, and its normal form are written in subsection 4.1.a. Another local bifurcation of codimension one that we encounter with the hysteresis phenomenon is the saddle-node one, which consists in the emergence of two fixed points of distinct stability.

2.2. Hysteresis curve. The nonlinear phenomenon of hysteresis is defined differently depending on the domain in which it appears. There are plenty of, such as electricity, economics, mechanics, physics of materials... In our context (that is the biology of cells), it is an irreversible one that concerns a system of which dynamics passes from one state to another one. As a threshold value is overtaken, the system jumps to the other state without residing between the two one's. Then, varying the parameter in the backward sense makes the system go back to its first state, but after having reached a threshold different from the first one (as seen in figure 2).

For our study we will need to precise mathematically this notion by defining a hysteresis curve (sometimes called switch):

Definition 2.1. By hysteresis we mean a connected curve in R 2 , defined by some equations of the form g (x, y) = 0, where g is a smooth function from R 2 to R, for which there exist two reals x * < x * * such that: (i)For every x < x * or x > x * * , there is only one z in R verifying g (x, z) = 0, and this unique critical point of the equation ẏ = g (x, y) is asymptotically stable.

(ii)For every x * < x < x * * , there are exactly three zeros (x, z i ) i=1,2,3 of g. Among them, two points say z 1 , z 3 are asymptotically stable critical points of ẏ = g (x, y) (and the other one is repulsive).

(iii)For x = x * or x = x * * there are exactly two zeros (x, z i ) i=1,2 of g. One of the z i 's is asymptotically stable and the other one is degenerate.

O x y 1 y 2 x * x * * Figure 1. A hysteresis curve.
Notation 2.2. In the rest of the text we will only use the notations (x * , y * ), (x * * , y * * ) so as to mention the first and second saddle-node bifurcation points of a given hysteresis.

Remark 2.3. Notice that the inferior and superior branches of a hysteresis defined by a function g are no more invariant when we perturb the equation ẏ = g(x, y) by adding an equation of the form ẋ = f(x, y), as we are going to do in the rest.

We will denote by G 2 the set of smooth real-valued functions g on R 2 defining such curves, and by G c,2 the subset of G 2 composed of hysteresis that can in addition be described as the graph of a (cubic-like) function in the variable y. (By a cubic function we mean a third degree polynomial function). We will also use hysteresis having flat inferior and superior branches and a linear unstable one. We will refer to them as hysteresis with flat branches. They are defined by their points (x * , y * ) and (x * * , y * * ), that we'll also call bifurcation points, So a hysteresis can be decomposed in three parts: two stable curves that we call the superior and inferior branches, separated by the unstable branch. If y * , y * * are the degenerate critical points of the vector fields g x * and g x * * let us remark that locally in (x * , y * ) appears a saddle-node bifurcation, since in a neighborhood of it we have a stable critical point and an unstable one that collapse and then disappear, as the parameter x varies. We call (x * , y * ) the first saddle-node bifurcation point. The same thing happens in the second saddle-node bifurcation point (x * * , y * * ). Notation 2.2. In the rest of the text we will only use the notations (x * , y * ), (x * * , y * * ) so as to mention the first and second saddle-node bifurcation points of a given hysteresis.

Remark 2.3. Notice that the inferior and superior branches of a hysteresis defined by a function g are no more invariant when we perturb the equation ẏ = g (x, y) by adding an equation of the form ẋ = f (x, y), as we are going to do in the following.

We will denote by G 2 the set of smooth real-valued functions g on R 2 defining such curves, and by G c,2 the subset of G 2 composed of hysteresis that can in addition be described as the graph of a (cubic-like) function in the variable y. (By a cubic function we mean a third degree polynomial function). We will also use hysteresis having flat inferior and superior branches and a linear unstable one (see figure 3), and will refer to them as hysteresis with flat branches. They are defined by their points (x * , y * ) and O x y 1 y 2

x * x * * (x * * , y * * ), that we'll also call bifurcation points, even though we cannot really say there is a bifurcation in this case. G p,2 will be the set of piecewise linear functions from R 2 to R, defining such curves. To construct an element h of G p,2 , one can for instance take a negative real a and define the real function e by:

∀x ∈ R, e (x) =    y * * if x ≤ a, (y * * -y * ) x/a + y * if x ∈ [a, 0]. y * if x ≥ 0
Then it suffices to bend the graph of e by considering the function h (x, y) = e (xby)y, for a convenient choice of the constant b, and to notice that the three constant branches of the hysteresis defined by h have the desired stability.

3. The V-system as the coupling of an oscillating system with a switch.

In this section, we present the analytical study of the V-system. We exhibit some conditions on the parameters for which its sub-systems V A1,A2 and V A3,A4 present the two desired dynamics.

3.1.

A local Hopf bifurcation in the sub-system V A1,A2 . The sub-system V A1,A2 is the one defined by the equations:

                           Ȧ1 = k 1 + k 11 A 1 j 11 2 + k 13 A 3 j 13 2 1 + A 1 j 11 2 + A 2 j 12 2 + A 3 j 13 2 -γ 1 A 1 Ȧ2 = k 21 A 1 j 21 2 1 + A 1 j 21 2 -γ 2 A 2 ,
where A 3 is considered here as a parameter.

For every A 3 , let F A3 (A 1 , A 2 ) be the vector field in R 2 associated to V A1,A2 . Making Ȧ2 = 0 and replacing the expression of A 2 in the first equation, we remark that the nullcline of this system, that is to say the set {(A 1 , A 2 , A 3 ) : F A3 (A 1 , A 2 ) = (0, 0)}, can be expressed as the graph of a function in the variable A 1 . The Jacobian matrix DF A3 (A 1 , A 2 ) of our vector field is:

           2A 1 (k 11 -γ 1 A 1 ) j 2 11 1 + A 1 j 11 2 + A 2 j 12 2 + A 3 j 13 2 -γ 1 -2γ 1 A 1 A 2 j 2 12 1 + A 1 j 11 2 + A 2 j 12 2 + A 3 j 13 2 2j 2 21 A 1 (k 21 -k 2 ) (j 2 21 + A 2 1 ) 2 -γ 2           
.

Naturally we are interested in the eigenvalues of this matrix for points (A 1 , A 2 , A 3 ) belonging to the nullcline. The calculations are inextricable, but using Mathematica, we find that for the following particular set of parameters:

     k 1 = 0.05, k 11 = 5, k 13 = 2.4, γ 1 = 0.1 j 11 = 2.5, j 12 = 0.5, j 13 = 2 k 2 = 0, k 21 = 0.3, j 21 = 17.5, γ 2 = 0.03
and for the values (A 1 , A 2 , A 3 ) = (16.5139, 4.71033, 3.30896), the jacobian matrix of the system V A1,A2 (associated to this set of parameters) has two pure imaginary eigenvalues, and so this system admits a local Hopf bifurcation.

3.2.

A hysteresis in the sub-system V A3,A4 . Now we study the second sub-system V A3,A4 defined by the two equations:

                     Ȧ3 = k 3 1 + A 4 j 34 2 -γ 3 A 3 Ȧ4 = k 4 + k 4 A 1 j 41 2 1 + A 1 j 41 2 + A 3 j 43 2 -γ 4 A 4 .
Making Ȧ4 = 0 and replacing the new expression of A 4 in the third equation, we get that the nullcline associated to this subsystem is:

              A 1 , A 3 , k 4 γ 4      1 + A 1 j 41 2 1 + A 1 j 41 2 + A 3 j 43 2           : (A 1 , A 3 ) ∈ C A1,A3          ,
where C A1,A3 designates the set:

             (A 1 , A 3 ) ∈ R 2 + : γ 3 A 3 = k 3 1 + A 1 j 41 2 + A 3 j 43 2 2 1 + A 1 j 41 2 + A 3 j 43 2 2 + k 4 2 γ 2 4 j 2 34 1 + A 1 j 41 2 2              .
We claim this last set is, under certain conditions on the parameters a hysteresis. To do this we prove the following property:

Property 3.1. Assuming we have k 4 γ 4 j 34 > 2, there exist two numbers 0 < A * 1 < A * * 1 such that: -∀A 1 < A * 1 or A 1 > A * * 1 , there is exactly one point in C A1,A3 -∀A 1 ∈ (A * 1 , A * * 1 ), there are exactly three points in C A1,A3 -If A 1 = A * 1 or A = A * *
1 , there are exactly two points in C A1,A3 . Proof. To reduce the expressions appeared above, we introduce a new variable:

X = 1 + A 1 j 41 2 .
As we are only interested in the positive values of our variables, we can use X instead of using A 1 in our calculations. Let f X (A 3 ) be the right-hand side of the equation defining the set C A1,A3 . A simple calculus shows that f X is strictly positive on R * + , tends to zero as A 3 tends to infinity, and that we have f X (0) = 0. Thus the graph of f X has a form given by the figure (4):

A 3 f X (A 3 ) O 1 Figure 4. Graph of the function f X .
Now the idea is to study the intersection of the previous graph with the line γ 3 A 3 when γ 3 varies. Let search on positive values A o 3 for which the tangent at the point

(A o 3 , f X (A o 3 )) passes through (0, 0). Such a point A o 3 verifies the equality f X (A o 3 ) A o 3 = f X (A o 3 )
, that is to say:

4 k 4 2 γ 2 4 j 2 34 A o 3 j 43 2 X 2 X + A o 3 j 43 2 2 + k 4 2 γ 2 4 j 2 34 X 2 = X + A o 3 jdc 2 ,
which is equivalent to:

C 0 3 + 3X 2 C 0 1 - k 4 2 γ 2 4 j 2 34 + 3XC 0 2 + X 3 1 + k 4 2 γ 2 4 j 2 34 = 0, where C 0 is equal to A o 3 j 43 2
. And since we have k 4 γ 4 j 34 > 2, there exist two strictly positive solutions (depending on the variable A 1 ) 0 < C 1 < C 2 of this equation, and so two values 0

< A 1 3 < A 2 3 for which the tangent at A i 3 , f X A i 3
passes through the origin. So for every positive A 1 , there exist two thresholds 0 < c 1 (A 1 ) < c 2 (A 1 ) (which are the two slopes of the tangents) such that:

-if

γ 3 < c 1 (A 1 ) or γ 3 > c 2 (A 1 ), there is only one point in C A1,A3 , -if γ 3 = c 1 (A 1 ) or γ 3 = c 2 (A 1
) there are two points in this curve, -if c 1 (A 1 ) < γ 3 < c 2 (A 1 ) there are three points in it. Considering again γ 3 as a fixed parameter, it remains to justify why we can express this last result in terms of the variable A 1 . To do this, we look at the monotony of c i . We have:

c i (A 1 ) = 4k 3 k 4 2 j 43 γ 2 4 j 2 34 X 2 √ C i (X + C i ) 3 1 + k 4 2 γ 2 4 j 2 34 X X + C i 2 2 , O γ 3 c 1 (A 1 ) c 2 (A 1 ) 1 Figure 5. The curve C A1,A3 for a fixed A 1 .
and because C i is a root of the equation written above, so:

(C i + X) 3 = k 2 4 γ 2 4 j 2 34 (3C i -X) X 2 .
From this last equality, we get two informations: the first one is that the quotient C i /X is a constant and thus does not depend on the variable A 1 , the second one is the following expression of the square root √ C i :

C i = γ 4 j 34 k 4 √ 3 C i + X (C i + X) 2 X 2 + k 4 2 γ 2 4 j 2 34 X C i + X .
Substituting in the expression of c i (A 1 ) we obtain the existence of two strictly positive constants (α i ) i=1,2 , such that for every positive A 1 , we have:

c i (A 1 ) = α i 1 + A 1 j 41 2 ,
which implies the strict monotony of the functions c i . Thus the numbers c -1 1 (γ 3 ) and c -1 2 (γ 3 ) are well defined and distinct because c 1 (A 1 ) = c 2 (A 1 ) for any positive A 1 . Setting A * 1 = c -1 1 (γ 3 ), and A * * 1 = c -1 2 (γ 3 ) we get the result. Property 3.2. Assuming the parameters satisfy the condition

k 4 γ 4 j 34 > 2, then C A1,A3 is a hysteresis defined by a function in G c,2 .
Proof. Since the set C A1,A3 can be expressed as the graph of a function in the variable A 3 , and thanks to the property (3.1), we are sure it has the form of a hysteresis defined by an element of the set G c,2 . To see that the stability conditions required in the definition 2.1 are satisfied, it suffices to consider the sign of the

derivatives ∂g ∂A 3 (A 1 , A 3 ) where g (A 1 , A 3 ) = f X (A 3 ) -γ 3 A 3 . For a point (A 1 , A 3 ) in C A1,A3
, this comes to looking at the sign of the following expression:

4k 3 γ 3 2 k 4 2 j 2 43 γ 4 2 j 2 34 1 + A 1 j 41 2 2 A 3 3 1 + A 1 j 41 2 + A 3 j 43 2 3 -γ 3 .
This expression tends to -γ 3 as A 1 tends to infinity: by continuity this proves that the inferior branch of the curve C A1,A3 is stable. Necessarily its other parts also have the desired stabilities.

Lastly, the isocline of V A3,A4 appears as the intersection of the cartesian product of C A1,A3 by R + with the set:

         (A 1 , A 3 , A 4 ) : A 4 = 1 γ 4      k 4 + k 4 A 1 j 41 2 1 + A 1 j 41 2 + A 3 j 43 2               .
This intersection is a smooth curve in the three-dimensional space (A 1 , A 3 , A 4 ) (see the figure 6).

O

A 4 A 1 A 3 Figure 6.
The isocline of the sub-system V A3,A4 .

3.3.

Using two scales of time to reduce the dimension of the system. As explained in the introduction, we aim at exhibiting a mechanism creating chaos in the V-system (and in other ones having the same dynamical features). This needs a geometrical analysis of its phase portrait, which is a difficult task since it is a four-dimensional system.

To reduce the dimension to three, we apply the Tychonoff theorem on slow-fast systems (see [START_REF] Wasow | Asymptotic expansions for ordinary differential equations[END_REF], [START_REF] Lobry | On Tykhonov's theorem for convergence of solutions of slow and fast systems[END_REF]), by considering the fourth equation as a fast equation. Indeed, this last equation can be rewritten under the form:

1 γ 4 Ȧ4 = k 4 γ 4 1 + A 1 j 41 2 1 + A 1 j 41 2 + A 3 j 43 2 -A 4 .
Thus, if the value γ 4 is very great and of the same order of magnitude as k 4 , the Tychonoff theorem tells us that the dynamics of the V-system is approximately the same as the following slow system S :

                                                 A 1 = k 1 + k 11 A 1 j 11 2 + k 13 A 3 j 13 2 1 + A 1 j 11 2 + A 2 j 12 2 + A 3 j 13 2 -γ 1 A 1 A 2 = k 2 + k 21 A 1 j 21 2 1 + A 1 j 21 2 -γ 2 A 2 A 3 = k 3 1 + A 1 j 41 2 + A 3 j 43 2 2 1 + A 1 j 41 2 + A 3 j 43 2 2 + k 4 2 γ 2 4 j 2 34 1 + A 1 j 41 2 2 -γ 3 A 3 ,
where stands for the derivative according to the fast time τ = γ 4 t. The assumptions needed to apply this result are very easy to verify here since the nullcline of the fourth equation in the variable A 4 is the graph of a function in the other variables and thus is a manifold. Moroever the use of two scales of time does not change the structure of the V-system: it only requires to choose the parameter γ 4 very great and k 4 such that the ratio k 4 /γ 4 be greater or equal to one.

Creating chaos from a generic family of vector fields admitting a Hopf bifurcation

From the analysis made in the part 3, we are led to investigate chaos in a general differential system obtained as the coupling of a Hopf bifurcation with a system of which nullcline is a hysteresis curve. This is the aim of part 4.

Let (f µ ) µ∈R be a generic family of smooth vector fields in R n , depending smoothly on the parameter µ, admitting a local Hopf bifurcation in a point x 0 of R n , for the value of parameter µ 0 . For each real µ, the associated vector field has a flow governed by the equation:

ẋ = f (x, µ) ,
where f is the smooth function naturally defined on R n+1 by the equality f (x, µ) = f µ (x). Our goal is to prove that under an assumption of contractility on the eigenvalues of the Jacobian at x 0 , and under a hypothesis on the curve of critical points of f, which is implicitely defined in a neighborhood of x 0 , we can construct a smooth function g such that the extended singularly perturbed system:

ẋ = f (x, µ) μ = g (x, µ) ,
where is a small positive number, is chaotic. 4.1. The result in dimension two. In all the section we deal with the two-dimensional version of the problem (i.e we set n = 2). 4.1.a. The hypotheses. Here we precise the assumptions we make on the function f. H1 There exists a value µ 0 and a critical point x 0 of the field f µ0 such that the Jacobian D x f µ0 x 0 has a pair of pure imaginary complex eigenvalues ±iβ with β > 0.

By the implicit function theorem, the curve of zeros of f is, in a neighborhood of x 0 , µ 0 , the graph of a smooth function in the variable µ, defined in a small open interval U. We denote by φ = (φ 1 , φ 2 ) this function.

H2 Let α (µ) ± iβ (µ) be the eigenvalues of the Jacobian D x f µ (φ (µ)) that are equal to ±iβ in µ 0 . We assume we have α (µ 0 ) = 0, which means that the two eigenvalues cross the pure imaginary axis with a non-zero velocity.

Under these two assumptions, there exists a change of variables for which the Taylor expansion of degree three of f is of the form:

ẋ1 = a 0 µ + b 0 x 2 1 + x 2 2 x 1 -ω + a 1 µ + b 1 x 2 1 + x 2 2 x 2 (1) ẋ2 = ω + a 1 µ + b 1 x 2 1 + x 2 2 x 1 + a 0 µ + b 0 x 2 1 + x 2 2 x 2 . (2)
Generically, the coefficients a 0 , b 0 are non null in the last equations: in this case, restricting the interval U if necessary, the Hopf bifurcation's theorem tells us that all the fixed points φ (µ) with µ smaller than µ 0 are focuses encircled by limit cycles (which disappear in the bifurcation value µ 0 ), and all the fixed points φ (µ) with µ greater than µ 0 are also focuses but with inverse stability. Recall that this Hopf bifurcation is supercritical if α (µ 0 ) < 0, which means that the cycles are stable so that the focuses in the planes µ > µ 0 are stable ones, or subcritical ifnot (in which case the stability is reverse). Remark 4.1. It is well known that Assumptions H1, H2 are not sufficient to get a Hopf bifurcation. This can be seen by looking at the system (written in complex form) ż = (µ + i) z. Our result will also work for such a degenerate case. If the function φ verifies (φ 1 (µ 0 ) , φ 2 (µ 0 )) = (0, 0), then there exists a smooth function h in G 2 , and a non zero real number a, such that for every sufficiently small positive number , the singularly perturbed system:

(3)    ẋ1 = f 1 (x 1 , x 2 , µ) ẋ2 = f 2 (x 1 , x 2 , µ) μ = h (x 1 + ax 2 , µ)
, taken in a neighborhood V of x 0 , µ 0 enough small, is chaotic.

Notice that the hypothesis on the curve of critical points of f is generic. In particular, it is satisfied by differential models of gene regulatory networks presenting a Hopf bifurcation, notably the V-system. Moreover this condition is the same as the one asked by E.Pécou, in her article [START_REF] Pécou | Desynchronization of one-parameter families of stable vector fields[END_REF]. If one wants to have a hysteresis that can be described as a cubic-like function in the variable µ, the analogous result is the following: Corollary 4.3. With the same assumptions as in Theorem (4.2), there exists a smooth function h in G c,2 , a non zero real a, and two small numbers 0 < 1 < 2 , such that for every in ] 1 , 2 [, the system (3), taken in a neighborhood V of x 0 , µ 0 enough small, is chaotic.

The proof of this corollary is the same as the one of (4.2), except for its last step (see subsection 4.1.d).

4.1.c. Proof of Theorem (4.2). All along the proof, we assume that we have α (µ 0 ) < 0 (in case of a Hopf bifurcation at µ 0 this means the cycles are stable ones). The case where α (µ 0 ) > 0 is totally similar (see the remark at the end of the subsection).

The idea of the demonstration is to take a function h defining a hysteresis curve of which stable part intersects the curve of fixed points of f in only one point, say the origin O that belongs to the Hopf bifurcation surface, so as to use the spiraling motion due to this bifurcation. There are infinitely many such smooth functions h: let's take any one in G c,2 , and denote by F the vector field defining system (3). Then, the Jacobian at the origin is:

DF (O) =            ∂f 1 ∂µ (O) Df(O) ∂f 2 ∂µ (O) 1 ∂h (0, 0) ∂x 1 a ∂h (0, 0) ∂x 1 1 ∂h (0, 0) ∂µ           
.

Because of the assumption on the curve of critical points of f, the terms on the top right of this matrix are non null, and so all the possible cases can happen concerning its spectrum when vary in ]0, 1[: it can, for instance, be only composed of real eigenvalues, in which case in the neighborhood of the origin the spiraling motion would not exist. So we cannot directly work with any smooth function h in G c,2 . On the other hand, if the function h is in G p,2 then the terms on the bottom left are equal to zero, and the dynamics locally in the origin is the one we want to have. More than that, the following proposition (of which proof is just an application of the implicit function theorem) tells us that given a very small , if the stable branches of the hysteresis have a slope of order , then the dynamics is the same as in this flat case. Proposition 4.4. Let B be a square matrix of size two, of which spectrum is (α ± iβ), with α, β > 0, and let A ( ) be the matrix of size three:

A ( ) =     a 1,3 B a 2,3 a 3,1 aa 3,1 -γ     ,
where γ is a strictly positive number. Then for any enough small, the spectrum S (A ( )) of A ( ) has the form: S (A ( )) = {α + z 1 ( ) ± i (β + z 2 ( )) , -γ/ + z 3 ( )}, where the z i ( ) are bounded functions of .

Thus we adopt the following strategy: we first describe the hysteresis curves with flat branches that are conveniently placed comparing with the curve of critical points φ in order to prove Theorem (4.2) in the (weaker) case when the function h is in G p,2 . And then we will transform the convenient function h into a smooth one.

Let us fix a value µ * * in U greater than µ 0 . Because of the generic assumption there exists a value µ * for which we have φ (µ * ) = φ (µ * * ). Without loss of generality, we suppose that the point (φ (µ * ) , µ * ) is the origin O. We want to work with a function g of which set of zeros is the cartesian product of the straight line (Oφ (µ * * )) with a hysteresis in the plane orthogonal to this line. To do this we begin by choosing two values

x * 1 , x * * 1 verifying x * 1 < x * * 1 < 0 if φ 2 (µ * * ) > 0 (or else x * 1 > x * * 1 > 0 if φ 2 (µ * * ) < 0)
and such that the hysteresis with flat branches defined by the two bifurcation points (x * 1 , 0) and (x * * 1 , µ * * ) does not intersect the graph of φ in another point than the origin. Besides we ask that the reals x * 1 , x * * 1 are close enough to zero, so that the points (x * 1 , 0) and (x * * 1 , 0) are inside the possible cycle C 0 belonging to the plane µ = 0 (see the figure 7).

In the rest of the text we will assume (without loss of generality) that φ 2 (µ * * ) > 0, and thus the values x * 1 , x * * 1 will be taken negative. Now that we have defined our function h, we set g (x 1 , x 2 , µ) = h (x 1 + ax 2 , µ), where x 1 + ax 2 = 0 is an equation of the straight line (Oφ (µ * * )). Our first goal is to prove that for a certain choice of the values x * 1 , x * * 1 and µ * * , such a function g verifies that the system:

   ẋ1 = f 1 (x 1 , x 2 , µ) ẋ2 = f 2 (x 1 , x 2 , µ) μ = g (x 1 , x 2 , µ) , ( 4
) 6 CAMILLE POIGNARD AND ELISABETH P ÉCOU µ x 1 x 2 (φ(µ * * ), µ * * ) (φ(µ 0 ), µ 0 ) O 1 x 1 x 2 C 0 φ(µ * * ) O x1 + a x2 = 0 Figure 2.
The Hopf bell-like surface and the surface defined by a piecewise-linear hysteresis.

2.3.a.

Using the singular perturbations theory to describe the flow of (2). Let's take a piecewise linear function g from R 3 to R of the same kind as above. By construction, the system (2) admits the origin as unique critical point. Remark that since the function x 1 → h (x 1 , 0) is constant, there exists locally in the origin, an invariant stable manifold of dimension one and an unstable one of dimension two. In fact, we even have that the inferior and superior half planes of the hysteresis, that are the cartesian products of the stable branches of h by (Oφ (µ * * )), are invariant except near the two straight lines ∆ * and ∆ * * (that we'll call the two fold lines of the hysteresis), respectively defined by the equations x 1 + ax 2 = x * 1 and x 1 + ax 2 = x * * 1 . This is clear by the Cauchy-Lipschitz theorem, which can be applied here because the function h is Lipschitz continuous (indeed, the function e defined at the end of 2.1.b is Lipschitz continuous). Using the singular perturbations theory to describe the flow of (4). Let's take a piecewise linear function g from R 3 to R of the same kind as above. By construction, the system (4) admits the origin as unique critical point. Remark that since the function x 1 → h (x 1 , 0) is constant, there exists locally in the origin, an invariant stable manifold of dimension one and an unstable one of dimension two. In fact, we even have that the inferior and superior half planes of the hysteresis, that are the cartesian products of the stable branches of h by the line (Oφ (µ * * )), are invariant except near the two straight lines ∆ * and ∆ * * (that we call the two fold lines of the hysteresis), respectively defined by the equations x 1 + ax 2 = x * 1 and x 1 + ax 2 = x * * 1 . This is clear by the Cauchy-Lipschitz theorem, which can be applied here because the function h is Lipschitz continuous (indeed, the function e defined at the end of 2.1 is Lipschitz continuous). Moroever, these invariant half planes can be described as (constant) graphs in the variables x 1 , x 2 : namely the graphs ξ -(x 1 , x 2 ) = 0 for (x 1 , x 2 ) above ∆ * and ξ + (x 1 , x 2 ) = µ * * for (x 1 , x 2 ) below ∆ * * . This fact allows us to use singular perturbation theory: Proposition 4.5. As > 0 tends to zero, the flow ϕ 2, of (4) is C 0 -approached by a flow ϕ 2 of which trajectories are successions of continuous arcs, each of them being the union of a segment of the form x 0 1 , x 0 2 , µ : 0 ≤ µ ≤ µ * * with a solution of one of the equations ẋ = f (x, ξ + (x)) and ẋ = f (x, ξ -(x)).

More precisely, we have:

∀M > 0, ∀p ∈ V, ∀t ∈ [0, M ], lim →0 ϕ 2, (p, t) = ϕ 2 (p, t) ,
where V is a neighborhood of x 0 , µ 0 enough small.

Proof. It suffices to apply again the Tychonoff theorem on slow-fast systems. The idea is the following: for > 0 enough small, any point which is not a zero of our function g will be carried vertically (that is to say along the µ-axis) by the flow ϕ 2, until it reaches a stable part of the hysteresis, in which case its motion will be defined by the flow reduced on this surface.

To precise this idea, let ϕ + be the reduced flow associated to the slow equation ẋ = f (x, ξ + (x)). By construction of our hysteresis defined by g, its critical point φ (µ * * ) (which is unique in a neighborhood of x 0 , µ 0 enough small) is stable. So, for initial conditions enough close to this point, the trajectories will hit the fold line ∆ * * in a finite time, afterwards they will not exist anymore. Moroever, asking that the value x * * 1 be closer to zero if necessary, we get that these trajectories are almost logarithmic spirals defined by a polar equation of the form:

ρ = ρ 0 e α (µ * * ) β (µ * * ) (θ-θ0)
, in the sense that locally in φ (µ * * ), the flow ϕ + is close to its linear part.

The same holds for the reduced flow ϕ -of the other slow equation ẋ = f (x, ξ -(x)). In this case, the trajectories are nearby repulsive logarithmic spirals, having a polar equation defined by the positive coefficient α (0) /β (0) (see figure 8). Now, let us denote by ϕ 2 the continuous flow, of which trajectories are those of the reduced flows ϕ + , ϕ - connected between them by vertical segments of the form

x * 1 , x 0 2 , µ : 0 ≤ µ ≤ µ * * or x * * 1 , x 0 2 , µ : 0 ≤ µ ≤ µ * * .
Tychonoff's theorem tells us that ϕ 2 is the limit of ϕ 2, as tends to zero, in the meaning of proposition (4.5). Definition 4.6. The flow ϕ 2 of proposition (4.5) is called the limit flow associated to the system (4).

A linear model creating chaos. Here we simplify again our problem by considering the case where the reduced flows ϕ -and ϕ + (of which orbits are repulsive spirals centered at the origin, and attractive ones centered at φ (µ * * )), are linear. The linear limit flow of proposition (4.5) obtained by this way is denoted by ϕ 3 . We claim we have:

x 1 x 2 ϕ0 ϕ1 O ∆ * ∆ * * x1 + ax2 = 0 1 Figure 9. The flow ϕ 3 .
Proposition 4.7. There exists a choice of the values x * 1 < x * * 1 < 0 and µ * * > 0 such that the linear limit flow ϕ 3 associated to the system (4) covers at least twice a segment I belonging to the plane µ = 0. More precisely, there exists a decomposition I = I 1 ∪ I 2 in two sub-intervals and a Poincaré return map P (associated to ϕ 3 ) defined on I, such that the images P 2 (I 1 ), P 2 (I 2 ) strictly contain I.

Proof. Let F 1 , F 2 be the linear flows defined in the plane µ = 0, of which orbits are respectively centered at the origin (0, 0) and at the point φ (µ * * ), and have the polar equations:

ρ = ρ 0 e α (0) β (0) (θ-θ0) and ρ = ρ 0 e α (µ * * ) β (µ * * ) (θ-θ0) .
We use these two flows so as to fix the position of the two fold lines ∆ * , ∆ * * of our hysteresis h (i.e to fix the values x * 1 , x * * 1 ). Recall that their slopes is -1/a. (i)We begin with the position of the fold ∆ * . Given an initial condition M 0 = (ρ 0 , θ 0 ) near the origin, any point M = (ρ, θ) belonging to the trajectory (F 1 (M 0 , t)) t∈R satisfies:

M =   ρ0e α (0) β (0) (θ-θ0) cos (θ) , ρ 0 e α (0) β (0) (θ-θ0) sin (θ)    ,
in cartesian coordinates. Derivating this expression in θ, we get that the locus of points at which the tangent of the flow F 1 has a slope equals to -1/a, is a straight line, having a constant angle with the half-axis Ox 1 equals to arctan

    α (0) β (0) a + 1 1 a - α (0) β (0)     ,
or π/2 in the case where a = β (0) /α (0) (see figure 10).

Denote by A the intersection of this straight line with one having an equation of the form x 1 + ax 2 = b where b is strictly negative. Call ∆ this last line, and consider the first return of the point A in ∆, that is to say the point B = F 1 (A, τ (A)), where τ (A) is the first strictly positive time necessary for A to hit this line. We obtain a straight line (OB) which is clearly the locus of first return of the points, in which the slope of the tangent to F 1 is -1/a, in this tangent. To see this, it suffices to observe that the image of a logarithmic spiral under a homotethy is still a spiral of same nature, defined by the same coefficient (here α (0) /β (0)). Remark that (OB) has a slope strictly smaller than a. From the construction of these two lines (OA) and (OB), we adjust the position of ∆ * by taking the circle centered at φ (µ * * ) of ray Oφ (µ * * ): it cuts (OB) in a point at which we draw ∆ * (figure 10). necessary for the point M to reach ∆ * * ), we get that the image J

= [A 1 B 1 ] of J under Ψ 1 is a segment of which extremity B
1 has a strictly positive x 2 -coordinate. Moreover, fixing the value x * * 1 so that the origin be enough close to ∆ * * , we obtain that J intersects the spiral S 0 in at least 3 points having a strictly negative x 2coordinate, and such that the extremity A 1 = Ψ 1 (A 1 ) does not belong to the spiral S 1 (see figure 6). (iii) Now that we have fixed our hysteresis h, let us prove that ϕ 3 is chaotic. By (ii) the segment J contains three points M 1 , M 2 , M 3 of the spiral S 0 (all having a strictly negative x 2 -coordinate) which are consecutive (i.e there exist two times t 1 , t 2 > 0, such that ϕ 3 (t 1 , M 3 ) = M 2 and ϕ 3 (t 2 , M 2 ) = M 1 ). Then let us consider the first return map Ψ in the axis (Ox 1 ) associated to the flow ϕ 3 (or equivalently F 1 ), and set:

I 1 = [Ψ (M 1 ) Ψ (M 2 )] I 2 = [Ψ (M 2 ) Ψ (M 3 )]
.

By construction we have: (ii) Then, we consider the flow F 2 , so as to position ∆ * * . Let us denote by J = [A 1 B 1 ] the segment of which extremities are the intersections of ∆ * with our two straight lines (OA) and (OB). We know there is a unique trajectory S 0 of the flow F 1 which is tangent to ∆ * , and that the tangency point is the boundary A 1 of this segment. The return of J in the line (Oφ (µ * * )) by rotations centered at φ (µ * * ) is a segment of the form [OC] with C having a strictly non zero x 2 -coordinate. As the value µ * * ≥ µ 0 was chosen arbitrarily in the construction of our hysteresis h, we can take it closer to µ 0 if necessary, so that we can make the spiraling motion of the flow F 2 be very close to a rotation one. Thus denoting by Ψ 1 the application of first return in ∆ * * (i.e defined by the equality Ψ 1 (M ) = F 2 (M, τ (M )), where τ (M ) is the time necessary for the point M to reach ∆ * * ), we get that the image J = [A 1 B 1 ] of J under Ψ 1 is a segment of which extremity B 1 has a strictly positive x 2 -coordinate. Moreover, fixing the value x * * 1 so that the origin be enough close to ∆ * * , we obtain that J intersects the spiral S 0 in at least 3 points having a strictly negative x 2 -coordinate, and such that the extremity A 1 = Ψ 1 (A 1 ) does not belong to the spiral S 1 (see figure 11).

Ψ(J ) = Ψ • Ψ 1 • Ψ 0 (I 1 ) I 1 ∪ I 2 Ψ • Ψ 1 • Ψ 0 (I 1 ) = Ψ • Ψ 1 • Ψ 0 (I 2 ) , which means that the map Ψ • Ψ 1 • Ψ 0 covers the interval I = I 1 ∪ I 2 twice.
(iii) Now that we have fixed our hysteresis h, let us prove that ϕ 3 is chaotic. By (ii) the segment J contains three points M 1 , M 2 , M 3 of the spiral S 0 (all having a strictly negative x 2 -coordinate) which are consecutive (i.e there exist two times t 1 , t 2 > 0, such that ϕ 3 (t 1 , M 3 ) = M 2 and ϕ 3 (t 2 , M 2 ) = M 1 ). Then let us consider the first return map Ψ in the axis (Ox 1 ) associated to the flow ϕ 3 (or equivalently F 1 ), and set:

I 1 = [Ψ (M 1 ) Ψ (M 2 )] I 2 = [Ψ (M 2 ) Ψ (M 3 )]
.

By construction we have:

Ψ(J ) = Ψ • Ψ 1 • Ψ 0 (I 1 ) I 1 ∪ I 2 Ψ • Ψ 1 • Ψ 0 (I 1 ) = Ψ • Ψ 1 • Ψ 0 (I 2 ) ,
which means that the map Ψ • Ψ 1 • Ψ 0 covers the interval I = I 1 ∪ I 2 twice.

CHAOS OBTAINED FROM A HOPF BIFURCATION 11 From all of this we conclude that the segment I is a Poincaré section for the linear limit flow ϕ 3 . Let P be the return map of ϕ 3 associated to this section. This From all of this we conclude that the segment I is a Poincaré section for the linear limit flow ϕ 3 . Let P be the return map of ϕ 3 associated to this section. This map does not cover I twice because we have P (I 2 ) = I 1 , but the map P 2 does (because P (I 1 ) = Ψ • Ψ 1 • Ψ 0 (I 1 ) contains I 2 ), which finishes the proof.

A1 B1 A 1 B 1 x 1 x 2 φ( µ * * ) S 0 J J ∆ * ∆ * * 1 M1 M2 M3 A1 B1 x 1 x 2 I 1 I 2 ∆ * * x1 + ax2 = 0
End of the proof of the result. Finally, we can transform the piecewise linear hysteresis h into a smooth one, by smoothing the rough edges near the points (x * 1 , µ * ) and (x * * 1 , µ * * ): this does not affect the dynamics of our system since the images (under the flow (φ t ) t∈R ) of the invariant Cantor set included in R (s 0 ) are all at a distance strictly positive of the flat branches, in particular of the rough edges of the hysteresis. The result is proved in the case α (µ 0 ) < 0. Remark 4.10. In the case where α (µ 0 ) > 0 (which corresponds to a subcritical Hopf bifurcation), it suffices to take the symmetric of the hysteresis we considered above with regard to the axis Ox 1 . The exactly same reasoning applies in this case. The proof is achieved.

4.1.d. Proof of Corollary (4.3). Let us consider again the set R (s 0 ) defined above. There exist two small numbers 0 < 1 < 2 such that for any in ] 1 , 2 [, both P 2 (R 1 (s 0 )) and P 2 (R 2 (s 0 )) will cover this set while being at a bounded distance of the plane µ = 0. We can therefore bend the inferior stable branch of the piecewise linear hysteresis h without modifying the dynamics established above (see Proposition (4.4) at the end of the remark 4.10 below), and also smooth the rough edge near the point (x * 1 , µ * ). We obtain an inferior branch that can be described as the graph of a smooth function in the variable µ. Applying the same for the superior branch of h, we get the result.

µ = 0 Γ(s 0 ) P 2 (R 1 (s 0 )) P 2 (R 2 (s 0 ))
Figure 12. The covering of the rectangle R (s 0 ) by map P 2 (for a fixed in ] 1 , 2 [).

4.2.

The result in the general case n ≥ 2. Theorem (4.2) can be extended to the case where (f µ ) µ∈R is a family of vector fields in R n , with n ≥ 2, provided we make an additional hypothesis of contractility on the other eigenvalues of the Jacobian at the bifurcation point. 4.2.a. The hypotheses. H1 There exists a value µ 0 and a critical point x 0 of the field f µ0 such that the Jacobian D x f µ0 x 0 has a pair of pure imaginary complex eigenvalues ±iβ with β > 0, and the other eigenvalues λ 3 , • • • , λ n have a strictly negative real part.

Here again, we denote by φ = (φ 1 , • • • , φ n ) the smooth function defined on a neighborhood U of the bifurcation value µ 0 , that locally defines the graph of fixed points associated to the function f. Hypothesis H2 on the cross of the pure imaginary axis with a non-zero velocity, has not changed. Remark that restricting U if necessary, we have that the eigenvalues λ 2 (µ) , • • • , λ n (µ) (which are equal to λ 2 , • • • , λ n in µ = µ 0 ) have also a strictly negative real part. For each µ in U, we denote by Π µ the plane span by the real and imaginary parts of the eigenvectors associated to α (µ) ± iβ (µ). Note that all those planes are almost parallel because these eigenvectors vary smoothly with the parameter µ.

Under these two assumptions, the center manifold theorem gives us the existence, for every µ in U, of a smooth manifold of dimension two W µ , tangent at the point (φ (µ) , µ) to the plane Π µ . This manifold is attracting for the flow of the equation ẋ = f (x, µ). Moroever, for every µ in U, the dynamics of this system restricted to W µ is given by equations of the form:

(5) u1 = a 0 µ + b 0 u 2 1 + u 2 2 u 1 -ω + a 1 µ + b 1 u 2 1 + u 2 2 u 2 u2 = ω + a 1 µ + b 0 u 2 1 + u 2 2 u 1 + a 0 µ + b 0 u 2 1 + u 2 2 u 2 .
4.2.b. The theorem and its proof. The n-dimensional version of our result is the following: Theorem 4.11. Let us suppose that Hypotheses H1 , H2 on the function f are satisfied. Assume moreover that the projection of the graph x = (φ 1 (µ) , • • • , φ n (µ)) on the plane Π µ0 is not reduced to the point (φ (µ 0 ) , µ 0 ). Then, there exist two distinct integers i, j in [1, n], two non zero real numbers a i , a j and a smooth function h in G 2 such that for every > 0 enough small, the singularly perturbed system:

(6) ẋ = f (x, µ) μ = h (a i x i + a j x j , µ) ,
taken in a neighborhood of x 0 , µ 0 enough small, is chaotic.

Proof. The exactly same reasoning as in the two-dimensional case can be applied.

(i)Indeed, let us first fix a value µ * * in U, strictly greater than µ 0 . By assumption on the graph x = φ (µ), there exists a value of the parameter (say 0) which is strictly smaller than µ 0 and such that the projection X 0 of the point (φ (µ * * ) , 0) on the plane Π 0 is not equal to the fixed point (φ (0) , 0). Without loss of generality, we assume that this last point is the origin O = (0,

• • • , 0) of R n+1 .
Then we consider two planes H * , H * * which are parallel to the plane containing the origin and the two points (φ (µ * * ) , 0) , X 0 . There exist two distinct integers i, j in [1, n], such that the equations of these planes are of the form:

a i x i + a j x j = x * a i x i + a j x j = x * * ,
where a i , a j and x * , x * * are four non zero numbers. The intersections H * ∩ Π 0 and H * * ∩ Π 0 are two lines that will be our two fold lines ∆ * , ∆ * * of subsection 4.1.c. To finish this first step we consider the element h in G p,2 defining the hysteresis with flat branches of which bifurcation points are (x * , 0) and (x * * , µ * * ), in order to prove the piecewise linear version of the result before smoothing this convenient hysteresis.

(ii) Without loss of generality, we can assume that the origin is very close to the two fold lines that intersect by the way a small neighborhood of the origin in R n+1 , in which the flow of the system (6) is very close to its linear part. Thanks to the hypothesis H1 , we can therefore construct a rectangular section I in the hyperplane µ = 0, that transversely cross the plane Π 0 , and which is covered twice by the limit flow associated to [START_REF] Bryan | Oscillatory behavior in enzymatic control processes[END_REF]. It remains to extend vertically this section I, that is to say to consider the sets

R (s) = {(x, µ) ∈ I × [0, s]} ,
for some small strictly positive s. There exist a value s 0 such that R (s 0 ) is a Horseshoe for the flow of the system (6). We conclude by taking off this rectangle from the hyperplane µ = 0, and by smoothing the hysteresis as in the two-dimensional case.

To finish section 4.2, we also have the similar corollary as in section 3:

Corollary 4.12. With the same assumptions as in Theorem (4.11), there exist two integers i, j in [1, n], two non zero real numbers a i , a j , a smooth function h in G c,2 and two small numbers 0 < 1 < 2 , such that for every in ] 1 , 2 [, the system (6), taken in a neighborhood of x 0 , µ 0 enough small, is chaotic.

5.

Numerical evidence of a chaotic motion in the V-system Finally, in order to apply Theorem (4.2) to find chaos in our V-system, it suffices to use three scales of time so that the slow system S written in section 3.3 be in fact a slow-fast one. This comes to considering the V-system under the form:

                                                           Ȧ1 = k 1 + k 11 A 1 j 11 2 + k 13 A 3 j 13 2 1 + A 1 j 11 2 + A 2 j 12 2 + A 3 j 13 2 -γ 1 A 1 Ȧ2 = k 2 + k 21 A 1 j 21 2 1 + A 1 j 21 2 -γ 2 A 2 1 γ 3 Ȧ3 = k 3 γ 3 1 1 + A 4 j 34 2 -A 3 1 γ 4 Ȧ4 = k 4 γ 4 1 + A 1 j 41 2 1 + A 1 j 41 2 + A 3 j 43 2 -A 4 ,
in which γ 3 is chosen very great but much less than γ 4 , and the ratios k 3 /γ 3 , k 4 /γ 4 are greater or equal to one. From this use of three scales of time and the study done in section 3, we are finally able to exhibit the key parameters (among the 17 ones) of which the chaotic dynamics in the V-system strongly depends on, by classifying all the parameters in three groups. The first group contains those associated to V A1,A2 (for which this sub-system must have a Hopf bifurcation), the second one is the subset {j 34 , k 3 , k 4 , γ 3 , γ 4 }, and the last one is composed of the parameters that remain. Recall that, the second group must satisfy the assumption of property 3.1 to have a hysteresis, besides the inequalities we have just mentioned above. We thus can enunciate the following:

Theorem 5.1. There exist some sets of parameters P, satisfying the conditions:

-the sub-systems V A1,A2 admit a Hopf bifurcation, -the parameters j 34 , k 3 , k 4 , γ 3 , γ 4 verify the inequalities:

1 γ 3 γ 4 , 2 < k 4 γ 4 j 34 , 1 ≤ k 3 γ 3 , 1 ≤ k 4 γ 4 ,
for which the associated V-systems V are chaotic.

We cannot say the two conditions of Theorem 5.1 are sufficient to get chaos, since we also need that the Hopf bifurcation surface be well positioned compared with the hysteresis surface (conformed with figure [START_REF] Guckenheimer | Nonlinear oscillations, dynamical systems, and bifurcations of vector fields[END_REF]) and that this last one be enough flat. But this condition, involving notably the parameters of the third group, cannot be quantified. From this hierarchy among the parameters, we are led to the following set of parameters .

The figure [START_REF] Lobry | On Tykhonov's theorem for convergence of solutions of slow and fast systems[END_REF] shows the bifurcation surfaces of the slow system S 1 , associated to the V-system V 1 . Numerical investigations realized with the software Xdim permitted us to find a Poincaré map that covers twice the section on which it is defined, in an exactly similar way as planned (see figures [START_REF] Lozi | Etude des motifs dans les oscillations de la réaction de Bélousov-Zhabotinsky(I): Existence d'un confineur principal pour un modèle contraint particulier[END_REF], [START_REF] Lu | Next-generation synthetic gene networks[END_REF], and ( 16)). We thus can conclude that V 1 is chaotic. is a hysteresis. The figure [START_REF] Kaufman | A new necessary condition on interaction graphs for multistationarity[END_REF] shows the bifurcation surfaces associated to the slow system S 1 .

Figure 11. The hysteresis and the curve of critical points associated to the Hopf bifurcation of the system S 1

With Xdim we found a Poincaré map that covers twice the section on which it is defined, in an exactly similar way as planned (see figures [START_REF] Klages | Microscopic chaos, fractals and transport in nonequilibrium statistical mechanics[END_REF] and ( 13)).

This convinces us that S 1 is chaotic. Such a section can be found for the Virginiator V 1 , by multiplying the parameters k sd , k sda , k dd by a very great number.

Conclusion

We have proved that there exist some values of the parameters for which the V-system admits chaos, and highlighted the key parameters responsible in a major part of this dynamics (see Theorem 5.1). This has been done by proving the existence of a horseshoe, which guarantees mathematically the robustness of chaos in our system: that is to say, there exists a small number η 0 > 0 such that for every 0 < η < η 0 , every vector field of the form F + ηG (where F denotes the one defining the V-system) is chaotic near the Hopf bifurcation point. This robustness is called structural stability in mathematics (see [START_REF] Guckenheimer | Nonlinear oscillations, dynamical systems, and bifurcations of vector fields[END_REF], [START_REF] Iooss | Comportement chaotique des systèmes déterministes[END_REF]), and must not be confused with the sensitive dependence on initial conditions of the system (indeed given a chaotic system, the trajectories are very sensitive to perturbations, but the system itself can still be chaotic under some small changes). In practice even a slight modification of the parameters in the V-system makes us leave the region on which we have the desired dynamics: as said in the introduction the chaos in it is weak even if it is theoretically robust. In this fact lies all the interest of the paper, for a purely numerical study would have attributed the variations in this system to noise.

Since a lot of gene networks present the same characteristics as the V-system, namely a positive feedback (induced here by a Hopf bifurcation) and a negative one (induced by a switch), we can hope to see chaotic variations in the dynamics of other such networks. Besides the interest mentioned in the introduction in controlling, modifying the activity of a given cell, the investigation of chaotic G.R.N.s, presenting a stronger and more visible chaos (for instance having strongly positive Lyapunov exponents) is also relevant because of the possibility of chaos as a source of diversity in the behavior of cells (see [START_REF] Munteanu | Phenotypic diversity and chaos in a minimal cell model[END_REF]). For instance, recent studies done on stem cells have pointed out that cells presenting chaotic oscillations in gene expression dynamics have the ability to differentiate into other cell types, leading to the hypothesis that chaos could be involved in differentiation (see [START_REF] Furusawa | Chaotic expression dynamics implies pluripotency: when theory and experiment meet[END_REF]). Of course, the variations in gene expression must be in a great part due to randomness but the idea that deterministic nonlinear effects (and in particular chaos) are also involved in them have been established experimentally in some recent cases (see [START_REF] Shimojo | Oscillations in notch signaling regulate maintenance of neural progenitors[END_REF] and [START_REF] Becks | Experimental demonstration of chaos in a microbial food web[END_REF]). These advances in favor of a chaotic gene expression have been done thanks to progress in experimental techniques, which may permit to corroborate or unvalidate in the future those new ideas.

Figure 2 .

 2 Figure 2. A hysteresis curve.

Figure 3 .

 3 Figure 3. A hysteresis with flat branches.

4. 1 Theorem 4 . 2 .

 142 .b. The Theorem. The two-dimensional version of our result is the following: Let us assume Hypotheses H1, H2 on f are satisfied.

Figure 7 .

 7 Figure 7. The Hopf bell-like surface and the surface defined by a piecewise-linear hysteresis.

Figure 3 .x1 + a x2 = 0 Figure 8 .

 308 Figure 3. The dynamics of the flow ϕ 2 .

1 Figure 5 .

 15 Figure 5. Construction of the fold line ∆ * .

Figure 10 .

 10 Figure 10. Construction of fold line ∆ * .

Figure 6 .

 6 Figure 6. The one-dimensional section I covered twice by the flow ϕ 3 .

Figure 11 .

 11 Figure 11. The one-dimensional section I covered twice by the flow ϕ 3 .

k 1

 1 = 0.05, k 11 = 5, k 13 = 2.4, γ 1 = 0.1 j 11 = 2.5, j 12 = 0.5, j 13 = 2 k 2 = 0, k 21 = 0.3, j 21 = 17.5, γ 2 = 0.03 k 3 = 330, γ 3 = 32, j 34 = 1.5 k 4 = 387, j 41 = 2.9119, j 43 = 0.75, γ 4 = 100.1

  cannot really increase the values because of numerical instabilities appeared in the software Xdim for such high values. Remark that the parameters of the system V 1A,B and V 0 A,B are the same, thus V 1 A,B still admits a Hopf bifurcation. Moroever we have k sd k dd j cd = 2.57742 and k sd = k sda thus by the proposition 3.2 the set C 1 A,C

Figure 13 .

 13 Figure 13. The hysteresis and the curve of critical points associated to the Hopf bifurcation of the system S 1 .

Figure 14 .

 14 Figure 14. A trajectory of the system V 1 .

Figure 15 .

 15 Figure 15. A Poincaré map for the flow of V 1 that covers twice the segment on which it is defined.

CAMILLE POIGNARD AND ELISABETH P ÉCOU

Acknowledgements

The author thanks Elisabeth Pécou for her help and her advices, especially on the biological aspects of this work.

Our application P of the previous proposition is piecewise continuous, because the boundary points of the segments I 1 , I 2 (denoted Ψ (M 1 ) , Ψ (M 2 ) , Ψ (M 3 ) in the previous proof) are discontinuous ones. It is well known that any piecewise continuous application of an interval covering at least twice this interval has strictly positive topological entropy. Thus, we have proved: Corollary 4.8. There exists a choice of the three values x * 1 < x * * 1 < 0 and µ * * > 0 such that the linear limit flow ϕ 3 associated to the system (4) is chaotic.

Proof of the result when the hysteresis is piecewise linear with flat branches. Now let us consider again the limit flow ϕ 2 of the system (4), but this time without assuming that ϕ + and ϕ -are linear. In this case we still have a Poincaré return map P associated to ϕ 2 and defined on a segment I = I 1 ∪ I 2 belonging to the plane µ = 0, such that the images P 2 (I 1 ) and P 2 (I 2 ) strictly contain I.

Indeed, making a Taylor development of the function x → f (x, 0) at the origin (0, 0), we get by the implicit function theorem, that the locus of points at which the tangent of the flow ϕ -is parallel to (Oφ (µ * * )) is (locally at (0, 0)) a curve which is tangent to the straight line (OA) of the construction we made above. Similarly the points of first return in these tangents form now a curve tangent (at the origin (0, 0)) to the straight line (OB) we had before. So the circle of ray (Oφ (µ * * )) still cuts this curve in a point at which we can draw the line ∆ * . Then we can go on the same construction as in the linear model.

Then, from this segment I constructed above, we can create a rectangular section transverse to the plane µ = 0, admitting a Horseshoe: Proposition 4.9. For the same choice of the values x * 1 , x * * 1 , µ * * taken in proposition (4.7), there exists a rectangle R transverse to the plane µ = 0, such that for any > 0 enough small, the flow ϕ 2, of the system (4) covers at least twice this section. More precisely, there exists a decomposition R = R 1 ∪ R 2 in two rectangles such that for any > 0 enough small, the flow ϕ 2, admits a Poincaré return map P defined on R, such that P 2 (R 1 ) , P 2 (R 2 ) are disjoint and contain both one rectangle that intersects R along all its length (see figure 12).

Proof. Let us fix a positive number s > 0. We consider the rectangle R 1 (s) which is above I 1 , that is to say defined by: R 1 (s) = {(x 1 , x 2 , µ) : (x 1 , x 2 ) ∈ I 1 and µ ∈ [0, s]} , and set R (s) = R 1 (s) ∪ R 2 (s), where R 2 (s) is the rectangle of height s, associated to the segment I 2 . As previously, for any enough small, there exists a Poincaré return map P associated to the flow ϕ 2, which is defined on the rectangle R (s). Each of the images P 2 (R 1 (s)) and P 2 (R 2 (s)) is a rectangle that crosses R (s) by covering all its lenght. Therefore, to obtain a horseshoe, it suffices to verify that these two images are disjoint.

We first remark that the images under P 2 of the interiors

and thus either z 1 or z 2 would return in R (s) in a time strictly smaller than the first return time in this section. Then, let Γ (s) be the edge common to our two rectangles, that is the segment above the point Ψ (M 2 ) belonging to both I 1 and I 2 . By construction of the return map P, the distance between P 2 (Ψ (M 2 )) and the boundary of I is strictly positive (see the proof of proposition (4.7)). Thus, as this distance does not depend on and by continuity of each flow ϕ 2, , there exists a value s 0 such that for any > 0 enough small, the image P 2 (Γ (s 0 )) does not intersect the section R (s 0 ). The assertion is proved.

The conclusion of all we have done in subsection 4.1.c, is that there exists a choice of the values x * 1 , x * * 1 , µ * * for which the associated function g constructed at the beginning of subsection 4.1.c verifies that (for every > 0 enough small) the system (4) is chaotic, which is the piecewise linear version of Theorem (4.2) .