Nathalie Corson
email: nathalie.corson@univ-lehavre.fr

M A Aziz-Alaoui
email: aziz.alaoui@univ-lehavre.fr

Complex emergent properties in synchronized neuronal oscillations

Keywords: Hindmarsh-Rose model, synchronization, complexity, emergent properties

This work adresses the dynamics and complexity of the Hindmarsh-Rose neuronal mathematical model. The general aim is the study of the asymptotic behaviour of neuron networks. In this paper, the analysis of these networks uses the synchronization theory via connections between neurons which give rise to emergent properties and self-organization. Our results lead to a classical law which describes many natural or artificial self-organized complex systems. This has been performed using numerical tools.

Introduction

Also called nerve cells, neurons are the most important cells of the nervous system. The Hindmarsh-Rose model (HR hereafter), see [START_REF] Hindmarsh | A model of the nerve impulse using two first-order differential equations[END_REF][START_REF] Hindmarsh | A model of neuronal bursting using three coupled first order differential equations[END_REF], describes the dynamics of the membrane potential in the axon of a neuron with a three dimensional system of nonlinear first order differential equations which read as,

   ẋ = y + ax 2 -x 3 -z + I ẏ = 1 -dx 2 -y ż = ǫ(b(x -x c ) -z) (1)
While x describes the membrane potential, y describes the exchange of ions accross the neuron membrane through fast ionic channels and z the exchange of ions through slow ionic channels. I is the applied current while ǫ is a recovery variable, which is very small. This last parameter controls the slow motion of the neuron activity. Parameters a, b, d are constants experimentally determined and x c is the equilibrium x-coordinate of the two dimensional Hindmarsh-Rose system composed of the first two equations of (1) with z = 0 and I = 0. In this paper, we mainly study synchronization motions of coupled HR systems which can be chaotic for some parameters choices. See for example [START_REF] Corson | Asymptotic dynamics of the slow-fast Hindmarsh-Rose neuronal system, To appear in DCDIS-B 6. J.M. Gonzalez-Miranda[END_REF]6,[START_REF] Innocenti | Dynamical phases of the Hindmarsh-Rose neuronal model : Studies of the transition from bursting to spiking chaos[END_REF][START_REF] Wang | Genesis of bursting oscillations in the Hindmarsh-Rose model and homoclinicity to a chaotic saddle[END_REF], or [START_REF] Mosekilde | Bifurcation structure of a model of bursting pancreatic cells[END_REF][START_REF] Terman | Chaotic spikes arising from a model of bursting in excitable membranes[END_REF][START_REF] Terman | The transition from bursting to continuous spiking in excitable membrane models[END_REF][START_REF] Yang | Characteristics of Period-Adding Bursting Bifurcation Without Chaos in the Chay Neuron Model[END_REF] in which the behaviour bursting, spiking chaos or adding bifurcation are studied. For the reader convenience, let us recall that synchronization is a phenomenon characteristic of many processes in natural systems and nonlinear science. It has remained an objective of intensive research and is today considered as one of the basic nonlinear phenomena studied in mathematics, physics, engineering or life science. Several different regimes of synchronization have been investigated by many authors, see for example [START_REF] Aziz-Alaoui | Synchronization of chaos[END_REF]4,[START_REF] Izhikevich | Dynamical systems in neuroscience -The geometry of excitability and bursting[END_REF][START_REF] Pecora | Master stability functions for synchronized coupled systems[END_REF][START_REF] Pikovsky | Synchronization, A universal concept in nonlinear sciences[END_REF]. This word has a greek root, syn = common and chronos = time, which means to share common time or to occur at the same time, that is correlation or agreement in time of different processes. Thus, synchronization of two dynamical systems generally means that one system somehow follows the motion of another. A lot of research has been carried out and, as a result, showed that even chaotic systems could synchronize when they are coupled. Many researchers have discussed the theory, the design or applications of synchronized motion in coupled chaotic systems. A broad variety of applications have emerged, for example to increase the power of lasers, to synchronize the output of electronic circuits, to control oscillations in chemical reactions or to encode electronic messages for secure communications. Here are some synchronization regimes,

• Identical (or complete) Synchronization, which is defined as the coincidence of states of interacting systems, • Generalized Synchronization, which extends the identical synchronization phenomenon and implies the presence of some functional relation between two coupled systems; if this relationship is the identity we recover the identical synchronization,

• Phase Synchronization, which means driving of phases of chaotic oscillators, whereas their amplitudes remain uncorrelated, • Lag Synchronization, which appears as a coincidence of shifted-in-time states of two systems.

Hereafter, for all numerical experiments, we use HR system with the following coordinates changes, see [START_REF] Belykh | Synchronization of Bursting Neurons: What matters in the Network Topology[END_REF],

y = 1-y, z = 1+I +z, d = a+α, c = -1-I -bx c .
From this changes of coordinates, we obtain,

   ẋ = ax 2 -x 3 -y -z ẏ = (a + α)x 2 -y ż = ǫ(bx + c -z) (2)
Let us consider a network composed by n HR neurons. These neurons are coupled by their first variable x i . A model of this network is given by,

   ẋi = ax 2 i -x 3 i + y i -z i -h(x i , x j ), i = j, i = 1, ..., n, j = 1, ..., n ẏi = (a + α)x 2 i -y i żi = ǫ(bx i + c -z i ) (3)
In the case of a network of n neurons linearly coupled, the coupling function h is chosen as,

h(x i , x j ) = k n n j=1 c ij Γ (x j ) (4)
in which the synaptic coupling Γ is a linear function, Γ (x j ) = x i -x j , for i = 1, 2, ..., n. Parameter k n represents the coupling strength and C n = (c ij ) is the n × n connectivity matrix, c ij = 1 if i and j are connected, i = 1, ...n, j = 1, ...n, i = j c ij = 0 if i and j are not connected.

This matrix C n can be symmetric or not, so that unidirectional or bidirectional coupling are possible.

In the case of a network of n neurons nonlinearly coupled, the coupling function h is given by [START_REF] Belykh | Synchronization of Bursting Neurons: What matters in the Network Topology[END_REF] and reads as,

h(x i , x j ) = (x i -V )k n n j=1 c ij Γ (x j ) (5)
in which the synaptic coupling Γ is a nonlinear function with a threshold as given,

Γ (x j ) = 1 1 + exp(-λ(x j -Θ)) ( 6
)
where Θ is the threshold reached by every action potential for a neuron and λ is defined in [START_REF] Belykh | Synchronization of Bursting Neurons: What matters in the Network Topology[END_REF]. As in the case of linear coupling, parameter k n corresponds to the synaptic coupling strength. Parameter V is the reversal potential and must be larger than x i (t) for all i and all t since synapses are supposed excitatory.

We have chosen the situation in which each neuron has the same number of inputs from other neurons. Indeed, according to [START_REF] Belykh | Synchronization of Bursting Neurons: What matters in the Network Topology[END_REF], it is a necessary condition for the synchronous solution to exist. Therefore, the synaptic connections we use have to be bidirectional. The C n matrix should be symmetric and in our case, C n = 1 n -Id n where 1 n = (1) nn and Id n = (Id) nn is the identical matrix.

The following parameters are fixed as follows throughout this paper,

a = 2.8, α = 1.6, c = 5, b = 9, ǫ = 0.001 (7) V = 2, λ = 10, Θ = -0.25 (8)
Interactions between nerve cells are possible thanks to synapses. The word synapse has a greek root, syn = common and haptein = to touch, which means connection. This is the functional contact part which exists between two neurons or between a neuron and another cell (muscular cell, sensorial receptor, ...). It operates the conversion of an action potential fired by the presynaptic neuron into a signal in the postsynaptic cell. Usually, two different types of synapses are discerned, which are,

• the chemical synapse, widely majoritary, which needs some neurotransmitters to transmit information, • the electrical synapse, in which the signal is transmitted electrically through gap-junctions.

The size of the synaptic cleft is characteristic of one or the other kind of synapse. In the case of electrical synapses, this synaptic cleft is about two nanometers, while it can reach from ten to forty nanometers in the case of chemical synapses.

This paper is organized as follows. In the next section, the coupling function between mathematical neurons is a linear function in order to represent electrical synaptic connections. In this section, we firstly assume that neurons are identical and then, we assume that they are slightly different. In section 3, the same steps are followed with a nonlinear function in order to study chemical synaptic connections.

Property emerging from synchronization with a linear coupling

In this paper, the asymptotic behaviour is studied only by numerical analysis in various regions of the parameter space commonly used in the litterature. We have done the numerical computations very accurately in double precision and for different variations of parameters. In order to obtain reliable numerical results, the step size has been chosen to be equal to 10 -4 and the first 10 6 steps are discarded to avoid the transient regime.

Coupling identical HR neurons

The first step is to consider the asymptotic behaviour of two neurons lineraly coupled. Therefore, we numerically study system (3) with coupling function (4) for n = 2. The connectivity matrix C 2 and the graph of the network are given in figure 2.

C2 = 12 -Id2 = 0 1 1 0 Fig. 2.
Bidirectional connection scheme for two neurons and the connectivity matrix C2.

