Guillemette Chapuisat
email: guillemette.chapuisat@univ-cezanne.fr

Romain Joly
email: romain.joly@ujf-grenoble.fr

Asymptotic profiles for a travelling front solution of a biological equation

Keywords: spreading depression, reaction-diffusion equation, travelling fronts, Sturm-Liouville theory AMS classification codes (2000): 34C10, 35B35, 35K57, 92C20

We are interested in the existence of depolarization waves in the human brain. These waves propagate in the grey matter and are absorbed in the white matter. We consider a two-dimensional model u t = ∆u + f (u)½ |y|≤R -αu½ |y|>R , with f a bistable nonlinearity taking effect only on the domain R × [-R, R], which represents the grey matter layer. We study the existence, the stability and the energy of nontrivial asymptotic profiles of possible travelling fronts. For this purpose, we present dynamical systems technics and graphic criteria based on Sturm-Liouville theory and apply them to the above equation. This yields three different behaviours of the solution u after stimulation, depending of the thickness R of the grey matter. This may partly explain the difficulties to observe depolarization waves in the human brain and the failure of several therapeutic trials.

Introduction

The propagation of depolarization waves, also called spreading depressions, may appear in a brain during strokes, migraines with aura or epilepsy. In rodent brains, the propagation of depolarization waves has been observed for more than fifty years [START_REF] Leão | Spreading depression of activity in the cerebral cortex[END_REF]. During stroke in rodent brain, they cause important damages and are therefore a therapeutic target. Pharmacological agents blocking the appearance of those waves have been studied and reduce strongly after-effects of stroke in the rodent brain [START_REF] Mies | Correlation between peri-infarct DC shifts and ischaemic neuronal damage in rat[END_REF][START_REF] Obeidat | Glutamate does not mediate acute neuronal damage after spreading depression induced by O2/glucose deprivation in the hippocampal slice[END_REF]. However propagation of depolarization waves during stroke in the human brain is still uncertain [START_REF] Aitken | Spreading depression in human hippocampal tissue in vitro[END_REF][START_REF] Mayevsky | Cortical spreading depression recorded from the human brain using a multiparametric monitoring system[END_REF][START_REF] Maclachlan | Spreading depression of Leao in rodent and human cortex[END_REF][START_REF] Somjen | Ions in the brain: Normal Function, Seizures, and Stroke[END_REF][START_REF] Sramka | Functional ablation by spreading depression: possible use in human stereptactic surgery[END_REF] and the pharmalogical agents used in rodent have seemed to have no effect on human stroke. Mathematical models of these depolarization waves may help to understand these points.

In [START_REF] Chapuisat | Discussion of a simple model of spreading depression[END_REF], the first author has build a mathematical model of these waves and has studied numerically how the morphology of the brain may influence their propagation. Simplifications of the mathematical model lead to the study of propagation phenomena for the following biological equation

∂u ∂t -△u = f (u)½ Ω -αu½ R N \Ω , t ∈ R, X ∈ R N , (1.1)
where f (u) = λ u(ua)(1u) is the usual bistable nonlinearity with a ∈]0, 1/2[ and λ > 0, and where α is a positive number. We denote by ½ Ω the caracteristic function of the domain Ω that is ½ Ω (X) = 1 if X ∈ Ω and ½ Ω (X) = 0 elsewhere. The function u represents the depolarization of the brain that is that if u(X) = 0 the brain is normally polarized at the point X whereas if u(X) = 1 the brain is totally depolarized. The domain Ω represents the grey matter of a brain, where the reaction-diffusion process that triggers depolarization waves takes place, whereas R N \ Ω contains the white matter, where the waves are absorbed. The problem is to understand the influence of the geometry of Ω on the propagation of waves. In particular, the layer of grey matter in the human brain is thiner than in other species and admits more circumvolutions. The influence of the circumvolutions of the grey matter is partially studied in [START_REF] Chapuisat | Existence and nonexistence of traveling wave solutions for a bistable reaction-diffusion equation in an infinite cylinder whose diameter is suddenly increased[END_REF].

We will focus here on the influence of the thickness of the layer of grey matter. We may assume that Ω is a straight cylinder of radius 2R. In [START_REF] Chapuisat | Discussion of a simple model of spreading depression[END_REF], numerical studies have shown that small values of R may prevent the spreading of the ionic waves. A partial theoretical study in any dimension N has been led in [START_REF] Chapuisat | Existence and non-existence of curved front solution of a biological equation[END_REF], where the first author has proved the non-existence of the depolarization waves if the thickness R is small enough and their existence if R is large enough.

The results of [START_REF] Chapuisat | Existence and non-existence of curved front solution of a biological equation[END_REF] are not completely satisfactory. First, they only deal with R small or large enough. Secondly, no complete description of the possible asymptotic profiles of the travelling waves, their stability, or their energy, has been pursued. As a consequence, even if the existence of waves is proved, nothing is known about their asymptotic profiles and their stability. The purpose of the present paper is to complete this study for any R, in dimension N = 2. More precisely, we consider the equation

∂u ∂t = △u + f (u)½ |y|≤R -αu½ |y|>R , t ∈ R, (x, y) ∈ R 2 . (1.2)
We are looking to solutions travelling in the x direction at speed c, that are solutions u of (1.2) which can be written u(x, y, t) = v(xct, y). Travelling fronts are solutions u(x, y, t) = v(xct, y) such that there are two asymptotic profiles

V -= V + with lim ξ→±∞ v(ξ, .)-V ± H 1 (R) = 0. Using standard elliptic estimates, a profile V is a solution in H 1 (R) of the equation V ′′ + f (V )½ |y|≤R -αV ½ |y|>R = 0, y ∈ R . (1.3)
Notice that a profile is trivially associated to an equilibrium point E(x, y) = V (y) that is a stationnary solution of (1.2). The trivial equilibrium point E ≡ 0 corresponds to the normal state of the brain, whereas a non-trivial equilibrium point corresponds to a depolarized state, deleterious for the brain. We are more particularly interested in the existence of travelling fronts with positive speed c connecting the profile V + ≡ 0 with a non-trivial profile V -. Such fronts correspond to the invasion of the equilibrium state 0 by a deleterious state. The purpose of this article is to achieve a complete study of the existence of the profiles, their stability and their energy (see the following sections for the definitions of stability and energy).

Main result 1.1. There exist two critical thicknesses R 1 > R 0 > 0 such that: i) if 0 < R < R 0 , there is no non-trivial profile, i.e. that V ≡ 0 is the only solution of

(1.3) in H 1 (R).
ii) If R 0 < R, there exist non-trivial profiles. One of them, denoted by V M , is larger than every other one. The largest profile V M is stable, and every other non-trivial profiles are unstable.

iii) The energy of the unstable profiles is always larger than the energy of the stable profiles 0 and

V M . If R 0 < R < R 1 , the energy of V M is larger than the energy of 0, whereas it is smaller if R > R 1 .
In this paper, we first adapt several graphic criteria for stability of profiles, based on Sturm-Liouville theory. Then, we present dynamical systems technics for studying energy of profiles. These criteria and technics are general ones and are valid for any nonlinearity f . We obtain the above results by applying these ideas to the particular nonlinearity f (u) = λ u(ua)(1u). The check on the graphic criteria sometimes relies on a numerical study of the phase plane.

Once the profiles and their stability completely described by the Main Result 1.1, one is able to obtain the existence of travelling fronts. Consequences 1.2. i) If 0 < R < R 0 then there is no travelling front for Equation (1.2). ii) If R 0 < R < R 1 there exists a globally stable travelling front u(x, y, t) = v(xct, y) with lim ξ→+∞ v(ξ, .) = 0 and lim ξ→-∞ v(ξ, .) = V M with a negative speed c < 0. ii) If R 1 < R there exists a globally stable travelling front u(x, y, t) = v(xct, y) with lim ξ→+∞ v(ξ, .) = 0 and lim ξ→-∞ v(ξ, .) = V M with a positive speed c > 0.

Assertion i) is a clear consequence of the Main Result 1.1. The proof that this Main Result also implies the other assertions is not the goal of this paper and will not be detailled. One way to prove it is to use a method recently introduced by Risler in [START_REF] Risler | Global convergence towards travelling fronts in nonlinear parabolic systems with a gradient structure[END_REF]. The fundamental idea is to use the existence of an energy functional for (1.2) in any Galilean frame travelling at constant speed s in the x direction. Keeping this basic idea, the original proof has been revised in [START_REF] Th | A variational proof of global stability for bistable travelling waves[END_REF] and [START_REF] Th | Global stability of travelling fronts for a damped wave equation with bistable nonlinearity[END_REF]. The adaptation of these technics to Equation (1.2) has been detailled in [START_REF] Chapuisat | Existence and non-existence of curved front solution of a biological equation[END_REF]. The original motivation of Risler's technics was to develop a method of proof of existence and stability of fronts in equations where no comparison principle is available. As suggested by the work of the first author in [START_REF] Chapuisat | Existence and non-existence of curved front solution of a biological equation[END_REF], the method of Risler is also usefull to get the existence and stability of fronts in equations as (1.2) for which even if there is a comparison principle, no particular positive subsolution is known.

The results of this article yield a good insight into the different behaviours of the depolarization in the grey matter after it has been stimulated. If the stimulation takes place in a part of a brain where the grey matter is thin (R < R 0 ), the neurons of the grey matter quickly repolarize, i.e. the solution of (1.2) goes back to zero, uniformly and exponentially fast. If the grey matter is slightly thicker (R 0 < R < R 1 ), the repolarization is slower and more progressive, the depolarized area disappearing by shrinking. Mathematically speaking, there exists a stable non-trivial profile but with an energy larger than zero. Thus, the equilibrium state zero invades the excited state by travelling fronts, reducing the excited area. Notice that it may take very long time to go back to rest if the initial excited area is large. Finally, if the grey matter is thick (R 1 < R), a depolarization wave propagates. Indeed, the excited state has a lower energy than zero. These three behaviours are illustrated in Figure 1. The fact that the behaviour depends on the width of the grey matter may explain why the depolarization waves have been observed or not in the human brain depending on the experiments, see [START_REF] Aitken | Spreading depression in human hippocampal tissue in vitro[END_REF][START_REF] Mayevsky | Cortical spreading depression recorded from the human brain using a multiparametric monitoring system[END_REF][START_REF] Maclachlan | Spreading depression of Leao in rodent and human cortex[END_REF] or [START_REF] Sramka | Functional ablation by spreading depression: possible use in human stereptactic surgery[END_REF] and the discussion of Section 5.

Figure 1: The three typical behaviours of the depolarization in the grey matter after stimulation. From left to right:

R < R 0 , R 0 < R < R 1 and R 1 < R.
The paper is splitted as follows. In Section 2 we explicit a relation between the profiles and the equilibrium points of a parabolic equation in (-R, R). In Section 3, we obtain graphic criteria of existence and stability of the profiles, mainly based on Sturm-Liouville theory. The energy of the different profiles is studied in Section 4, using technics coming from the infinite dimensional dynamical systems theory. Finally, in Section 5, we discuss the relations between our mathematical results and the biological phenomena of spreading waves in the brains.

