Branda Goncalves

AN INTERACTING NEURONAL NETWORK WITH INHIBITION: THEORETICAL ANALYSIS AND PERFECT SIMULATION

Keywords: spiking rate, interacting neurons, perfect simulation algorithm, classical contour technique

We study a purely inhibitory neural network model where neurons are represented by their state of inhibition. The study we present here is partially based on the work of Cottrell [3] and Fricker et al. [5]. The spiking rate of a neuron depends only on its state of inhibition. When a neuron spikes, its state is replaced by a random new state, independently of anything else and the inhibition state of the other neurons increase by a positive value. Using the Perron-Frobenius theorem, we show the existence of a Lyapunov function for the process. Furthermore, we prove a local Doeblin condition which implies the existence of an invariant measure for the process. Finally, we extend our model to the case where the neurons are indexed by Z. We construct a perfect simulation algorithm to show the recurrence of the process under certain conditions. To do this, we rely on the classical contour technique used in the study of contact processes, and assuming that the spiking rate lies on the interval [β * , β * ], we show that there is a critical threshold for the ratio δ = β * β * -β * over which the process is ergodic.

Introduction

For the operation of a neural network, neurons excite and or inhibit each other. Here, we study a model of a purely inhibitory neural network where neurons are represented by their inhibitory state. The study we present is partially based on the work of Cottrell [START_REF] Cottrell | Mathematical analysis of a neural network with inhibitory coupling[END_REF]. Her model consists of considering N interacting neurons described their state of inhibition. In her work, a neuron spikes when its state touches the value 0. When a neuron spikes, the state of inhibition of the other neurons increase by a non-negative deterministic constant θ. The spiking neuron immediately receives a random inhibition independent of anything else. In Cottrell's work the state of inhibition is just the waiting time until the next spike.

In the present work we generalize Cottrell's model in several natural ways. Actually, in Cottrell's model, the next spiking time in the neural net is deterministic and we will lift this assumption. A random spiking time is more realistic than deterministic one since stochasticity is present all over in the brain functioning. Secondly, to allow formal general models we allow the state of inhibition to decrease at a general rate in between the successive spikes of the network while in Cottrell's work the drift of flow is equal to -1.

In the first part of this paper, we consider systems of N interacting neurons, in which any neuron can spike at any time. The spiking neuron takes a new random state of inhibition, and the others increase their inhibitory state by a deterministic quantity that we will call the inhibition weight, which depends on the distance between the spiking neuron and the "receiving" neuron, so that a neuron located far away of the spiking neuron is not impacted by the spike. The model thus presented obviously extends Cottrell [START_REF] Cottrell | Mathematical analysis of a neural network with inhibitory coupling[END_REF] and Fricker et al. [START_REF] Fricker | Analysis of a Network Model[END_REF] in two ways: the spiking time is no more deterministic but it is random; the dynamic of the process is no more constant. Firstly, we show the existence of a Lyapunov function that allows us to formulate a sufficient condition of non-evanescence of the process in the sense of Meyn and Tweedie [START_REF] Meyn | Stability of Markovian processes III: Foster-Lyapunov criteria for continuous-time processes[END_REF], i.e. a condition ensuring that the process does not escape at infinity. To do so, we introduce a reproduction matrix H and we suppose the spectral radius of H is lower than 1. The eigenvector associated with the spectral radius of H allows us to find a Lyapunov function for the process.

Secondly, we study the recurrence of the process relying on Doeblin conditions which we establish for the embedded chain sampled at the jump times. We show the existence of an invariant probability measure for the process. We do this in the case the distribution of the new states has an absolutely continuous density and the jump rate is bounded.

In a second part, we consider the case where we have an infinite number of neurons indexed by Z (see Comets et al. [START_REF] Comets | Processes with long memory: Regenerative construction and perfect simulation[END_REF], Galves and Löcherbach [START_REF] Galves | Infinite systems of interacting chains with memory of variable length -a stochastic model for biological neural nets[END_REF] and Galves et al. [START_REF] Galves | Perfect simulation of infinite range Gibbs measures and coupling with their finite range approximations[END_REF]). In the work of Ferrari et al. [START_REF] Ferrari | Phase Transition for Infinite Systems of Spiking Neurons[END_REF], considering an infinite system of interacting point processes with memory of variable length, the authors investigated the conditions for the existence of a phase transition using the classical contour technique, based on the classical work of Griffeath [START_REF] Griffeath | The basic contact process[END_REF] on a contact process. Following the idea of Ferrari et al. [START_REF] Ferrari | Phase Transition for Infinite Systems of Spiking Neurons[END_REF] and Griffeath [START_REF] Griffeath | The basic contact process[END_REF], we construct a perfect simulation algorithm that allows us to show the recurrence of the process. Assuming that the spiking rate takes values in the interval [β * , β * ], we show that there is a critical threshold for the ratio δ = β * β * -β * over which the process is ergodic. This paper is organized as follows. In section 2 we describe the model and study the law of the first jump time of the process. The Foster-Lyapunov and Doeblin conditions are discussed to find nonevanescence criteria and to show the existence of the invariant probability measure of the process in section 3 which is our first main result. Finally, in Section 4, we present a perfect simulation algorithm and we simulate the law of the state of inhibition of a given neuron in its invariant regime.

The model

2.1. Description of the model. In our paper, let us consider we have N neurons that are related to each other. For all i ∈ {1, • • • , N }, X i,N t describes the state of inhibition of neuron i at time t. When the neuron i ∈ {1, • • • , N } spikes,

• The current state of inhibition of neuron i is replaced by a new value Y i independently of anything else with distribution F i . • The state of inhibition of any neuron j = i is increased by a positive value W i→j at time t.

In between successive jumps of the system, each neuron i follows the deterministic dynamic

. x i t = -α i x i t , x i 0 = x i , with α i x i continuous on [0, ∞), positive on (0, ∞) and non-negative on [0, ∞) and x = (x 1 , • • • x N ).
Let β i x i be a continuous positive and decreasing rate function on [0, ∞). We have taken β i to be decreasing so that the larger x i t is, the lower its probability of spiking and the smaller x i t is, the higher its probability of spiking.

