BENO Jeongho Kim

AND Î T Perthame

Delphine Salort

FAST VOLTAGE DYNAMICS OF VOLTAGE-CONDUCTANCE MODELS FOR NEURAL NETWORKS

Keywords: 2010 Mathematics Subject Classification. 35Q92, 35Q84, 92B20 Voltage-conductance model, Integrate-and-Fire, Asymptotic behavior, Neuron assemblies

We present the conductance limit of the voltage-conductance model with random firing voltage when conductance dynamics are slower than the voltage dynamics. The result of the limiting procedure is a transport/Fokker-Planck equation for conductance variable with a non-linear drift which depends on the total firing rate. We analyze the asymptotic behavior of the limit equation under two possible rescalings which relate the voltage scale, the conductance scale and the firing rate. We provide the sufficient framework in which the limiting procedure can be rigorously justified. Moreover, we also suggest a sufficient condition on the parameters and firing distribution in the limiting conductance equation under which we are able to obtain a unique stationary state and its asymptotic stability. Finally, we provide several numerical illustrations supporting the analytic results.

Introduction

Dynamics of neurons is characterized by the membrane potential and the synaptic or ionic channels conductance. When considering an assembly of neurons, a mean-field equation describes the neural network, called the voltage-conductance model which has been introduced in the neuroscience literature, see [START_REF] Cai | Kinetic theory for neuronal network dynamics[END_REF][START_REF] Cai | An effective kinetic representation of fluctuation-driven neuronal networks with application to simple and complex cells in visual cortex[END_REF][START_REF] Rangan | Kinetic theory for neuronal networks with fast and slow excitatory conductances driven by the same spike train[END_REF], and studied mathematically in [START_REF] Cáceres | A numerical solver for a nonlinear fokker-planck equation representation of neuronal network dynamics[END_REF][START_REF] Perthame | On a voltage-conductance kinetic system for integrate & fire neural networks[END_REF][START_REF]Derivation of an integrate and fire equation for neural networks from a voltage-conductance kinetic model[END_REF]. This model enters a larger class of equations that are used to describe neural nets as the time-elapse model [START_REF] Cañizo | Asymptotic behaviour of neuron population models structured by elapsedtime[END_REF][START_REF] Mischler | Relaxation in time elapsed neuron network models in the weak connectivity regime[END_REF][START_REF] Chevallier | Microscopic approach of a time elapsed neural model[END_REF][START_REF] Pakdaman | Adaptation and fatigue model for neuron networks and large time asymptotics in a nonlinear fragmentation equation[END_REF] and the Integrate and Fire model [START_REF] Brette | Adaptive exponential integrate-and-fire model as an effective description of neural activity[END_REF][START_REF] Cáceres | Analysis of nonlinear noisy integrate & fire neuron models: blow-up and steady states[END_REF][START_REF] Carrillo | Qualitative properties of solutions for the noisy integrate and fire model in computational neuroscience[END_REF][START_REF] Fourcaud | Dynamics of the firing probability of noisy integrate-and-fire neurons[END_REF] (and the references therein). Mathematically, these models are focused on a single variable, either time from the last firing or voltage itself, assuming equilibrium for the others.

The voltage-conductance model describes the dynamics of a probability density function p(t, v, g), the probability of finding neuron with membrane potential v and conductance g at time t. Here we consider that the domain of voltage variable is v ∈ (V R , V E ), where V R and V E denote the rest potential and excitatory reversal potential. The conductance variable g varies over nonnegative real number, i.e., g ∈ [0, ∞). The prototype of voltage-conductance model reads as

∂ t p + ∂ v (g L (V L -v) + g(V E -v))p + ∂ g g in (t) -g σ E p - a(t) σ E ∂ 2 gg p = 0.
Here, g L > 0 and V L , with V R < V L < V E , stand for leak conductance and leak potential respectively. For the Integrate-and-Fire model, a firing potential V F is introduced which determines the boundary conditions in v. Due to this difficulty, the original voltageconductance model was simplified in [START_REF]Derivation of an integrate and fire equation for neural networks from a voltage-conductance kinetic model[END_REF] assuming a distribution of V F which leads to the equation

∂ t p + ∂ v [(g L (V L -v) + g(V E -v))p] + ∂ g [(G eq (v, bN (t)) -g)p]
-a∂ 2 gg p + φ F (v)p = 0, t ≥ 0, V R < v < V E , g ≥ 0,

(1.1)

with the no-flux boundary conditions at g = 0 and v = V E , and entering flux at

v = V R (because V R < V L < V E ),
   (G eq (v, bN ) -g)p(t, v, g) -a∂ g p(t, v, g) = 0, for g = 0,

[g L (V L -V R ) + g(V E -V R )] p(t, V R , g) = N (t, g), p(t, V E , g) = 0.
(1.

2)

The individual firing distribution of neurons with potential v is given by φ F (v), which is generally a non-negative increasing function. The function N (t, g) is the network firing rate of neurons with conductance g and N (t) denotes the total firing rate of neurons, which are defined as

(1.3) N (t, g) := V E V R φ F (v)p dv, N (t) := ∞ 0 N (t, g) dg.
Moreover, G eq (v, •) ≥ 0 is the conductance equilibrium at voltage v when there is no noise in conductance variable. Finally, the positive constants a and b denote synaptic noise and synaptic strength of network coupling respectively. The flux in v-direction vanishes at a single point which appears often in the present analysis, therefore we introduce the notation

(1.4) 0 ≤ V L < V * (g) := g L V L + gV E g + g L < V E .
We complete the equation with an initial data p 0 satisfying p 0 (v, g) ≥ 0,

V E V R ∞ 0 p 0 (v, g) dv dg = 1, (1.5) ∞ 0 V E V R g k p 0 (v, g) dv dg < ∞, ∀k ≥ 0.
Note that the no-flux boundary conditions, together with sufficiently fast decay of p at g = ∞ implies the conservation of total number of neurons:

V E V R ∞ 0 p(t, v, g) dv dg = 1.
In the present paper, we are interested in limiting procedure when voltage dynamics is fast compared to the conductance dynamics. More precisely, we first consider the scaled equation of (1.1) given by (1.6)

           ∂ t p ε + 1 ε ∂ v g L (V L -v) + g(V E -v) p ε + ∂ g G eq (v, bN ε (t)) -g p ε -a∂ 2 gg p ε + φ F (v)p ε = 0, t ≥ 0, V R < v < V E , g ≥ 0, 1 ε g L (V L -V R ) + g(V E -V R ) p ε (t, V R , g) = N ε (t, g) := V E V R φ F (v)p ε dv.
Among other possible scales, we also mention the rescaled version of (1.1) with fast firing regime (1.7)

             ∂ t p ε + 1 ε ∂ v g L (V L -v) + g(V E -v) p ε + ∂ g G eq (v, bN ε (t)) -g p ε -a∂ 2 gg p ε + 1 ε φ F (v)p ε = 0, g L (V L -V R ) + g(V E -V R ) p ε (t, V R , g) = N ε (t, g) := V E V R φ F (v)p ε dv.
From now on, to simplify our settings, we only consider the case when the conductance equilibrium G = G eq is bounded above, depends only on the total firing rate N but not on the voltage variable v. Under this assumption, the heuristic limiting equation is obtained by integration in v and should be written as

(1.8) ∂ t n + ∂ g G(bN (t)) -g n -a∂ 2 gg n = 0, t ≥ 0, g ≥ 0,
together with boundary condition at g = 0. We will specify the boundary conditions for the case a = 0 and a > 0 respectively in the next sections as well as the definition of N (t).

Notice that, in [START_REF]Derivation of an integrate and fire equation for neural networks from a voltage-conductance kinetic model[END_REF], the authors considered the fast conductance limit, i.e., following scaled version of equation (1.1) in which fast relaxation of the conductance variable g was considered:

∂ t p ε + ∂ v g L (V L -v) + g(V E -v) p ε + 1 ε ∂ g G eq (v, bN ε (t)) -g p ε , - a ε ∂ 2 gg p ε + φ F (v)p ε =
0, and they derived the following Integrate-and-Fire type equation by letting the scaling parameter ε → 0:

           ∂ t n + ∂ v G(v, bN (t)) V (bN (t)) -v n + φ F (v)n = 0, t ≥ 0, V R < v < V E , N (t) := V E V R φ F (v)n(t, v) dv, G(V R , bN (t)) V(bN (t)) -V R n(t, V R ) = N (t), n(t, V E ) = 0.
(1.9)

A major issue in neural networks is the appearance of synchronisation. This was observed in the voltage-conductance model by [START_REF] Cáceres | A numerical solver for a nonlinear fokker-planck equation representation of neuronal network dynamics[END_REF], in the Integrate-and-Fire model by [START_REF] Brunel | Dynamics of sparsely connected networks of excitatory and inhibitory spiking networks[END_REF][START_REF] Brunel | Fast global oscillations in networks of integrate-and-fire neurons with long firing rates[END_REF][START_REF] Fourcaud | Dynamics of the firing probability of noisy integrate-and-fire neurons[END_REF][START_REF] Sirovich | Dynamics of neural populations: Stability and synchrony[END_REF] and is mathematically related to blow-up in such equation [START_REF] Cáceres | Analysis of nonlinear noisy integrate & fire neuron models: blow-up and steady states[END_REF][START_REF] Dumont | Synchronization of an excitatory integrate-and-fire neural network[END_REF]. Synchronisation also occurs for the equation (1.9). A major issue is therefore to know if this phenomena also persists for one of the versions of equation (1.8) derived from the two proposed rescalings.