2 Relations with a parabolic equation on (-R, R)

In this preliminary section, we enhance an obvious relation between the profiles and the equilibrium points of a parabolic equation on the segment (-R, R). We also recall some basic properties of the dynamics of this parabolic equation. In fact, in this paper, we could perform all the arguments with the parabolic equation on the whole domain R. However, this would bring useless technicalities.

Back to a bounded interval

The profiles V are the solutions of Equation (1.3). In other word, they are the stationnary solutions of the evolution equation

∂u ∂t = ∂ 2 u ∂y 2 + g(y, u) (2.1)
where g(y, u) = f (u)½ |y|≤R -αu½ |y|>R . The energy corresponding to this reaction-diffusion equation is the functional E : H 1 (R) -→ R defined by

E(u) = 1 2 R du dy (y) 2 dy - R G(y, u(y))dy (2.2)
where G(y, u) = u 0 g(y, v)dv is a primitive of g.

We notice that any profile

V satisfies V (y) = c + e √ αy + c -e -√ αy outside (-R, R). Since V belongs to H 2 (R), one has c + = 0 for y > R and c -= 0 for y < R. Since H 2 (R) ⊂ C 1 (R), V is continuous and solving (1.3) is equivalent to solving the following problem    V ′′ (y) + f (V (y)) = 0 y ∈ (-R, R) V ′ (-R) = √ α V (-R) V ′ (R) = - √ α V (R) .
(2.3)

Equation (2.
3) caracterises the equilibrium points of a parabolic equation on (-R, R) with Robin boundary conditions

             ∂u ∂t (y, t) = ∂ 2 u ∂y 2 (y, t) + f (u(y, t)) (y, t) ∈ (-R, R) × R + ∂u ∂y (-R, t) = √ α u(-R, t) t > 0 ∂u ∂y (R, t) = - √ α u(R, t) t > 0 (2.4)
The energy corresponding to Equation (2.4) is given by

E(u) = 1 2 R -R du dy (y) 2 dy - R -R F (u(y))dy + √ α 2 |u(R)| 2 + |u(-R)| 2 (2.5)
where

F (u) = u 0 f (v)dv is a primitive of f . It is straightforward to verify that if V is a profile, then E(V ) = E(V |(-R,R)
). As a conclusion, the problems of the existence and the energy of the profiles are equivalent to the study of the equilibrium points of the parabolic equation (2.4) and their energy. We will often abusively denote V |(R,R) by simply V . However, notice that the dynamics of the equation (2.1) on the whole real line and the ones of (2.4) on (-R, R) are different. Moreover, the spectrum of the linearization of (2.1) at a profile V is different from the spectrum of the linearization of (2.4) at V |(-R,R) , even if we prove in Section 3.2 that the number of positive eigenvalues is the same.

Basic facts about the parabolic equation on (-R, R)

The theory of the parabolic equation (2.4) with Robin boundary conditions is very wellknown, see for example [START_REF] Henry | Geometric theory of semilinear parabolic equations[END_REF]. We recall here some basic dynamical properties, which will be used in this paper.

First, (2.4) admits a comparison principle.

Lemma 2.1. Let u(t) and v(t) be two solutions of (2.4) with respective initial data u 0 and v 0 belonging to H 1 (-R, R). Assume that u 0 (y) ≤ v 0 (y) for all y ∈ (-R, R). Then, for all (y, t)

∈ (-R, R) × R + , u(y, t) ≤ v(y, t).
Secondly, the parabolic type of (2.4) and the fact that f (u)u < 0 for large u implies that any solution u(t) with initial data

u 0 ∈ H 1 (-R, R) is bounded in H k (-R, R)
uniformly for t ∈ [t 0 , +∞) for any t 0 > 0 and for any k ∈ N. This yields a compacity result.

Lemma 2.2. Let (t n ) n∈N be a sequence of times such that lim inf t n > 0, let (u 0 n ) n∈N ⊂ H 1 (-R, R) be a sequence of initial data and let u n (t) be the solutions of (2.4) with u n (0) = u 0 n . Then, there exist an extraction ϕ and a function

u ∞ ∈ H 1 (-R, R) such that u ϕ(n) (t ϕ(n) )- u ∞ H 1 (-R,R) -→ 0 when n goes to +∞.
The energy E defined by (2.5) is a strict Lyapounov functional for the dynamics of (2.4): for each u(t) ∈ H 1 (-R, R) solution of (2.4), the function E(u(t)) is decreasing in time, except if u(t) is constant, that is except if u is an equilibrium point, and thus a profile. Indeed, it is straightforward to get

∂ t E(u(t)) = - R -R | ∂u ∂t (y, t)| 2 dt.
The existence of a strict Lyapounov function yields what is called gradient dynamics. In addition of the compacity property of Lemma 2.2, it is very classical to obtain Lasalle principle, see [START_REF] Hale | Asymptotic behavior of dissipative systems[END_REF].

Lemma 2.3. Let u(t), t ≥ 0, be a solution in H 1 (-R, R) of (2.4). Let ω(u) be the ω-limit set of u, that is ω(u) = v ∈ H 1 (-R, R) / ∃(t n ) n∈N , t n ----→ n→+∞ +∞ and u(t n ) -v H 1 (-R,R) ----→ n→+∞ 0 .
Then, ω(u) is a non-empty connected compact set which only consists in equilibrium points, i.e. in profiles. Similarly, if u(t) is defined for all t ∈ R and uniformly bounded, the α-limit set

α(u) = v ∈ H 1 (-R, R) / ∃(t n ) n∈N , t n ----→ n→+∞ -∞ and u(t n ) -v H 1 (-R,R) ----→ n→+∞ 0 is a non-empty connected compact set which only consists in profiles. Moreover, if V -∈ α(u) and V + ∈ ω(u) then E(V -) > E(V + ) and in particular V -≡ V + .

Existence and stability of profiles

We are interested in the asymptotic profiles of the possible travelling fronts of Equation (1.2). We recall that a profile V ∈ H 1 (R) is solution of (1.3). We wonder if there exist non-trivial profiles (i.e. profiles V ≡ 0) and if they are stable. A profile V is stable if the linear operator

L V : H 2 (R) -→ L 2 (R) defined by L V ϕ = ϕ ′′ + f ′ (V )½ |y|≤R -α½ |y|>R ϕ (3.1)
has no spectrum in the half-plane {z ∈ C, Re(z) ≥ 0}.

In this section, we give graphic criteria for existence and stability of profiles and we apply them to Equation (1.3) to obtain the following result. Proposition 3.1. For any given α, λ and a, there exists a positive number R 0 such that: i) if R < R 0 , then there does not exist any profile different from 0. ii) if R ≥ R 0 , then there exist non-trivial profiles. One of them, denoted by V M is strictly larger than every others, smaller than 1 and is even. Moreover, if R > R 0 then V M and 0 are stable profiles and are the only ones to be stable. iii) if α ≥ aλ then if R > R 0 , there exist exactly two non-trivial profiles.

A graphic criterium for existence of profiles

We recall here the simple graphic arguments already mentionned in [START_REF] Chapuisat | Existence and non-existence of curved front solution of a biological equation[END_REF]. As already noticed, looking for profiles V solutions of (1.3) is equivalent to looking for solutions V of the equation (2.3). In the phase plane, the flow Φ y 0 : R 2 -→ R 2 associated to the differential equation V ′′ + f (V ) = 0 between -R and y 0 is given by

Φ y 0 (V 0 , V 1 ) = (V (y 0 ), V ′ (y 0 )) where V solves V ′′ (y) + f (V (y)) = 0 , y ∈ (-R, y 0 ) (V (-R), V ′ (-R)) = (V 0 , V 1 ) Let D and D ′ ⊂ R 2 be the straight lines R.(1, √ α) and R.(1, - √ α). Notice that the conditions V ′ (-R) = √ α V (-R) and V ′ (R) = - √ α V (R) correspond to (V, V ′ )(-R) ∈ D and (V, V ′ )(R) ∈ D ′ respectively.
Hence, the graphic interpretation of (2.3) is the following.

Graphic criterium 3.2. The profiles, i.e. the solutions of (2.3), correspond to the trajectories of the flow Φ y that start on D when y = -R and finish on D ′ when y = R. Therefore, they correspond to the intersections of the curve Φ R D with the straight line D ′ .

The results of existence stated in Proposition 3.1 follow from Graphic Criterium 3.7. The critical radius R 0 is the first (positive) R such that Φ R D intersects D ′ elsewhere than at (0, 0). In Figures 2 and3, we illustrate the phase plane, the flow Φ y and the curve Φ R D for different cases. Figure 4 shows the graph of four profiles V (y) in [-R, R] and the corresponding intersections in the phase plane.

V ′ 0 D ′ D a Φy 1 1 V V ′ 0 D ′ D Φy 2 a 1 V V ′ 0 D ′ D ΦR a 1 V Figure 2:
The phase plane for α ≥ aλ and Φ R D for three increasing values of R: two cases with no non-trivial profiles and the last with two non-trivial profiles.

Figures 2 and3 show that there are two cases. If D is always outside the domain delimited by the homoclinic orbit to zero, then there cannot exist more than two nontrivial profiles. If D intersects the homoclinic orbit to zero, then for R large enough

V 1 0 V ′ D ′ D a ΦRD V 1 0 V ′ D ′ D a ΦRD V 1 0 V ′ D ′ D a ΦRD Figure 3:
The phase plane for α < aλ and Φ R D for three increasing values of R: no non-trivial profiles, two ones and four ones. there exist as much non-trivial profiles as wanted. The critical value of α is obtained by computing the unstable direction of 0 in the phase plane. This direction is given by the eigenvector (1, √ aλ) related to the eigenvalue √ aλ of the matrix 0 1 aλ 0 , which corresponds to the linearization of the flow at the equilibrium point (0, 0). This shows that D intersects the homoclinic orbit to zero if and only if α < aλ.

4 2 1 3 V 1 0 V ′ D ′ D a Φ R D -4 -3 -2 -1 0 1 2 3
The fact that there exists a profile V M larger than the other ones is classical, see for example [START_REF] Rodríguez-Bernal | Extremal equilibria for reaction-diffusion equations in bounded domains and applications[END_REF]. It is a consequence of Lemmas 2.1 and 2.3: if u(t) is a solution of (2.4) with u(0) larger than any profile, then the ω-limit set ω(u(0)) consists in a profile V M larger than any profile. The fact that V M is even follows from the symmetry of the problem since V M (-.) is also a profile. Notice that this symmetry corresponds to the symmetry of the phase plane and to the fact that the trajectory y → Φ y (V M (-R), V ′ M (-R)) intersect D ′ for the first time at y = R.

A graphic criterium for stability of the profiles

Let V be a profile solution of (1.3). We consider the linearization L V of the flow along the profile V . We recall that L V is given by (3.1). The operator L V is self-adjoint and is a relatively compact disturbance of ∆α and therefore has the same essential spectrum, see [START_REF] Henry | Geometric theory of semilinear parabolic equations[END_REF]. This yields the following characterisation of the stability of a profile V , see [START_REF] Smoller | Shock waves and reaction-diffusion equations, 2nd Edition. Grundlehren der Mathematischen Wissenschaften[END_REF].