This gives the six-dimensional system representing two bidirectionaly coupled HR neurons given in [START_REF] Mosekilde | Bifurcation structure of a model of bursting pancreatic cells[END_REF],

                 ẋ1 = ax 2 1 -x 3 1 + y 1 -z 1 -k 2 (x 1 -x 2 ) ẏ1 = (a + α)x 2 1 -y 1 ż1 = ǫ(bx 1 + c -z 1 ) ẋ2 = ax 2 2 -x 3 2 + y 2 -z 2 -k 2 (x 2 -x 1 ) ẏ2 = (a + α)x 2 2 -y 2 ż2 = ǫ(bx 2 + c -z 2 ) (9)
Figure 3 shows the evolution of the (x 1 , x 2 ) phase portrait of system (9) while making the coupling strength k 2 increase after a certain period of time. For low values of this coupling strength, the behaviour of x 1 is really different from the one of x 2 , whereas when this coupling strength is large enough, x 1 and x 2 behaviours become synchronous. We observe that synchronization motion between the first variable x 1 of one neuron and the first variable x 2 of the other, i.e :

x 1 = f (x 2 ) (here, f (x) = Id(x))
, appears for a coupling strength k 2 ≥ 0.376. For every coupling strength larger that this value of k 2 , x 1 and x 2 remain synchronous. This phenomenon is also observed for variables y and z, for the same values of parameters. Therefore, this synchronization is total.

In the case of three neurons, the connectivity matric C 3 and the graph of the network are given in figure 4. System (3) with coupling function (4) for n = 3 reads as,

C3 = 13 -Id3 = 0 1 1 1 0 1 1 1 0
                                 ẋ1 = ax 2 1 -x 3 1 + y 1 -z 1 -k 3 (x 1 -x 2 ) -k 3 (x 1 -x 3 ) ẏ1 = (a + α)x 2 1 -y 1 ż1 = ǫ(bx 1 + c -z 1 ) ẋ2 = ax 2 2 -x 3 2 + y 2 -z 2 -k 3 (x 2 -x 1 ) -k 3 (x 2 -x 3 ) ẏ2 = (a + α)x 2 2 -y 2 ż2 = ǫ(bx 2 + c -z 2 ) ẋ3 = ax 2 3 -x 3 3 + y 3 -z 3 -k 3 (x 3 -x 1 ) -k 3 (x 3 -x 2 ) ẏ3 = (a + α)x 2 3 -y 3 ż3 = ǫ(bx 3 + c -z 3 ) (10)
As shown in the case of two neurons, figure 5 shows the evolution of x 1 according to x 2 for system [START_REF] Innocenti | Dynamical phases of the Hindmarsh-Rose neuronal model : Studies of the transition from bursting to spiking chaos[END_REF] while making the coupling strength increase.

For low values of this coupling strength, the behaviour of x 1 is really different from the one of x 2 , whereas when this coupling strength is large enough, x 1 and x 2 behaviours become identical. If figure 5 only shows the (x 1 , x 2 ) view of the phase portraits of system (10), a similar phenomenon is observed for between x 1 and x 3 . Therefore, the minimum coupling strength k 3 to exhibit synchronization between the first variables x 1 , x 2 and x 3 of three HR neurons is about k 3 = 0.255. Once again, we observe also this phenomenon for the variables y and z for the same numerical values of parameters. This is a synchronization threshold.

There are different possible ways for coupling four, five or more neurons. However, as we already said, we consider the case in which each neuron is connected to all the others. Therefore, there is a unique possible connection scheme.

Here, we consider system (3) with coupling function (4) for n = 4, n = 5, n = 6, n = 7, n = 8. Figures 6, 8, 10, 12 and 14 show the connection graphs of the network and the connectivity matrix C n respectively in the case of 4, 5, 6, 7 and 8 neurons coupled all together. Fig. 6. Bidirectional connection scheme for four neurons. The connectivity matrix C4 is C4 = 14 -Id4. Figures 7,[START_REF] Mosekilde | Bifurcation structure of a model of bursting pancreatic cells[END_REF][START_REF] Izhikevich | Dynamical systems in neuroscience -The geometry of excitability and bursting[END_REF][START_REF] Pikovsky | Synchronization, A universal concept in nonlinear sciences[END_REF][START_REF] Terman | The transition from bursting to continuous spiking in excitable membrane models[END_REF] show the evolution of the (x 1 , x 2 ) phase portrait while making the coupling strength increase in the case of 4, 5, 6, 7 and 8 neurons coupled all together using the coupling function given in (4). These phase portraits are plotted after a period of time. These figures show x 2 according to x 1 , but the same motion is observed comparing the behaviour of the first variable of the other corresponding variables. 1. Tabular summarizing minimal coupling strength kn to observe synchronous motion of n neurons, with n = 2, ..., 8 when coupling identical HR neurons using the linear function (4). Kn is the theoretical value of coupling strength obtained with the heuristic law given in [START_REF] Izhikevich | Dynamical systems in neuroscience -The geometry of excitability and bursting[END_REF].

HR synchronization and complex emergent properties

In this tabular, the first column shows the number n of neurons in the network.

The second column shows the numerical results we obtained while the last column shows the theoretical results obtained thanks to the heuristic law we observed,

K n = k 2 0.62 * (n -1) (11)
Figure 16 (a) uses table 1 to show k n versus n (dots) and K n versus n (solid lines). Figure 16 (b) shows the curve obtained plotting Log(K n ) according to Log(n). The curve of figure 16 (b) corresponds to a classical law (for n large enough) which can describe many self-organized complex systems, like earthquakes (for 1000 earthquakes of magnitude 4 on Richter scale, there are only 100 of magnitude 5 and 10 of magnitude 6), linguistics (for 1000 occurrences of "the" in an english text, there are only 100 occurrences of "I" and 10 of "say"), urban systems (big cities are rare and small ones are frequent in an exponential way).

Coupling non-identical HR neurons

In this subsection, the same coupling function ( 4) is used to study the asymptotic behaviour of networks composed of HR neuron which parameters are slightly different from one another. Indeed, in nature, it is not realistic to assume that a neuron network is composed of exactly identical neurons. We then make all parameters vary slightly from one neuron to another. Let p i be a generic way of writing every parameters of neuron i (i.e. a i , α i , b i , c i ). The variation of parameters from one neuron to another is defined as in [START_REF] Pecora | Master stability functions for synchronized coupled systems[END_REF]. It has been choosen around e = 10 -4 ,

     p i = p 1 + (1 + i 2 10 -1
) × e, i = 2j, j = 1, 2, 3, 4

p i = p 1 + (1 - i -1 2 10 -1 ) × e, i = 2j + 1, j = 1, 2, 3 (12)
With these variations of parameters from a neuron to another, we obtain the following tabular (Table 2) of minimal coupling strength to obtain synchronization, and figure 17 showing a similar law as given in figure 16. Then, in order to make the difference between neurons bigger, the variation is chosen as done in ( 12) with e = 10 -3 . With these variations of parameters from a neuron to another, we obtain the following tabular (Table 3) of synchronization thresholds, and figure 18 showing a similar law as previously presented. 3. Tabular summarizing minimal coupling strength to observe synchronous motion of n neurons, with n = 2, ..., 8 when coupling slightly distinct HR neurons using linear function (4). The difference betwenn neurons is around 10 -3 .

As the previous numerical results show, the synchronization threshold follows the heuristic law,

K n = k 2 /(n -1) (15)
This result has been first presented in [START_REF] Belykh | Synchronization of Bursting Neurons: What matters in the Network Topology[END_REF].

Coupling non-identical HR neurons

In this subsection, the coupling function [START_REF] Corson | Asymptotic dynamics of the slow-fast Hindmarsh-Rose neuronal system, To appear in DCDIS-B 6. J.M. Gonzalez-Miranda[END_REF], with Γ defined as in (6), is used to study the asymptotic behaviour of networks composed of coupled HR neurons, the parameters of which are slightly different from one another, as in subsection 2.2. The variation of parameters from one neuron to another is defined as in ( 12) with e = 10 -4 . With these different neurons, we obtain the following tabular (Table 5) which gives the synchronization thresholds from two to eight coupled neurons. The obtained results are similar to those we observed in the previous section. They show the same heuristic law given by equation [START_REF] Terman | The transition from bursting to continuous spiking in excitable membrane models[END_REF]. These results plotted figure 20. Emergence and complexity refer to the appearance of higher-level properties and behaviours of a system that obviously comes from the collective dynamics of that system's components, see [START_REF] Aziz-Alaoui | Complex emergent properties and chaos (De) synchronization[END_REF]. These properties are not directly deductable from the lower-level motion of that system. Emergent properties are properties of the 'whole' that are not possessed by any of the individual parts making up that whole. For example, an air molecule is not a cyclone, an isolated species does not form a food chain and an isolated neuron is not conscious : emergent behaviours are typically novel and unanticipated. Thus, if a synchronization phenomenon is exhibited in a neuron network, it may develop some kind of consciousness. This is an emergent property which comes from the collective dynamics of n neurons. Moreover, as given in figure (19), as the number of neurons n increases, the synchronization threshold k n decreases.