Thus, the goal of this paper is twofold. Firstly, we rigorously show the limiting procedure from (1.6) to (1.8); we derive a uniform-in-ε estimate for the moments of probability distribution p ε as well as the firing rate N ε , from which the weak convergence of p ε can be obtained. Secondly, we study the asymptotic behaviors of limiting equation (1.8), considering the cases a = 0 and a > 0 separately. The asymptotic distributions are the Dirac mass and Gaussian distribution for each case respectively.

The paper is organized as follows. In Section 2, we show the slow-fast limit from the voltage-conductance equation (1.6) to the conductance equation. Section 3 contains several properties and asymptotic behavior of the limiting conductance equation (1.8). In Section 4, we consider the other scaling (1.7) where the firing rate φ F has also a fast dynamics. In Section 5, we provide several numerical tests for the models, supporting our analytical results in the previous sections. Finally, Section 6 is devoted to concluding remarks and discussion.

From voltage conductance model to conductance model

To begin with, from equation (1.6), we derive the limiting equation (1.8) at the formal level, and then we provide the rigorous convergence result. We start with the case without a conductance noise, i.e., a = 0 and apply a similar methodology to the case of a > 0.

2.1. Hyperbolic equation. When the synaptic noise a is neglected, the scaled equation (1.6) becomes

(2.1) ∂ t p ε + 1 ε ∂ v g L (V L -v) + g(V E -v) p ε + ∂ g G(bN ε (t)) -g p ε + φ F (v)p ε = 0
and the corresponding no-flux boundary condition becomes

(2.2)      1 ε g L (V L -V R ) + g(V E -V R ) p ε (t, V R , g) = N ε (t, g) := V E V R φ F (v)p ε dv, p ε (t, v, 0) = 0.
To derive the conductance-only equation, we consider a v-marginal of p ε defined as

n ε (t, g) := V E V R p ε (t, v, g) dv.
Next, we integrate (2.1) with respect to v-variable and utilizing the boundary condition (2.2) to obtain

∂ t n ε + ∂ g G(bN ε (t)) -g n ε = 0, where the total firing rate N ε is still N ε (t) := ∞ 0 V E V R φ F (v)p ε (t, v, g) dv dg.
In order to close the equation, we need to express the total firing rate N ε in terms of n ε . However, due to the scaling, the formal limit density p(t, v, g) should satisfy

∂ v g L (V L -v) + g(V E -v) p = 0.
This implies that there exists some C = C(t, g) such that

g L (V L -v) + g(V E -v) p = C(t, g).
However, since p ε (t, V E , g) = 0 for all ε > 0, we have C(g) = 0 and this implies that the limit density p should be concentrated on v = V * (g) (using the notation (1.4)), i.e., (2.3) p(t, v, g) = n(t, g)δ V * (g) (v).

Therefore, at the formal level, the limiting conductance equation should be

∂ t n + ∂ g G(bN (t)) -g n = 0, t ≥ 0, g ≥ 0,
subjected to the zero-flux boundary condition n(t, 0) = 0 and

N (t) = ∞ 0 V E V R φ F (v)n(t, g)δ V * (g) (v) dv dg = ∞ 0 φ F (V * (g)) n(t, g) dg.
2.2. Diffusive equation. Similar to the hyperbolic equation, the limit probability density p(t, v, g) should be concentrated according to (2.3) and hence, the reduced equation becomes,

∂ t n + ∂ g G(bN (t)) -g n -a∂ 2 gg n = 0, N (t) = ∞ 0 φ F (V * (g)) n(t, g) dg.
However, in the diffusive case, the zero-flux boundary condition is given as

G(bN (t))n(t, 0) = a∂ g n(t, 0).
Of course, this boundary condition can include the case when a = 0. However, the analysis, and more precisely the proof of compactness of N ε (t), is much harder when a > 0. Therefore, we separate the two cases when a = 0 and a > 0. In the next section, we prove the slow-fast limit for the hyperbolic case first. The diffusive case is discussed later.

2.3.

Uniform-in-ε estimates on moments and firing rate, a = 0. In order to rigorously derive the limit equation, we follow the strategy in [START_REF]Derivation of an integrate and fire equation for neural networks from a voltage-conductance kinetic model[END_REF]. The first argument we need is the uniform-in-ε boundedness of moments of p so as to ensure that there is no mass loss at infinity in g. with finite g-moments, i.e., (1.5) holds uniformly in ε. Then, for all k ≥ 0, there exists C(k) independent of ε satisfying

M (k) g,ε (t) := ∞ 0 V E V R g k p ε (t, v, g) dv dg ≤ C(k), ∀t ≥ 0.
Proof. Since the case when k = 0 is nothing but the conservation of total mass, we only consider the case when k ≥ 1. We multiply (2.1) by g k , for k ≥ 1, and use integration by parts to find, thanks to the boundary condition at

V R in (2.2), dM (k) g,ε dt - ∞ 0 g k N (t, g)dg -k ∞ 0 V E V R g k-1 (G(bN ε (t)) -g)p ε dv dg + ∞ 0 V E V R g k φ F (v)p ε (t, v, g) dv dg = 0.
Since G is bounded, we have

dM (k) g,ε dt + kM (k) g,ε ≤ CkM (k-1)
g,ε .

However, for any R > 0,

M (k-1) g,ε (t) = V F V R R 0 g k-1 p ε (t, v, g) dv dg + V F V R ∞ R g k-1 p ε (t, v, g) dv dg ≤ R k-1 + 1 R V F V R ∞ 0 g k p ε (t, v, g) dv dg = R k-1 + 1 R M (k) g,ε (t).
Therefore, dM

(k) g,ε dt + kM (k) g,ε ≤ CM (k-1) g,ε ≤ CkR k-1 + Ck R M (k) g,ε .
By taking R sufficiently large, we have

dM (k) g,ε dt + M (k) g,ε ≤ CR k-1 , i.e., M (k) g,ε (t) ≤ M (k) g,ε (0)e -t + CR k-1
1 -e -t which implies the uniform boundedness of any k-th g-moment.

Next, we study the uniform-in-ε boundedness of total firing rate N ε (t). Lemma 2.2 (Compactness of N ε (t), a = 0). Let {p ε } ε>0 be a family of solution to (2.1) with initial data p 0 ε ε>0 with finite g-moments, i.e., (1.5) holds uniformly in ε for 0 < ε < 1. Then, we have the uniform bound

N ε L ∞ (0,∞) ≤ φ F L ∞
Moreover, after initial time layer τ of order √ ε, we also have the Lipschitz bound

dN ε dt L ∞ (τ,∞) ≤ C(τ ).
Proof. The uniform boundedness of N ε follows immediately from its definition (1.3). Next, we estimate the derivative of N ε as

dN ε dt = ∞ 0 V E V R φ F (v) - 1 ε ∂ v g L (V L -v) + g(V E -v) p ε -∂ g G(bN ε (t)) -g p ε -φ F (v)p ε dv dg = ∞ 0 1 ε φ F (V R ) g L (V L -V R ) + g(V E -V R ) p ε (t, V R , g) dg + 1 ε ∞ 0 V E V R φ F (v) g L (V L -v) + g(V E -v) p ε dv dg - ∞ 0 V E V R φ F (v) 2 p ε dv dg.
Using the boundary condition, this yields

dN ε dt = ∞ 0 V E V R φ F (V R ) -φ F (v) φ F (v) p ε dv dg + 1 ε ∞ 0 V E V R φ F (v) g L (V L -v) + g(V E -v) p ε dv dg.
(2.4)

In order to show that N ε (t) is uniformly Lipschitz continuous in ε, we need to control the second term. To do so, we define

Q ε (t) := ∞ 0 V E V R g L (V L -v) + g(V E -v) p ε dv dg. To estimate Q ε , we multiply (2.1) by |g L (V L -v) + g(V E -v)| and use integration by parts to obtain d dt Q ε (t) ≤ 1 ε ∞ 0 g L (V L -V R ) + g(V E -V R ) g L (V L -V R ) + g(V E -V R ) p ε (V R , g) dg - 1 ε ∞ 0 V E V R (g L + g) g L (V L -v) + g(V E -v) p ε dv dg + ∞ 0 V E V R ∂ g g L (V L -v) + g(V E -v) G(bN ε (t)) -g p ε dv dg - ∞ 0 V E V R g L (V L -v) + g(V E -v) φ F (v)p ε dv dg.
Here, we have used the relation

∂ v |g L (V L -v) + g(V E -v)| g L (V L -v) + g(V E -v) = -(g + g L ) g L (V L -v) + g(V E -v) .
Using the boundary condition, this implies

d dt Q ε (t) + 1 ε ∞ 0 V E V R (g L + g) g L (V L -v) + g(V E -v) p ε dv dg ≤ ∞ 0 g L (V L -V R ) + g(V E -V R ) -g L (V L -v) + g(V E -v) N (t, g) dg + ∞ 0 V E V R (V E -v) G(bN ε (t)) -g p ε dv dg.
Since the first g-moment of p ε is uniformly bounded by Lemma 2.1, we have

(2.5) d dt Q ε (t) + g L ε Q ε (t) ≤ C,
which, together with the Grönwall inequality, implies

Q ε (t) ≤ Q ε (0)e -g L ε t + Cε g L .
Therefore, for τ = O( √ ε), we have

Q ε (τ ) ≤ Q ε (0)e -C √ ε + Cε g L ≤ Ce -C √ ε + Cε g L ≤ Cε
for sufficiently small ε. Therefore, we conclude that for ε 1, using again (2.5),