Proposition 3.3. The spectrum of L V consists in ] -∞, -α] and a finite number of isolated real eigenvalues µ i > -α of finite multiplicity. As a consequence, the profile V is asymptotically stable if and only if L V has no nonnegative eigenvalue.

We are looking to nonnegative eigenvalues of L V that is to λ ≥ 0 such that there exists ϕ ∈ H 2 (R) such that ϕ ≡ 0 and

ϕ ′′ (y) + f ′ (V (y))½ |y|≤R -α½ |y|>R ϕ(y) = µϕ(y) . (3.2)
By the same argument used for V in Section 2, ϕ is explicit outside (-R, R) and it is equivalent to look to ϕ satisfying

   ϕ ′′ (y) + f ′ (V (y))ϕ(y) = µϕ(y) y ∈ (-R, R) ϕ ′ (-R) = √ µ + αϕ(-R) ϕ ′ (R) = - √ µ + αϕ(R) (3.3)
In this section, we adapt the standard Sturm-Liouville theory to the eigenvalue problem (3.3). Keeping in mind the geometric interpretation of Sturm-Liouville arguments, we will obtain a graphic criterium to count the number of nonnegative eigenvalues of L V . The idea of using such graphic arguments for studying the one-dimensional parabolic equation is not new, see [START_REF] Fusco | A permutation related to the dynamics of a scalar parabolic PDE[END_REF]. In fact, one can understand the whole dynamics of a one-dimensional parabolic equation by similar arguments as shown in [START_REF] Fiedler | Heteroclinic orbits of semilinear parabolic equations[END_REF]. See for example [START_REF] Fiedler | Spatio-temporal dynamics of reaction-diffusion patterns[END_REF] for a review on this subject.

For any µ ≥ 0 we define θ µ by

θ ′ µ = -sin 2 θ µ + (µ -f ′ (V (y))) cos 2 θ µ , y ∈ (-R, R) θ µ (-R) = arctan( √ µ + α) (3.4)
and we set for µ > -α

h(µ) = θ µ (R) + arctan( √ µ + α) . (3.5)
Lemma 3.4. A nonnegative number µ is a nonnegative eigenvalue of L V if and only if h(µ) belongs to πZ.

Proof: Let µ ≥ 0 and ϕ ∈ C 2 ([-R, R]), ϕ ≡ 0, satisfying (3.3). We introduce two functions of class C 2 , ρ > 0 and θ ∈ R, such that (ϕ, ϕ ′ )(y) = ρ(y)(cos θ(y), sin θ(y)).

Notice that these functions exist and are regular since the vector (ϕ, ϕ ′ ) cannot vanish because ϕ is a non-trivial solution of a linear second-order differential equation. Up to the change of the sign of ϕ and to the addition of a constant in πZ to θ, ϕ

′ (-R) = √ µ + αϕ(-R) implies θ(-R) = arctan( √ µ + α). Moreover, it follows from (3.3) that, for y ∈ (-R, R), θ ′ = -sin 2 θ + (µ -f ′ (V (y))) cos 2 θ ρ ′ = ρ sin θ cos θ (1 + µ -f ′ (V (y))) (3.6) Hence, θ = θ µ and ϕ ′ (R) = - √ µ + αϕ(R) is equivalent to θ µ = -arctan( √ µ + α) + kπ with k ∈ Z.
To show the other implication, it is sufficient to follow the previous arguments in the opposite way: we start with a function θ µ , we construct ρ by (3.6) and we check that ϕ = ρ cos θ µ is a solution of (3.3). Lemma 3.5. Let V be a profile solution of (1.3)

. i) If µ > f ′ (V ) ∞ then, for all y ∈ [-R, R], θ µ (y) ∈ (0, π
2 ) and h(µ) ∈ (0, π). ii) If µ 1 > µ 2 then, for all y ∈ [-R, R], θ µ 1 (y) > θ µ 2 (y) and h is an increasing function.

Proof: We prove assertion i) by contradiction. Let

y 0 = inf{y ∈ [-R, R]/θ µ (y) ∈ (0, π/2)}. Notice that y 0 > -R since θ µ (-R) ∈ (0, π/2). If θ µ (y 0 ) = 0, then θ ′ µ (y 0 ) = µ -f ′ (V (y 0 )) > 0 which is absurd since θ µ (y) > 0 for y < y 0 by definition of y 0 . The contradiction is similar if θ µ (y 0 ) = π/2. The proof of ii) is very similar: let y 0 = inf{y ∈ [-R, R]/θ µ 1 (y) ≤ θ µ 2 (y)}. We have θ µ 1 (y 0 ) = θ µ 2 (y 0 ). If θ µ 1 (y 0 ) ∈ π/2 + πZ, then θ ′ µ 1 (y 0 ) -θ ′ µ 2 (y 0 ) = (µ 1 -µ 2 ) cos 2 θ µ 1 (y 0 ) > 0 and it contradicts the definition of y 0 . If θ µ i (y 0 ) ∈ π/2 + πZ, then θ ′ µ 1 (y 0 ) = θ ′ µ 2 (y 0 ), θ ′′ µ 1 (y 0 ) = θ ′′ µ 2 (y 0 ) and θ ′′′ µ 1 (y 0 ) = 2 + µ 1 -f ′ (V (y 0 )) > θ ′′′ µ 2 (y 0 ).
In both cases, we obtain the desired contradiction.

As a consequence of both preceding lemmas, we obtain the following Sturm-Liouvilletype result. Proposition 3.6. Let V be a profile solution of (1.3). Let L V be the operator defined by (3.1) and let h be the function defined by (3.5). i) k ∈ N is such that h(0) ∈ (-kπ, (1k)π] if and only if the operator L V has k nonnegative eigenvalues. ii) if L V has a nonnegative eigenvalue, there exists a positive eigenfunction associated to the largest eigenvalue of L V . iii) the nonnegative eigenvalues of L V are simple.

Proof: By Lemma 3.4, the nonnegative eigenvalues correspond to the values µ for which h(µ) ∈ Z. Since, as shown in Lemma 3.5, µ → h(µ) is increasing and belongs to (0, π) for large µ, then Assertion i) follows from the continuity of h. Moreover, if it exists, the largest eigenvalue µ is necessary such that h(µ) = 0. Thus θ µ (R) ∈ (-π/2, π/2). Since θ µ (-R) ∈ (-π/2, π/2) and since θ ′ µ (y) < 0 as soon as θ µ (y) = ±π/2, we must have θ µ (y) ∈ (-π/2, π/2) for all y ∈ [-R, R]. Following the arguments in the proof of Lemma 3.4, we construct an associated eigenfunction ϕ = ρ cos θ µ which is positive on [-R, R]. It is positive on R by the extension imposed by (3.2). Finally, the eigenspace of a nonnegative eigenvalue is one-dimensional since any eigenfunction must satisfy (3.3). Since L V is self-adjoint, any nonnegative eigenvalue is simple.

The interest of the first assertion of Proposition 3.6 is to have a graphic interpretation. Indeed, let ϕ be the solution of (ϕ, ϕ ′ )(-R) = (1, √ α) and ϕ ′′ + f ′ (V )ϕ = 0. In other words, y → (ϕ, ϕ ′ )(y) is a trajectory of the linearization of the flow D V Φ y along the profile V such that (ϕ, ϕ ′ )(-R) is in the tangent space of D at (V, V ′ )(-R). For each y ∈ (-R, R), (ϕ, ϕ ′ )(y) is therefore in the tangent space of Φ y D at (V, V ′ )(y). By the same arguments as in the proof of Lemma 3.4, θ 0 (y) defined by (3.4) is exactly the argument of the vector (ϕ, ϕ ′ )(y). Therefore, one can observe the angle θ 0 (y) by looking at the argument of the vector tangent to Φ y D at (V, V ′ )(y). Thus, we can graphically interpret Assertion i) of Proposition 3.6 in the following way. Graphic criterium 3.7. Let V be a profile, solution of (1.3) and let L V be the operator defined by (3.1).

Let y ∈ [-R, R] -→ A V (y) ∈ R 2 be a continuous function such that for all y ∈ [-R, R], A V (y) is a unit vector tangent to Φ y D at (V, V ′ )(y).
Then, the number of times A V (y) crosses the direction of D ′ when y describes [-R, R] is exactly the number of nonnegative eigenvalues of L V . The crossings have to be counted in a algebraic way: positively in the clockwise sense and negatively otherwise.

The previous graphic criterium is illustrated in Figure 5. This criterium requires the knowledge of the curve Φ y D for each y ∈ (-R, R). Therefore, it is difficult to apply as it stands. Fortunately, there is a way to count the number of times A V (y) crosses the direction of D ′ with the knowledge of the curve Φ R D only. Indeed, in the particular case of this article, the phase plane associated to

f is such that the curve Φ R D is above D ′ in {(u, v) ∈ R 2 , u ≥ 1} and below D ′ in {(u, v) ∈ R 2 , u < 0}.
For this reason, using homotopy and topological arguments, the vectors A 0 (y) and A V M (y) corresponding to the extremal profiles 0 and V M must have a trivial number of crossings with the direction of D ′ . Moreover, for another profile V , we can compute the number of crossings of A V (y) with D ′ as follows.

1 0 D ′ a Φ R D Φ y4 D Φ y3 D Φ y2 D Φ y1 D V ′ V D = Φ -R D 1 0 D ′ a Φ R D Φ y4 D Φ y3 D Φ y2 D Φ y1 D V ′ V D = Φ -R D Figure 5
: an example of application of the graphic criterium 3.7. One counts along a profile V the number of crossings between the tangent vector to the curve Φ y D and the direction of D ′ . This gives the number of nonnegative eigenvalues of L V . Left, the largest profile has no nonnegative eigenvalue. Right, the other non-trivial profile has one nonnegative eigenvalue. Graphic criterium 3.8. Let V be a given profile. We follow the curve Φ R D from 0 (or from V M ) to V . We count in a algebraic way the number of times the unit tangent vector to Φ R D crosses the direction of D ′ during this course. The resulting number is exactly the number of nonnegative eigenvalues corresponding to the profile V .

The equivalence of the graphic criteria 3.7 and 3.8 can be easily "seen on the figure". However, the rigourous proof uses homotopies and the topology of the plane, similarly to Jordan theorem. This proof is not the subject of this paper and will be omitted.

We can apply the above graphic criterium to the particular phase plane of this article. Because the flow of V ′′ + f (V ) is turning around the equilibrium (a, 0), every profile V , which is not one of the extremal profiles 0 or V M , is unstable, see Figure 6.