Consciousness is more important when the number of neurons is larger?

Fig. 1 .

 1 Fig. 1. System (1) for a = 3, b = 4, d = 5, xc = -1 2 (1 + √ 5), I = 3.25 and ǫ = 0.001. (a) Time series (t, x), (b) (x, y, z) view.

Fig. 3 .

 3 Fig.3. Phases portraits for system[START_REF] Mosekilde | Bifurcation structure of a model of bursting pancreatic cells[END_REF] showing x2 according to x1 for the coupling strength (a) k2 = 0.2, (b) k2 = 0.25, (c) k2 = 0.36, (d) k2 = 0.376. For two neurons linearly coupled, the synchronization threshold for the coupling strength is about k2 = 0.376.

Fig. 4 .

 4 Fig. 4. Bidirectional connection scheme for three neurons and the connectivity matrix C3.

Fig. 5 .

 5 Fig. 5. Phases portraits for system (10) showing x2 according to x1 for the coupling strength (a) k3 = 0.1, (b) k3 = 0.2, (c) k3 = 0.24, (d) k3 = 0.255. For three linearly coupled neurons, the synchronization threshold of the coupling strength is around k3 = 0.255.

Fig. 7 .

 7 Fig. 7. x2 according to x1 when four neurons are connected for the coupling strength (a) k4 = 0.1, (b) k4 = 0.15, (c) k4 = 0.18, (d) k4 = 0.191. Minimum coupling strength to observe synchronization between variables xi and xj, yi and yj , zi and zj (i = 1, ..., 4, j = 1, ..., 4, i = j) when coupling four neurons is k4 = 0.191.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Coupling five neurons. The connectivity matric C5 is C5 = 15 -Id5.

Fig. 10 .

 10 Fig. 10. Coupling six neurons. The connectivity matrix C6 is C6 = 16 -Id6.

Fig. 11 .

 11 Fig. 11. x2 according to x1 when six neurons are connected for the coupling strength (a) k6 = 0.03, (b) k6 = 0.05, (c) k6 = 0.1, (d) k6 = 0.122. Minimum coupling strength to observe synchronization between variables xi and xj, yi and yj , zi and zj (i = 1, ..., 6, j = 1, ..., 6, i = j) when coupling six neurons is k6 = 0.122.

Fig. 12 .

 12 Fig. 12. Coupling seven neurons. The connectivity matrix C7 is C7 = 17 -Id7.

Fig. 13 .

 13 Fig. 13. x2 according to x1 when seven neurons are connected for the coupling strength (a) k7 = 0.025, (b) k7 = 0.05, (c) k7 = 0.1, (d) k7 = 0.104. Minimum coupling strength to observe synchronization between variables xi and xj , yi and yj, zi and zj (i = 1, ..., 7, j = 1, ..., 7, i = j) when coupling seven neurons is k7 = 0.104.

Fig. 14 .Fig. 15 .

 1415 Fig. 14. Coupling eight neurons. The connectivity matrix C8 is C8 = 18 -Id8.

Fig. 16 .

 16 Fig. 16. (a) Minimum coupling strength kn to observe a synchronous behaviour of n neurons, according to the number of neurons in the network (dots) plotted together with Kn = k2/(0.62(n -1)) (solid line). (b) Log(Kn) according to Log(n).

Table 2 .Fig. 17 .

 217 Fig. 17. (a) Minimum coupling strength kn to observe a synchronous behaviour of n neurons according to the number of neurons in the network n (dots) plotted together with Kn = k2/(0.62(n -1)) (solid line). (b) Log(Kn) according to Log(n).

Table 5 .

 5 Tabular summarizing minimal coupling strength kn to observe synchronous motion of n neurons, with n = 2, ..., 8 when coupling slightly distinct HR neurons using the nonlinear function (5) Γ defined as in (6). The third column shows the coupling strength obtained with the heuristic law given in[START_REF] Terman | The transition from bursting to continuous spiking in excitable membrane models[END_REF].

Fig. 20 .

 20 Fig. 20. (a) Minimum coupling strength kn to observe a synchronous behaviour of n neurons according to the number of neurons in the network (dots) plotted together with Kn = k2/(n -1) (solid line). (b) Log(Kn) according to Log(n).

  Synchronization thresholds of the coupling strength we obtained are summarized in table 1.

	n	2	3	4	5	6	7	8
	kn 0.376 0.255 0.191 0.152 0.122 0.104 0.088
	Kn		0.303 0.202 0.151 0.121 0.101 0.086
	Table

3 Property emerging from synchronization with a nonlinear coupling

Coupling identical HR neurons

For n = 2, system (3) with the sigmoïdal coupling function h defined as in [START_REF] Corson | Asymptotic dynamics of the slow-fast Hindmarsh-Rose neuronal system, To appear in DCDIS-B 6. J.M. Gonzalez-Miranda[END_REF] and Γ defined as in (6) reads as,

For n = 3, system (3) with the coupling function h defined as in [START_REF] Corson | Asymptotic dynamics of the slow-fast Hindmarsh-Rose neuronal system, To appear in DCDIS-B 6. J.M. Gonzalez-Miranda[END_REF] and Γ defined as in (6) reads as,

The same numerical experiments as in the previous section lead us to table 4, in which theorical coupling strength K n given by the heuristic law ( 14) is also presented,


Nathalie Corson
email: nathalie.corson@univ-lehavre.fr

M A Aziz-Alaoui

ASYMPTOTIC DYNAMICS OF THE SLOW-FAST HINDMARSH-ROSE NEURONAL SYSTEM

Keywords: neuron model, asymptotic dynamics, bifurcation, chaos AMS (MOS) subject classification: 34C15, 37D45, 37G35, 37N25, 65P20, 92C20

This work adresses the asymptotic dynamics of a neuronal mathematical model. The aim is first the understanding of the biological meaning of existing mathematical systems concerning neurons such as Hodgkin-Huxley or Hindmarsh-Rose models. The local stability and the numerical asymptotic analysis of Hindmarsh-Rose model are then developed in order to comprehend bifurcations and dynamics evolution of a single Hindmarsh-Rose neuron. This has been performed using numerical tools borrowed from the nonlinear dynamical system theory.

Introduction and neurophysiology

Also called nerve cells, neurons are the most important cells of the nervous system. They are composed of a cell body, or soma, extended by an axon, an axon terminal and some dendrites (Fig. 1). These extensions are useful not only for the transmission of information through a neuron but also for the transmission from a neuron to another. Connections between neurons are possible thanks to synapses. Ionic channels can be seen as macrocellular pores in the neuron membrane. They enable molecules to pass through the membrane. They are the link between intra-cellular and extra-cellular space. There are many types of ionic channels. For example, we may cite those which are always open, those which are voltage-dependent or those which select the molecule allowed to cross the membrane. Sodium channels are so called because they are specific to sodium ions. They can be in an active state or in an inactive one. Potassium channels open and close with delay. Leak channel are always open.

Neurons are enclosed by a membrane which separates the interior of the cell from the extracellular space. On the inside of the cell, the concentration of ions is different than in the surrounding liquid. The difference of concentration generates an electrical potential which plays an important role in neuronal dynamics. The last is the membrane potential.

When not sending any signal, a neuron is at rest. At rest, there are relatively more sodium ions (Na+) on the outside of the neuron and more potassium (K+) ions on the inside. The inside of the neuron is negative compared to the outside. The difference in the voltage between the inside and the outside of the neuron gives the resting potential, the value of which is in general about -70 mV.

The equilibrium potential is predictable with the Nernst equation. An action potential, also called spike or impulse (Fig. 2) occurs when a neuron sends information along an axon away from the cell body. The action potential is an explosion of electrical activity that is created by a depolarizing current. A stimulus makes the resting potential move up towards positive values. If the potential reaches a certain threshold (about -55 mV), the neuron fires an action potential the amplitude of which is always the same. If the potential does not reach this threshold value, no action potential fires. Therefore, a neuron respects the all or none principle. At that time, sodium channels start to close. This induces the repolarization. Indeed, it makes the membrane potential go back towards the resting potential. Then, there is the hyperpolarization of the neuron, since the potential goes past the resting potential because the potassium channels close with delay. Gradually, the ion concentrations go back to resting levels and the cell returns to -70 mV (Fig. 2). For more details, see for example [START_REF] Hille | Ionic channels of excitable membranes[END_REF].

A burst is a group of at least two action potentials that occur close together in time, separated from other action potentials by large time intervals which are called quiescent or silent phases.

Mathematical models

Two neurophysiologists, Alan Lloyd Hodgkin and Andrew Fielding Huxley, have developed an empirical kinetic description of ionic mechanisms in a neuron, see [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]. This description was not only simple enough to make practical computation of electrical responses but also sufficiently precise so as to predict conduction and major features of excitability. Their model comprises mathematical equations and suggests the main behaviour of the gating mechanisms, [START_REF] Hille | Ionic channels of excitable membranes[END_REF][START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]. For shortness, HH is a reference to the Hodgkin Huxley model.