Q ε (t) ≤ Q ε (τ )e -g L ε (t-τ ) + Cε g L ≤ Cε, for t ≥ τ.
Thus, we can control the quantity Q ε after initial layer, and therefore, we can substitute this estimate to (2.4) to obtain

dN ε dt = ∞ 0 V E V R φ F (V R ) -φ F (v))(φ F (v) p ε dv dg + 1 ε ∞ 0 V E V R φ F (v) g L (V L -v) + g(V E -v) p ε dv dg ≤ 2 φ F 2 L ∞ + φ F L ∞ ε Q ε (t) ≤ C, for t ≥ τ.
Remark 2.3. One can get rid of the initial layer region if we consider a well-prepared initial data. More precisely, if we consider the initial data p 0 ε with Q ε (0) < Cε, then, it automatically follows from (2.5) that Q ε (t) ≤ Cε, for t ≥ 0.

2.4. Uniform-in-ε estimates on moments and firing rate, a > 0. Next, we prove the compactness of the total firing rate when there is a synaptic noise. The major difficulty lies in estimating the Lipschitz norm of N ε and this requires more step than in the case a = 0.

Again, we start with an estimate for moments.

Lemma 2.4. Let {p ε } ε>0 be a family of solution to (1.6) with initial data {p 0 ε } ε>0 with finite g-moments, i.e., (1.5) holds uniformly in ε. Then, for all k ≥ 0, there exists C(k) independent of ε satisfying

M (k) g,ε (t) := ∞ 0 V E V R g k p ε (t, v, g) dv dg ≤ C(k), ∀t ≥ 0.
Proof. Again, we only consider the case when k ≥ 1. After exactly same procedure with Lemma 2.1 with k ≥ 2, we obtain

dM (k) g,ε dt + kM (k) g,ε ≤ C M (k-1) g,ε + M (k-2) g,ε
.

However, for any R > 0, we have,

M (k-1) g,ε ≤ R k-1 + 1 R M (k) g,ε , M (k-2) g,ε ≤ R k-2 + 1 R 2 M (k) g,ε ,
and by taking R sufficiently large, it holds

dM (k) g,ε dt + M (k) g,ε ≤ C(R k-1 + R k-2 ), k ≥ 2, which implies a uniform-in-ε bound of M (k)
g,ε for k ≥ 2. The case of k = 1 immediately follows by using the Cauchy-Schwarz inequality, together with mass conservation.

As before, the last and main step is to estimate the total firing rate. Lemma 2.5 (Compactness of N ε (t), a > 0). Let {p ε } ε>0 be a family of solution to (1.6) with initial data {p 0 ε } ε>0 satisfying (1.5) uniformly in ε for 0 < ε < 1 and

(2.6) φ F (V L ) = φ F (V L ) = 0.
Then, we have the estimate

N ε L ∞ (0,∞) ≤ φ F L ∞ .
Moreover, after initial time layer τ of order √ ε, we also have the Lipschitz bound

dN ε dt L ∞ (τ,∞) ≤ C(τ, T ).
Proof. Since the estimate of N ε is identical to previous case, we focus on estimate of dNε dt . Considering boundary condition, the equality (2.4) stills holds, and we need to estimate the term

R ε (t) = ∞ 0 V E V R φ F (v) g L (V L -v) + g(V E -v) p ε dv dg = O(ε)
The previous strategy based on the quantity Q ε (t) is not compatible with the diffusion in g. Therefore we introduce another approach and set

R 1 ε (t) = ∞ 0 V E V R φ F (v) -φ F (V * (g)) g L (V L -v) + g(V E -v) p ε dv dg, R 2 ε (t) = ∞ 0 V E V R φ F (V * (g)) g L (V L -v) + g(V E -v) p ε dv dg.
We are going to prove that both terms are of order ε for t ≥ τ .

Step 1. The term R 1 ε (t). This term can be treated as Q ε (t) in the case a = 0 because, here, the multiplier is smooth. To begin with, notice that, because of linearity

g L (V L -v) + g(V E -v) = -(g L + g)(v -V * (g)),
we can write using that φ F is Lipschitz

|R 1 ε (t)| ≤ CR 1 ε (t), R 1 ε (t) := ∞ 0 V E V R (g L + g) v -V * (g) 2 p ε dv dg, It remains to prove that R 1 ε (t) is of order ε. We estimate R 1 ε by using equation (1.6
) and integration by parts as

dR 1 ε (t) dt = - 2 ε ∞ 0 V E V R (g L + g) 2 v -V * (g) 2 p ε dv dg + ∞ 0 V E V R ∂ g (g L + g)(v -V * (g)) (G -g)p ε -a∂ g p ε dv dg + O(1)
where the terms O(1) are computed exactly as in the case a = 0. Only the control of the term with the derivative in g is new. We write, after one more integration by parts ,

dR 1 ε (t) dt = - 2 ε ∞ 0 V E V R (g L + g) 2 v -V * (g) 2 p ε dv dg + a V E V R ∂ g (g L + g)(v -V * (g)) p ε dv g=0 + O(1)
Using the formula (1.4), it remains to notice that

∂ g (g L + g)(v -V * (g)) = v -V E ≤ 0.
This allows to conclude, as in Lemma 2.2, that |R 1 ε (t)| ≤ Cε for t ≥ τ and conclude the estimate on the term R 1 ε .

Step 2. The term R 2 ε (t). For this term, we define

S(t, g) := V E V R g L (V L -v) + g(V E -v) p ε dv, R 2 ε (t) = ∞ 0 φ F (V * (g))S(t, g) dg.
By using (1.6), one can compute the equations for S(t, g) as

∂ t S(t, g) + g L + g ε S = (g + g L ) V E V R (v -V R )φ F (v)p ε dv - V E V R g L (V L -v) + g(V E -v) ∂ g G(bN ε (t)) -g p ε -a∂ g p ε dv.
We multiply e g L +g ε t on the both sides and integrate in time to derive

S(t,g) = S(0, g)e -g L +g ε t + t 0 e g L +g ε (s-t) (g + g L ) V E V R (v -V R )φ F (v)p ε (s, v, g)dvds - t 0 e g L +g ε (s-t) V E V R g L (V L -v) + g(V E -v) ∂ g G(bN ε (t)) -g p ε -a∂ g p ε dv ds
Therefore, we can recover that

R 2 ε (t) = ∞ 0 φ F (V * (g))S(0, g)e -g L +g ε t dg + I + II,
where

I = ∞ 0 φ F (V * (g)) t 0 e g L +g ε (s-t) (g + g L ) V E V R (v -V R )φ F (v)p ε (s, v, g)dv ds dg
which, using the moment bounds in Lemma 2.4, is estimated as

I ≤ φ F ∞ φ F ∞ (V E -V R ) ∞ 0 t 0 e -g L ε (t-s) (g + g L ) V E V R p ε (s, v, g) dv ds dg ≤ φ F ∞ φ F ∞ (V E -V R ) t 0 e -g L ε (t-s) ds sup 0≤s≤t ∞ 0 V E V R (g + g L )p ε (s, v, g) dv dg ≤ Cε.
And the other term is written, thanks to the zero flux boundary condition at g = 0,

II = - ∞ 0 φ F (V * (g)) t 0 e g L +g ε (s-t) × V E V R g L (V L -v) + g(V E -v) ∂ g G(bN ε (t)) -g p ε -a∂ g p ε dv ds dg = ∞ 0 V E V R t 0 ∂ g φ F (V * (g))e g L +g ε (s-t) g L (V L -v) + g(V E -v) × G(bN ε (t)) -g p ε -a∂ g p ε dv ds dg = O(ε) -a ∞ 0 V E V R t 0 ∂ g φ F (V * (g))e g L +g ε (s-t) g L (V L -v) + g(V E -v) ∂ g p ε dv ds dg = O(ε) + a ∞ 0 V E V R t 0 ∂ 2 gg φ F (V * (g))e g L +g ε (s-t) g L (V L -v) + g(V E -v) p ε dv ds dg
and the boundary term at g = 0 vanishes thanks to the assumption (2.6). Using again Lemma 2.4, this provides a control

R 2 ε (t) ≤ Ce -g L t ε + Cε.
The remaining procedure is identical to that of Lemma 2.2.

2.5. Weak convergence. We have gathered the material to present our first main theorem, which holds both for the hyperbolic and diffusive regimes (noiseless or noisy models)

Theorem 2.6. Let {p ε } ε>0 be a family of solution to (2.1) with initial data {p 0 ε } ε>0 with finite g-moments, i.e., (1.5) holds uniformly in ε. Then, as ε → 0,

p ε (t, v, g) n(t, g)δ V * (g) (v), t ≥ τ, τ = O( √ ε)
(weakly in the sense of probability measures), up to a subsequence. the probability density n(t, g) satisfies

(2.7) ∂ t n + ∂ g ((G(bN (t)) -g)n) -a∂ 2 gg n = 0, t ≥ 0, g ≥ 0, (2.8) n 0 (g) := V E V R p 0 (v, g) dv, N (t) := ∞ 0 φ F V * (g) n(t, g) dg,
with the no-flux boundary condition

G(bN (t))n(t, 0) = a∂ g n(t, 0).
Proof. Since N L ∞ (0,∞) and dN dt L ∞ (τ,∞) are uniformly bounded, by the Arzela-Ascoli theorem, for any T > 0, there exists a Lipschitz function

N : (τ, T ) → [0, ∞) such that lim ε→0 N ε -N L ∞ (τ,T ) = 0.
Moreover, since {p ε } ε≥0 is a family of probability measures which moments are uniformly bounded, there exists a probability measure p such that, in the weak sense of bounded measures, and up to extraction of a subsequence p ε p, gp ε gp.

Then, after applying the arguments used in deriving the formal limit equation in Section 2.1 and Section 2.2, we conclude that the limit probability measure p is of the form

p(t, v, g) = n(t, g)δ V * (g) (v), t ≥ τ,
with n satisfying (2.7)-(2.8).

Asymptotic stability of conductance equation

Since the conductance equation (2.7) is established, we now look for the behavior of solutions. We first establish the existence of a stationary solution, and then we study its asymptotic stability. Again, we first study the hyperbolic equation.

∂ t n + ∂ g G(bN (t)) -g n = 0, N (t) := ∞ 0 φ F (V * (g)) n(t, g) dg
with no-flux boundary condition n(t, 0) = 0. From now on, we consider the linear case when G(bN (t)) = bN (t) for analytic simplicity. Then, a stationary solution n ∞ should satisfy

∂ g bN -g n ∞ = 0, N = ∞ 0 φ F (V * (g)) n ∞ (g) dg.
Therefore, considering the boundary condition, the stationary solution n ∞ has to be n ∞ (g) = δ bN (g).

We again substitute this form of stationary solution to the equation for N which yields

N = φ F (V * (bN )) = φ F g L V L + bN V E bN + g L .
Therefore, finding a stationary solution is equivalent to finding a fixed point of the map

Φ : N → φ F g L V L + bN V E bN + g L =: φF (bN ).
Here, we defined the auxiliary function φ for notational convenience.

Proposition 3.1. Suppose that b is sufficiently small so that the following relation is satisfied:

(3.2) b < g L (V E -V L ) φ F L ∞ .
Then, the map Φ has the unique fixed point N ∞ . Thus, there is a unique stationary state of (3.1).

Proof. First of all, let us briefly study the property of φF . Since the function φ F is increasing and the map (using the notation (1.4))

g → V * (g) = V E - g L (V E -V L ) g + g L
is increasing with respect to g, so is Φ, the composition of two increasing maps when b > 0.

Since Φ(0) = φF (0) = φ F (V L ) > 0, Φ has a unique fixed point N ∞ if Φ < 1. However, this condition is equivalent to

dΦ dN = d dN φF (bN ) = b φ F (bN ) = bφ F g L V L + bN V E bN + g L g L (V E -V L ) (bN + g L ) 2 < 1.
Indeed, assuming the condition (3.2) on b, we have

dΦ dN = bφ F g L V L + bN V E bN + g L g L (V E -V L ) (bN + g L ) 2 < g L bN + g L 2 < 1
and we conclude that Φ has a unique fixed point if (3.2) holds.

From now on, we assume (3.2) so that we have the unique fixed point N ∞ . Our next goal is to study the long term convergence of solution n(t, g) toward the stationary solution n ∞ (g) = δ bN∞ (g). In order to show the convergence, we consider the Lyapunov functional

L(t) := ∞ 0 (g -bN ∞ ) 2 n dg.
Note that this Lyapunov functional can be understood as a measure of how far the solution n is from the stationary solution δ bN∞ (g). The following theorem shows that the condition (3.2) which guarantees the existence and uniqueness of fixed point N ∞ is also sufficient condition to exponential decay of L. Theorem 3.2. Let n = n(t, g) satisfy (3.1) and assume the smallness condition (3.2) on b. Let N ∞ be the unique fixed point of Φ obtained in Proposition 3.1. Then, the Lyapunov functional L decays exponentially towards 0, namely

L(t) ≤ L(0)e -2(1-b φ F L ∞ )t .
Let us mention that the method based on Monge-Kantorowich distance (see [START_REF] Fournier | Monge-Kantorovich distance for PDEs: the coupling method[END_REF]) allows to prove that solutions will converge to a steady state, even when it is not unique. However, the norm is weaker and should be adapted which is beyond the scope of the present analysis.

Proof. We directly take a derivative to the definition of L to obtain

dL dt = ∞ 0 (g -bN ∞ ) 2 -∂ g (bN (t) -g)n dg = 2 ∞ 0 (g -bN ∞ )(bN (t) -g)n dg = -2 ∞ 0 (bN ∞ -g) 2 n dg + 2b(N (t) -N ∞ ) ∞ 0 (g -bN ∞ )n dg ≤ -2L + 2b|N (t) -N ∞ | √ L (3.3)
where we have used the Cauchy-Schwarz inequality together with a mass conservation law

∞ 0 |g -bN ∞ |n dg ≤ ∞ 0 n dg 1/2 ∞ 0 |g -bN ∞ | 2 n dg 1/2 = √ L.
On the other hand, since N ∞ is the unique fixed point of Φ, it satisfies N ∞ = φF (bN ∞ ). Therefore, we have

|N (t) -N ∞ | = ∞ 0 φF (g) -φF (bN ∞ ) n dg ≤ ∞ 0 φF (g) -φF (bN ∞ ) n dg ≤ ∞ 0 φ F L ∞ |g -bN ∞ |n dg ≤ φ F L ∞ √ L.
(3.4)

Now, we substitute (3.4) into (3.3) to obtain

dL dt ≤ -2(1 -b φ F L ∞ )L.
Therefore, if we have b φ F L ∞ < 1, then we have a desired exponential decay of L. Note that this condition is exactly same condition with (3.2).

Diffusive equation.

When considering the noisy conductance equation, the limit equation is given by

(3.5) ∂ t n + ∂ g bN (t) -g n -a∂ 2
gg n = 0 with the zero-flux boundary condition bN (t)n(t, 0) = a∂ g n(t, 0).

As in the hyperbolic case, we first search for a stationary solution n ∞ which should satisfy

∂ g bN ∞ -g n ∞ = a∂ 2
gg n ∞ . Thanks to the boundary condition, we have

(bN ∞ -g)n ∞ = a∂ g n ∞ or equivalently n ∞ (g) = 1 Z N∞ exp - (g -bN ∞ ) 2 2a
where N ∞ is given by the implicit equation and normalizing constant

N ∞ = ∞ 0 φF (g)n ∞ (g) dg, Z N∞ := ∞ 0 exp - (g -bN ∞ ) 2 2a dg.
Based on our experience in hyperbolic case, finding a stationary state is equivalent to finding a fixed point of the following map:

Ψ : N → ∞ 0 φF (g) 1 Z N exp - (g -bN ) 2 2a dg.
The following proposition provides a sufficient condition so that the map Ψ has a unique fixed point. Proof. First of all, we note that

Ψ(0) = 1 Z 0 ∞ 0 φF (g) exp - g 2 2a dg > 0.
Therefore, if we have Ψ (N ) < 1, then we have the existence and uniqueness of fixed point. Hence we first investigate Ψ (N ). In fact, we have

d dN Z N = b a ∞ 0 (g -bN ) exp - (g -bN ) 2 2a dg = b a ∞ -bN g exp - g 2 2a dg = b a -a exp - g 2 2a ∞ bN = b exp - b 2 N 2 2a . (3.6)
We take a derivative to Ψ to obtain

dΨ dN = - d dN Z 2 N ∞ 0 φF (g) exp - (g -bN ) 2 2a dg + 1 Z N ∞ 0 φF (g) b a (g -bN ) exp - (g -bN ) 2 2a dg.
(3.7)

However, the integration by parts yields

1 a ∞ 0 φF (g)(g -bN ) exp - (g -bN ) 2 2a dg = - ∞ 0 φF (g)∂ g exp - (g -bN ) 2 2a dg = φF (0) exp - b 2 N 2 2a + ∞ 0 φ F (g) exp - (g -bN ) 2 2a dg,
which, combined to (3.7), gives

dΨ dN = - d dN Z N Z 2 N ∞ 0 φF (g) exp - (g -bN ) 2 2a dg + b Z N φF (0) exp - b 2 N 2 2a + b Z N ∞ 0 φ F (g) exp - (g -bN ) 2 2a dg = b Z N ∞ 0 φ F (g) exp - (g -bN ) 2 2a dg - d dN Z N Z 2 N ∞ 0 φF (g) -φF (0) exp - (g -bN ) 2 2a dg = 1 Z N ∞ 0 b φ F (g) - d dN Z N Z N ( φF (g) -φF (0)) exp - (g -bN ) 2 2a dg.
Therefore, if the network connectivity b satisfies

(3.8) b φ F (g) - d dN Z N Z N ( φF (g) -φF (0)) < 1
for any g and N , we have

dΨ dN < 1 Z N ∞ 0 exp - (g -bN ) 2 2a dg = 1,
which implies the existence and uniqueness of fixed point. However, the condition (3.8) is easily satisfied when b φ F L ∞ < 1 since φF (g) ≥ φF (0) and d dN Z N > 0 and therefore,

b φ F (g) - d dN Z N Z N φF (g) -φF (0) ≤ b φ F L ∞ < 1.
Remark 3.4. We investigate the condition (3.8) more precisely as follows.

1. We can estimate Z N also as follows:

Z N := ∞ -bN exp - g 2 2a dg = aπ 2 + bN 0 exp - g 2 2a = aπ 2 + √ 2a bN √ 2a 0 exp(-x 2 ) dx = aπ 2 1 + erf bN √ 2a = aπ 2 2 -erfc bN √ 2a
However, the complementary error function erfc(x) can be estimated as in [START_REF] Chiani | New exponential bounds and approximations for the computation of error probability in fading channels[END_REF], erfc(x) ≤ e -x 2 , which implies

aπ 2 2 -exp - b 2 N 2 2a ≤ Z N ≤ √ 2aπ.
2. The condition (3.8) can be rewritten as

(3.9) b φ F (g) < 1 + d dN Z N Z N φF (g) -φF (0) ,
for any g and N . On the other hand, according to (3.6) and estimate on Z N in (1),

d dN Z N Z N ≤ b exp -b 2 N 2 2a aπ 2 2 -exp -b 2 N 2 2a → 0, as N → ∞.
Since (3.9) should be satisfied as N → ∞, it becomes b φ F (g) < 1 for all g ≥ 0. Therefore, the condition b φ F L ∞ < 1 is indeed equivalent to (3.8). We now study the asymptotic behavior of (3.5). Basically, we use a relative entropy method. To this end, we define h(t, g) := n(t, g) n ∞ (g) .

Lemma 3.5. For any convex function S, we have

d dt ∞ 0 n ∞ S(h(g, t)) dg = -a ∞ 0 n ∞ S (h)(∂ g h) 2 dg+b(N ∞ -N ) ∞ 0 n ∞ ∂ g (S(h)-S (h)h) dg
Proof. Since n = hn ∞ , we have

∂ g n = n ∞ ∂ g h + h∂ g n ∞ , and ∂ 2 gg n = n ∞ ∂ 2 g h + 2∂ g n ∞ ∂ g h + ∂ 2 gg n ∞ h. Therefore, ∂ t h = 1 n ∞ (-∂ g ((bN (t) -g)n) + a∂ 2 gg n) = 1 n ∞ (n + (g -bN (t))∂ g n + a∂ 2 gg n) = n n ∞ + (g -bN ) ∂ g h + ∂ g n ∞ n ∞ h + a∂ 2 gg h + 2a ∂ g n ∞ n ∞ ∂ g h + a ∂ 2 gg n ∞ n ∞ h = g -bN + 2a ∂ g n ∞ n ∞ ∂ g h + a∂ 2 gg h + n n 2 ∞ n ∞ + (g -bN )∂ g n ∞ + a∂ 2 gg n ∞ = g -bN + 2a ∂ g n ∞ n ∞ ∂ g h + a∂ 2 gg h + b(N ∞ -N ) n n 2 ∞ ∂ g n ∞
where we used the equation for n ∞ :

n ∞ + (g -bN ∞ )∂ g n ∞ + a∂ 2 gg n ∞ = -∂ g ((bN ∞ -g)n ∞ ) + a∂ 2 gg n ∞ = 0. Therefore, ∂ t S(h) = S (h)∂ t h = S (h) g -bN + 2a ∂ g n ∞ n ∞ ∂ g h + a∂ 2 gg h + b(N ∞ -N ) n n 2 ∞ ∂ g n ∞ = g -bN + 2a ∂ g n ∞ n ∞ ∂ g S(h) + a ∂ 2 gg S(h) -S (h)(∂ g h) 2 + b(N ∞ -N ) n n 2 ∞ S (h)∂ g n ∞
and, consequently,

∂ t (n ∞ S(h)) = (g -bN )n ∞ ∂ g S(h) -a∂ 2 gg n ∞ S(h) + a∂ 2 gg (n ∞ S(h)) -an ∞ S (h)(∂ g h) 2 + b(N ∞ -N )hS (h)∂ g n ∞ = ∂ g ((g -bN ∞ )n ∞ S(h)) + a∂ 2 gg (n ∞ S(h)) -an ∞ S (h)(∂ g h) 2 + b(N ∞ -N )(S (h)hn ∞ + n ∞ ∂ g S(h)).
Hence, we have

d dt ∞ 0 n ∞ S(h) dg = bN ∞ n ∞ (0)S(h(0)) -a∂ g (n ∞ S(h))| g=0 -a ∞ 0 n ∞ S (h)(∂ g h) 2 dg + b(N ∞ -N ) ∞ 0 S (h)hn ∞ + n ∞ ∂ g S(h) dg = -an ∞ ∂ g S(h)| g=0 -a ∞ 0 n ∞ S (h)(∂ g h) 2 dg + b(N ∞ -N ) ∞ 0 n ∞ (S (h)h -S(h)) dg -b(N ∞ -N )n ∞ (0)S(h(0)) (3.10)
However, since we have

∂ g S(h) = S (h) (∂ g n)n ∞ -n(∂ g n ∞ ) n 2 ∞ and boundary conditions bN (t)n(0) = a∂ g n(0), bN ∞ n ∞ (0) = a∂ g n ∞ (0), we obtain (3.11)
-an ∞ ∂ g S(h)| g=0 = b(N ∞ -N )n(0)S (h(0)).
Therefore, we substitute (3.11) into (3.10) to derive

d dt ∞ 0 n ∞ S(h) dg = -a ∞ 0 n ∞ S (h)(∂ g h) 2 dg + b(N ∞ -N ) n(0)S (h(0)) -n ∞ (0)S(h(0)) + ∞ 0 n ∞ (S (h)h -S(h)) dg = -a ∞ 0 n ∞ S (h)(∂ g h) 2 dg + b(N ∞ -N ) ∞ 0 n ∞ ∂ g (S(h) -S (h)h) dg.
Now, we take a specific form of convex function S(x) = (x -1) 2 and we define a relative entropy Lyapunov functional E as

E(t) := ∞ 0 n ∞ S(h) dg = ∞ 0 n ∞ (h -1) 2 dg,
and its dissipation D as

D(t) := 1 2 ∞ 0 n ∞ S (h)(∂ g h) 2 dg = ∞ 0 n ∞ (∂ g h) 2 dg.
Then the entropy equality in Lemma 3.5 becomes (3.12)

dE dt = -2aD + b(N ∞ -N ) ∞ 0 n ∞ ∂ g (1 -h 2 ) dg.
In order to estimate further, we use the following Poincaré-type inequality [START_REF] Cáceres | Analysis of nonlinear noisy integrate & fire neuron models: blow-up and steady states[END_REF][START_REF] Carrillo | Qualitative properties of solutions for the noisy integrate and fire model in computational neuroscience[END_REF]: there exists a positive constant γ such that for any function h satisfying

∞ 0 n ∞ (g)h(g) dg = 1, γ ∞ 0 n ∞ (g)(h(g) -1) 2 dg ≤ ∞ 0 n ∞ (g)(∂ g h) 2 (g) dg,
or in our notation, γE ≤ D. Below, we provide a sufficient condition on parameters and initial data which guarantee the exponential decay of relative entropy E. Theorem 3.6. Suppose n = n(t, g) satisfies (3.5). Moreover, we assume the following conditions on parameters

b φ F L ∞ < 1, b 1 + 2 φ F 2 L ∞ γ < 2a
and initial data

2bE(0) < 2a -b 1 + 2 φ F 2 L ∞ γ .
Let N ∞ be the unique fixed point of Ψ obtained in Proposition 3.3. Then, the relative entropy functional E decays exponentially towards 0: there exists a positive constant µ > 0 such that

E(t) ≤ E(0)e -µt .
Proof. From equation (3.12), we have

dE dt = -2aD + 2b(N -N ∞ ) ∞ 0 n ∞ ((h -1)∂ g h + ∂ g h) dg = -2aD + 2b(N -N ∞ ) ∞ 0 n ∞ ((h -1)∂ g h) dg + 2b(N -N ∞ ) ∞ 0 n ∞ ∂ g h dg ≤ -2aD + b (N -N ∞ ) 2 + ∞ 0 n ∞ ((h -1)∂ g h) dg 2 + b (N -N ∞ ) 2 + ∞ 0 n ∞ ∂ g h dg 2 ≤ -2aD + b((N -N ∞ ) 2 + ED) + b((N -N ∞ ) 2 + D) = -D(2a -b -bE) + 2b(N -N ∞ ) 2 (3.13)
However, from the definition of N ,

(N -N ∞ ) 2 = ∞ 0 φF (g)(n -n ∞ ) dg 2 ≤ ∞ 0 φF (g)n ∞ |h -1| dg 2 ≤ φ F 2 L ∞ ∞ 0 n ∞ dg ∞ 0 n ∞ (h -1) 2 dg = φ F 2 L ∞ E
Therefore, we substitute the estimate of (N -N ∞ ) 2 in (3.13) to obtain

dE dt ≤ -D(2a -b -bE) + 2b φ F 2 L ∞ E.
Now, by the assumption on initial data, 2a -b -bE(0) > 0 and hence, by the Poincaré-type inequality, we have,

dE dt t=0 ≤ -γE 2a -b -bE(0) - 2b φ F 2 L ∞ γ ≤ 0.
Therefore, E decays in times and hence

2bE(t) < 2a -b 1 + 2 φ F 2 L ∞ γ .
This implies the following exponential decay of E with decay rate µ = γ 2 2a -b -

2b φ F 2 L ∞ γ : dE dt ≤ -γE 2a -b -bE - 2b φ F 2 L ∞ γ ≤ - γ 2 2a -b - 2b φ F 2 L ∞ γ E.

Slow-fast limit with fast firing regime

We now consider the case when both the voltage dynamics and the firing rate have fast dynamics. A scaling of interest, because it leads to a different analysis than in the previous section, is when both terms have the same strength. More precisely, we consider the scaled equation, where the factor 1 ε also appears in front of firing rate φ F , (

∂ t p ε + 1 ε ∂ v g L (V L -v) + g(V E -v) p ε + ∂ g G eq (v, bN ε (t)) -g p ε -a∂ 2 gg p ε + 1 ε φ F (v)p ε = 0, t ≥ 0, v ∈ (V L , V E ), g ≥ 0, 4.1)
with the corresponding no-flux boundary conditions (4.2)

                 g L (V L -V R ) + g(V E -V R ) p ε (t, V R , g) = N ε (t, g) := V E V R φ F (v)p ε dv, p ε (t, V E , g) = 0, G eq (v, bN ε (t)) -g p ε -a∂ g p ε (t, v, 0) = 0, N ε (t) := ∞ 0 N ε (t, g)dg = ∞ 0 V E V R φ F (v)p ε dv dg.
We first identify the limiting density as ε → 0. Formally, the probability density function p = lim ε→0 p ε should satisfy

∂ v g L (V L -v) + g(V E -v) p + φ F (v)p = 0, which is equivalent to (4.3) ∂ v (log p) = g L + g -φ F (v) g L (V L -v) + g(V E -v)
.

When V R ≤ v < V * (g) (see (1.4)), we integrate (4.3) from V R to v to obtain p(t, v, g) = p(t, V R , g) exp v V R g L + g -φ F (w) g L (V L -w) + g(V E -w) dw , V R ≤ v < V * (g). For V * (g) < v ≤ V E , we integrate (4.3) from v to V E to get p(t, v, g) = p(t, V E , g) exp - V E v g L + g -φ F (w) g L (V L -w) + g(V E -w) dw = 0.
In summary, this case is very close to that studied in [START_REF]Derivation of an integrate and fire equation for neural networks from a voltage-conductance kinetic model[END_REF], and the limit density p is more complex, given as

(4.4) p(t, v, g) =    p(t, V R , g) exp v V R g L + g -φ F (w) g L (V L -w) + g(V E -w) dw V R ≤ v < V * (g), 0 V * (g) < v < V E .
For the special case when φ F is constant, i.e., φ F ≡ φ, we can explicitly calculate the limiting distribution p. First of all, the exponent in (4.4) can be calculated as

v V R g L + g -φ g L (V L -w) + g(V E -w) dw = - (g L + g -φ) g L + g log((g L V L + gV E ) -(g L + g)w) v V R = φ -g L -g g L + g log g L V L + gV E -(g L + g)v g L V L + gV E -(g L + g)V R .
Thus, the density p is given by

p(t, v, g) =p(t, V R , g) exp φ -g L -g g L + g log g L V L + gV E -(g L + g)v g L V L + gV E -(g L + g)V R = p(t, V R , g) g L V L + gV E -(g L + g)v g L V L + gV E -(g L + g)V R φ-g L -g g L +g
, which behaves like x φ-g L -g g L +g near x = 0 as v → V * (g). Since φ-g L -g g L +g > -1, the distribution is integrable with respect to v, although it blows-up when g > φ -g L .

For the sake of completeness, we mention the a priori bounds for the g-moments Lemma 4.1. Let {p ε } ε>0 be a family of solution to (4.1), (4.2) with initial data p 0 ε ε>0 with finite g-moments, i.e., (1.5) holds uniformly in ε. Then, for all k ≥ 0, there exists C(k) independent of ε satisfying

M (k) g,ε (t) := ∞ 0 V E V R g k p ε (t, v, g) dv dg ≤ C(k), ∀t ≥ 0.
We omit the proof which is similar to those of Lemma 2.1 for a = 0 or Lemma 2.4 for a > 0.

Even though the convergence theory seems much harder than with the previous scalings, the behavior as ε → 0 is tested numerically in Subsection 5.1 and the formal asymptotic is confirmed numerically.

Numerical simulation

We now illustrate our theoretical results with numerical simulations. First, we test the convergence of solution to (1.6) and (1.7) as ε → 0. In particular, we observe that, for (1.6), the solution converges to the Dirac delta distribution concentrated on the 1-dimensional manifold (v, g) : v = V * (g)} and that, for (1.7), we obtain that the formal limit distribution for fast firing regime has algebraic order near the critical voltage v = V * (g). Secondly, we move on to the numerical simulation for limit equations (3.1) and (3.5), for both a = 0 and a > 0. We provide the dynamics of N as well as the dynamics of n(t, g) and the numerical results indicate that the asymptotic convergence analysis to a stationary state on the limit model is also observable in the numerical simulation. Moreover, with specific choice of φ, we can numerically observe the emergence of a bifurcation behavior with respect to the parameter b.