We are thus led to consider the piecewise deterministic Markov process (PDMP)

X N t = (X 1,N t , • • • , X N,N t ) ∈ R N + . For i ∈ {1, • • • , N }, the dynamic of X i,N
t is given by:

(1)

dX i,N t = -α i X i,N t- dt + ∞ 0 ∞ 0 (y i -X i,N t-)1 {r≤β i (X i,N t-)} M i dt, dr, dy i + j =i W j→i ∞ 0 ∞ 0 1 {r≤β j (X j,N t-)} M j dt, dr, dy j ,
where M i is a random Poisson measure with intensity dtdrF i (dy) and for all i, the M i are all independent. This model extends that of Goncalves et al. [START_REF] Goncalves | On decay-surge population models[END_REF] in the multidimensional case.

Remark 1. For all i ∈ {1, • • • N }, X i,N
t can be interpreted as the inhibition state of the neuron i at time t and W j→i as the inhibition weight of the neuron j on the neuron i. When W i→j ≤ 0, we say that the neuron i is excitatory for the neuron j and when W i→j ≥ 0, we say that the neuron i is inhibitory for the neuron j. In our paper we are interested in the case where neuron i is inhibitory for neuron j i.e., W i→j ≥ 0.

Remark 2. The formula (1) is well-posed in the sense that there is non explosion of the process . Since β i (X i,N s ) ≤ β i (0) for all i we deduce that t 0 β i (X i,N s )ds < ∞ whence the non explosion, that is, almost surely, the process has only a finite number of jumps within each finite time interval.

The infinitesimal generator associated with this model is given by:

(2) G N V (x) = - N i=1 α i (x i ) ∂ ∂x i V (x) + N i=1 β i (x i ) ∞ 0 F i (dy i )[V (x + e i y i -e i x i + j =i e j W i→j ) -V (x)]
where V is a smooth function and e i is the i -th unit vector.

In other words, at each jump of the process, a single neuron spikes. If it is neuron i then its state is replaced by Y i and all other neurons receive the inhibition weight W i→j ≥ 0 for any j = i.

First jump time.

Let N i t be the counting process of successive jumps of neuron i, that is,

N i t = t 0 R+ R+ 1 {r≤β i (Xs(x i ))} M i (ds, dr, dy i )
and S i 1 the first jump time of neuron i, so we have

S i 1 = inf{t > 0|N i t = 1} and P(S i 1 > t) = e -t 0 β i (x i s (x i ))ds .
Let S 1 be the first jump time of the first neuron to jump, that is, S 1 = min i S i 1 . For all t > 0,

(3)

P(S 1 > t) = P(min i S i 1 > t) = N i=1 P(S i 1 > t) = N i=1 e -t 0 β i (x i s (x i ))ds .
Moreover, if t < min i t 0 (x i ) where

t 0 (x i ) := x i 0 dy α i (y)
is the time for the neuron i hit 0 starting from x i , we can write by making a change of variables that is no longer valid after touching 0, that

P(S 1 > t) = N i=1 e -[Γi(x i )-Γi(x i t (x i ))] ,
with Γ i (x i ) :=

x i γ i (y)dy and

γ i (x i ) = β i (x i )/α i (x i ). Assumption 1. Γ i (0) = -∞ for all 1 ≤ i ≤ N .
Proposition 1.

1. Suppose Assumption 1 holds. Then S 1 < ∞ almost surely. 2. Suppose Assumption 1 does not hold.

-If t 0 (x i ) < ∞ and α i (0) = 0 then S 1 < ∞ almost surely if and only if β i (0) > 0 for all i.

-If t 0 (x i ) = ∞ then P(S 1 = ∞) > 0 i.e. with a positive probability the first jump time is infinite.

Proof. Let N be fixed and suppose Assumption 1 holds.

If t 0 (x i ) = ∞ and letting t tend to ∞ in (3) we have

P(S 1 = ∞) = N i=1 e -[Γi(x i )-Γi(x i ∞ (x i ))] = N i=1 e -[Γi(x i )-Γi(0)] , since x i ∞ (x i ) = 0. Then P(S 1 = ∞) = 0 that is S 1 < ∞ almost surely. If t 0 (x i ) < ∞
and letting t ↑ min i t 0 (x i ) in (3), we obtain

P(S 1 ≥ min i t 0 (x i )) = lim t↑mini t0(x i ) P(S 1 > t) = N i=1
e -[Γi(x i )-Γi(0)] = 0 implying that S 1 < ∞ almost surely.

Suppose Assumption 1 does not hold. If α i (0) = 0 (this means that the flow of the process brings us to 0 at most ) the time for the neuron i hit 0 starting from x i is finite i.e t 0 (x i ) < ∞ then it is obvious (by definition of t 0 (x i )) to see that it is enough that β i (0) > 0 to have S 1 < ∞ almost surely.

If Assumption 1 does not hold and t 0 (x i ) = ∞ then by making t → ∞ in (3) we have P(S 1 = ∞) > 0 that is S 1 = ∞ with a positive probability.

We finish this section with a simulation of the process starting from some fixed initial configuration (x 1 0 , • • • , x N 0 ). For this, we assume that for all i the jump rate β i (x i ) is bounded and lower bounded, that is,

β i (x i ) ∈ [β * , β * ] for all x i > 0, where 0 < β * < β * < ∞.
The following variables will be used to write our simulation algorithm.

• T is the time vector • L is the label associated with T. It will be {sure} or {uncertain} • P = (P 1 , • • • , P N ) is the vector of states of the N neurons at a fixed instant • I is the vector which represents the number of the neuron which spikes.