1 2 V ′ 0 D ′ D Φ R D a V 1 4 2 1 3 0 V ′ D ′ D a Φ R D 1 V
Figure 6: examples of application of the graphic criterium 3.8. One follows the curve Φ R D from 0 to a profile V and counts in a algebraic way the number of crossings between the tangent vector to the curve and the direction of D ′ . One then obtains the number of nonnegative eigenvalues of L V . At the left, the largest profile [START_REF] Chapuisat | Discussion of a simple model of spreading depression[END_REF] has no nonnegative eigenvalue whereas the other non-trivial profile [START_REF] Aitken | Spreading depression in human hippocampal tissue in vitro[END_REF] has one nonnegative eigenvalue. At the right, the four non-trivial profiles have respectively one, two, one and zero nonnegative eigenvalues.

Energy of the profiles

We recall that the relevant energy for the profiles is the function E defined by (2.2). We also recall that the energy E(V ) of a profile V is equal to the energy E(V |(-R,R) ) of the restriction of V to (-R, R), see Section 2. Let R 0 be the critical thickness introduced in Proposition 3.1. We know that for R > R 0 there exists a unique non-trivial profile V M which is stable and that V M is the largest profile. This section is devoted to the following result. Proposition 4.1. i) If V is an unstable profile, then E(V ) is larger than E(0) = 0 and E(V M ). ii) There exists R 2 ≥ R 0 such that on (R 0 , R 2 ) (this set being possibly empty, for example if α ≥ aλ) the energy E(V M ) of the largest profile is an increasing function of the radius R and on (R 2 , +∞), E(V M ) is a decreasing function of the radius R.

iii) E(V M ) is positive for R > R 0 close to R 0 and converges to -∞ when R goes to +∞. As a consequence, there exists a radius

R 1 > R 0 such that E(V M ) is positive for R 0 < R ≤ R 1 and negative for R > R 1 .
The proof of Proposition 4.1 relies on technics which are more general than the framework of this paper. However, to obtain Proposition 4.1 we have to use properties of the phase plane corresponding to (2.3). These particular properties are observed on the phase plane and on the numerical simulations.

Energy of the unstable profiles

We prove here the first assertion of Proposition 4.1. The proof relies on the following lemma.

Lemma 4.2. Let V be an unstable profile. Let µ > 0 be the first eigenvalue of L V and let ϕ be a positive associated eigenfunction, the existence of which is stated in Proposition 3.6. Then, there exist two globally bounded solutions u -(t) and u + (t) of (2.4) with the following prescribed asymptotic behaviour in H 1 (-R, R):

u ± (t) = V ± e µt ϕ + o e µt when t -→ -∞ . (4.1)
As a consequence, for all (y, t)

∈ (-R, R) × R, u -(y, t) ≤ V (y) ≤ u + (y, t).
Proof: The lemma is a classical application of the theory of stable and unstable manifolds near an equilibrium point. We refer to [START_REF] Chen | Invariant foliations for C 1 semigroups in Banach spaces[END_REF] and [START_REF] Henry | Geometric theory of semilinear parabolic equations[END_REF]. Indeed, we can split the spectrum of the linearisation L V into {µ} and the spectrum contained in the half-plane {z ∈ C, Re(z) ≤ µ -ε} with ε > 0 small enough. As shown in Proposition 3.6, µ is simple and admits a positive eigenfunction ϕ. We defined the strongly unstable set by

W uu (V ) = {u 0 ∈ H 1 (-R, R), ∃ u(t) global solution of (2.4) such that lim t→-∞ e (-µ+ε/2)t u(t) -V H 1 (-R,R) = 0 } .
The theory of invariant manifolds near an equilibrium point shows that W uu (V ) is an invariant one-dimensional manifold which is tangent at V to the line V + Rϕ. The manifold W uu (V ) consists in V and two globally defined trajectories u + and u -satisfying the asymptotic behaviour (4.1). Of course, the last assertion is a direct consequence of this asymptotic behaviour, of the positivity of ϕ and of Lemma 2.1.

The first assertion of Proposition 4.1 is deduced from Lemma 4.2 as follows. Let V be an unstable profile. We know from Section 3 that V lies between 0 and the largest profile V M , both being stable. Let u -(t) be the solution given by Lemma 4.2 for the profile V . By Lemma 2.2, the ω-limit set of u -is non-empty and consists in profiles. We wonder which profiles may belong to this ω-limit set. Assume that there is an unstable profile Ṽ in ω(u -). This profile must satisfy Ṽ ≤ V since u -(t) ≤ V for all t. Moreover, by Lemma 2.3, Ṽ ≡ V and thus, by uniqueness of ODE solutions, Ṽ (y) < V (y) for all y ∈ R. Let ũ+ (t) be the solution given by Lemma 4.2 for the profile Ṽ . For t close to -∞, ũ+ (t) < u -(t). Using Lemma 2.1, Ṽ ≤ ũ+ (t) ≤ u -(t) ≤ V for all times. Since Lemma 2.3 prevents ũ+ (t) to go back to Ṽ , u -(t) cannot converge to Ṽ , which yields a contradiction. Therefore, since ω(u -) is non-empty, the only possibility is ω(u

-) = {0}. Lemma 2.3 then shows that E(V ) > E(0) = 0. 0 < E(u n -(t)) < E(V Rn )
and thus E(u ∞ ) = 0. Let u ∞ (t) be the solution of (2.4) for R = R 0 , with initial data u ∞ . This solution is not a profile, and so its energy E(u(t)) decreases and is negative for t > 0. However, u ∞ (t) must converge to a profile when t goes to +∞ due to Lemma 2.3. Since all the profiles at R = R 0 have an energy equal to 0 by assumption, this is impossible and we get a contradiction. Therefore,

E(V R 0 ) = E(V R 0 M ) = E(V R 0 M
) must be positive.

Discussion

In this paper, we have proved that there exist two critical thicknesses R 1 > R 0 > 0 such that if 0 < R < R 0 , there is no non-trivial profile solution of (1.3). If R 0 < R, there exist non-trivial profiles. One of them is larger than every other one and is stable, whereas every other non-trivial profiles are unstable. Finally, the energy of the unstable profiles is always larger than the energy of the stable profiles. If R 0 < R < R 1 , the energy of the largest profile is larger than the energy of 0, whereas it is smaller if R > R 1 . These results give us informations on the propagation of travelling fronts solution of equation (1.1) as stated in Consequence 1.2. We recall that in equation (1.1), u represents the depolarization of the brain so if u(X) = 0 the brain is normally polarized at the point X, and if u(X) = 1 the brain is totally depolarized. A depolarization wave in the brain corresponds to a travelling front solution of (1.1) where a non-trivial stable profile invades the zero state. Thus, the above mathematical analysis of the stability and energy of the asymptotic profiles of equation (1.3) yields informations on the propagation of depolarization waves in the human brain during stroke. If the width of the grey matter is smaller than R 1 no depolarization wave can propagate through the brain whereas if the width of the grey matter is larger than R 1 , there exist depolarization waves. This may explain why attempts to observe depolarization waves in the human brain have received opposite conclusions in different studies [START_REF] Aitken | Spreading depression in human hippocampal tissue in vitro[END_REF][START_REF] Mayevsky | Cortical spreading depression recorded from the human brain using a multiparametric monitoring system[END_REF][START_REF] Maclachlan | Spreading depression of Leao in rodent and human cortex[END_REF][START_REF] Somjen | Ions in the brain: Normal Function, Seizures, and Stroke[END_REF][START_REF] Sramka | Functional ablation by spreading depression: possible use in human stereptactic surgery[END_REF]. In the human brain, the grey matter is particularly thin and its thickness may vary a lot. The difficulties to observe depolarization waves in the human brain could thus be explained simply by the morphology of the brain. If the initial depolarization takes place in a part of the brain where the grey matter is large enough, then the depolarization can spread through a part of the grey matter. But if the initial depolarization takes place in a part of the brain where the grey matter is very thin then no propagation will occur. Hence it is possible to observe depolarization waves in the human brain [START_REF] Aitken | Spreading depression in human hippocampal tissue in vitro[END_REF][START_REF] Mayevsky | Cortical spreading depression recorded from the human brain using a multiparametric monitoring system[END_REF][START_REF] Maclachlan | Spreading depression of Leao in rodent and human cortex[END_REF], but they will not appear in all the cases [START_REF] Somjen | Ions in the brain: Normal Function, Seizures, and Stroke[END_REF][START_REF] Sramka | Functional ablation by spreading depression: possible use in human stereptactic surgery[END_REF]. Moreover the depolarizations would not travel over large distances because they will stop as soon as the grey matter becomes too thin. For example, they will stop at the bottom of the large sulkus, where the grey matter becomes thiner. This has already be observed in the case of the migraine with aura. The aura may be due to a depolarization wave [START_REF] Gorji | Spreading depression: a review of the clinical relevance[END_REF][START_REF] James | Cortical spreading depression and migraine: new insights from imaging?[END_REF] and for most of the patients the aura stops at the bottom of the Rolando sulkus.

This paper also gives informations on how the initial depolarization is erased as illustrated in Figure 1. If the grey matter is very thin, then the depolarization is quickly absorbed uniformly in the excited area. If the grey matter is a little bit larger, the depolarization occurs in a different way. The depolarized area shrinks progressively while the cells stay totally depolarized as long as they are in the middle of the excited area. In this case, the repolarization is due to a travelling front where the normally polarized state invades the depolarized state. The neurons in the center of the excited area may stay depolarized a very long time, which may cause local damages even if no spreading wave is observed. To our knowledge this behaviour has never been observed in experiments.

In order to verify biologically these results and to validate the model, it would be interesting to estimate numerically the values of the thresholds R 0 and R 1 . They depend on the values of the parameters of equation (1.1) that are hard to compute due to the few possibilities of quantitative mesures in the brain. Moreover these parameter values may vary a lot from one person to another and from a species to another. Estimations of this parameter values for the rodent and for the human may give us more informations on the difficulties to observe depolarization waves in the human brain. Indeed even if the mecanisms are the same in the rodent brain and in the human brain, the densities of cells are totally different and this must influence the values of the parameters of this model and thus of the thresholds. For the moment we are just trying to understand numerically the influence of each parameters on these thresholds.

Figure 4 :

 4 Figure 4: The intersections of Φ R D and D ′ in the phase plane and the corresponding profiles in the (y, V ) plane

Acknowledgements: the authors would like to thank Thierry Gallay for several fruitful discussions.

By using u + (t), we prove similarly that E(V ) > E(V M ). Notice that the above arguments also prove that every unstable profile intersects the other unstable ones as illustrated in Figure 4.

Energy of the stable profile

The proof of the second and third assertions of Proposition 4.1 is splitted in several parts. We use in this section the notations of Section 3.1 with obvious changes to add the dependence of the different objects with respect to R.

M , we invoke the graphic interpretation of the profiles developped in Section 3.1. By symmetry of the largest profile, R)) belongs to D ′ and at this point of the phase plane, the flow Φ y crosses D ′ transversally from the top to the bottom. Thus, for small ε > R)) is below D ′ and by continuity, there exists an intersection of Φ R+ε D and

M piles up to the singularity (1, 0) of the phase plane.

By a standard implicit functions argument, one shows

We admit that it is always the case when R > R 0 .