5.1. Convergence of solution as ε → 0. We first simulate equation (1.6) and (1.7) for various values of ε and observe the behavior of the solutions. We take a set of parameters as follows:

g L = 3, g M ax = 20, V R = 0, V L = 0.1, V E = 1, a = 1, b = 1,
and the conductance equilibrium function G and firing function φ F (v) as

G(v, bN ) = bN , φ F (v) = 1 + v.
The initial data is chosen as a Gaussian distribution given as

p 0 (v, g) := 1 Z exp -30(v -0.5) 2 -0.5(g -5) 2 ,
where Z is a normalization constant making p 0 as a probability density. The numerical solutions are calculated for ε = 0.1, 0.15, 0.2, 0.3, 0.5 and 1 until t = 0.05. The results are shown in Figure 1. Here the red line is a reference line for the 1-dimensional manifold (v, g)

: v = V * (g) := g L V L +gV E g L +g
which carries the limiting Dirac distribution as ε → 0. As we proved it, the solutions to (1.6) indeed converges to Dirac delta distribution in v variable as ε → 0.

Next, we discuss the dynamics of the total firing rate N ε (t) for different values of ε and we compare their dynamics, still for equation (1.6). The dynamics of N ε of the equation (1.6) for the same set of ε as before is depicted in Figure 2. As we expect, the upper bound for N ε is around 1.65, which is always less than the φ F L ∞ = 2. Moreover, the time derivative of N ε (t) is also bounded after the initial layer 0 ≤ t ≤ τ , where τ decreases as ε decays.