Algorithm (1) We set T 1 ∼ exp(β * N ) -L 1 = {uncertain} with probability β * -β * β * -L 1 = {sure} with probability β * β *
(2) We initialize the vector P with the values (

x 1 0 , • • • , x N 0 ) (3) We choose I 1 = k with probability 1 N -If L 1 = {sure}, (a) P k ∼ F k , P i ← x i T1 (P i ) + W k→i -If L 1 = {uncertain} we accept the jump with probability p = β k (x k T1 (P k )) -β * β * -β * ,
and we apply (a).

-Else P j ← x j T1 (P j ), ∀j ∈ {1, • • • , N }. (4) We update the vector P and start the procedure again from (1).

We stop the procedure after a fixed finite number n of iterations.

We plot in the following figure a typical trajectory of X i,N t with N = 2 neurons. In both figures N = 2 neurons and n = 50 iterations. α(x) = x, β(x) = 3 + 1(x ≤ 2). In the figure on the left, W j→i = i/N for all j = i and in the figure of right, W j→i = 1/2 for all j = i.

Foster-Lyapunov and Doeblin conditions

In this section, we want to find conditions of non-evanescence of the process and show the existence of an invariant probability measure of the process.

3.1. Foster-Lyapunov condition. We suppose that γ i := β i /α i is bounded and we define W the matrix of inhibition weight by W ij := W j→i , i = j and W ii = 0.

It is further assumed that the matrix W is irreducible in the sense that there exists an integer p > 0 such that W p > 0. We introduce the reproduction matrix

H ij = W j→i γ i ∞ , i = j, H ii = γ i ∞ ∞ 0 y i F i (dy i )
which is also irreducible.

Suppose that

ρ(H) < 1
where ρ(H) is the largest eigenvalue of H that is the spectral radius of H. Then (Perron Frobenius) there exists a left eigenvector κ associated to this eigenvalue ρ, that is, for all i,

j κ j H ji = ρκ i .
On the other hand, put

m i = κ i γ i ∞ and et V (x) = m i x i .
Finally, let V : R N + → R such that V (x) = N i=1 m i x i and we recall that the infinitesimal generator is given by:

G N V (x) = - N i=1 α i (x i ) ∂ ∂x i V (x) + N i=1 β i (x i ) ∞ 0 F i (dy i )[V (x + e i y i -e i x i + j =i e j W i→j ) -V (x)]
So by replacing V by its expression in the infinitesimal generator G N V (x) we have:

G N V (x) = - N i=1 α i (x i )m i + N i=1 β i (x i ) ∞ 0 dF i (y i )[ N j=1,j =i (W i→j + x j )m j + y i m i - N j=1 x j m j ] = - N i=1 α i (x i )m i + N i=1 β i (x i )(m i ∞ 0 y i F i (dy i ) + j =i W i→j m j ) - N i=1 β i (x i )x i m i Then, since -β i (x i )x i ≤ 0, G N V (x) ≤ - i α i (x i )m i + i β i (x i )(m i ∞ 0 y i F i (dy i ) + j =i W i→j m j ) = - i α i (x i )   m i -γ i (x i )   m i c i H ii + j =i c j H ji m j     = - i α i (x i )   m i -γ i (x i )   κ i H ii + j =i κ j H ji     = - i α i (x i ) γ i ∞ κ i -γ i (x i )ρκ i = - i α i (x i ) γ i ∞ κ i 1 - γ i (x i ) γ i ∞ ρ .
Definition 1. We call the process non evanescent if there exists a compact K such that for all x, P xalmost surely,

lim sup t 1 K (X t ) = 1. Proposition 2. If ρ < 1, then the process is non-evanescent. Proof. V (x) defined in G N V (x) above is a norm-like function and G N V (x) < 0. The theorem (CD0)
of Meyn and Tweedie [START_REF] Meyn | Stability of Markovian processes III: Foster-Lyapunov criteria for continuous-time processes[END_REF] implies the result.

Example 1. (Mean-field interaction) Suppose we have N neurons. We suppose also that γ i = γ, where γ is bounded and F i = F, W j→i = θ for all i. In this case the reproduction matrix is

H ij = θ γ ∞ , i = j, H ii = γ ∞ E(Y ). Suppose ρ(H) is the spectral radius of H. Then, ρ(H) = γ ∞ (E(Y ) + (N -1)θ) and its associated eigenvector is κ = (1, • • • , 1). The condition ρ(H) < 1 is therefore equivalent to γ ∞ (E(Y ) + (N - 1)θ) < 1. For θ = ε/N with ε > 0 the condition ρ(H) < 1 becomes γ ∞ (E(Y ) + ε) < 1.
Example 2. (Torus) Suppose we have N ≥ 3 neurons such that each neuron interacts with its two nearest neighbors (its left and right neighbors). Neuron 1 interacts with neuron 2 and neuron N .

Neuron N interacts with neuron N -1 and neuron 1, so we have a torus.

We suppose for all i, we have γ i = γ, bounded and

F i = F. W j→i = θ for all j ∈ {i + 1, i -1} and W j→i = 0 if j = {i + 1, i -1}.
In this case the reproduction matrix is

H ij =      θ γ ∞ , if i = j, j ∈ {i + 1, i -1} 0 , if i = j, j = {i + 1, i -1} γ ∞ E(Y ) , if i = j. If ρ(H) is the spectral radius of H then ρ(H) = γ ∞ (E(Y ) + 2θ) and its associated eigenvector is κ = (1, • • • , 1). The condition ρ(H) < 1 is equivalent to γ ∞ (E(Y ) + 2θ) < 1. 3.2. Doeblin condition. Let S 0 < S 1 < • • • < S n < • • •
be the instants of successive jumps of the N neurons. It is obvious that the embedded chain Z n := X Sn is a Markov chain. Let I n be the index of the neuron which jumps at time S n . Proposition 3. Suppose that the assumptions of proposition 1 hold. Then, (Z n , I n ) is a Markov chain and its transition Q(x, dy) is given by:

(4) P(Z n ∈ dy, I n = j|Z n-1 = x, I n-1 = i) = ∞ 0 dse -s 0 dl N i=1 β i (x i l (x i )) β j (x j s (x j )) × F j (du)δ (x 1 s (x 1 )+Wj→1,••• ,x j-1 s (x j-1 )+Wj→j-1,u, x j+1 s (x j+1 )+Wj→j+1,••• ,x N s (x N )+W j→N ) (dy).
Theorem 4. Suppose for all 1 ≤ i ≤ N, α i ∈ C 1 and there exists a compact set

K ⊂]0, ∞[ N such that for all x ∈ K, for all 1 ≤ i ≤ N, β i (x i + i-1 j=1 W j→i ) > 0.
Moreover we suppose that F i (dy) is absolutely continuous and β i ∞ < ∞ for all i. Then there exist d ∈ (0, 1) and a probability measure

ν on (IR N + , B(IR N + )), such that (5) Q N (x, dy) ≥ d1 K (x)ν(dy)
where Q is the transition operator of embedded chain Z n = X Sn and Q N is its N -th iterate.

To prove the above result we fix any deterministic sequence

s 1 < • • • < s N .
In the sequel we shall work on the event

S 1 = s 1 , • • • , S N = s N , I 1 = 1, • • • , I N = N and Y 1 = y 1 , • • • , Y N = y N .
This means that the jumps are ordered such that neuron 1 jumps before neuron 2 and etc. Let y = (y 1 , • • • , y N ) where y i is the new state of inhibition of neuron i after the spike.

Let

t k = s k -s k-1 for all 1 ≤ k ≤ N the inter jump times of the N neurons which implies that s k = t 1 + • • • + t k .
Conditionally on this event, let Ψ s N be the vector of states of the process at time s N . We can define Ψ s N as a function of the states y 1 , • • • , y N such that Ψ s N : R N → R N is given by:

Ψ k s N (y) = ψ k,N t N • • • • • ψ k,k+1 t k+1 (y k ) , if 1 ≤ k < N y N , if k = N
where for all l = k,

(6) ψ k,l s (u) = x k s (u) + W l→k and x k s (u) means the solution of the deterministic dynamic . x k s = -α i (x k s ), x k 0 = u. Remark 3.
In the definition of Ψ k s N (y), we note that it depends only on y k . Therefore we have for all i = j, ∂Ψ i s N ∂y j = 0.

Proposition 5. For all 1 ≤ k ≤ N let α k be a globally Lipschitz function. For all y ∈ R N + , there exists an open neighborhood B of y such that

Ψ s N : B → Ψ s N (B) is a local diffeomorphism.
Proof. Let J Ψs N (y) be the Jacobian matrix of Ψ s N (y). Using the remark 3 we have :

det(J Ψs N (y) ) = det          ∂Ψ 1 s N (y) ∂y 1 • • • ∂Ψ 1 s N (y) ∂y N . . . . . . . . . ∂Ψ N s N (y) ∂y 1 • • • ∂Ψ N s N (y) ∂y N          = det         ∂Ψ 1 s N (y) ∂y 1 0 • • • 0 0 . . . . . . . . . . . . 0 0 • • • 0 ∂Ψ N s N (y) ∂y N         .
We obtain det(J Ψs N (y) ) = 0 if and only if

N j=1 ∂Ψ j s N (y) ∂y j = 0 that is ∂Ψ j s N (y) ∂y j = 0, ∀1 ≤ j ≤ N. It is obvious to see that ∂Ψ N s N (y) ∂y N = 1.
For all 1 ≤ j ≤ N -1, we have:

(7) ∂Ψ j s N ∂y j (y) = N -j i=1 exp - s N -(i-1) s N -i α j x j s Ψ j s N -i (y) ds = 0. It means that | det(J Ψs N (y) )| = 0 then Ψ s N (y) is a local diffeomorphism.
Localizing, we may therefore conclude that for each y there exists B such that Ψ s N :

B → Ψ S N (B) is a diffeomorphism.
Proof of theorem 4. Let ε > 0 fixed. We will work on the event

E = {S 1 ≤ ε, • • • , S n+1 -S n ≤ ε, ∀n < N : (I 1 , • • • , I N ) = (1, • • • , N )}.
In particular, on E, the index I n of the n-th neuron is equal to n for all n ∈ {1, • • • , N }.

Knowing that the first jump takes place at time S 1 = s 1 , the probability that the index I 1 of the first jump is equal to 1 is given by:

P(I 1 = 1|S 1 = s 1 ) = P(S 1 1 < S j 1 , ∀j = 1) = β 1 (x 1 s1 (x 1 )) N j=1 β j (x j s1 (x j ))
.

We want to compute, P(

I 1 = 1, • • • , I N = N |S 1 = s 1 , S 2 = s 2 , • • • , S N = s N ).
To obtain a compact formula, using formula (6) we define

φ k j (x k , y k , s 1 , • • • , s N ) = ψ k,j tj • • • • • ψ k,k+1 t k+1 (y k ), if 1 ≤ k ≤ j -1 ψ k,j tj • • • • • ψ k,1 t1 (x k ), if j ≤ k ≤ N
giving the states of neuron k at time S j depending on whether neuron k jumped before or after time

S j . Let x k j = x k tj (φ k j-1 (x k , y k , s 1 , • • • , s N
) be the state of neuron k before the j -th jump. We know that as long as neuron k has not yet jumped, it receives each time a quantity W j→k , ∀j = k from the other neurons that jumped before it. So knowing all the jump times where other neurons jumped, we have:

P(I 1 = 1, • • • , I N = N |S 1 = s 1 , S 2 = s 2 , • • • , S N = s N ) = β 1 (x 1 s1 (x 1 )) N i=1 β i (x i s1 (x i )) R N -1 + N i=2 β i (x i i ) N i=2 ( N k=1 β k (x k i )) N -1 k=1 P(Y k ∈ dy k ).
For any Borel subset B of R N we have

Q N (x, B) ≥ P x (Z N ∈ B, E) = [0,ε] N dt 1 • • • dt N R N F 1 (dy 1 ) • • • F N (dy N ) × ( N k=1 β k (x k k ))e -s N 0 β k (Ψ k t (y))dt 1 B (Ψ s N (y)).
Remark that on the event E,

x k k ≤ x k + k-1 j=1 W j→k . Recall β k is decreasing function and let µ k = inf k {β k (x k + k-1 j=1 W j→k ) : x ∈ K} the lower-bound on K of β k (x k + k-1 j=1 W j→k ). Using the fact that β i ∞ < ∞ for all i, let c = ( N k=1 µ k )e -N β i ∞N ε . Then we have (8) Q N (x, B) ≥ c [0,ε] N dt 1 • • • dt N R N F 1 (dy 1 ) • • • F N (dy N )1 B (Ψ s N (y)).
Following the arguments of Benaïm et al. [START_REF] Bénaïm | Qualitative properties of certain piecewise deterministic Markov processes[END_REF], for any t * ≤ N ε, there exists a ball B r (t * ) of radius r, of center t * and an open set I ⊂ R N such that we can find for all s N ∈ B r (t * ), an open set W s N ⊂ R N :

Ψs N : W s N → I y → Ψ s N (y)
is a diffeomorphism (see Benaïm et al. [START_REF] Bénaïm | Qualitative properties of certain piecewise deterministic Markov processes[END_REF], Lemma 6.2). In the above formula, Ψs N denotes the restriction of Ψ s N to W s N . This allows us to apply the theorem of a change of variables in the inequality [START_REF] Goncalves | On decay-surge population models[END_REF].

α j x j s Ψ j s N -i (y)
is upper bounded since α j is a global Lipschitz function. Then, for all 1 ≤ j ≤ N -1 we obtain:

∂Ψ j s N ∂y j (y) ≤ exp(N -j)ε α j ∞ .
Then, ∀y ∈ W S N , c | det(J Ψs N (y) )| -1 ≥ c > 0 and the inequality (8) becomes :

Q N (x, B) ≥ c [0,ε] N dt 1 • • • dt N B F 1 (dy 1 ) • • • F N (dy N )1 B (Ψ s N (y)) ≥ c Br(t * ) dt 1 • • • dt N Ws N ∩B dy1 B ( Ψs N (y))| det(J Ψs N (y) )| ≥ c λ(B r (t * )) I 1 B (x)dx = d1 B (x)ν(I)
where d = c λ(B r (t * )) with λ(B r (t * )) the Lebesgue measure of the ball B r (t * ) and ν(I) the uniform measure of I.

Corollary 6. If for all k ≤ N, β k is strictly lower-bounded and bounded, then the process is recurrent.

Remark 4. When β k is strictly lower-bounded and bounded, we can notice that the lower bound obtained in theorem 4 holds on the whole state space R + , that is, without 1 K . This allows us to have the global lower bound Q N (x, dy) ≥ dν(dy) and thus the uniform ergodicity of the process.

Perfect simulation

In this section, we consider a framework with an infinity of neurons indexed by Z. We want to build a perfect simulation algorithm to show in another way the recurrence of our process under certain conditions. Let V .→i = {j : W j→i = 0} and V i→. = {j : W i→j = 0} be the incoming and out-coming neighborhoods of the neuron i (see Comets et al. [START_REF] Comets | Processes with long memory: Regenerative construction and perfect simulation[END_REF] and Galves and Löcherbach [START_REF] Galves | Infinite systems of interacting chains with memory of variable length -a stochastic model for biological neural nets[END_REF]).

We consider the case where each neuron has a finite number of neighbors.

We assume that for all i the jump rate

β i (x i ) is bounded, that is, β i (x i ) ∈ [β * , β * ] for all x i > 0, where 0 < β * < β * < ∞.
The following variables will be used to write the perfect simulation algorithm:

-T is the time vector -P is the matrix of states where each row of this matrix represents the different states of the N neurons at a fixed instant -I is the vector which represents the index of the neuron which spikes.

We fix a neuron i ∈ Z and in what follows we are interested in finding the state of i at time 0 in the stationary regime, that is, assuming that the process starts from -∞. To do so we explore the past in order to determine all sets of indices and times which affect the value of neuron i at time 0. The set of all such couples (j, s) will be called the clan of ancestors of neuron i (see Galves and Löcherbach [START_REF] Galves | Infinite systems of interacting chains with memory of variable length -a stochastic model for biological neural nets[END_REF], Galves et al. [START_REF] Galves | Perfect simulation of infinite range Gibbs measures and coupling with their finite range approximations[END_REF]). The clan of ancestors is a process that evolves in time by successive jumps. We start with C i 0 = {i} and in the following we will define the updates of this process at the time of the jumps. More precisely we do the following:

-We simulate , ∀ l ∈ Z, N l,s t and N l,p t two Poisson processes with respective intensities β * and β * -β * . The jump times of N l,s t and N l,p t are respectively T l,s n and T l,p n for the neuron l after n jumps.

-Let i ∈ Z be fixed and T 1 = inf{T j,s 1 , T j,p 1 > 0 : j ∈ V .→i , T i,s 1 > 0} where V .→i is the incoming neighborhood of i.

-

If T 1 = T j,s 1 , we set C i T1 = {i} and I 1 = j. -If T 1 = T j,p
1 , we set C i T1 = {i, j} and we set I 1 = j. -If T 1 = T i,s 1 , we set C i T1 = ∅ and we stop the algorithm. In this case we set I 1 = i.

-Suppose T n is the n -th jump time of C i Tn . Then,

T n+1 = inf{T j,s m , T j,p m > T n : ∃l ∈ C i Tn , j ∈ V .→l , T k,s m > T n , k ∈ C i Tn }. -If T n+1 = T j,s
m we set I n+1 = j and then

C i Tn+1 = C i Tn -If T n+1 = T j,p m we set I n+1 = j and then C i Tn+1 = C i Tn ∪ {j} -If T n+1 = T k,s m we set I n+1 = k and then C i Tn+1 = C i Tn \ {k} where k ∈ C i Tn .
We stop the procedure at time T i stop = inf{t : C i t = ∅}. To ensure that the algorithm stops it will be necessary to find a criterion so that T i stop < ∞. This will be done in Theorem 8 below. The above algorithm is called the backward procedure.

In the following we will write a forward procedure of the process in case where each neuron has a finite number of neighbors. For this we define:

N i stop = inf{n > 0 : C i Tn = ∅}, Ci = ∪ N i stop n=0 C i Tn and ∂ ext (C i t ) = {j / ∈ C i t : ∃k ∈ C i t , j → k} where N i
stop is the number of steps of the backward procedure, Ci is the union of all clans of ancestors up to N i stop and ∂ ext (C i t ) is the set of neurons not belonging to the clan of ancestor of neuron i but having an interaction with at least one neuron in the ancestor clan of neuron i.

In this algorithm, we will rely on the a priori realizations of the processes N i,s t and N i,p t . Algorithm (forward procedure)

(1) We initialize the set of sites for which the decision to accept can be made by

S i = {(I m , T m ) ∈ Ci × R + , C Im Tm = ∅} For n = N i stop we have P In n ∼ F In . Starting from n → n -1 : (2) If (I n-1 , T n-1 ) ∈ S i then P In-1 n-1 ∼ F In-1 .
-If for j ∈ V In-1→. , we have j ∈ C i Tn-1 then

P j n-1 = x j Tn-Tn-1 (P j n ) + W In-1→j -If for j / ∈ V In-1→.
, we have j ∈ C i Tn-1 then

P j n-1 = x j Tn-Tn-1 (P j n ) (3) If k := I n-1 ∈ ∂ ext (C i
Tn-1 ), we have

P l n-1 = x l Tn-Tn-1 (P l n ) + W k→l where there exists l such that k → l ∈ C i Tn-1 (4) If (I n-1 , T n-1 ) ∈ ( Ci × R + ) \ S i then:
-We decide according to the probabilities

p = β(x
In-1

Tn-Tn-1 (P

In-1 n )) -β * β * -β *
to accept the presence of a spike of neuron I n-1 .

We update

S i ← S i ∪ {(I m , T m ) ∈ Ci × R + , C Im
Tm ⊂ S i } and go back to step 2.

-Else with the probabilities 1 -p we reject the presence of a spike of neuron I n-1 and P

In-1 n-1 = x

In-1 Tn-Tn-1 (P

In-1 n

). We consider all the elements of S i and we always start with the last element to get out of the clan. The update of S i allows us to start the procedure again.

We stop the procedure when all the elements of Ci are filled.

Remark 5. For any site (i, t) ∈ Z × R + , C i t is a Markov jump process taking values in the finite subset of Z (see Galves et al. [START_REF] Galves | Perfect simulation of infinite range Gibbs measures and coupling with their finite range approximations[END_REF]) and its infinitesimal generator is given by

A clan g(C) = j∈C β * [g(C \ {j}) -g(C)] + j∈∂ext(C) (β * -β * )[g(C ∪ {j}) -g(C)]
where g is a test function.

Proposition 7. Let d j = min x j β j (x j ), d j = max x j β j (x j ) and b j = k→j (d k -d k ). If b j < d j for all j then T i stop is finite almost surely i.e. the process is subcritical.

Proof. We shall construct a process Z n such that for all n, |C i Tn | ≤ Z n and such that Z n evolves as follows: with probability j∈Zn dj j∈Zn (bj +dj ) we have Z n+1 = Z n -1 and with probability j∈Zn bj j∈Zn (bj +dj ) we have

Z n+1 = Z n + 1.
In this general case where a neuron has a finite number of neighbors (more than two neighbors) with which it interacts, we can say no more than proposition 7. Thus, in the following, we put ourselves in the case where each neuron i has exactly two neighbors so that the neuron i interacts only with the neurons i + 1 and i -1. In other words, the incoming neighborhood of i is V .→i = {i + 1, i -1}.

Algorithm (backward procedure)

(1) We simulate , ∀ l ∈ Z, N l,s t and N l,p t two Poisson processes with respective intensities β * and β * -β * . The jump times of N l,s t and N l,p t are respectively T l,s n and T l,p n for the neuron l after n jumps. The jump times T l,s n will be considered as times of sure jumps (counted by the process N l,s t ) and the jump times T l,p n will be considered as times of possible jumps (counted by the process N l,p t )

(2) Let i ∈ Z fix and

T 1 = inf{T i±1,s 1 , T i±1,p 1 , T i,s 1 }. We set I 1 = i ± 1 -If T 1 = T i±1,s 1 , we set C i T1 = i and we put I 1 = i ± 1. -If T 1 = T i±1,p 1 , we set C i T1 = {i, i ± 1}. We put I 1 = i ± 1 -If T 1 = T i,s 1 ,
we set C i T1 = ∅ and we stop the algorithm. We put I 1 = i.

(3) Suppose T n is the n-th jump time of C i Tn . We have:

T n+1 = inf{T j,s m , T j,p m > T n : |j -C i Tn | ≤ 1, T k,s m > T n , k ∈ C i Tn }. -If T n+1 = T j,p m we set: If j ∈ C i Tn , C i Tn+1 = C i Tn If j / ∈ C i Tn , C i Tn+1 = C i Tn ∪ {j} -If T n+1 = T k,s m we set: If k ∈ C i Tn , C i Tn+1 = C i Tn \ {k} If k / ∈ C i Tn , C i Tn+1 = C i Tn
We update C i t and start the procedure again. We stop the procedure at time T i stop where T i stop = inf{t : C i t = ∅}. Indeed, the whole procedure makes sense only if T i stop < ∞ almost surely. Remark 6. The forward procedure is the same as in the first case where each neuron has a finite number of neighbors.

The following theorem gives conditions on the finitude of the extinction time.

Theorem 8. We set δ = β * β * -β *
. There exists a critical value 0 < δ c < ∞ such that:

-if δ > δ c , then the extinction time is finite almost surely that is, P(∀i, T i stop < ∞) = 1 -if δ < δ c , then the extinction time is infinite with a positive probability that is, P(∀i, T i stop = ∞) > 0.

Proof. We first show that T i stop < +∞ almost surely for sufficiently large δ. We observe that we can upper bound |C i t | (where |C i t | is the cardinal of C i t ) by Z t almost surely for all t ≥ 0 where Z 0 = 1 and (Z t ) t≥0 is a branching process. With a rate n(β * -β * ) the transition from Z t is from n to n + 1 and with a rate nβ * this transition is from n to n -1.

We can therefore define for any bounded test function f, the associated infinitesimal generator of (Z t ) t≥0 as follows :

Af (n) = n[(β * -β * ) (f (n + 1) -f (n)) + β * (f (n -1) -f (n))].
Take f (n) = n, we obtain :

Af = f [(β * -β * ) -β * ] = f (β * -β * )(1 -δ).
Then, for δ > 1, we have

Af (n) = -cf (n) where -c = (β * -β * )(1 -δ). Assuming x t = E(f (Z t )
) and using the Itô formula, we have:

x t = x 0 + E t 0 Af (Z s )ds = x 0 -c t 0 x s ds = x 0 e -ct .
Therefore, when t → ∞, we have x t → x 0 = 1. Which implies that if δ > 1, P(T i stop < ∞) ≥ P(lim t→∞ Z t = 0) = 1 thus ensuring that δ c < 1.

We now show that for all δ < δ c , T i stop = +∞ with positive probability. For this proof, we will use the classical graphical construction of C i t (see Ferrari et al [START_REF] Ferrari | Phase Transition for Infinite Systems of Spiking Neurons[END_REF], Griffeath [START_REF] Griffeath | The basic contact process[END_REF]). We work within the space-time diagram Z × [0, ∞[. For each i ∈ Z, we consider N i,s t and N i,p t two independent Poisson processes with respective intensities β * and β * -β * . The jump times of N i,s t and N i,p t are respectively T i,s n and T i,p n for the neuron i after n jumps. For each i ∈ Z, we draw graphical sequences as follows. First draw arrows pointing from (i -1, T i,p n ) to (i, T i,p n ) and from (i + 1, T i,p n ) to (i, T i,p n ) for all n ≥ 1, i ∈ Z. Second, δ's at all (i, T i,s n ), for all n ≥ 1, i ∈ Z. We also suppose that time is going up which implies that we thus obtain a random graph P. Let us say that there is a chain of vertical upward and horizontal directed edges in the random graph that leads from (i, 0) to (j, t) ( with j ∈ {i + 1, i -1}) without passing through a δ. Notice that C i t is the set of the clan of ancestors of site (i, t), that is C i t = {j : there is a path from (i, 0) to (j, t) for j = i ± 1}. It is obvious to notice that T i stop < ∞ if and only if Ci = ∪ t≥0 C i t is a finite set. We will therefore show that P(T i stop < ∞) = P(| Ci | < ∞) < 1 for sufficiently small values of δ using classical contour techniques. (see Griffeath [START_REF] Griffeath | The basic contact process[END_REF].) For this, on | Ci | < ∞, we draw the contour of Ci as follow.

Starting from (i-1 2 , 0). Let Γ be a possible path of the graph P. Γ consists of 4n alternating vertical and horizontal edges for some n ≥ 1 which we encode as a succession of direction vectors (D 1 , • • • , D 2n ). Each of the D i can be one the seven triplets dld, drd, dru, ulu, uru, urd, dlu, where d, u, l and r stand for down, up, left and right, respectively. Note that uld cannot occur in a possible path Γ because the direction of uld is counter-clockwise. We start at (i - (for more details, see Ferrari et al. [START_REF] Ferrari | Phase Transition for Infinite Systems of Spiking Neurons[END_REF].)

We first observe that the occurrence of either uru, urd, or drd can be upper bounded by δ. This is due the fact that the probability associated with uru or drd is δ 1+2δ and that of urd is δ 2+δ . In the same way, we observe that the occurrence of either dld, ulu or dlu can be upper bounded by 1. Indeed, the associated probability with its directions is 1 1+δ . Therefore we obtain the following list of upper bounds uru occurs with probability at most δ urd occurs with probability at most δ drd occurs with probability at most δ dru occurs with probability at most 1 dld occurs with probability at most 1 ulu occurs with probability at most 1 dlu occurs with probability at most 1.

In the above list, we have upper bounded the probability associated with dru which is given by δ 3δ = 1 3 , by 1. For a given contour having 4n edges, with n ≥ 2, its probability is therefore upper bounded by δ N (drd)+N (uru)+N (urd) = δ n-N (dru) ≤ δ n/2 . Indeed, for each triplet we have 4 possible choices. The first entry of a given triplet is always fixed by the previous triplet in the sequence, and for the first triplet D 1 the first entry is always u.

Then, for n = 1, the probability of appearance of a contour of length 4 is equal to P(D 1 = urd) = δ 2+δ ≤ δ. We also have, for n = 2, the probability of appearance of a contour of length 8 is equal to

P(D 1 = ulu, D 2 = urd, D 3 = drd) + P(D 1 = ulu, D 2 = uru, D 3 = urd)+ P(D 1 = uru, D 2 = urd, D 3 = dld) + P(D 1 = urd, D 2 = drd, D 3 = dld) ≤ 4δ 2 .
Remark 7. In the above probabilities, we have not put the direction D 4 = dlu because it is a certain direction. It is common to all possible paths and its probability of occurrence is 1.

Therefore, a very approximate upper bound on the total number of possible triplets (D 1 , • • • , D 2n ) is given by 4 2n = 16 n . We get for all δ < 1 (16) 2 ,

P(T i stop < ∞) ≤ δ + 4δ 2 + n≥3 (16) n δ n/2 = δ + 4δ 2 + (16 √ δ) 3 1 -16 √ δ . We set φ : δ → φ(δ) = δ + 4δ 2 + (16 √ δ) 3 1-16 √ δ . Then, P(T i stop < ∞) ≤ φ(δ).
As δ → 0, φ(δ) → 0 which implies that there exists δ c such that φ(δ c ) = 1. As a consequence, P(T i stop < ∞) < 1, ∀ 0 < δ < δ c . We therefore conclude that δ c exists and 0 < δ c < 1. 4.1. Some simulations. We simulate the state X 0 (i) in the stationary regime for a fixed neuron i ∈ Z at time 0 and estimate its density. The main purpose of this simulation is to have an idea about the theoretical distribution of X 0 (i) in its stationary regime and whether this distribution is impacted by the specification of F i .

We denote by D the set of neurons which belong to a clan of ancestors of neuron i at a time t or to its neighborhood.

To do this, we apply the following algorithm:

(1) Initialize the family V .→i of non empty neighborhoods of the neuron i (2) Initialize C i 0 = i the clan of ancestors of neuron i at time t = 0.

(

i tnext = C i t \ {i} -If j ∈ C i t and the jump is possible C i tnext = C i t -If j ∈ V (C i t ) (where V (C i t ) := ∪ j∈C i t V .→j ) and the jump is sure, then C i tnext = C i t -If j ∈ V (C i t ) and the jump is possible C i tnext = C i t ∪ {j} -We update t ← t next end While.
(5) We determine the chronological list of the different jump times from 0 to the last time which makes the clan empty.

-For each of these jump times, we indicate the associated neuron and the nature of the jump.

-If the jump is sure, we simulate a random state following a distribution F i at the neuron associated with this jump time. -Let dt = T r-s+1 -T r-s and j the neuron associated with the jump time T r-s . The state of k at time T r-s is x k -dt (x) + W j→k * 1 {sure jump of j at Tr-s} with W j→k the inhibition weight of j on k.

-We determine rather the occurence is effective or not of the jump of k at time T m thanks to its state at time T m+1 . 1-If the jump is effective, we simulate a random state for k at time T m following a distribution F i . Otherwise, we determine the state of k at time T m as x Tm-Tm+1 (x) where x = state of k at time T m+1 . Let m be the new rank of the last possible jump time of D and repeat the procedure.

end While.

Remark 8. After this step, we know the exact nature of all jumps.

(8) Determine for neuron i its first safe jump time T n where n is the rank of this time in the chronology of jump times.

(9) The state of neuron i is determined recursively from its state at time T n to T 0 as follows:

-For s ∈ {1, • • • , n -1} let x = state of neuron i at time T n-s+1 .

-Let dt = T n-s+1 -T n-s and j the neuron associated with the jump time T n-s . The state of neuron i at time T n-s is x i -dt (x) + W j→i * 1 {sure jump of j at Tr-i} with W j→i the inhibition weight of j on i.

Remark 9. The last value determined is the initial state of neuron i.

In the three following examples we consider α i (x) = x, β i (x) = 3 + 1 {x≤2} , W i→j = 1. To verify if the distribution of inhibition state depends on the distribution F i , we consider three different distributions for F i that are E(1), E(10) and 0.5δ 1 +0.5δ 2 . We simulate, with the algorithm described above N = 1000 values for the inhibition state . We then estimate non parametrically the distribution of the inhibition state in these three cases of distribution F i and we compare them.

The stationary distribution of the process in the three following cases seems to be continuous. We do not provide a proof here, this is outside the scope of this paper. We can remark that the distribution of state of inhibition X 0 (i) in stationary regime is concentrated in the interval (0, 4) when F i = E(1) whereas this distribution is rather concentrated on the interval (20000, 40000) when F i = E [START_REF] Meyn | Stability of Markovian processes III: Foster-Lyapunov criteria for continuous-time processes[END_REF]. This shows that these two distributions of state X 0 (i) are different. In this example, the distribution of the state of inhibition X 0 (i) in stationary regime seems to be continuous although F i is discrete. We do not provide a proof here, this is outside the scope of this paper. We observe two local extrema at 1 and 2 which are linked to the jumps because of the Dirac. These extrema suggest that jumps are very frequent in this process.

1 2 ,

 2 0) and move clockwise around the curve. The two figures below show examples of possible paths for n = 3 and n = 4. Figure.1 shows a possible path with n = 3 and in this case we have Γ : ulu, ulu, urd, drd, drd, dlu.

For n = 4 ,

 4 Figure.2 gives Γ : ulu, ulu, urd, drd, dru, urd, dld, dlu. Writing N (dld), N (drd), • • • for the number of appearances of the different direction vectors, we have that N (dlu) = 1 (dlu is the last triplet of which appears exactly one single time) and N (dru) = N (urd) -1 ≤ n/2, N (drd) + N (dru) + N (uru) + N (urd) = n.

( 6 )

 6 We set m = ∞. While m > 0 do -Let m be the rank of the last possible jump time T m of D in the chronology of jump times. Let k be the neuron associated with this jump.[START_REF] Galves | Perfect simulation of infinite range Gibbs measures and coupling with their finite range approximations[END_REF] We determine the rank r of the last certain jump time T r > T m of k in the chronology of jump times. The state of k is determined recursively from its state at time T r to its state at time T m as follows:-For s ∈ {1, • • • , r -m -1} let x = state of k at time T r-s+1 .

  ) For all time t > 0 we let C i t the clan of ancestors of neuron i at time t (4) While |C i t | > 0 (where |C i t | denotes the cardinality of C i t ) do -Determine the next jump time t next > t in the clan of ancestors of neuron i at time t next and in ∂ ext (clan), the correspondant neuron j and the nature of jump -If neuron j ∈ C i t and the jump is sure, then C

The jump occurs with a Bernoulli distribution with parameter (β(x k -dt (x)) -β * )/(β * -β * )

Acknowledgments:

The author thanks Eva Löcherbach for the many discussions that led to this version of the paper. This research was conducted within the part of the Labex MME-DII(ANR11-LBX-0023-01) project and the CY Initiative of Excellence (grant "Investissements d'Avenir" ANR-16-IDEX-0008), Project EcoDep PSI-AAP 202-00000013