The model is based on Sodium, Potassium and leakage ion flow. HH model has separated equations for sodium conductance g N a and potassium conductance g K . They are expressed as maximum conductances ḡNa and ḡK multiplied by coefficients representing the fraction of the maximum conductances actually expressed. The multiplying coefficients are numbers varying between 0 and 1. In the model, the conductance changes depend only on voltage and time. Because of their observations, Hodgkin and Huxley noted that such kinetics would be obtained if the opening of a K channel is controlled by several independent "particles". Suppose that there are four identical particles, each one with a probability n of being in the correct position to set up an open channel. The probability that all four particles are correctly placed is n 4 . Mathematically, I K is represented in the HH model by,

I K = n 4 ḡK (V -E(K))
The voltage and time dependent changes of n are given by the reaction, 1 -n ⇋ αn βn n where the gating particles make transitions between the permissive and nonpermissive forms with voltage-dependent rate constants α n and β n . If the initial value of the probability n is known, subsequent values can be computed by solving the simple differential equation,

dn dt = α n (1 -n) -β n n
The HH model uses a similar formalism to describe I N a , with four hypothetical gating particles making independent first order transitions between permissive and non permissive state to control the channel. However, because there are two opposing gating processes, activation and inactivation, there must be two kinds of gating particles, say h and m. Three m particles control activation and one h particle controls inactivation. Therefore, the probability that they are all in permissive position is m 3 h. I N a is then represented by,

I N a = m 3 hḡ N a (V -E(N a))
As for the n parameter of K channels, m and h are assumed to undergo transitions between permissive and non-permissive forms with rates satisfying differential equations,

dm dt = α m (1 -m) -β m m dh dt = α h (1 -h) -β h h
HH model for the squid giant axon describes ionic current accross the membrane in terms of three components,

I N a = m 3 hḡ N a (V -E(N a)) , I K = n 4 ḡK (V -E(K)) , I L = ḡL (V -E(L))
where ḡL is a fixed background leakage conductance.

Therefore, the HH model reads as follows,

                       -C dV dt = m 3 hḡ N a (V -E(N a)) + n 4 ḡK (V -E(K)) + ḡL (V -E(L)) -I dn dt = α n (1 -n) -β n n dm dt = α m (1 -m) -β m m dh dt = α h (1 -h) -β h h
in which C represents the membrane capacity, V the total membrane potential, m the N a activation variable, h the N a inactivation variable, ḡNa the maximum sodium conductance, E(N a) the N a equilibrium potential, n the K activation variable, ḡK the maximum potassium conductance, E(K) the K equilibrium potential, ḡL the maximum leakage conductance, E(L) the leakage equilibrium potential, I the external current (applied during an experiment), α i the gate inactivation rate (i = m, n, h), β i the gate activation rate (i = m, n, h).

In 1982, J.L. Hindmarsh and R.M. Rose simplified the Hodgkin-Huxley model, see [START_REF] Hindmarsh | A model of the nerve impulse using two first-order differential equations[END_REF]. The model they obtained is of Fitzhugh-Nagumo type. Indeed, they observed that some variables could be replaced by constants and they found relations between different variables. Therefore, they simplified the Hodgkin-Huxley system into a more simple one with two euqations, as follows,

ẋ = y + ax 2 -x 3 ẏ = 1 -dx 2 -y (1)
Two years later, Hindmarsh and Rose decided to add a third equation to their model, so that the motion of their model could be closer to the motion of a real neuron, see [START_REF] Hindmarsh | A model of neuronal bursting using three coupled first order differential equations[END_REF]. This Hindmarsh-Rose model (HR) describes the dynamics of the membrane potential x in the axon of a neuron. It is a three dimensional system of non linear first order differential equations, which reads as,

   ẋ = y + ax 2 -x 3 -z + I ẏ = 1 -dx 2 -y ż = ǫ(b(x -x c ) -z) (2)
While x describes the membrane potential, y describes the exchange of ions accross the neuron membrane through fast ionic channels and z the exchange of ions through slow ionic channels. Parameters a, b, d are constants experimentally determined and

x c = - 1 2 (1 + √ 5
) is the equilibrium x-coordinate of system (1) for the following parameters,

a = 3, b = 4, d = 5, x c = - 1 2 (1 + √ 5) (3)
In this paper, I and ǫ are chosen as bifurcation parameters due to their biological meaning. Indeed, I is the applied current while ǫ is a recovery variable, which is very small. This last parameter controls the slow motion of the neuron activity. This system, as well as HH model, has been studied in different ways, see for example [START_REF] Françoise | Oscillations en Biologie : Analyse qualitative et Modèles, Collection Mathématiques et Applications[END_REF][START_REF] Izhikevich | Dynamical systems in neuroscience -The geometry of excitability and bursting[END_REF] or [START_REF] Ginoux | Slow manifold of a neuronal bursting model, Emergent Properties in Natural and Artificial Dynamical Systems[END_REF][START_REF] Shilnikov | Methods of the qualitative theory for the Hindmarsh-Rose model : a case study[END_REF][START_REF] Vidal | Periodic orbits of tritrophic slow-fast system and double homoclinic bifurcations[END_REF] and the references therein cited. Interesting studies are also presented in [START_REF] Belykh | Synchronization of Bursting Neurons: What matters in the Network Topology[END_REF][START_REF] Corson | Complex emergent properties in synchronized neuronal oscillations, From System Complexity to Emergent Properties[END_REF][START_REF] Gonzalez-Miranda | Complex bifurcation structures in the Hindmarsh-Rose neuron model[END_REF][START_REF] Guckenheimer | Bifurcations in the fast dynamics of neurons: implications for bursting, Bursting, The Genesis of Rhythm in the Nervous System[END_REF][START_REF] Innocenti | Dynamical phases of the Hindmarsh-Rose neuronal model : Studies of the transition from bursting to spiking chaos[END_REF].

3 Bifurcations and Chaos

Dissipativity and existence of attractors

Theorem 1 If 0 < a < √ 3(1 + ǫ) 1 2 then system (2) is dissipative.
Proof :

The variation of the volume V (t) of a small element δΩ(t) = δxδyδz in the phase space is determined by the divergence of the flow,

∇V = ∂ ẋ ∂x + ∂ ẏ ∂y + ∂ ż ∂z System (2) leads to, ∇V = -3x 2 + 2ax -1 -ǫ ∇V < 0 ⇔ -3x 2 + 2ax -1 -ǫ < 0 Therefore, if -3x 2 + 2ax -1 -ǫ < 0, then ∇V < 0. So that with 0 < a < √ 3(1 + ǫ) 1 2 , system (2)
is dissipative (obviously, one can prove the existence of a larger parameters range for which this system remains dissipative).

Besides, since -3x 2 + 2ax -1 -ǫ < a 2 3 -1 -ǫ < 0, system (2) has, for the volume element δΩ(t) = δxδyδz, an exponential contraction rate satisfying,

δΩ(t) ≤ e a 2
That is, a volume element V 0 is contracted by the flow into a volume element less than V 0 e ( a 2 3 -1-ǫ)t in time t. In other words, each volume containing the system trajectory shrinks to zero as t → ∞ at an expontential rate less than a 2 3 -1 -ǫ, which is independent of x, y and z. Thus, all trajectories are ultimately confined to a specific subset having zero volume (which, of course, can not ensure the existence of bounded attractor), this has been confirmed by our computer simulations.

Chaos dans HR84

The chaotic behaviour of the Hindmarsh-Rose system is shown by numerical analysis in different regions of the parameters space, by bifurcation diagrams, Lorenz plots, Lyapunov exponents, time series and phase portraits. Since ǫ << 1, the first two equations of HR system correspond to the fast dynamics while the last equation controls the slow dynamics. In figure 3 (e), the fast dynamics is the spiking part of time series while slow dynamics is the resting part. In figure 3 (a),(b),(c),(d), the fast dynamics corresponds to the right part of the attractor, the part where we observe spirals, while slow dynamics is on the left of the figures. An important phenomenon in neuron activity is the transition between spiking and bursting. Spiking is represented by a generation of action potentials, while bursting is represented by a membrane potential changing from resting to repetitive firing state. We can see also these phenomena looking at time series and phase portraits, but bifurcation diagrams that are shown in this section do not give any information on this "slow-fast" motion. Indeed, in order to obtain these bifurcation diagrams, the chosen Poincaré section crosses the phase portrait within the hole of the attractor which is on the right part on figures 3 (b,c,d,e) and corresponds to the fast dynamics. Therefore, the following diagrams give only information about the fast dynamics of HR system. In the case of neuron dynamics, the location of bifurcation values is important to determine the transition between a quiescent state and an oscillatory one, and also between different kinds of oscillatory motion. Therefore, a bifurcation corresponds to a qualitative change of the information transmitted through the axon of the neuron.

In this paper, codimension-1 bifurcations (that is the number of independent conditions determining this bifurcation is one) are numerically computed. We firstly make the slow parameter ǫ be the control parameter and then, the injected current I is chosen for this assignment.

Dynamics with respect to ǫ