As a final numerical test of this subsection for (1.6), we compare the case of the slow firing scaling (1.6) and the fast firing scaling (1.7). In Figure 3, we plot the section of the distributions at g = 5, where ε = 0.1 and t = 0.05. As we derived in Section 4, the limit distribution p(t, v, g) := lim ε→0 p ε (t, v, g) for fast firing regime has algebraic order near the critical voltage v = V * (g), while the limit distribution for slow firing regime is just a delta distribution. 5.2. Long term behaviors of the limit equations (3.1) and (3.5). We now move to the numerical simulations for asymptotic behaviors of the limit equations (3.1) and (3.5). We first consider the equation (3.1) which is a hyperbolic case, i.e., a = 0. For numerical simulation, we use b = 1 and the initial data

n 0 (g) = 1 Z exp -0.2(g -10) 2 .
Figure 4 shows the graph of φ

F g L V L +bN V E bN +g L
with respect to N . The intersection gives the stationary state. Since

b = 1 < 10 3 = 3 (1 -0.1)1 = g L (V E -V L ) φ F L ∞ , the condition (3.
2) is satisfied and we expect that a fixed point of map Φ is unique. As we can find in Figure 4, the map Φ has a unique fixed point N ∞ ≈ 1.384. Moreover, Figure 5 and Figure 6 show that the total firing rate N (t) converges to unique fixed point N ∞ exhibited in Figure 4 and the distribution n(t, g) converges to the Dirac delta distribution n ∞ (g) = δ bN∞ (g).