• E(V R M ) converges to -∞ when R goes to +∞. This fact has already been proved in [START_REF] Chapuisat | Existence and non-existence of curved front solution of a biological equation[END_REF] by the first author. We only give here brief arguments. When R increases, the trajectory y ∈ R)) goes to the right and V R M (0) converges to 1. For large R, the trajectory y ∈

M (y) is differentiable. Let E ρ be the energy functional defined by (2.2) with R = ρ. Let V be a given fixed profile. We consider the variation of ρ → E ρ (V ) at ρ = R.

Since V is an equilibrium of (2.1), it is a critical point for the energy E R . In other words, for any

It remains to notice that R → V R M (0) is increasing and converges to 1 when R goes to +∞. Moreover, there exists β > 0 such that F is negative on (0, β) and positive on (β, 1]. In fact, (β, 0) is the point of the phase plane where the homoclinic orbit to 0 intersects the horizontal axis, see Figure 7.

When R passes the value R 0 , the creation of both non trivial profiles occurs through a saddle-node bifurcation: two equilibrium states appear at the same point, one stable and one unstable. The stable profile is V R M , the largest one, and we denote by V R the unstable profile which lies between 0 and V R M . From now on, we work in (-R, R) by using the analogy presented in Section 2. We want to show that E(V R M ) is positive for R close to R 0 . When R decreases to R 0 : the profiles V R M and V R collide. As E(V R ) > 0 for R > R 0 as shown in Section 4.1, we must have at the limit

R 0 be a sequence of thickness decreasing to R 0 and let V Rn be the associated unstable profiles. Let u n -(t) be the sequence of solutions of (2.4) given by Lemma 4.2 applied to V Rn . We know by Section 4.1 that u n -(t) converges to 0 when t goes to +∞. Let K = lim n→+∞ V Rn H 1 (-R,R) . Notice that K is positive since for each n there exists y ∈ (-R, R) such that V Rn (y) > a. We set t n to be a time such that u n -(t n ) H 1 (-R,R) = K/2. By a compacity argument similar to Lemma 2.2 (R is moving, but nothing singular happens), one can extract a subsequence such that u ϕ(n) -(t ϕ(n) ) converges in H 1 (-R, R) to a function u ∞ . Notice that by construction, u ∞ is neither 0 nor V R 0 . The gradient structure of (2.4) shows that, for all n and t,


A Garenne

J Henry

C O Tarniceriu
email: oana.tarniceriu@inria.fr

Analysis of Synchronization in a Neural Population by a Population Density Approach

Keywords: single neuron model, population density approach, synchronization AMS subject classification: 92C20, 92D25, 34D10, 35L65 1. Introduction

In this paper we deal with a model describing the evolution in time of the density of a neural population in a state space, where the state is given by Izhikevich's two -dimensional single neuron model. The main goal is to mathematically describe the occurrence of a significant phenomenon observed in neurons populations, the synchronization. To this end, we are making the transition to phase density population, and use Malkin theorem to calculate the phase deviations of a weakly coupled population model.

It is usually a very complex and tedious task to figure out the population coding and information processing modalities of neuronal populations in experimental contexts involving living models. Thanks to the existing electrophysiological tools and setups, we now have reliable and real-time access to various parts of the nervous system signalling activity like isolated spikes or local field potentials recordings. These signals reveal significant but minor fractions of the whole brain (or brain regions) electrical activity and most of the information stream dynamics which are yet essential to a global understanding of its intrinsic properties remain out of sight. Among these network-scale properties, the synchronization mechanisms in large neuronal assemblies is one of the most widely investigated by the neuroscientists since it is known to play key roles in physiological as well as in pathological processes [START_REF] Dejean | Dynamic changes in the cortex-basal ganglia network after dopamine depletion in the rat[END_REF][START_REF] Sheeba | Neuronal synchrony during anesthesia: a thalamocortical model[END_REF]. A widespread method to deal with this problem, and more generally to investigate complex biological structures and having a full access to all their components, consists in building their most reliable possible computational model. Its qualitative and quantitative properties are thus made available through the adequate numerical implementation. The validated computer model can then be used by the experimenter as a comprehension and prediction tool and exploited to corroborate, invalidate or foresee working hypotheses. In practical terms, the implementation of these mathematical neuronal models mainly confronts two major issues: (i) the computation power required by neuronal modelling usually involves a drastic trade-off between the network size and model refinement degree [START_REF] Meunier | Playing the devil's advocate: is the Hodgkin-Huxley model useful?[END_REF][START_REF] Izhikevich | Which model to use for cortical spiking neurons?[END_REF]) and (ii) the mathematical formalism potentialities of the equations are lost and scaled-down to their numerical schemes [START_REF] Chauvet | The use of representation and formalism in a theoretical approach to integrative neuroscience[END_REF][START_REF] Bennani | A hierarchical modeling approach of hippocampus local circuit[END_REF]. To overcome these two difficulties and to provide innovative tools and formalism to neuroscientists, a population density approach is presented here.

The population density approach has been successfully used until now to describe the evolution of structured populations in the context of demography, epidemiology, economy etc. In these cases, the theoretical framework has been very well developed and certain properties of the solutionssuch as stability for example -have been intensively studied. Since we devoted our attention to the study of a population of neurons having an internal structure given by their physiological properties, we find appropriate to use the same method to describe the evolution of neural populations. The main advantage of it lies in the fact that it allows us to obtain an analytical expression that describes the synchronization of the neurons regardless the number of neurons in the population, and consequently, offering us a way of understanding the mechanisms that leads to the occurrence of the phenomenon, and, perhaps, a way of controlling it. The long time behavior of the solution expressing the synchronized state and the conditions under which the stability of such a solution can be achieved represents our next point of interest. In particular, we are interested to highlight the effect of different type of synaptic couplings on the stability of the synchronized solution.

Synchronization of neurons has made the object of investigations, both analytically and numerically, of a number of papers, as [START_REF] Faugeras | Abolute stability and complete synchronization in a class of neural fields models[END_REF], [START_REF] Izhikevich | Phase equations for relaxation oscillators[END_REF], [START_REF] Medvedev | Synchronization and transient dynamics in the chains of electrically coupled Fitzhugh-Nagumo oscillators[END_REF], [START_REF] Koppel | Mechanisms of phase-locking and frequency control in pairs of coupled neural oscillators[END_REF], to remind just few of them.

We start by considering a neural population composed by a large number of neurons, each of them being characterized by a state vector denoted by w = (v, u). The evolution of the state of a single neuron is described by:

       v(t) = F v (v, u), u(t) = F u (v, u), (1.1)
where v denotes the membrane potential and u is the so-called recovery variable. It is assumed that every time membrane's potential v reaches a critical value denoted by v f , a spike is emitted, v is reset to a reset potential value v r and u is increased by a fixed value d (v r , v f and d are known). Therefore, the reset mechanism is formally described by:

v > v f ⇒ v r ← v and u + d ← u. (1.2)
For the sake of simplicity we are using general notations for the functions that give the evolution of the two state variables, and we will assume regularity properties such that certain behaviors to be obtained. In particular, choosing

       F v (v, u) = 0.04v 2 + 5v + 140 -u + I(t), F u (v, u) = a(bv -u),
we arrive to the Izhikevich model, where {a, b, v r , d} are parameters and I is an applied current.