Parameter ǫ is the ratio of time scales between spiking (fast dynamics) and resting (slow dynamics). Therefore, it controls the difference between the slow and the fast dynamics of HR model corresponding to the difference between fast fluxes of ions accross the membrane to slow ones. Thereby, making this parameter evolving to observe the neuron reaction is really interesting, see [START_REF] Ginoux | Slow manifold of a neuronal bursting model, Emergent Properties in Natural and Artificial Dynamical Systems[END_REF].

A bifurcations diagram show the evolution of the aysmptotic behaviour of solutions according to one parameter. The successive local maxima map (first return map, Lorenz plot) is a reliable numerical method that gives information about the presence (or not) of chaos for some fixed parameters. This map gives the local maxima of z against the immediatly preceding local maxima. Looking at phase portraits, it is justified to wonder if, for some parameters, HR model exhibits a chaotic dynamics.

The Lorenz plot gives a good idea of this possible chaotic behaviour, when the resulting plot is unimodal with a shape similar to a skewed Hénon map. Lyapunov exponents measure the sensibility to initial conditions and are numerical tools to have an idea of the presence of chaos. If all Lyapunov exponents are negative, two initialy close trajectories remain close to each other, while a positive Lyapunov exponent means that two initialy close trajectories diverge. Therefore, plotting the larger Lyapunov exponent gives a reliable evidence of the presence of chaos. shows, among other things, that there is an ǫ 1 ∈ [0.00041, 0.00049] for which the neuron behaviour changes abruptly. Indeed, ∀ǫ < ǫ 1 , the neuron exhibits a tonic spiking motion and, ∀ǫ > ǫ 1 , the neuron exhibits a bursting motion. Moreover, this figure shows that system

(2) with parameters given in (3) and I = 3.25 does not exhibit chaotic behaviour for ǫ ∈ [ǫ 1 , 0.002]. On figure 5(a) one can see that this model does not exhibit any chaotic behaviour for ǫ = 0.0005. Figure 5(b) confirms the fact that there is no presence of chaos for ǫ ∈ [ǫ 1 , 0.002], since the larger Lyapunov exponent remains really close to zero. Time series and (x, y, z) phase portraits of figure 6 show how the variation of parameter ǫ makes the dynamics of system (2) change from spiking to periodical bursting. These figures confirm the non-chaotic behaviour of this system for ǫ ∈ [ǫ 1 , 0.002].

The enlargement of figure 4(a) for ǫ ∈ [0.005, 0.015] shown in figure 4(c) exhibits not only inverse period doubling cascades starting with period 3, period 4 or period 5 but also some dark parts, which is a numerical sign of chaotic motion. Some other numerical evidence can be seen, for ǫ = 0.008, on the Lorenz plot of figure 7

Dynamics with respect to I

Parameter I corresponds to the current which is injected in the neuron. Thus, it can be controlled during experiments and can then play the role of bifurcation parameter. Time series and phase portraits given in figure 13 show a periodic behaviour of system (2) for the chosen values of parameter I. Moreover, for each value of I, there is the same number of branches within the bifurcation diagram (figure 11), the same number of spikes within bursts in time series and the same number of laps in the phase portraits (figure 13). This number correponds to the period of the associated limit cycle.

Conclusion

In this paper, an accurate numerical analysis of the asymptotic dynamics of Hindmarsh-Rose model has been done. We have discussed the chaotic dynamics of this system by use of different nonlinear dynamics numerical tools such as bifurcation diagrams, Lorenz plots and Lyapunov exponents for different parameter ranges. Two of the HR parameters, ǫ the time scale ratio between fast and slow dynamics and I the applied current, have been chosen as bifurcation parameters due to their biological meaning. It is obvious that the slow-fast motions due to the small parameter ǫ makes the dynamics even more interesting.

Figure 1 :

 1 Figure 1: Structure of a typical neuron.

  It is given by, E(S) = RT ZF ln [S] out [S] in , where E(S) represents the Nernst potential for ion S (measured as for membrane potential, inside with respect to outside), [S] out represents the concentration of S outside the cell, [S] in the concentration of S inside the cell, R is the gas constant, T represents the absolute temperature, Z represents the valence of ion S : Z(K+) = +1, Z(N a+) = -1, and F is Faraday's constant.

Figure 2 :

 2 Figure 2: Action potential (spike or impulse).

Figure 3 :

 3 Figure 3: Numerical integration of system (2) for parameters given in (3) with I = 3.25 and ǫ = 0.001. (a) : (t, x) times series showing a sequence of action potentials. (b), (c), (d), (e) : phase portraits respectively in x-y projection, x-z projection, y-z projection and (x, y, z) tridimensional view.

Figure 4 (

 4 Figure 4(a) gives the bifurcation diagram with respect to the control parameter ǫ in the range [0, 0.05].In order to have a more accurate analysis of the dynamics of system (2), we present in figure4(b),(c),(d) enlargements of figure4(a). Figure4(b) shows, among other things, that there is an ǫ 1 ∈ [0.00041, 0.00049] for which the neuron behaviour changes abruptly. Indeed, ∀ǫ < ǫ 1 , the neuron exhibits a tonic spiking motion and, ∀ǫ > ǫ 1 , the neuron exhibits a bursting motion. Moreover, this figure shows that system

Figure 4 (

 4 Figure 4(a) gives the bifurcation diagram with respect to the control parameter ǫ in the range [0, 0.05].In order to have a more accurate analysis of the dynamics of system (2), we present in figure4(b),(c),(d) enlargements of figure4(a). Figure4(b) shows, among other things, that there is an ǫ 1 ∈ [0.00041, 0.00049] for which the neuron behaviour changes abruptly. Indeed, ∀ǫ < ǫ 1 , the neuron exhibits a tonic spiking motion and, ∀ǫ > ǫ 1 , the neuron exhibits a bursting motion. Moreover, this figure shows that system

  (a) which gives a typical unimodal map corresponding to a chaotic behaviour. Moreover, either the positive larger Lyapunov exponent of figure 7(b) or the time series and (x, y, z) phase portraits of figure 8 point out the chaotic behaviour of the model for some ranges of ǫ ∈ [0.005, 0015]. The enlargement of figure 4(c) for ǫ ∈ [0.0138, 0.0148] shown in figure4(d) also exhibits a chaotic behaviour of system (2). This is confirmed not only by the Lorenz plot for ǫ = 0.0145 (figure9(a)) but also by the positive larger Lyapunov exponent in the left part of figure9(b) and by the time series and (x, y, z) phase portrait presented in figure10(a 1,2 ). The right part of figure4(a) exhibits a reverse period doubling cascade. As ǫ becomes larger, the number of spikes within a burst decreases until the bursting motion of the neuron disappears to let the spiking motion arise. This phenomenum is also obvious while looking at figure 10(b 1,2 ),(c 1,2 ),(d 1,2 ).

Figure 4 :Figure 5 :

 45 Figure 4: Bifurcation diagrams in (ǫ, y) plane for system (2) with parameters given in (3) and with I = 3.25. (a) An inverse period doubling cascade is observed for ǫ ∈ [0, 0.05]. (b) Enlargement of figure (a) for ǫ ∈ ]0; 0.002]. (c) Enlargement of figure (a) for ǫ ∈ [0.005; 0.015]. (d) Enlargement of figure (c) for ǫ ∈ [0.0138; 0.0147].

Figure 6 :

 6 Figure 6: Time series and (x, y, z) phase portraits of system (2) with parameters given in (3) and I = 3.25 for (a 1 , a 2 ) ǫ = 0.0002, (b 1 , b2) ǫ = 0.0005, (c 1 , c 2 ) ǫ = 0.001, (d 1 , d 2 ) ǫ = 0.002. According to the chosen parameter ǫ, the neuron behaviour changes from tonic spiking (a 1 , a 2 ) to bursting ((b1, b 2 ), (c 1 , c 2 ), (d 1 , d 2 )).

Figure 7 :

 7 Figure 7: System (2) with parameter given in (3) and I = 3.25. (a) Successive local maxima in z (i.e. 'Lorenz plot') for ǫ = 0.008. It produces a typical tent first-return map, which is a numerical signature of chaos. (b) Larger Lyapunov exponent for ǫ ∈ [0.004; 0.016].

Figure 8 :

 8 Figure 8: Time series and (x, y, z) phase portraits of system (2) with parameters given in (3) and with I = 3.25, for (a 1 , a 2 ) ǫ = 0.008, (b 1 , b 2 ) ǫ = 0.0095, (c 1 , c 2 ) ǫ = 0.01, (d 1 , d 2 ) ǫ = 0.011. We observe a period doubling cascade starting with a period three solution.

Figure 9 :

 9 Figure 9: System (2) with parameters given in (3) and I = 3.25. (a) First return map to maxima in z for ǫ = 0.0145. It produces a typical tent first-return map, which is a chaos numerical signature, (b) Larger Lyapunov exponent for ǫ ∈ [0; 0.05]. Since it is positive for ǫ = 0.0145, it also confirms numerically the presence of chaos for these values of parameters.