For the case of diffusive equation (3.5), we choose the diffusion coefficient a = 1 and the same initial data as in the case of hyperbolic case. Here, we chose parameter sets which may not satisfy the condition in Theorem 3.6, since we don't know the exact value of the Poincaré constant γ. However, even in this case, we are able to observe that the total firing rate N converges to equilibrium N ∞ (Figure 7) and the distribution n(t, g) converges to the Maxwellian equilibrium n ∞ (g) (Figure 8). 5.3. Bifurcating behavior of equilibrium. As a last numerical simulation, we present the behavior of equation (3.1) when φ F is not smooth and there are more than one steady state, i.e., a fixed point of φ F . As an example, we consider the following choice of φ F :

φ F (x) =      0, 0 < x < 0.3, 1, 0.3 < x < 0.8, 2,
x > 0.8.

We observe the firing rate dynamics for several values of b in Figure 9. The result of the numerical simulation implies that the equilibrium of the firing rate N (t) suddenly changed 10), this bifurcation is not due to the structure of φ F which does not vary qualitatively.

Conclusion

For a voltage-conductance nonlinear model (1.6) describing neural assemblies coupled through the total activity (firing rate) N of the network. We have established rigorously the conductance limit using rescalings of the equation. To do so, we obtain the uniform boundedness of moments of the solution of (1.6), together with the W 1,∞ boundedness of the total firing rate after the initial layer.

We have also established the asymptotic stability analysis of the limit conductance equation (1.8) for both cases a = 0 (hyperbolic model when noise is neglected) and a > 0 (cases including Gaussian noise). For this stability analysis, we provide a sufficient condition for the existence of a unique stationary solution to (1.8) and provide the asymptotic convergence of the solution toward this unique stationary solution. The method uses Lyapunov functionals measuring the distance of the solution from the stationary solution. All the results in the paper are supported by numerical simulations. In particular, we have focused on showing the stability of the solution to the limit equation. No oscillations or synchronisation effect occurs in such a simple conductance model, neither in the voltage-conductance with the regime ε 1, which does not contradict the periodic regime exhibited in [START_REF] Cáceres | A numerical solver for a nonlinear fokker-planck equation representation of neuronal network dynamics[END_REF], neither that in the voltage only model in [START_REF]Derivation of an integrate and fire equation for neural networks from a voltage-conductance kinetic model[END_REF]. However, there are still several remaining questions. A first issue would be to make more precise the convergence to the Dirac limiting solution and establish convergence rates. Secondly, the compactness in the case of fast firing, equation (1.7) is left open. Thirdly, as pointed out in Section 5.3, if there exist multiple fixed points of φ F , it is not proved that the solution of (3.1) will converge to Dirac mass on one of the fixed points, let alone specifying the fixed point on which the solution will concentrated. A hint towards this goal is a result using Monge-Kantorovich distance in [START_REF] Fournier | Monge-Kantorovich distance for PDEs: the coupling method[END_REF]. Furthermore, all the numerical solutions converge to its equilibrium and it does not show any oscillatory or periodic behavior. Thus, it will be also interesting if one can understand more physically based models where the conductance-only limit also shows periodic phenomena.