For more details about this model see [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The geometry of excitability and bursting[END_REF].

The Izhikevich model has been proven to admit different behaviors for different choice of parameters, in particular the existence of a periodic limit cycle attractor, which is of a particular interest for us.

The method we are using in this paper (Population Density Approach) is based on the following assumptions: the number of neurons in the population is considered large and all the neurons have the same dynamics given by (1.1) and the same coupling pattern. Then, the evolution of the population is described by the dynamics of the population density function, p(t, v, u), which satisfies:

p(t, v, u) dv du = N,
where N is the number of neurons in the population. We are dealing therefore with a population structured by the state variables (v, u). The evolution of p(t, v, u) is given by the conservation law:

∂ ∂t p(t, v, u) = -divJ(t, v, u), (1.3)
where J is the neural flux flowing through the state (v, u). It is considered that the flux is divided in two components, a streaming flux, J s (t, v, u), and an interaction flux, J i (t, v, u). The streaming flux is thought as the flux of p(t, v, u) neurons flowing through the state (v, u) at time t, with the velocity F (w) = ẇ(t), where we have denoted F = (F v , F u ) and w = (v, u). On the other hand, the interaction flux stands for the interactions between a neuron and the rest of the population, and depends on the synaptic afferents from other neurons in the population. We consider it in the form

J i (t, v, u) = σ(t)ē v v v-
p(t, ṽ, u) dṽ, (1.4) where σ represents the average individual spike reception rate of a single neuron in the population, ēv is the unit vector in the direction v, and is the amplitude of a jump in potential of a neuron when receiving a spike. The variable σ can have various forms, depending on the parameters related to the neural connectivity which can be taken into account. If synaptic delays are considered, then one way to express σ is

σ(t) = G N τ + τ - r(t -τ )α(τ ) dτ, (1.5)
where G is the average number of afferents, N is the number of neurons in the population, τ τ -/+ are the extremum values for spike conduction delays, and α is a spike conduction delay kernel. In the following we assume that there is no conduction delay. The firing rate r(t) is given by the flux in the direction ēv at the boundary v = v f , i.e.

r(t) = umax u min J(t, v f , u) du ēv . (1.6)
The reset mechanism (1.2) is naturally translated here to the following boundary condition

J v (t, v + r , u + d) = J v (t, v - r , u + d) + J v (t, v f , u), (1.7)
where J v denotes the component of the flux in the direction ēv . The model described above has been introduced by J. Modolo in his Ph.D thesis; a full description of the model as well as extensions of it and numerical results can be found in [START_REF] Modolo | Modélisation et analyse mathématique des effets de la stimulation cérébrale profonde dans la maladie de Parkinson[END_REF], [START_REF] Modolo | Development and validation of a neural population model based on the dynamics of discontinuous membrane potential neuron model[END_REF], [START_REF] Modolo | Dynamics of the subthalamo-pallidal complex in Parkinson's Disease during deep brain stimulation[END_REF], [START_REF] Modolo | New insights offered by a computational model of deep brain stimulation[END_REF].

The weakly-coupled neurons case

In what follows we will assume that the amplitude of the jump in the potential of a neuron which receives a spike satisfy 1,

and we approximate v v-p(t, ṽ, u) dṽ ≈ p(t, v, u) by ignoring the second order terms in . Then, the conservation law simplifies to

∂ ∂t p(t, v, u) = -div [(F (v, u) + σ(t)ē v ) p(t, v, u)] . (2.1)
As we argued at the beginning, our main goal is to derive the phase density model in the case of weakly coupled neurons via specific tools of weakly perturbed systems theory. Therefore we start by considering the uncoupled case of the population density approach model, i.e. the conservation law is given by:

∂ ∂t p(t, v, u) = -div [F (v, u)p(t, v, u)] , (2.2)
which can be written equivalently along the trajectories as

d dt p(t, v(t), u(t)) = -div [F (v(t), u(t))] p(t, v(t), u(t)), (2.3)
where the evolution of (v, u) is given by (1.1) together with initial conditions and the reset mechanism described by (1.2).

We assume in the following the existence of a stable limit cycle attractor C for (1.1) and that it has an attractive normally hyperbolic compact invariant manifold. In the Izhikevich model's case, it has been shown (see [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The geometry of excitability and bursting[END_REF]) that there exists I such that the model admits a stable limit cycle corresponding to periodic spiking. Then we suppose that the neural population is distributed along this periodic trajectory, and consider next the weakly coupled case. Note that the evolution of weakly-coupled population densities is also given by (2.3) but now along the weakly perturbed trajectories described by

       v(t) = F v (v(t), u(t)) + σ(t), u(t) = F u (v(t), u(t)).
(2.4) Then (see for example [START_REF] Hoppensteadt | Weakly connected neural networks[END_REF], Theorem 4.1), the perturbed system has an attractive normally hyperbolic invariant manifold ω which is an -perturbation of the original one. For small enough we take a stable neighborhood W of ω such that it includes C in it, and then fix the domain

D = [v min , v f ] × [u min , u max ] such that W ⊂ D.
Therefore, considering the dynamics of (2.3) in the domain D, we can see that the active part of the boundary ∂D is ({v f } × (u min , u max )) W , and, we can assume without loss of generality that p = 0 for the rest of the boundary.

Then, the evolution of the weakly coupled population is given by:

                               ∂ ∂t p(t, v, u) = -div [(F (v, u) + σ(t))p(t, v, u)] , t > 0, (v, u) ∈ D, J v (t, v + r , u + d) = J v (t, v - r , u + d) + J v (t, v f , u), t ≥ 0 p(0, v, u) = p 0 (v, u), (v, u) ∈ W with p 0 (u, v) = 0 (v, u) ∈ D \ W p(t, v, u) = 0, t ≥ 0, (v, u) ∈ ∂D \ Γ l , (2.5)
with

Γ l = ({v f } × (u min , u max )) W .
Under these conditions the support of p(t, u, v) remains in W . Correspondingly, the evolution of an uncoupled population is given by

                               ∂ ∂t p(t, v, u) = -div [F (v, u)p(t, v, u)] , t > 0, (v, u) ∈ D, J v (t, v + r , u + d) = J v (t, v - r , u + d) + J v (t, v f , u), t ≥ 0 p(0, v, u) = p 0 (v, u), (v, u) ∈ W with p 0 (u, v) = 0 (v, u) ∈ D \ W p(t, v, u) = 0, t ≥ 0, (v, u) ∈ ∂D \ Γ l , (2.6)
Let us look more closely at the problems we defined above. Since we consider in the uncoupled case (2.6) the neural population distributed along the limit cycle, solving the problem means integrating the differential equation (2.3) on the limit cycle given by (1.1). Note that a discontinuity occurs at the time of the spike t s , but, the problem is well defined also after the time t s since we can solve the same equation starting with an initial condition at time t + s , initial condition which is given by (2.6) 2 . Since the trajectory we integrate on is periodic, then (v(t), u(t)) remain bounded for all t ≥ 0, which makes the problem to be well-posed for all t ≥ 0. In the weakly coupled case, the equation (2.5) 1 can be written also in the form (2.3), with the only difference that the trajectories on which we integrate are given now by (2.4). So, we can argue in the same way the extension of the solution after the times of spikes, but, in this case, the integration along the trajectories is restricted for the case when the trajectories defined by (2.4) do not cross.

Phase equations. Malkin theorem.

The first step we take in this section is to derive a phase density model corresponding to the uncoupled population in the specific case we mentioned, i.e. when the population is distributed along the limit cycle. We consider next the weakly perturbed trajectories with respect to the limit cycle, which will correspond to the weakly connected population, and see how the model changes in this case. We mention that in the case of discrete neural networks, an adapted model of Malkin theorem has been given in [START_REF] Hoppensteadt | Weakly connected neural networks[END_REF].

Let us consider thus first the case of a uncoupled neural population, therefore the conservation law is given by:

∂ ∂t p(t, v, u) + div (F (v, u)p(t, v, u)) = 0, t > 0, (v, u) ∈ D, (3.1)
or, equivalently (2.3), and the evolution of the state of a single neuron is given by the system (1.1) introduced in the first section. In order to avoid confusion, we will denote by (ṽ, ũ) the stable limit cycle attractor. Obviously, (ṽ, ũ) is a periodic solution to the system (1.1), and we denote by T its period. We want to describe next the activity on the limit cycle in terms of its phase of oscillation.

To this end, let us denote by θ(t) = θ 0 + t, t ≥ 0, the natural phase of the limit cycle, where θ 0 ∈ [0, T ) denotes the initial phase; then

θ = 1.
Without loss of generality, we can choose θ such that the next spike arrival corresponds to the phase θ = T . The parametrization of the curve is given by the mapping (ṽ, ũ) :

[0, T ) → R 2 , (ṽ, ũ)(θ(t)) = (v(θ(t)), u(θ(t))),
and (ṽ, ũ)(θ(t)) is a solution to:

dṽ(θ(t)) dt = F v (ṽ, ũ) θ, t ∈ [0, t s ), dũ(θ(t)) dt = F u (ṽ, ũ) θ, t ∈ [0, t s ).
The reset mechanism will take place at θ = T (= θ(t s )) and is described by

(ṽ(T ), ũ(T )) → (v r , ũ(T ) + d) = (v r , ũ0
), due to the periodicity of (ṽ, ũ).

The next step is to make the transition from the population densities in the state space, to phase densities populations. Let us denote by p(t, s) the distribution at moment t along the periodic trajectory, i.e., denoting by C = trace(γ), p = p(t, s)δ C , where δ is the Dirac function and s is the arclength. Since p satisfies (3.1), it follows that, for

ζ ∈ {f ∈ C ∞ ([0, ∞) × [0, T ] × [0, T ]); f -periodic
function in the last two arguments and with compact support w.r.t. first argument},

T 0 D ∂ ∂t p(t, v, u) + ∂ ∂u (F u (v, u)p(t, v, u)) + ∂ ∂v (F v (v, u)p(t, v, u)) ζ(t, v, u) du dv dt = - T 0 D p(t, v, u) ∂ ∂t ζ(t, v, u) + F (v, u) • ∇ζ du dv dt = - T 0 C p(t, s) ∂ζ ∂t + F • ∇ζ ds dt.
We choose the test function ζ such that:

ζ(v, u) = ψ(s), ∇ζ = dψ ds (τ v , τ u ),
where τ = (τ v , τ u ) is the tangent unit vector. Then

T 0 C p(t, s) ∂ζ ∂t + F • ∇ζ ds dt = = T 0 C ψ(s) ∂ p(t, s) ∂t + ∂ ∂s [(F • τ )p(t, s)] ds dt. Taking τ = dṽ ds , dũ ds = dṽ/ dt ds/ dt , dũ/ dt ds/ dt = F v ṡ , F u ṡ , then F • τ = F 2 ṡ = ṡ.
Therefore, on the periodic trajectory (ṽ, ũ), the evolution of population density is given by the following equation:

∂ ∂t p(t, s) + ∂ ∂s (p(t, s) ṡ) = 0, t ∈ [0, t s ), (3.2)
and since the population is distributed on the limit cycle, we can extend it by periodicity to t > 0 (we use the same argument in the following computations) . We define the phase density q(t, θ) such that the following relation holds:

q(t, θ) dθ = p(t, s) ds, (3.3)
which implies that q(t, θ) = p(t, s) ds dθ = p(t, s) ds dt dt dθ .

Let us compute now the total derivative of q defined as above with respect to t:

d dt q(t, θ(t)) = d dt (p(t, s) ṡ) = ṡ ∂ ∂t p(t, s) + sp(t, s) + ṡ2 ∂ ∂s p(t, s) = ṡ ∂ ∂t p(t, s) + s ṡ p(t, s) + ṡ ∂ ∂s p(t, s) = ṡ ∂ ∂t p(t, s) + ∂ ∂s (p(t, s) ṡ) = 0.
We obtained thus d dt q(t, θ(t)) = 0, t > 0, which is equivalent to

∂ ∂t q(t, θ) + θ ∂ ∂θ q(t, θ) = 0, t > 0, θ ∈ (0, T ), (3.4)
and, since θ = 1, ∂ ∂t q(t, θ) + ∂ ∂θ q(t, θ) = 0, t > 0, θ ∈ (0, T ).

(3.5)

The boundary condition reads: q(t, 0) = q(t, T ), t ≥ 0.

Let us consider next the weakly connected system, i.e. (2.5), where

σ(t) = c 0 umax u min J v (t, v f , u) du,
( c 0 is a constant), and in our case:

J v (t, v, u) = (F v (v, u) + σ(t))p(t, v, u).
We remind that, in this case, the evolution of a single neuron is given by:

       v(t) = F v (v, u) + σ(t), t ∈ [0, t s ), u(t) = F u (v, u), [0, t s ), (3.6)
and, at the spike moment t s the reset condition takes place. The solution to the above system will be denoted by (v , u ).

We consider next the corresponding phase density to the perturbed system, q (t, τ, θ ), i.e. the population density having at time t the phase θ , where by θ we denote the perturbed phase. The dependence of q on the variable τ will be explained later. The goal is to derive the phase density equation for q , as in the uncoupled case. For that, we shall look first at the changes produced in the evolution of the perturbed phase. For every θ 0 ∈ [0, T ) an initial phase, the corresponding perturbed phase at time t will be given by θ = θ 0 + t + ϕ(τ, θ 0 ), with ϕ slow phase deviations due to the weak connections, and τ the slow time.

Then, for any t ∈ [0, t s ) such that θ < T , the corresponding perturbed trajectory can be written in terms of phase deviation and small perturbation with respect to the unperturbed trajectory as:

       v (t) = ṽ(θ 0 + t + ϕ(τ, θ 0 )) + y 1 (θ 0 + t + ϕ(τ, θ 0 )), u (t) = ũ(θ 0 + t + ϕ(τ, θ 0 )) + y 2 (θ 0 + t + ϕ(τ, θ 0 )). (3.7)
In the above expression, (v , u ) are evaluated at time t which corresponds to an unperturbed state (ṽ, ũ) on the orbit C with phase θ 0 + t + ϕ(τ, θ 0 ). We choose the functions (y 1 , y 2 ) periodic with the same period as (ṽ, ũ), and since ṽ, ũ satisfy the reset relations, (u , v ) will satisfy the reset relations (up to an O( ) correction.)

The following theorem is the main result of this section.

Theorem 1. (Malkin)

Consider the uncoupled system (2.6) and suppose that the corresponding single neuron dynamics given by (1.1) have an exponentially orbitally stable periodic solution (ṽ, ũ). Let us consider next the weakly coupled system (2.5) with the corresponding single neuron dynamics given by (3.6), and denote τ = t the slow time, and ϕ the slow phase deviations from the natural oscillation (ṽ(t), ũ(t)). Then ϕ is a solution to

dϕ(τ, θ 0 ) dτ = H(ϕ), ∀θ 0 ∈ (0, T ), ∀τ > 0, (3.8)
where

H(ϕ) = c 0 ṡ(T ) T T 0 q (t -θ 0 -ϕ(τ, θ 0 ), τ, T )λ 1 (t) dt, (3.9)
and λ is the unique nontrivial periodic solution to the adjoint system

       d dt λ 1 (t, ϕ) = - ∂F v ∂v (ṽ, ũ)λ 1 (t, ϕ) - ∂F u ∂v (ṽ, ũ)λ 2 (t, ϕ), t ∈ [0, T ), d dt λ 2 (t, ϕ) = - ∂F v ∂u (ṽ, ũ)λ 1 (t, ϕ) - ∂F u ∂u (ṽ, ũ)λ 2 (t, ϕ), t ∈ [0, T ), (3.10)
satisfying the normalization condition λ(0) • F (ṽ(0), ũ(0)) = 1.

Proof. In the following, wherever there are missing arguments of a variable, it should be understood that it is calculated in (θ 0 + t + ϕ(τ, θ 0 )); this in done in order to avoid overnotation.

We linearize first the functions in (3.6) with respect to (ṽ, ũ), and obtain:

       d dt v (t) = F v (ṽ, ũ) + ∇F v (ṽ, ũ) • y + σ(t), t ∈ [0, T ), d dt u (t) = F u (ṽ, ũ) + ∇F u (ṽ, ũ) • y, t ∈ [0, T ), (3.11)
where we denoted y = (y 1 , y 2 ). We remind that the expression of σ is given by:

σ(t) = c 0 umax u min J v (t, v f , u) du,
where J v was the flux in the direction v calculated at the firing threshold v f ; we want thus to express σ in terms of phase densities and to approximate it with respect to the periodic trajectory (ṽ, ũ), as we have done with the rest of the functions. First notice that, for each trajectory, the firing threshold v f corresponds to the phase T . Then we can write the flux at the firing threshold as the speed on each perturbed trajectory calculated in the phase T multiplied by the phase densities population on the corresponding trajectory. Linearizing it with respect to the periodic trajectory, and keeping into account that, since the perturbation is small (of the order ), then the perturbation of the periods of the perturbed trajectories are of the same order, therefore we can approximate the speed on each trajectory at phase T with the speed on the periodic trajectory at the same phase, by neglecting the terms of second order in . Then the total flux is given by:

σ(t) = c 0 ṡ(T )q (t, τ, T ), (3.12)
where ṡ(T ) = ((F v (ṽ(T ), ũ(T ))) 2 + (F u (ṽ(T ), ũ(T )) 2 )) 1/2 .

Differentiating now (3.7) with respect to t, we obtain:

       dv dt (t) = dṽ dt (1 + d dτ ϕ(τ, θ 0 )) + dy 1 dt , t ∈ [0, T ), du dt (t) = dũ dt (1 + d dτ ϕ(τ, θ 0 )) + dy 2 dt t ∈ [0, T ). (3.13)
Here and in the following ϕ is taken as a parameter.

Taking into consideration (3.11) and (3.13) it follows (after dividing by ):

       d dt y 1 = ∇F v (ṽ, ũ) • y + b 1 , d dt y 2 = ∇F u (ṽ, ũ) • y + b 2 .
(3.14)

We are looking for (y 1 , y 2 ) periodic of period T satisfying these equations. In the above expressions we have denoted:

b 1 = c 0 ṡ(T )q (t, τ, T ) -F v (ṽ, ũ) d dτ ϕ(τ, θ 0 ), and b 2 = -F u (ṽ, ũ) d dτ ϕ(τ, θ 0 ).
The above system is a linear nonhomogeneous one, and to study the existence and uniqueness of the solutions one must consider the adjoint linear homogeneous system for the adjoint variable λ(θ 0 + t + ϕ) := λ(t, ϕ), with λ = (λ 1 , λ 2 ):

       d dt λ 1 (t, ϕ) = - ∂F v ∂v λ 1 (t, ϕ) - ∂F u ∂v λ 2 (t, ϕ), t ∈ [0, T ), d dt λ 2 (t, ϕ) = - ∂F v ∂u λ 1 (t, ϕ) - ∂F u ∂u λ 2 (t, ϕ), t ∈ [0, T ), (3.15)
where the functions with missing arguments are calculated in (ṽ(θ 0 + t + ϕ(τ, θ 0 )), ũ(θ 0 + t + ϕ(τ, θ 0 ))).

We take the following normalization condition:

1 T T 0 (F v (ṽ, ũ)λ 1 (t, ϕ) + F u (ṽ, ũ)λ 2 (t, ϕ)) dt = 1.
We impose the condition λ(T ) = λ(0), and using Fredholm's alternative, the system (3.14) has a unique periodic solution if and only if the orthogonality condition

1 T T 0 λ(t, ϕ(τ, θ 0 )) • b(θ 0 + t + ϕ(τ, θ 0 )) dt = 0 (3.16)
is satisfied. Replacing b in the above expression and taking into account the normalization condition, we get:

dϕ(τ, θ 0 ) dτ = c 0 ṡ(T ) T T 0 λ 1 (t, ϕ(τ, θ 0 ))q (t, τ, T ) dt.
The above result was obtain for the time interval [0, T ), and allows us to compute the evolution in slow time of the phase deviations during an inter-spike interval. The results can be extended for all intervals [nT, (n+1)T ), n ∈ N by taking the initial values λ(nT

+ ) = λ(0), y(nT + ) = y(0), n ∈ N * .
Notice now that the solutions to the adjoint system satisfy λ(θ 0 + t + ϕ, 0) = λ(t, ϕ), therefore it suffices to find a solution to the adjoint system corresponding to θ 0 + ϕ = 0.

Thus we can write dϕ(τ, θ 0 ) dτ = c 0 ṡ(T )

T T 0 λ 1 (θ 0 + t + ϕ(τ, θ 0 ), 0)q (t, τ, T ) dt = c 0 ṡ(T ) T T 0 λ 1 (s, 0)q (s -θ 0 -ϕ(τ, θ 0 ), τ, T ) ds := H(ϕ(τ, θ 0 )).
Moreover, it is obvious that F (ṽ, ũ) is a solution to the homogeneous linear problem

ż(t) = ∇F (ṽ(t), ũ(t))z(t),
and λ(t, 0) is a solution to the adjoint linear system λ = -(∇F (ṽ(t), ũ(t))) λ.

Any such solution satisfies λ(t, 0) • z(t) = const, ∀t ≥ 0, and from the normalization condition we get that

λ(t, 0) • z(t) = λ(t, 0) • F (ṽ(t), ũ(t)) = 1, ∀t ≥ 0,
and in particular for t = 0, which allows us to rewrite the normalization condition as

λ(0) • F (ṽ(0), ũ(0)) = 1.
We remind that θ = θ 0 + t + ϕ(τ, θ 0 ), and the phase on the limit cycle was given by

θ = θ 0 + t.
Since the density along the perturbed trajectory should satisfy q(t, θ) dθ = q (t, τ, θ ) dθ , it follows that q(t, θ) = q (t, τ, θ ) dθ dθ .

Differentiating with respect to t the above equality, and taking into account that d dt q(t, θ) = 0, we obtain that the following relation holds:

d dt q (t, τ, θ ) dθ dθ = 0. Writing dθ dθ = dθ dt dt dθ = dθ dt ,
and since dθ dt := g(θ ), it follows that

d dt q (t, τ, θ ) dθ dθ = g(θ ) d dt q (t, τ, θ ) + q (t, τ, θ ) d dt g(θ ) = g(θ ) d dt q (t, τ, θ ) + g (θ )q (t, τ, θ ) = 0.
Dividing the last expression by g, we get

d dt q (t, τ, θ ) + g (θ )q (t, τ, θ ) = 0, or, equivalently ∂ ∂t q (t, τ, θ ) + ∂ ∂θ ( θ q (t, τ, θ )) = 0.
Having in mind the expression of θ , we get:

dθ dt = 1 + dϕ dτ = 1 + H(ϕ),
where H(ϕ) is given by (3.9). Therefore we arrived to the following system of equations for the phase densities of the perturbed system:

               ∂ ∂t q (t, τ, θ ) + ∂ ∂θ [(1 + H(ϕ(τ, θ 0 )))q (t, τ, θ )] = 0, t ≥ 0, θ ∈ (0, T ), q (t, τ, 0) = q (t, τ, T ), t ≥ 0, q (0, τ, θ ) = q 0 (θ ), θ ∈ [0, T ). (3.17)
Actually, since q evolves in two different time-scales, the slow and fast times, we are looking for an expression of q in the following form:

q (t, τ, θ ) = q (τ, θ ) + q (t, τ, θ ), (3.18)
where q is the mean part of q and q is its oscillating part. For the sake of generality we choose in the following t ∈ [nT, (n + 1)T ), i.e. an inter-spike interval. We make the choice of oscillating and mean part such that, the evolution of the oscillation part, q , to describe the periodic behavior of the phase-densities in every inter-spike interval, therefore we will impose the periodicity in t (the fast time), but also the periodicity in θ which is due to (3.17). On the other hand, q will describe the phase density evolution in slow time, and again, we have the periodicity of q with respect to θ given by the original system. We must choose then q such that to be periodic with respect to t and to have the mean equal to zero, i.e.

(n+1)T nT q(t, τ, θ ) dt = 0. (3.19)
If we replace the expression of q (3.18) in (3.17) and regroup the terms with , we obtain two equations for q and q respectively. The equation for q is

∂ ∂t q (t, τ, θ ) + ∂ ∂θ q (t, τ, θ ) = - ∂ ∂θ q (τ, θ ), (3.20)
Let us assume more that q(nT, τ, θ ) = 0 for n ∈ N. Then, integrating (3.20) along the characteristics θ = t + const., we obtain for 0 < t -nT < θ that q (t, τ, θ ) = q (τ, θ -t) -q (τ, θ ).

Actually, using the periodicity in θ of q , the same expression holds true for θ < t -nT < T . Taking in the above expression t = (n + 1)T and using again the periodicity in θ of q , we obtain that q ((n + 1)T, τ, θ ) = 0, therefore the function defined in this way satisfies the periodicity condition in t and one can verify directly that (3.19) is satisfied.

Therefore it remains to analyze the evolution of q . Note that, by finding the evolution of q , the evolution of q , and thus that of q , is fully determined. Taking the terms in obtained by replacing (3.18) in (3.17), we get:

∂ ∂τ q (τ, θ ) + ∂ ∂θ (H(ϕ)(q (τ, θ ) -q (t, τ, θ ))) = 0, (3.21)
and, taking into account the expression of q

∂ ∂τ q (τ, θ ) + ∂ ∂θ (H(ϕ)q (τ, θ -t)) = 0.
Since here t is considered fixed, making the change of variable θ → x := θ -t we obtain

∂ ∂τ q (τ, x) + ∂ ∂x (H(ϕ)q (τ, x)) = 0. (3.22)
The periodicity condition in θ implies the periodicity in x:

q (τ, 0) = q (τ, T ). Now, due to the expression of q given by (3.18) and the expression of x, we can rewrite H in terms of the new independent variables (τ, x) as follows:

H(ϕ(τ, θ 0 )) = c 0 ṡ(T ) T T 0 λ 1 (t)q (t -θ 0 -ϕ(τ, θ 0 ), τ, T ) ds = c 0 ṡ(T ) T T 0 λ 1 (t)q (τ, T -t + x) dt.
Therefore, we arrived to the following integro-differential system:

                       ∂ ∂τ q (τ, x) + ∂ ∂x (H(τ, x)q (τ, x)) = 0, τ > 0, x ∈ (0, T ), H(τ, x) = c 0 ṡ(T ) T T 0 q (τ, x -t)λ 1 (t) dt, τ ≥ 0, x ∈ [0, T ),
q (τ, 0) = q (τ, T ), τ > 0, q (0, x) = q 0 (x),

x ∈ [0, T ).

(3.23)

4. Special solutions. Some stability results.

We remind that the method we have used to derive the form of the system describing the evolution in time of phase densities is based on the assumption that the characteristics given by θ (t), and consequently x(τ ), are not crossing. Actually we will show in the following that this is not happening in finite time.

We define the characteristic lines in the following way: for every τ fixed, ∀x ∈ [0, T ), we denote by c[x 0 ](τ ) ( to avoid confusion) the characteristic line starting from x 0 that satisfies:

c[x 0 ](τ ) = x,
and is solution to

       d dτ c[x 0 ](τ ) = H(τ, c[x 0 ](τ )), τ > 0, c[x 0 ](0) = x 0 . (4.1)
Actually for τ fixed, ∀x ∈ [0, T ), the characteristic line which passes trough (τ, x) start either from a point (0, x 0 ) in the down boundary (putting x on the abscissa axis), or from a point (τ , 0) (or (τ , T ), depending on the sign of H); but due to the periodicity of x the second case can be extended to the first one. Then, along the characteristic lines defined by (4.1), the first equation in (3.23) writes as

d dτ q (τ, c[x 0 ](τ )) = - ∂H ∂c (τ, c[x 0 ](τ ))q (τ, c[x 0 ](τ )). (4.2)
Or, equivalently:

q (τ, x) = q (τ, c[x 0 ](τ )) = q 0 (x 0 ) exp{- τ 0 ∂H ∂c (τ , c[x 0 ](τ )) dτ }.
Here, as before, the solution on a characteristic that starts from the left or right boundary can be written as above due to the periodicity condition q(τ, 0) = q(τ, T ).

We want to see the way the characteristic lines depend on their initial state; to do that we calculate dc[x 0 ](τ )

dx 0 = 1 + τ 0 ∂H ∂c dc[x 0 ](τ ) dx 0 dτ , which yields dc[x 0 ](τ ) dx 0 = exp{ τ 0 ∂H ∂c (τ , c[x 0 ](τ )) dτ }.
Let us look closely to the expression of H; due to the commutativity of the convolution product, we can write:

H(τ, x) = c 0 ṡ(T ) T T 0 λ 1 (x -ξ)q (τ, ξ) dξ = c 0 ṡ(T ) T T 0 λ 1 (x -c[x 0 ](τ ))q (τ, c[x 0 ](τ ) dc = c 0 ṡ(T ) T T 0 λ 1 (x -c[x 0 ](τ ))q 0 (x 0 ) exp{- τ 0 H x (τ , c[x 0 ](τ )) dτ } dc
and due to the expression of dc[x 0 ](τ ) dx 0 , we obtain:

H(τ, x) = c 0 ṡ(T ) T T 0 λ 1 (x -c[x 0 ](τ ))q 0 (x 0 ) dx 0 = c 0 ṡ(T ) T T 0 λ 1 (x -c[x 0 ](τ ))q 0 (x 0 ) dx 0 .
Then, the following relation holds:

∂H ∂x (τ, x) = c 0 ṡ(T ) T T 0 ∂λ 1 ∂x (x -c[x 0 ](τ ))q 0 (x 0 ) dx 0 .
Now, since λ 1 is given by the adjoint linear system (3.15) and depends intrinsically on the functions of Izhikevich model computed along the limit cycle, we have that the derivative of λ 1 in the above expression exists and is bounded. Therefore the derivative of H with respect to x stays bounded and it implies that the characteristics starting from different initial points will not cross. In this way the solution defined along the characteristics is well defined for every τ .

The stationary solution.

Let us consider now the stationary problem:

               ∂ ∂x H 0 (x)q(x) = 0, x ∈ (0, T ) H 0 (x) = c 0 ṡ(T ) T T 0 λ 1 (t, 0)q(x -t) dt, x ∈ [0, T ) q(0) = q(T ). (4.3)
From the first equation of (4.3) we get that H 0 (x)q(x) = constant. (4.4)

In the case when H = 0, the obvious solution to the the above equation is given by q(x) = q, ∀x ∈ [0, T ] where q = constant. Replacing in H 0 , we obtain

H 0 = k 0 T 0 λ 1 (t, 0) dt,
where k 0 = c 0 ṡ(T )q. This solution corresponds to the a constant repartition of neurons with respect to their phase. If we denote by q := q -q, with q solution to (3.23) and q solution to (4.3), we get that q is a solution to the following problem:

               ∂ ∂τ q (τ, x) + ∂ ∂x (H(τ, x)q (τ, x)) = - ∂H ∂x (τ, x)q, τ > 0, x ∈ (0, T ), q (τ, 0) = q (τ, T ), τ > 0, q (0, x 0 ) = q 0 (x 0 ) -q, x 0 ∈ [0, T ). (4.5)
Let us assume now that q 0 has the support included in the region of [0, T ) where λ 1 is increasing. We can write the problem (4.5) as before along the characteristic line defined by (4.2):

               d dτ q (τ, c[x 0 ](τ )) = - ∂H ∂c (τ, c[x 0 ](τ )) (q (τ, c[x 0 ](τ )) + q) , τ > 0, x 0 ∈ [0, T ), q (τ, 0) = q (τ, T ), τ > 0, q (0, x 0 ) = q 0 (x 0 ) -q, x 0 ∈ [0, T ), (4.6
) where the expression of H is given by:

∂H ∂c (τ, c[x 0 ](τ )) = c 0 ṡ(T ) T supp q 0 ∂λ 1 ∂x (c[x 0 ](τ ) -c[ξ 0 ](τ ))q 0 (ξ 0 ) dξ 0 .
Under our hypothesis it follows that ∂H ∂c (τ, c[x 0 ](τ )) ≥ 0. Taking now into account that q ≥ 0 (being a density of population) and using a comparison result, we obtain that the solution q ≤ z, with z solution to the following problem:

               d dτ z(τ, c[x 0 ](τ )) = - ∂H ∂c (τ, c[x 0 ](τ ))z(τ, c[x 0 ](τ )), τ > 0, x 0 ∈ [0, T ), z(τ, 0) = z(τ, T ), τ > 0, z(0, x 0 ) = q 0 (x 0 ) -q, x ∈ [0, T ). (4.7) In fact z(τ, c[x 0 ](τ )) = z(0, x) exp{- τ 0 ∂H ∂c (τ , c[x 0 ](τ )) dτ },
and, making τ → +∞, we get that, in this case z(τ, x) → 0, as τ → +∞, which implies that q (τ, x) → 0, as τ → +∞, or, equivalently, q (τ, x) → q, as τ → +∞.

Then, it follows also that H(τ, x) → H 0 , as τ → +∞ (and, therefore c[x 0 ](τ ) → x 0 + H 0 τ as τ → +∞), therefore the solution of our problem will converge toward a constant repartition under the special hypothesis that we reminded earlier.

The synchronized solution.

As we stressed at the beginning, we have a special interest in the existence of a solution to (3.23) that expresses the synchronized state of the neural population. We will look for the solution in the following form q (τ, x) = q 0 δ(x -hτ -x 0 ) (

where δ the Dirac function and h will be determined later. Such a solution will express the fact that all the neurons in the population will be "concentrated" in the phase state x = hτ +x 0 , therefore the synchronization of the population. Let us check first that q(τ, x) defined by (4.8) satisfies (3.23). First of all, the corresponding expression of H is given by H(τ, x) = c 0 ṡ(T ) T T 0 λ 1 (t)q (τ, x -t) dt = c 0 ṡ(T ) T q 0 λ 1 (x -hτ -x 0 ).

Using the fact that λ 1 (x -hτ -x 0 )q (τ, x) = q 0 λ 1 (0)δ(x -hτ -x 0 ) one can check directly that, indeed, q satisfies ∂ ∂τ q (τ, x) + ∂ ∂x (H(τ, x)q (τ, x)) = 0, with H(τ, hτ + x 0 ) = h = c 0 ṡ(T ) T λ 1 (0)q 0 . We consider that the Dirac mass is concentrated at the initial moment in the initial state ξ 0 , and suppose that q 0 has the support included in the interval [ξ 0 -β, ξ 0 + β]. We take next two initial points such that ξ 0 -β ≤ x 1 0 ≤ x 2 0 ≤ ξ 0 + β; we are interested to see how it evolves in time the distance between the characteristic lines that start from x 1 0 respectively x 2 0 . We will show that if λ 1 is decreasing in the neighborhood of the origin, we have the convergence of the distance between the two characteristics to zero when τ tends to infinity. First of all, due to the fact that dc[x 0 ](τ ) dx 0 > 0 , we have that q (τ, c[x 0 ](τ )) → ∞, Therefore, if there exists a β > 0 such that λ 1 (t) > λ 1 (0) for T -β < t < T, λ 1 (t) < λ 1 (0) for 0 < t < β, the synchronized solution is stable, that is any initial condition with support in that phase interval will converge to a synchronized solution for some x 0 chosen as above. So the quantity λ 1 (0) characterizes the speed at which the phase is changed or in other words λ 1 (0)T is the variation of period due to the coupling.

The mechanisms leading to occurrence of such a solution as well as a more general result of the stability of it makes the object of further investigations in our future work.

c[x 1 0

 1 ](τ ) < c[x 2 0 ](τ ).If we denote ζ(τ ) = c[x 2 0 ](τ ) -c[x 1 0 ](τ ), ζ is a solution to dζ dτ = H(τ, c[x 2 0 ](τ )) -H(τ, c[x 1 0 ](τ ), ) ζ(0) = x 2 0 -x 1 0 .Writing H explicitly as before, it follows that:dζ dτ = H(τ, c[x 2 0 ](τ )) -H(τ, c[x 1 0 ](τ )) = k 0 ξ 0 +β ξ 0 -β λ 1 (c[x 1 0 ](τ ) -c[x 0 ](τ )) -λ 1 (c[x 2 0 ](τ ) -c[x 0 ](τ ))q 0 (x 0 ) dx 0 ,and one can see that under our assumptions we have that λ1 (c[x 1 0 ](τ ) -c[x 0 ](τ )) -λ 1 (c[x 2 0 ](τ )c[x 0 ](τ )) ≤ 0, therefore ζ(τ ) → 0 as τ → ∞.Moreover, from the expression of the solution along the characteristics it follows also:

Acknowledgements: We thank T. Boraud for his guidance in neuroscience. We acknowledge the support of CNRS by "programme neuroinformatique".