Figure 10 :

 10 Figure 10: Phase portraits illustrating the bursting bifurcation phenomena in the HR model with parameters given in (3) and I = 3.25. (a 1 , a 2 ) ǫ = 0.0145, (b 1 , b 2 ) ǫ = 0.017, (c 1 , c 2 ) ǫ = 0.03, (d 1 , d 2 ) ǫ = 0.05. As this slow parameter increases, the period, which corresponds to the number of spikes per burst, decreases. A classical period doubling cascade is observed.

Figure 11 (

 11 a) shows the bifurcation diagram according to the bifurcation parameter I for I ∈ [1.25, 4] while figure 11(b) is an enlargement of figure 11(a) for I ∈ [3.25, 3.3]. Figure 12(a) shows that the motion of this system is not chaotic for I = 3.25. This observation is confirmed by the larger Lyapunov exponent presented in figure 12(b) which remains close to zero for I ∈ [3, 3.5].

Figure 11 :

 11 Figure 11: (a) Bifurcation diagram of the HR model for parameters (3) and ǫ = 0.001. As the magnitude of injected current I increases, the number of branches on the diagram also increases. Biologicaly, the fast dynamics of the neuron is evolving. (b) Enlargements of (a) for I ∈ [3.25; 3.3].

Figure 12 :

 12 Figure 12: System (2) with parameters given in (3) and ǫ = 0.001. (a) Successive local maxima in z (i.e. 'Lorenz plot') for I = 3.25. Since this Lorenz plot has not the classical shape of unimodal chaotic map, for these parameters, HR model does not exhibit chaotic motion, (b) Larger Lyapunov exponent for I ∈ [3; 3.5]. Obviously, for I ≈ 3.25, the larger Lyapunov exponent is close to zero.

Figure 13 :

 13 Figure 13: Time series and (x, y, z) phase portraits of system (2) with parameters given in (3), ǫ = 0.001 and (a 1 , a 2 ) I = 1.5, (b 1 , b 2 ) I = 2, (c 1 , c 2 ) I = 3, (d 1 , d 2 ) I = 3.25. The larger the injected current I is (I ∈ [1.4; 3.25]), the larger the number of spikes within a burst is.


Nathalie Corson

Bifurcations of a neuron oscillator

Keywords: Neuronal model, asymptotic dynamics, Hopf bifurcation

This work adresses the study of the threedimensional autonomous ordinary differential equations Hindmarsh-Rose neuronal model. General bifurcation diagrams are first given after a brief presentation of the model. Then, the existence of a Hopf bifurcation according to a small parameter which corresponds to the ratio of time scales between the fast and the slow dynamics is proved. Using the Hassard method we show that, under some conditions, a Hopf bifurcation occurs for a critical value of this parameter. The direction, stability and period of this bifurcation are also discussed. Numerical simulations are done to observe this bifurcation and to illustrate theoretical results.

I. INTRODUCTION

In 1952, a mathematical model that describes neuron activity has been given by two neurophysiologists, A.L. Hodgkin and A.F.Huxley, see [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]. Different neuron models have been then developped and studied, see for example [START_REF] Izhikevich | Which model to use for cortical spiking neurons[END_REF], [START_REF] Izhikevich | Dynamical systems in Neuroscience[END_REF], [START_REF] Morris | Voltage oscillations in the barnacle giant muscle fiber[END_REF] and references therein cited. In this paper, we focus on one of them, the Hindmarsh-Rose model (HR), which results from a simplification and a generalization of the Hodgkin-Huxley model, see [START_REF] Hindmarsh | A model of the nerve impulse using two first-order differential equations[END_REF], [START_REF] Hindmarsh | A model of neuronal bursting using three coupled first order differential equations[END_REF]. As observed in various biology systems, neuron activity presents different time scales. This can be explicitly observed in HR, which is a slow-fast autonomous three ordinary differential equations. The two first equations control the fast dynamics while the third one controls the slow one. Besides, periodic phenomena or oscillations are observed as in many natural systems such as neuron models. Those phenomena can be closely related to Hopf bifurcation.

The HR model reads as follows, (HR)

   ẋ = y + ax 2 -x 3 -z + I ẏ = 1 -dx 2 -y ż = ǫ(b(x -c x ) -z) (1)
Parameters a, b and d are experimentally determined, c x is the equilibrium x-coordinate of the two-dimensional system given by the two first equations of (1) when I = 0 and z = 0 and parameter I corresponds to the applied current. It is easy to experimentaly change its value and it is therefore often used as the bifurcation parameter. Indeed, in the next part, bifrucation diagrams according to I are presented. Finally, parameter ǫ represents the ratio of time scales between fast and slow fluxes accross the membrane of a neuron and, therefore, plays a very special role in neuron activity. It is chosen, in nathalie.corson@univ-lehavre.fr Laboratoire de Mathématiques Appliquées du Havre, 25 rue Philippe Lebon, BP 540, 76058 Le Havre Cedex, France this paper, as the bifurcation parameter, as in [START_REF] Corson | Asymptotic dynamics of the slow-fast Hindmarsh-Rose neuronal system[END_REF] or in [START_REF] Gonzalez-Miranda | Complex bifurcation structures in the Hindmarsh-Rose neuron model[END_REF], in which numerical simulations are done, among other, to study this system according to parameter ǫ. In the last section of this paper, parameters a, b, d and c x are fixed as follows,

a = 3, b = 4, d = 5, c x = - 1 2 (1 + √ 5). (2)
Equilibria are given by ẋ = ẏ = ż = 0, that is to say by,

x 3 + (d -a)x 2 + bx -bc x -I -1 = 0 (3)
Let us denote,

x = ξ + a -d 3 , p = b - (a -d) 2 3 q = - 2(a -d) 3 27 + b(a -d) 3 -bc x -I -1 (4)
Then, (3) reads as, ξ 3 + pξ -q = 0. Solving this equation gives the equilibria of system (1).

Proposition 1.

With notations (4), if 4p 3 + 27q 2 > 0, then system (1) has a unique equilibrium S e = (x e , y e , z e ) given by,

                   x e = - q 2 + q 2 4 + p 3 27
1 2 1 3 + - q 2 - q 2 4 + p 3 27 1 2 1 3 + a -d 3 y e = 1 -dx e z e = b(x e -c x ) (5)
In the next section, a presentation of some bifurcation diagrams of the Hindmarsh-Rose model according to parameter I and parameter ǫ is done. Then the existence of a Hopf bifurcation according to parameter ǫ is studied. Indeed, even if this HR model dates from 1984 and has been widely numerically studied, see for example [START_REF] Arena | Locally active Hindmarsh-Rose neurons[END_REF], [START_REF] Corson | Asymptotic dynamics of the slow-fast Hindmarsh-Rose neuronal system[END_REF], [START_REF] Duan | Codimension-two bifurcation analysis in Hindmarsh-Rose model with two parameters[END_REF], [START_REF] Gonzalez-Miranda | Complex bifurcation structures in the Hindmarsh-Rose neuron model[END_REF], [START_REF] Innocenti | Dynamical phases of the Hindmarsh-Rose neuronal model: Studies of the transition from bursting to spiking chaos[END_REF], no theoretical proof has ever been published as far as we know.

II. BIFURCATION DIAGRAMS

A bifurcation diagram shows the evolution of the aysmptotic behaviour of solutions according to one parameter.

Parameter I corresponds to the current which is injected in the neuron. Thus, it can be controlled during experiments and can then play the role of bifurcation parameter. 2(c) exhibits not only inverse period doubling cascades starting with period 3, period 4 or period 5 but also some dark parts, which is a numerical sign of chaotic motion. Of course, this argument is not sufficient to clame that this system is chaotic for some given ranges of parameters. A more acurate study is done, for example, in [START_REF] Corson | Asymptotic dynamics of the slow-fast Hindmarsh-Rose neuronal system[END_REF]. The enlargement of figure 2(c) for ǫ ∈ [0.0138, 0.0148] shown in figure 2(d) also exhibits a chaotic behaviour of system (1). The right part of figure 2(a) exhibits a reverse period doubling cascade. As ǫ becomes larger, the number of spikes within a burst decreases until the bursting motion of the neuron disappears to let the spiking motion arises.

III. EXISTENCE, DIRECTION, STABILITY AND PERIOD OF A HOPF BIFURCATION ACCORDING TO ǫ

In this section, the existence, direction, stability and period of a Hopf bifurcation according to ǫ is studied, see [START_REF] Corson | Hopf bifurcation in the slowfast Hindmarsh-Rose neuronal system[END_REF]. Under the coordinate transformation, x 1 = x-x e , y 1 = y -y e and z 1 = z -z e , system (1) becomes,

       ẋ1 = (2ax e -3x 2 e )x 1 + y 1 -z 1 + F1 (x 1 , y 1 , z 1 ) ẏ1 = -2dx e x 1 -y 1 + F2 (x 1 , y 1 , z 1 ) ż1 = ǫbx 1 -ǫz 1 + F3 (x 1 , y 1 , z 1 ) (6)
where Fj (x 1 , y 1 , z 1 ), j = 1, 2, 3 are the nonlinear terms.

The Poincaré-Andronov-Hopf theorem applied to system (1) leads to the following proposition, Proposition 2. With notations ( 12) and ( 13), if the two following conditions hold,