Lemma 2 . 1 .

 21 Let {p ε } ε>0 be a family of solution to (2.1), (2.2) with initial data p 0 ε ε>0

3. 1 .

 1 Hyperbolic equation. The limiting equation for the hyperbolic model is (3.1)

Proposition 3 . 3 .

 33 Suppose that b satisfies condition (3.2) in Proposition 3.1. Then, the map Ψ has a unique fixed point N ∞ .

Figure 1 .

 1 Figure 1. Numerical results for (1.6) at t = 0.05 with various values of ε illustrating the convergence to a Dirac distribution.

Figure 2 .

 2 Figure 2. The dynamics of N (t) for various ε in equation (1.6)

Figure 3 .

 3 Figure 3. Comparison of the section of distributions p ε (v, 5) with ε = 0.1 for slow and fast firing regimes.

Figure 4 .

 4 Figure 4. The graph of N and φ F g L V L +bN V E bN +g L .

Figure 5 .Figure 6 .

 56 Figure 5. The dynamics of N for a = 0.

Figure 7 .

 7 Figure 7. The dynamics of N for a = 1.

Figure 8 .

 8 Figure 8. The dynamics of n(t, g) when a = 1.

Figure 9 .Figure 10 .

 910 Figure 9. Bifurcating behavior of equilibrium

Acknowledgment. BP has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 740623).


Élodie Puybareau

Edwin Carlinet

Alessandro Benfenati

Hugues Talbot

Spherical fluorescent particle segmentation and tracking in 3D confocal microscopy

Keywords:

Spherical fluorescent particle are micrometer-scale spherical beads used in various areas of physics, chemistry or biology as markers associated with local physical media. They are useful for example in fluid dynamics to characterize flows, diffusion coefficients, viscosity or temperature; they are used in cells dynamics to estimate mechanical strain and stress at the micrometer scale. In order to estimate these physical measurements, tracking these particles is necessary. Numerous approaches and existing packages, both open-source and proprietary are available to achieve tracking with a high degree of precision in 2D. However, little such software is available to achieve tracking in 3D. One major difficulty is that 3D confocal microscopy acquisition is not typically fast enough to assume that the beads are stationary during the whole 3D scan. As a result, beads may move between planar scans. Classical approaches to 3D segmentation may yield objects are not spherical. In this article, we propose a 3D bead segmentation that deals with this situation.

Introduction

The use of fluorescent beads as physical markers is common in many research problems and applications related to fluid mechanics [START_REF] Park | Temperature measurement for a nanoparticle suspension by detecting the brownian motion using optical serial sectioning microscopy (ossm)[END_REF], rheology [START_REF] De Gennes | Scaling Concepts in Polymer Physics[END_REF], and cell biology [START_REF] Wang | Mechanotransduction across the cell surface and through the cytoskeleton[END_REF].

Spherical bead tracking in confocal microscopy is a classical and well-studied problem in image analysis. Recently, an ISBI challenge was organized to compare state-of-the-art bead segmentation and tracking methods [START_REF] Chenouard | Objective comparison of particle tracking methods[END_REF]. However, this very challenge has exhibited a number of shortcomings.

First among them is the common assumption that beads are very small, similar in diameter to the optical limit of resolution of the microscope, i.e. 0.2-0.4µm. In reality, beads can be much larger than this. 1-2µm diameters are commonly used. This is due to the fact that larger beads are less expensive, more visible and better adapted to some problems, for example when the surface of the bead needs to be activated. In particular, tracking of very small beads is difficult when there are many beads. Due to their larger dimension, a more precise location of their center is also possible, whereas smaller beads are typically only made of a few pixels and so a precise location of their center is imprecise.

A second assumption is that since beads are expected to be very small, they are located in a single plane and are so acquired in a single pass. Unfortunately, larger beads span several planes and need to be acquired over several scans, which takes time.

Thirdly, a common assumption is that scanning time is negligible compared with the studied motion. This is often true in mechanical cell studies. However, in fluid analysis, beads are subjected to Brownian motion [START_REF] Brown | A brief account of microscopical observations made in the months of june, july and august, 1827, on the particles contained in the pollen of plants; and on the general existence of active molecules in organic and inorganic bodies[END_REF] and can move significantly even in the absence of flow. This is a problem even during a single 3D acquisition since a bead can visible move as the confocal scan is taking place.

In this article, we study and propose solutions for comparatively large bead segmentation and tracking. Due to Brownian motion, tracking must take place both at the intra and inter-volume level. We present an application where we measure the diffusion coefficient of the Brownian motion. This is useful for measuring viscosity or temperature in fluid media [START_REF] Puybareau | Morphological analysis of brownian motion for physical measurements[END_REF].

Particle tracking

Our objective is to track relatively large beads (0.8-2µm in diameter) in 3D confocal microscopy. We do not assume that the diameter is known precisely. Typically most bead production processes allow a standard deviation of ± 10%. Traditional confocal microscopy proceeds by scanning a 3D volume slice-by-slice. Scanning a small area of a single slice, such as the few lines that may include a bead, is fast enough that the intersection of a scanning plane with any bead can be assumed to be a disk. However, scanning a whole slice can be slow, and so the particle may move slightly in between slice acquisition.

We do not assume that the particle forms a perfect sphere. We do assume that deviations from the sphere are relatively small. We also assume that the slice separation is small enough so that each bead is scanned multiple times. This is typically an adjustable parameter of the microscope.

Fig 1 summarize the pipeline we developed. We named "2D" the slices as shown in Fig. 2. These 5 slices composed the 3D volume for a fixed time.

2D segmentation

The first step of our pipeline is a 2D segmentation of the beads and uses both standard and modern mathematical morphology techniques [START_REF]Mathematical Morphology: from theory to applications[END_REF]. The objective is to obtain a binary image of the beads suitable for masking and labelling. The steps are the following:

1. Supression of small local minimas: this is performed using a closing by reconstruction with an initial dilation by a disk of radius 4. This yields image g. 2. Calculation of the min-tree of g.

3.

Filtering of nodes according to multiple criteria: size must be smaller than 500, the grey-level intensity must be lower than 100, and the contrast of the component greater than 80. We keep the maximal components meeting all these criteria. 4. Regularisation of the selected blobs: holes are filled, and the components are dilated with a 4-connected structuring element and opened with a disk of radius 5.

An example of results of this procedure is shown in Fig. 3 2

.2 3D regularization

For each acquisition time, we stack the corresponding 2D segmented slices. As the segmentation is performed in 2D, and the beads may move significantly between slices, some irregularities may appear at the bead boundaries, especially near the bead poles. The first step is hence a regularization in the z axis to reconnect these incomplete beads by applying a 2D closing with a disk of radius 1 in all the xz planes.

3D individualization

As some beads can touch, we separated them using a classical watershed procedure. After noise removal with an opening on the binary image, we extract beads markers using the distance transform with a threshold and the background marker. A watershed is then applied. The result is cleaned by removing small objects using a volume criterion. All the steps are illustrated in Fig. 4. Morphological parameters were empirically optimized on a small subset of the data.

Temporal identification

To label the beads at time t, we propagate the labels of the beads at the previous time t -1 by finding the intersection between the beads at t and t -1. These intersections are the markers used for the reconstruction by dilation. We hence obtain for time t labels for the beads also present at time t -1. The remaining unlabeled beads are given new labels. Overall, this procedure is simple but effective, as we will see in the result section.

In the previous section, 2D bead segmentation is performed and label clustered in 3D. We consequently assume that a uniquely labeled object in 3D constitutes a single bead.

Estimating the 3D position of each bead from their 2D slice is a challenge. We propose to estimate 2D positions first, and derive the 3D position an estimation of the radius of the bead.

Centroid estimation in 2D

Various methods can be used to estimate the centroid of a disk. A least-square best-fit disk can for example be used [START_REF] Gander | Least-squares fitting of circles and ellipses[END_REF][START_REF] Fitzgibbon | Direct least square fitting of ellipses[END_REF]. However these methods are more appropriate when only a subset of the disk contour is known. Since the beads are small, we can assume that most of them do not lie near the border of the field of view, and so are totally visible. In this case, a moments-based method seems more appropriate [START_REF] Hu | Visual pattern recognition by moment invariants[END_REF]. Let M pq be a moment of order p + q defined as follows:

M pq = x y x p y q I(x, y), (1)
where I(x, y) is the intensity of the image at position x, y. If we restrict I(x, y) to be equal to 1 for the segmented bead with label l and 0 elsewhere, then

xl = M 10 M 00 ; ȳl = M 01 M 00 , ( 2
)
yield the coordinates of the centroid of bead l. Moments of order 0 and 1 are not affected by isotropic blur or partial volume effects and are relatively resistant to noise, and so this estimation is rather robust, as we will see later in simulations. Moreover, we can also estimate the radius via the second order centered moments

µ 20 = M 20 -xM 10 ; µ 02 = M 02 -ȳM 01 , (3)
and then for a disk µ 20 = µ 02 = π 8 r 4 , or more simply from the disk area equation r = M 00 /π.

Centroid estimation in 3D with known bead radius

If we know the radius R of a bead to a good precision, we can estimate the depth p R by noting that