4r 3 + 27s 2 > 0 (7) 2 3 (a -d) < x e < 0 (8)
then, when parameter ǫ passes the value ǫ c , system (1) undergoes a Hopf bifuration at the equilibrium S e , where, Proof: The existence of a Hopf bifurcation point in system (1) is studied using the linearized system (6) at S e . First of all, its jacobian matrix M (ǫ) is,

ǫ c = -(1 -m 11 ) 2 -m 11 b + ∆ 1 2 2(1 -m 11 + b) (9
M (ǫ) = (m ij ) 1≤i,j≤3 (10)
The corresponding characteristic equation is,

f (λ(ǫ)) = λ 3 (ǫ) + P (ǫ)λ 2 (ǫ) + Q(ǫ)λ(ǫ) + R(ǫ) (11)
where,

P (ǫ) = 1 -m 11 + ǫ Q(ǫ) = (1 -m 11 + b)ǫ -m 11 -m 21 R(ǫ) = ǫ(b -m 11 -m 21 ) (12)
Setting,

λ(ǫ) = ν(ǫ) -P (ǫ) 3 r(ǫ) = Q(ǫ) - P 2 (ǫ) 3 s(ǫ) = 2P 3 (ǫ) 27 - P (ǫ)Q(ǫ) 3 + R(ǫ) (13)
equation ( 11) reads as,

ν 3 (ǫ) + r(ǫ)ν(ǫ) + s(ǫ) = 0,
which is the equation giving M (ǫ) eigenvalues.

The sign of 4r 3 (ǫ) + 27s 2 (ǫ) provides the number of real and complex eigenvalues of this matrix. Indeed, if 4r 3 (ǫ) + 27s 2 (ǫ) > 0, that is if condition [START_REF] Hassard | Theory and Applications of Hopf bifurcation[END_REF] holds, then M (ǫ) has two complex eigenvalues λ 1,2 (ǫ) = α(ǫ) + iω(ǫ) and one real, λ 3 (ǫ).

Now, let us study the existence of a critical value ǫ c of parameter ǫ. From ( 10), [START_REF] Innocenti | Dynamical phases of the Hindmarsh-Rose neuronal model: Studies of the transition from bursting to spiking chaos[END_REF] and [START_REF] Izhikevich | Which model to use for cortical spiking neurons[END_REF], polynomial rules lead to the existence of,

ǫ c = -(1 -m 11 ) 2 -m 11 b ± ∆ 1 2 2(1 -m 11 + b) .
Algebraic computations show that under condition (8), ǫ c > 0.

Moreover, since x e < 0, it is obvious that m 11 < 0 and thus, P (ǫ c ) > 0. Therefore, λ 3 (ǫ) < 0.

The derivative according to ǫ of the characteristic equation given in [START_REF] Innocenti | Dynamical phases of the Hindmarsh-Rose neuronal model: Studies of the transition from bursting to spiking chaos[END_REF]

is, ∂f (ǫ) ∂ǫ = 3λ 2 (ǫ) ∂λ(ǫ) ∂ǫ + ∂P (ǫ) ∂ǫ λ 2 (ǫ) +2P (ǫ)λ(ǫ) ∂λ(ǫ) ∂ǫ + ∂Q(ǫ) ∂ǫ λ(ǫ) +Q(ǫ) ∂λ(ǫ) ∂ǫ + ∂R(ǫ) ∂ǫ . (14)
Therefore, solving ∂f ∂ǫ (ǫ c ) = 0 and separating imaginary and real parts, we obtain,

∂α ∂ǫ (ǫ c ) = ∂R ∂ǫ (ǫ c ) - ∂P ∂ǫ (ǫ c )Q(ǫ c ) -P (ǫ c ) ∂Q ∂ǫ (ǫ c ) 2Q(ǫ c ) + 2P (ǫ c ) 2 . Since 2Q(ǫ c ) + 2P (ǫ c ) 2 > 0, ∂α ∂ǫ (ǫ c ) < 0.
Finally, if 4r 3 (ǫ) + 27s 2 (ǫ) > 0 and 2 3 (a -d) < x e < 0, then all the conditions of the Poincaré-Andronov-Hopf theorem hold and (S e , ǫ c ) is a Hopf bifurcation point of system [START_REF] Arena | Locally active Hindmarsh-Rose neurons[END_REF].

Let us now study direction, stability and period of this Hopf bifurcation occuring at ǫ c using Hassard method, see [START_REF] Hassard | Theory and Applications of Hopf bifurcation[END_REF] and see also [START_REF] Hassard | Bifurcation formulae derived from center manifold theory[END_REF], [START_REF] Lü | Local bifurcation of the Chen system[END_REF], [START_REF] Zhou | Hopf bifurcation of the Liu system[END_REF] .

Let us denote by ω 0 the value ω(ǫ c ) > 0 and let v j , j = 1, 2, 3, be the eigenvectors of the matrix M (ǫ c ), given in [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF], corresponding to the eigenvalues λ j . We have, λ j = ±iω 0 = ±iQ 1/2 (ǫ c ), j = 1, 2, and λ 3 = -P (ǫ c ).

The eigenvector v 1 associated with λ 1 = iω 0 is,

v 1 = 1, m 21 (1 -iω 0 ) 1 + ω 2 0 , ǫ c b(ǫ c -iω 0 ) ǫ 2 c + ω 2 0 T ,
and the eigenvector v 3 associated with λ 3 is,

v 3 = 1, m 21 m 11 -ǫ c , ǫ c b m 11 -1 T . Let us define P such that (x 1 , y 1 , z 1 ) T = [P (x 2 , y 2 , z 2 )] T , P = (Re(v 1 ), -Im(v 1 ), v 3 ) = (p ij ) 1≤i,j≤3 (15)
The inverse matrix is given by P

-1 = (p -1 ij ) 1≤i,j≤3 . Thus,    ẋ2 = ω 0 y 2 + F 1 (x 2 , y 2 , z 2 ), ẏ2 = -ω 0 x 2 + F 2 (x 2 , y 2 , z 2 ), ż2 = λ 3 z 2 + F 3 (x 2 , y 2 , z 2 ),
where F 1 , F 2 and F 3 are the nonlinear terms, satisfying F i (x 2 , y 2 , z 2 ) = P -1 ( Fi (x 1 , y 1 , z 1 )), Procedures proposed by Hassard et al. [START_REF] Hassard | Theory and Applications of Hopf bifurcation[END_REF] are used to calculate the following quantities, evaluated at ǫ = ǫ c .

g 11 = 1 2 p -1 11 (a -3x e ) -p -1 12 d + i 2 p -1 21 (a -3x e ) -p -1 22 d g 02 = g 11 g 20 = g 11 G 21 = - 3 4 p -1 11 + ip -1 21 .
Moreover, let us calculate the quantities,

h 11 = 1 2 p -1 31 (a -3x e ) -p -1 32 d h 20 = h 11 .
Then, solving the two equations,

λ 3 w 1 = -h 11 (λ 3 -2iω 0 )w 20 = -h 20
gives,

w 11 = p -1 31 (a -3x e ) -p -1 32 d 2(1 -m 11 + ǫ c ) , w 20 = 1 2 p -1 31 (a -3x e ) -p -1 32 d (1 -m 11 + ǫ) 2 + 4ω 2 0 . (1 -m 11 + ǫ c + 2iω 0 )
Furthermore, calculating the quantities,

G 110 = 1 4 p -1 11 (a -3x e ) -p -1 12 d + i 4 p -1 21 (a -3x e ) -p -1 22 d G 101 = G 110 g 21 = G 21 + (2G 110 w 11 + G 101 w 20 )
and by setting,

c 1 = i 2ω 0 g 20 g 11 -2|g 11 | 2 - 1 3 |g 02 | 2 + 1 2 g 21
we can give the main result below in which,

µ 2 = - Re(c 1 ) ∂α ∂ǫ (ǫ c ) , τ 2 = - Im(c 1 ) + µ 2 ∂ω ∂ǫ (ǫ c ) ω 0 , β 2 = 2Re(c 1 ).
Theorem 1. Under the hypothesis of proposition 2, system (1) undergoes a Hopf bifurcation at the equilibrium point (x e , y e , z e ) as ǫ passes through ǫ c with the following properties.

1) If µ 2 < 0 (reps. µ 2 > 0) , then the direction of bifurcation is ǫ < ǫ c (resp. ǫ > ǫ c ) and the bifurcation is supercritical (resp. subcritical), 2) If β 2 < 0 (resp. β 2 > 0), the bifurcating periodic solutions are orbitally stable (resp. unstable), 3) If τ 2 > 0 (resp. τ 2 < 0), the period of bifurcating periodic solutions increases (resp. decreases).

The period and characteristic exponents are given by,

T = 2π ω 0 (1 + τ 2 E 2 + O(E 4 )) β = β 2 E 2 + O(E 4 )
Where,

E 2 = ǫ -ǫ c µ 2 + O(ǫ -ǫ c ) 2 (provided µ 2 = 0).
The periodic solutions themselves are,

  x y z   =   x e y e z e   + P.   u 1 u 2 u 3   (16)
where,

u 1 = Re(ζ) , u 2 = Im(ζ) u 3 = w 11 |ζ| 2 + Re(w 20 ζ 2 ) + O(|ζ| 3 ) and ζ = Ee 2iπt/T + iE 2 6ω 0 (g 02 e -4iπt/T
-3g 20 e 4iπt/T + 6g 11 ) + O(E 3 )

Now, numerical computations are done to illustrate these theoretical results. Thereafter, we consider system (1), with parameters a, b, c x and d fixed as in (2) and I = 3.25. Equilibria of this system are studied as presented in the first section of this paper, and using notation of (4), we obtain,

4p 3 + 27q 2 ≈ 76.443755 > 0
Therefore, proposition 1 leads to the existence an uniqueness of system (1) equilibrium (x e , y e , z e ), given by ( 5),

x e ≈ -0.722126 , y e ≈ -1.607329 , z e ≈ 3.583632

Let us verify if proposition 2 can be applied to system [START_REF] Arena | Locally active Hindmarsh-Rose neurons[END_REF] with the fixed values of parameters [START_REF] Corson | Asymptotic dynamics of the slow-fast Hindmarsh-Rose neuronal system[END_REF]. The bifurcation value of ǫ, given by equation ( 9) of proposition 2, is, ǫ c ≈ 0.125912.

This value of ǫ c is really close to the one we observe on the bifurcation diagram given in figure 3.

For this value of ǫ c , we have, 4r 3 (ǫ c ) + 27s 2 (ǫ c ) ≈ 443.299666 > 0 and condition [START_REF] Hassard | Theory and Applications of Hopf bifurcation[END_REF] holds. Therefore, the jacobian matrix M (ǫ c ) has one real eigenvalue λ 3 (ǫ c ) and two complex ones λ 1,2 = α(ǫ c ) ± iω(ǫ c ).

Since 2(a -d)/3 = -4/3, the condition 2(a -d)/3 < x e < 0 is verified. Furthermore, ∂α ∂ǫ (ǫ c ) ≈ -0.748444 = 0, and λ 3 (ǫ c ) ≈ -7.025406 < 0.

Thus, thanks to proposition 2, ǫ c ≈ 0.125912 is the Hopf bifurcation value of parameter ǫ for system (1) at (-0.722126, -1.607329, 3.583632). with parameter fixed as in (2) and ǫ = 0.12 < ǫc, the asymptotic solution is a stable limit cycle. (b) (x, y, z) view of the phase portrait and time series of system (1) with parameter fixed as in (2) and ǫ = 0.13 > ǫc, the asymptotic solution is stable focus.

The computation of matrix P and its inverse matrice P -1 gives, = -0.063437 -0.009879i

P =   1
Finally, computations give, µ 2 = -0.084758 < 0 , β 2 = -0.126873 < 0 , τ 2 = 0.384677 > 0.

According to theorem 1, the Hopf bifurcation occuring at ǫ c is supercritical and the direction of bifurcation is ǫ < ǫ c . Moreover, the bifurcating periodic solutions are asymptotically orbitally stable and the period of bifurcating periodic solutions increases.

The period of the solution is given by, T = 28.686363 -130.193932(ǫ -0.125912) +O((ǫ -0.125912) 2 ) and this period increases as ǫ decreases. The characteristic exponents is given by,

IV. CONCLUSION

In this paper, after a presentation of the three-dimensional autonomous ordinary differential equations Hindmarsh-Rose neuronal model, bifurcation diagrams according to parameters I and ǫ are presented. Then, the existence of a Hopf bifurcation according to parameter ǫ in this model is discussed. Indeed, for a critical value ǫ c of this parameter, a Hopf bifurcation occurs under some conditions. Using Hassard algorithm, the direction, stability and period of this bifurcation are then studied. Finally, numerical simulations are done to observe this bifurcation and to illustrate theoretical results.

Figure 1 .

 1 Figure 1. (a) Bifurcation diagram of the HR model for parameters (2) and ǫ = 0.001. As the magnitude of injected current I increases, the number of branches on the diagram also increases. Biologicaly, the fast dynamics of the neuron is evolving. (b) Enlargements of (a) for I ∈ [3.25; 3.3].

Figure 2 (

 2 Figure 2(a) gives the bifurcation diagram with respect to the control parameter ǫ in the range [0, 0.05]. In order to have a more accurate analysis of the dynamics of system (1), we present in figure 2(b),(c),(d) enlargements of figure 2(a).Figure 2(b) shows, among other things, that there is an ǫ 1 ∈ [0.00041, 0.00049] for which the neuron behaviour changes abruptly. Indeed, ∀ǫ < ǫ 1 , the neuron exhibits a tonic spiking motion and, ∀ǫ > ǫ 1 , the neuron exhibits a bursting motion. Moreover, this figure shows that system (1) with parameters given in (2) and I = 3.25 does not exhibit chaotic behaviour for ǫ ∈ [ǫ 1 , 0.002]. The enlargement of figure 2(a) for ǫ ∈ [0.005, 0.015] shown in figure2(c) exhibits not only inverse period doubling cascades starting with period 3, period 4 or period 5 but also some dark parts, which is a numerical sign of chaotic motion. Of course, this argument is not sufficient to clame that this system is chaotic for some given ranges of parameters. A more acurate study is done, for example, in[START_REF] Corson | Asymptotic dynamics of the slow-fast Hindmarsh-Rose neuronal system[END_REF]. The enlargement of figure2(c) for ǫ ∈ [0.0138, 0.0148] shown in figure2(d) also exhibits a chaotic behaviour of system (1). The right part of figure2(a) exhibits a reverse period doubling cascade. As ǫ becomes larger, the number of spikes within a burst decreases until the bursting motion of the neuron disappears to let the spiking motion arises.

  Figure 2(a) gives the bifurcation diagram with respect to the control parameter ǫ in the range [0, 0.05]. In order to have a more accurate analysis of the dynamics of system (1), we present in figure 2(b),(c),(d) enlargements of figure 2(a).Figure 2(b) shows, among other things, that there is an ǫ 1 ∈ [0.00041, 0.00049] for which the neuron behaviour changes abruptly. Indeed, ∀ǫ < ǫ 1 , the neuron exhibits a tonic spiking motion and, ∀ǫ > ǫ 1 , the neuron exhibits a bursting motion. Moreover, this figure shows that system (1) with parameters given in (2) and I = 3.25 does not exhibit chaotic behaviour for ǫ ∈ [ǫ 1 , 0.002]. The enlargement of figure 2(a) for ǫ ∈ [0.005, 0.015] shown in figure2(c) exhibits not only inverse period doubling cascades starting with period 3, period 4 or period 5 but also some dark parts, which is a numerical sign of chaotic motion. Of course, this argument is not sufficient to clame that this system is chaotic for some given ranges of parameters. A more acurate study is done, for example, in[START_REF] Corson | Asymptotic dynamics of the slow-fast Hindmarsh-Rose neuronal system[END_REF]. The enlargement of figure2(c) for ǫ ∈ [0.0138, 0.0148] shown in figure2(d) also exhibits a chaotic behaviour of system (1). The right part of figure2(a) exhibits a reverse period doubling cascade. As ǫ becomes larger, the number of spikes within a burst decreases until the bursting motion of the neuron disappears to let the spiking motion arises.

Figure 2 .

 2 Figure 2. Bifurcation diagrams in (ǫ, y) plane for system (1) with parameters given in (2) and with I = 3.25. (a) An inverse period doubling cascade is observed for ǫ ∈ [0, 0.05]. (b) Enlargement of figure (a) for ǫ ∈ ]0; 0.002]. (c) Enlargement of figure (a) for ǫ ∈ [0.005; 0.015]. (d) Enlargement of figure (c) for ǫ ∈ [0.0138; 0.0147].

) and m 11 =

 11 2ax e -3x 2 e , m 21 = -2dx e and ∆ = [(1 -m 11 ) 2 -m 11 b] 2 + 4(1 -m 11 + b)(m 11 + m 21 )(1 -m 11 ).

Figure 3 .

 3 Figure 3. Bifurcation diagram of system (1) with parameters given in (2) according to parameter ǫ.

Figure 4 .

 4 Figure 4. (a) (x, y, z) view of the phase portrait and time series of system (1) with parameter fixed as in (2) and ǫ = 0.12 < ǫc, the asymptotic solution is a stable limit cycle. (b) (x, y, z) view of the phase portrait and time series of system (1) with parameter fixed as in (2) and ǫ = 0.13 > ǫc, the asymptotic solution is stable focus.

β = 1 . 2 )u 1 =

 121 496888(ǫ -0.125912) + 0((ǫ -0.125912) The periodic solutions are, Re(ζ), u 2 = Im(ζ), u 3 = 0.359215|ζ| 2 + Re((0.357823 + 0.022319i)ζ 2 ) + O(|zη| 3 ) and, ζ = Ee 2iπt/T + iE 2 6ω 0 (0.060399 + 0.071870i)e -4iπt/T -3(0.060399 + 0.071870i)e 4iπt/T +6(0.060399 + 0.071870i) + O(E 3 )