p R = ± R 2 -r 2 , ( 4
)
with r the estimated radius of the disk, as seen on Fig. 5 Unfortunately there is a sign uncertainty with this estimation, so a correct estimation must include information from the other slices associated with the same bead. Fig. 6 illustrates what happens, in a 2D x + z simplified setting. Each intersection between the bead and a scanning plane yields a disk (a line in the drawing), from which two depths can be estimated, only one of which is correct. The beads is moving between acquisitions, so there are two depths estimations per slice. We assume the correct ones are those that minimize the overall distance between them. 1 2 3 Fig. 6. When estimating the depth of the bead centroid over multiple slice acquisition (labeled 1, 2 and 3 in the drawing), as the bead is moving during the acquisition, the centroid depth varies and can be one of two possibilities. We select the depth estimations that minimizes the overall Euclidean distance between them (i.e. corresponding to the perimeter of the yellow triangle).

More formally, for every slice s, for every bead with label l we estimate two centroid positions

c l,0 s = (x l , ȳl , d(s) -p R (s)) (5) c l,1 s = (x l , ȳl , d(s) + p R (s)), (6)
with d(s) the depth of slice s and p given by (4). For every slice s we have two choices. If v ∈ {0, 1} ns is a binary vector of length n s , the total square distance is

d 2 v,R = ns-1 s=0 c l,v[s] s -c l,v[(s+1) mod ns] s 2 , ( 7
)
where a mod b is the remainder of the integer division of a by b, and assuming without loss of generality that slices are numbered 0 to n s -1.

We propose to select the vector v corresponding to the combination that minimizes d v . Consequently, if there are n s slices for a bead, we need to test 2 ns combinations. In our experiments, we have 3-6 slices for each bead, meaning up to 64 tests per bead, which is manageable. If we had more slices, this could be a problem. In this case, a simplification would be to select the middle slice, and for both centroid estimations related to this slice, compute the distance to the two centroid position estimations to every other slides. In other words we now minimize d 2 v,R defined as follows:

d 2 v,R = ns-1 s=0 c l,v[s] s -c l,v[ ns/2 ] s 2 , ( 8
)
where a is the integer part of a. The difference here is that each term of this sum can be minimized independently, and so there are only 4n s combinations to test, which is beneficial as soon as n s ≥ 4.

Centroid estimation in 3D with imprecise bead radius

If the radius is not known precisely and varies from bead to bead, which is often the case, the problem is more difficult. We assume that the bead radius for each bead is a deviate from a known Gaussian distribution of mean R and standard deviation σ R . Typically, σ r ≈ 0.1 R.

R l ∼ G( R, σ R ) (9)
Given the known observations for each bead, we also know that r l is greater than the largest measured radius r max = max{r s , s ∈ 0..n s -1 }. We solve the following mixed problem:

R * l , v * = arg min R,v d 2 v,R + 1 2λ R -R 2 (10) s.t. R ≥ r max ., ( 11
)
where λ is a Lagrangian parameter, and the second term of ( 10) is the loglikelihood of a Gaussian distribution for R, thus we expect λ to be proportional to σ 2 R . The tighter the distribution, the more important this term becomes. For a given v binary vector, ( 4) is convex due to the constraint [START_REF] Tao | The dc (difference of convex functions) programming and dca revisited with dc models of real world nonconvex optimization problems[END_REF], the constraint itself is linear since r max is an observation, but both ( 7) and ( 8) are smooth differences of convex functions, so this problem cannot be solved by standard constrained convex optimisation methods [START_REF] Boyd | Convex optimization[END_REF]. However, it can be tackled by DC methods [START_REF] Tao | The dc (difference of convex functions) programming and dca revisited with dc models of real world nonconvex optimization problems[END_REF]. In addition, the problem remains combinatorial in v, as seen above, so need to run several such DC problems, up to 2 ns or 4n s depending on whether we choose to optimise [START_REF] Hu | Visual pattern recognition by moment invariants[END_REF] or [START_REF]Mathematical Morphology: from theory to applications[END_REF].

Since DC method solvers are complex to handle, we propose to solve (10) using gradient descent with the d measure of ( 8), but taking as reference the slice with the largest measured radius r max instead of the middle slice. We set max{r max , R} as the initialization for R and λ = σ R . In practice, the true R * is very close to the largest observed radius.

Results on simulations

To test our method, we simulate N discrete sphere beads at high resolution with a random radius sampled from a Gaussian distribution as described in the previous section. The beads are subjected to isotropic Brownian motion and depth-ordered slice acquisition. We down-sample the result to ensure sub-pixel accurate positions. We estimate the radius of the intersection disk between each slice and the bead via the area method, and the x, ȳ position via moments. The z position and the bead radius is estimated as described above using the d method. For this simulation, we simulate high-resolution beads 80 pixels in radius, in a 350 3 cube, acquired through 7 layers and subjected to Brownian motion with an isotropic diffusion coefficient of 0.1 relative to the size of the cube, or 35 pixels at high resolution. Low resolution images were down-sampled by a factor of 10 with nearest-neighbor interpolation. The diffusion coefficient of the Brownian motion can be interpreted as the standard deviation of the the components of the position of the center of the bead from one slice to the next. With Brownian motion, the average deviation is zero.

Figure 7 shows the intersection of a single bead through successive slices. The Brownian motion influence is significant and clearly seen. We simulated 100 beads, which generated 540 slices. Via our optimisation procedure, each slice provides a 3D center estimate from which various measurements can be derived. In particular, we compare the estimation of the center position to the simulated ground truth. The results are shown in Table 1 On this table, we see that the average Euclidean position error in x and y is very low. Indeed, a pixel error in 350 (the high-resolution size of the cube) would translate into 2.8 10 -3 . However, the z estimation error is 5× greater, which is still sub-pixel accurate at the low sub-sampled resolution, and so still quite acceptable. We note that these estimates are the same as our estimate of the average Brownian deviation, which should be zero. On Fig. 8, we show the estimation of the relative Brownian diffusion coefficient as a function of the sample size. We expect a relative coefficient of 10 -1 , i.e. the one we simulated. The x and y are very close to this estimate, but the z diffusion estimation appears biased with a value about 70% as much as it should be. All three estimators converge quickly and can be deemed reliable from about 50 slices, i.e. 10 beads or so, with a standard deviation of the estimate of diffusion coefficient close to 10 -2 , ie. 10% error. With 40 beads, or 200 slices, the standard deviation of the diffusion coefficient drops down to 5 10 -3 or 5% relative error.

Currently, this means that our optimisation method is indeed good enough for positional estimation in 3D, but that for more subtle measurements, like diffusion coefficient estimations, 2D measurements are sufficient. We note that we are able to measure this diffusion coefficient from single 3D acquisition only.

Results on real data

We applied our method on a video of 3µm beads, with ∆t = 0.5s. Resolution is 0.125µm per pixel in x and y, and only 0.5µm per slide in z. Comparison with manual labelling results show that we manage to correctly detect and label all beads, as well as track each bead individually with zero label error between 3D volumes. However, we overestimate beads presence in the z-axis due to the point spread function (PSF) of the microscope. In practice this means that beads faintly appear in slices above and below their true position. This overestimation is depicted in Fig. 9(c). On slice 0, we tend to find beads that are not supposed to be segmented (false positives), i.e.the central bead of the illustration. This "oversegmentation" is nonetheless acceptable as it remains coherent in 3D and does not imply bias in 2D. In 3D, this can cause the bead radius to be overestimated, and will need to be corrected in future work. On this sample, using our 2D method, we estimate an average diffusion coefficient of 4.3 10 -2 relative to the beads average diameter. Currently we do not have the ground truth for this measurement.

Conclusion

In this article we tacked the problem of tracking fluorescent beads in confocal microscopy, when we do not assume that beads are small enough to be assimilated to point sources, but instead span several slices; when they move between slices due to Brownian motion; and when their radius is not precisely known. Starting from 2D-based segmentation and 3D clustering using shape-based mathematical morphology, we have shown that estimating bead radius and true 3D positions from a small number of 2D slices is a difficult, but tractable problem. We have also shown that tracking these beads between 3D acquisitions is possible with a high degree of accuracy. One benefit of our approach is that we can obtain good estimates of the diffusion coefficient of Brownian motion from single 3D acquisition, which is useful for various physical measurements.

In future work, we plan to improve the quality and speed of our optimisation method, to make better use of prior knowledge in 3D particularly for initial segmentation, and we will release our software as a free, open-source package.

Fig. 1 .

 1 Fig. 1. Flowchart of our pipeline

Fig. 2 .Fig. 3 .

 23 Fig. 2. Consecutive depth (z) slices corresponding to a single time.

Fig. 4 .

 4 Fig. 4. Results of the successive steps: (a) initial image, (b) inital segmentation, (c) 3d regularization, (d) bead separation, (e) individualization and labeling.

Fig. 5 .

 5 Fig. 5. Bead depth estimation: R is the real radius, r the observed radius, and p the estimated depth.

Fig. 7 .

 7 Fig. 7. A single simulated bead acquired over several layers subjected to Brownian motion.

Fig. 8 .

 8 Fig. 8. Convergence of the diffusion estimator to a stable value (a) and its standard deviation as a function of the number of slices (b).

Fig. 9 .

 9 Fig. 9. Comparison of segmentation results: (a) original frame, (b) results of the first step of the 2D segmentation procedure before regularization, and (c) after 3D regularization. Top images are the entire xy slices, and the bottom is a xz view along the z-axis

Table 1 .

 1 Positional RMS relative error

	x	y	z
	Relative error 1.6 10 -3 1.6 10 -3 8.2 10 -3