Susanne Ditlevsen

Adeline Samson

ESTIMATION IN THE PARTIALLY OBSERVED STOCHASTIC MORRIS-LECAR NEURONAL MODEL WITH PARTICLE FILTER AND STOCHASTIC APPROXIMATION METHODS

Keywords: Sequential Monte Carlo, diffusions, pseudo likelihood, Stochastic Approximation Expectation Maximization, motoneurons, conductance-based neuron models, membrane potential

Parameter estimation in multi-dimensional diffusion models with only one coordinate observed is highly relevant in many biological applications, but a statistically difficult problem. In neuroscience, the membrane potential evolution in single neurons can be measured at high frequency, but biophysical realistic models have to include the unobserved dynamics of ion channels. One such model is the stochastic Morris-Lecar model, defined by a non-linear two-dimensional stochastic differential equation. The coordinates are coupled, i.e. the unobserved coordinate is non-autonomous, the model exhibits oscillations to mimick the spiking behavior, which means it is not of gradient-type, and the measurement noise from intra-cellular recordings is typically negligible. Therefore the hidden Markov model framework is degenerate, and available methods break down. The main contributions of this paper are an approach to estimate in this ill-posed situation, and non-asymptotic convergence results for the method. Specifically, we propose a sequential Monte Carlo particle filter algorithm to impute the unobserved coordinate, and then estimate parameters maximizing a pseudo-likelihood through a stochastic version of the Expectation-Maximization algorithm. It turns out that even the rate scaling parameter governing the opening and closing of ion channels of the unobserved coordinate can be reasonably estimated. An experimental data set of intracellular recordings of the membrane potential of a spinal motoneuron of a red-eared turtle is analyzed, and the performance is further evaluated in a simulation study.

Introduction.

In neuroscience, it is of major interest to understand the principles of information processing in the nervous system, and a basic step is to understand signal processing and transmission in single neurons. Therefore, there is a growing demand for robust methods to estimate biophysical relevant parameters from partially observed detailed models. Statistical inference from experimental data in biophysically detailed models of single neurons is difficult. Often these models are compared to experimental data by hand-tuning to reproduce the qualitative behaviors observed in experimental data, but without any formal statistical analysis. It is of particular interest to estimate conductances, which reflect the synaptic input from the surrounding network. These can be estimated from intracellular recordings, where the neuronal membrane potential is recorded at high frequency, and are typically done using only subthreshold fluctuations, ignoring the dynamics during action potentials [START_REF] Berg | Synaptic inhibition and excitation estimated via the time constant of membrane potential fluctuations[END_REF][START_REF] Berg | Balanced inhibition and excitation drive spike activity in spinal halfcenters[END_REF][START_REF] Borg-Graham | Visual input evokes transient and strong shunting inhibition in visual cortical neurons[END_REF][START_REF] Monier | In vitro and in vivo measures of evoked excitatory and inhibitory conductance dynamics in sensory cortices[END_REF][START_REF] Pospischil | Extracting synaptic conductances from single membrane potential traces[END_REF][START_REF] Rudolph | A method to estimate synaptic conductances from membrane potential fluctuations[END_REF]. The aim of this article is to estimate such biophysical parameters during the dynamics of spiking from intra-cellular data.

The Morris-Lecar model [START_REF] Morris | Voltage oscillations in the barnacle giant muscle fiber[END_REF] is a simple biophysical model, and a prototype for a wide variety of neurons. It is a conductance-based model [START_REF] Gerstner | Spiking neuron models -Single neurons, populations, plasticity[END_REF], introduced to explain the dynamics of the barnacle muscle fiber. It is given by two coupled first order differential equations, the first modeling the membrane potential evolution, and the second the activation of potassium current. If both current and conductance noise should be taken into account, the stochastic Morris-Lecar model arises, where diffusion terms have been added on both coordinates. If one of these noise sources are zero, a hypoelliptic diffusion arises leading to singular transition densities and particular statistical challenges [START_REF] Pokern | Parameter estimation for partially observed hypoelliptic diffusions[END_REF][START_REF] Samson | A contrast estimator for completely or partially observed hypoelliptic diffusion[END_REF]. Typically, the membrane potential will be measured discretely at high frequency, whereas the second variable cannot be observed. Our goal is to estimate model parameters from discrete observations of the first coordinate in the non-singular case of non-negligible noise on both coordinates. This includes estimation of a central rate parameter characterizing the channel kinetics of the unobserved component, which we believe has not been done before.

Estimation in these conductance-based models is not straightforward. Because of the coupling between the coordinates of the stochastic differential equation (SDE), the unobserved coordinate is non-autonomous, and the model does not fit into the (non-degenerate) Hidden Markov Model (HMM) framework, as explained in Section 2.3. Furthermore, the diffusion is not time reversible and the likelihood is generally not tractable. Thus, the problem of inference is complex. The literature contains various methodologies when all the coordinates are observed [START_REF] Aït-Sahalia | Maximum likelihood estimation of discretely sampled diffusions: a closed-form approximation approach[END_REF][START_REF] Beskos | Exact and computationally efficient likelihood-based estimation for discretely observed diffusion processes[END_REF][START_REF] Durham | Numerical techniques for maximum likelihood estimation of continuous-time diffusion processes[END_REF][START_REF] Jensen | Markov chain Monte Carlo approach to parameter estimation in the FitzHugh-Nagumo model[END_REF][START_REF] Pedersen | A new approach to maximum likelihood estimation for stochastic differential equations based on discrete observations[END_REF][START_REF] Sørensen | Parametric inference for diffusion processes observed at discrete points in time: a survey[END_REF][START_REF] Sørensen | Statistical methods for stochastic differential equations, chapter Estimating functions for diffusion-type processes[END_REF] or the hidden state is Markovian [START_REF] Ionides | Iterated filtering[END_REF]. They strongly rely on the Markov property and are hard to generalize to the non-Markovian case we are studying. In the non-Markovian case, methods are mainly based on data augmentation. The idea is that the likelihood can be approximated given the entire path or a sufficient partition of it. Therefore the unobserved coordinates are treated as missing data and are imputed. Most methods propose to approximate the transition density by the Euler-Maruyama scheme and consider a Bayesian point of view to estimate the posterior distribution of the parameters [START_REF] Elerian | Likelihood inference for discretely observed nonlinear diffusions[END_REF][START_REF] Eraker | MCMC analysis of diffusion models with application to finance[END_REF][START_REF] Golightly | Bayesian sequential inference for nonlinear multivariate diffusions[END_REF][START_REF] Golightly | Bayesian inference for nonlinear multivariate diffusion models observed with error[END_REF]. [START_REF] Golightly | Bayesian sequential inference for nonlinear multivariate diffusions[END_REF] study a model similar to us but with low frequency data. So they need to impute data between observations, which is computationally costly. Furthermore, there exists a strong dependence between the imputed sample paths and the diffusion coefficient and it is not possible to estimate the diffusion parameter with this kind of approach. An alternative is reparametrisation of the diffusion but is limited to scalar diffusions [START_REF] Roberts | On inference for partially observed nonlinear diffusion models using the Metropolis-Hastings algorithm[END_REF] or an autonomous hidden coordinate [START_REF] Kalogeropoulos | Likelihood-based inference for a class of multivariate diffusions with unobserved paths[END_REF].

In this paper, we propose to estimate the parameters with a maximum likelihood approach. We approximate the SDE through an Euler-Maruyama scheme to obtain a tractable pseudo-likelihood. Then we consider the statistical model as an incomplete data model, and maximize the pseudo-likelihood through a stochastic Expectation-Maximization (EM) algorithm, where the unobserved data are imputed at each iteration of the algorithm. We are in the setting of high frequency data so we do not need to impute data between observations, but our approach could be extended to that type of data as well. A similar but different method has been proposed by [START_REF] Huys | Efficient estimation of detailed single-neuron models[END_REF], where up to 10 4 parameters are estimated in a detailed multi-compartmental single neuron model. However, only parameters entering linearly in the loss function are considered, and channel kinetics are assumed known. It is a quadratic optimization problem solved by least squares, and shown to work well for low noise and high frequency sampling. When either the discretization step or the noise increase, a bias is introduced. In [START_REF] Huys | Smoothing of, and Parameter Estimation from, Noisy Biophysical Recordings[END_REF] they extend the estimation to allow for measurement noise, first smoothing the data by a particle filter, and then maximizing the likelihood through a Monte Carlo EM-algorithm. Because of the measurement noise, the model fits into the HMM framework and they can use a standard particle filter. But again, only parameters entering linearly in the pseudo-likelihood are considered. In particular, all parameters of the hidden coordinate are assumed known.

Here, we also want to estimate parameters from the hidden coordi- [START_REF] Fearnhead | Particle filters for partially observed diffusions[END_REF] generalises this particle filter to a non-autonomous hidden path but with drift of gradient type. In the ergodic case this corresponds to a time reversible diffusion. In particular, models exhibiting oscillations are not covered, which is the case of any realistic neuronal model. These algorithms cannot be directly applied because we are studying a multi-dimensional coupled SDE that is not of gradient type. Thus, we consider the SMC algorithm proposed by [START_REF] Doucet | An introduction to sequential Monte Carlo methods[END_REF] for more general dynamic models than HMM. As we combine this SMC with the Stochastic Approximation Expectation-Maximization (SAEM) algorithm which maximizes the pseudo-likelihood based on an Euler-Maruyama approximation of the SDE defining the model, we need non-asymptotic convergence results for the SMC to obtain the convergence of the SAEM-SMC. Non-asymptotic results for SMC, such as deviation inequalities, have been proposed in the literature only in the HMM framework [START_REF] Del Moral | Branching and interacting particle systems. approximations of Feynman-Kac formulae with applications to non-linear filtering[END_REF][START_REF] Del Moral | The Monte-Carlo method for filtering with discrete-time observations[END_REF][START_REF] Douc | Sequential Monte Carlo smoothing for general state space hidden Markov models[END_REF][START_REF]Recursive Monte Carlo filters: algorithms and theoretical analysis[END_REF], and the Markovian structure of the hidden path is a key element in the proofs. A major contribution here is that we are able to extend this result to a SMC for a non-Markovian hidden path. Then we prove that the estimator obtained from this combined SAEM-SMC algorithm converges with probability one to a local maximum of the pseudo-likelihood. We also prove that the pseudo-likelihood converges to the true likelihood as the time step between observations go to zero.

The paper is organized as follows: In Section 2 the model is presented, the noise structure is motivated, and the pseudo likelihood arising from the Euler-Maruyama approximation is found. In Section 3, the filtering problem is presented, as well as the SMC algorithm and deviation inequalities. In Section 4 we present the estimation procedure and the assumptions needed for the convergence results to hold. In Section 5 we apply the method on an experimental data set of intracellular record-ings of the membrane potential of a motoneuron of a turtle, and in Section 6 we conduct a simulation study to document the performance of the method. Proofs and technical results can be found in the Appendix.

Stochastic Morris

dV t = f (V t , U t )dt + γd Bt , dU t = b(V t , U t )dt + σ(V t , U t )dB t ,
where

f (V t , U t ) = 1 C (-g Ca m ∞ (V t )(V t -V Ca ) -g K U t (V t -V K ) -g L (V t -V L ) + I) , b(V t , U t ) = (α(V t )(1 -U t ) -β(V t )U t ) , m ∞ (v) = 1 2 1 + tanh v -V 1 V 2 , α(v) = 1 2 φ cosh v -V 3 2V 4 1 + tanh v -V 3 V 4 , β(v) = 1 2 φ cosh v -V 3 2V 4 1 -tanh v -V 3 V 4 ,
and the initial condition (V 0 , U 0 ) is random with density p(V 0 , U 0 Various noise sources are present in single neurons, and they act on many different spatial and temporal scales [START_REF] Gerstner | Spiking neuron models -Single neurons, populations, plasticity[END_REF][START_REF] Longtin | Neuronal noise[END_REF]. A main component arises from the synaptic bombardment from other neurons in the network, and in the diffusion limit appears as an additive noise on the current equation. Parameter γ scales this current noise. Conductance fluctuations caused by random opening and closing of ion channels leads to multiplicative noise on the conductance equation. Function σ(V t , U t ) models this channel or conductance noise. We consider the following function that ensures that U t stays bounded in the unit interval if σ ≤ 1

[13]: σ(V t , U t ) = σ 2 α(Vt)β(Vt) α(Vt)+β(Vt) U t (1 -U t ).
A trajectory of the model is simulated in Fig. 1. The peaks of (V t ) correspond to spikes of the neuron.

2.2.

Observations and approximate model. Data are discrete measurements of (V t ) while (U t ) is not measured. We denote t 0 ≤ t 1 ≤ • • • ≤ t n the discrete observation times. We denote V i = V t i the observation at time t i and V 0:n = (V t 0 , . . . , V tn ) the vector of all the observed data. Let θ ∈ Θ ⊆ R p be the vector of parameters to be estimated. We consider estimation of all identifiable parameters of the observed coordinate, and the rate parameter of the unobserved channel dynamics θ = (g Ca , g K , g L , V Ca , V K , I, γ, φ). Note that C is a proportionality factor of the conductance parameters and thus unidentifiable, as well as the constant level in f (•) is given by g L V L + I, and thus V L (or I) is unidentifiable. We conjecture that the information about σ in the observed coordinate is close to zero, and thus, in practice also σ is unidentifiable from observations of V 0:n only, at least for any finite sample size. This happens because σ is mainly shaping the dynamics of U t between spikes, while the dynamics during spikes resemble deterministic behavior, and the influence of U t on V t is only strong during spikes. This is confirmed in Sections 5 and 6 where misspecification of σ is shown not to deteriorate the estimation of θ. Finally, we assume the scaling parameters V 1 -V 4 known because otherwise the model does not belong to an exponential family, as required by assumption (M1) below. This could be solved by introducing an extra optimization step in the EM-algorithm at the cost of precision and computer time. It is not pursued further in this work.

The aim is to estimate θ by maximum likelihood. However, this likelihood is intractable, as the transition density of model ( 1) is not explicit. Let ∆ denote the step size between two observation times, which we for simplicity assume does not depend on i. The extension to unequally spaced observation times is straightforward. The Euler-Maruyama approximation of model ( 1) leads to a discretized model defined as follows

V i+1 = V i + ∆f (V i , U i ) + √ ∆ γ ηi , (2)
U i+1 = U i + ∆b(V i , U i ) + √ ∆σ(V i , U i )η i ,
where (η i ) and (η i ) are independent centered Gaussian variables. To ease readability the same notation (V i , U i ) is used for the original and the approximated processes. This should not lead to confusion, as long as the transition densities are distinguished, as done below.

2.3.

Property of the observation model. The observation model is a degenerate HMM. Let us recall the definition proposed by [START_REF] Cappé | Inference in Hidden Markov Models[END_REF]: A HMM with not countable state space is defined as a bivariate Markov chain (X i , Y i ) with only partial observations Y i , whose transition kernel has a special structure: both the joint process (X i , Y i ) and the marginal hidden chain (X i ) are Markovian.

In our model,

(U i ) is not Markovian, only (V i , U i ) is Markovian. So set X i = (V i , U i ), with Markov kernel R(X i-1 , dX i ) = p ∆ (dV i , dU i |V i-1 , U i-1 ),
the transition density of model ( 2), and

Y i = X (1) i , the first coordinate of X i with transition kernel F (X, dY ) = 1 {Y =X (1) } .
Here, 1 x is the Dirac measure in x. Thus, the kernel F is zero almost everywhere and the HMM is degenerate. This leads to an intrinsic degeneracy of the particle filter used in the standard HMM toolbox, as explained below.

Therefore we consider the observation model as a bivariate Markov chain (V i , U i ) with only partial observations V i whose hidden coordinate U i is not Markovian. It is not a HMM but a general dynamic model as considered by [START_REF] Andrieu | Sequential Monte Carlo methods for optimal filtering[END_REF]. The hidden process U i is distributed as

U 0 ∼ µ(dU 0 ), U i |(U 0:i-1 , V 0:i-1 ) ∼ K(dU i |U 0:i-1 , V 0:i-1 )
for some conditional distribution function K and the observed process V i is distributed as

V i |(U 0:i , V 0:i-1 ) ∼ G(dV i |U 0:i , V 0:i-1 )
for some distribution function G. Given the Markovian structure of the pair

(V i , U i ), we have K(dU i |U 0:i-1 , V 0:i-1 ) = K(dU i |U i-1 , V i-1 ) and G(dV i | U 0:i , V 0:i-1 ) = G(dV i |U i-1:i , V i-1 ).
To simplify, we use the same notation for random variables and their realizations and assume that G(

dV i |U 0:i , V 0:i-1 ) = G(V i |U 0:i , V 0:i-1 )dV i .
2.4. Likelihood function. We want to estimate the parameter θ by maximum likelihood of the approximate model, with likelihood

(3) p ∆ (V 0:n ; θ) = p(V 0 , U 0 ; θ) n i=1 p ∆ (V i , U i |V i-1 , U i-1 ; θ)dU 0:n .
It corresponds to a pseudo-likelihood for the exact diffusion. The multiple integrals of equation ( 3) are difficult to handle and it is not possible to maximize the pseudo-likelihood directly.

A solution is to consider the statistical model as an incomplete data model. The observable vector V 0:n is then part of a so-called complete vector (V 0:n , U 0:n ), where U 0:n has to be imputed. To maximize the likelihood of the complete data vector (V 0:n , U 0:n ), we propose to use a stochastic version of the EM algorithm, namely the SAEM algorithm [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF]. Simulation under the smoothing distribution p ∆ (U 0:n |V 0:n ; θ) is likely to be difficult, and direct simulation of the non-observed data (U 0:n ) is not possible. A SMC algorithm, also known as Particle Filtering, provides a way to approximate this distribution [START_REF] Doucet | An introduction to sequential Monte Carlo methods[END_REF]. We have adapted this algorithm to handle a coupled two-dimensional SDE, i.e. the unobserved coordinate is non-autonomous and non-Markovian. Then, we combine the SAEM algorithm with the SMC algorithm, where the unobserved data are filtered at each iteration step, to estimate the parameters of model [START_REF] Andrieu | Sequential Monte Carlo methods for optimal filtering[END_REF]. Details on the filtering are given in Section 3, and the SAEM algorithm is presented in Section 4.1. To prove the convergence of this new SAEM-SMC algorithm, a non-asymptotic deviation inequality is required for the SMC algorithm. Then we derive the convergence of the SAEM-SMC algorithm to a maximum of the likelihood.

Filtering.

3.1. The filtering problem and the SMC algorithm. For any bounded Borel function f : R → R, we denote π n,θ f = E ∆ (f (U n )|V 0:n ; θ), the conditional expectation under the exact smoothing distribution p ∆ (U 0:n |V 0:n ; θ) of the approximate model. The aim is to approximate this distribution for a fixed value of θ. When included in the stochastic EM algorithm, this value will be the current value θ m at the given iteration. For notational simplicity, θ is omitted in the rest of this Section.

We now argue why the HMM point of view is ill-posed for the filtering problem. Considering the model as a HMM,

X i = (V i , U i ) is the hidden Markov chain and Y i = X (1) i . But then the filtering problem π n f is the ratio of µ(dU 0 )R(X 0 , dX 1 )F (X 0 ; Y 1 ) • • • R(X n-1 , dX n )F (X n-1 ; Y n )f (X n ) and µ(dU 0 ) R(X 0 , dX 1 ) F (X 0 ; Y 1 ) • • • R(X n-1 , dX n ) F (X n-1 ; Y n ). Since F (X n-1 ; Y n ) = 1 {Yn=X (1)
n-1 } and the state space is continuous, the denominator is zero almost surely and the filtering problem is ill-posed. Now consider the model in a more general framework with the hidden state U i not Markovian, and introduce for i = 1, . . . , n the kernels H i from R into itself by [START_REF] Berg | Balanced inhibition and excitation drive spike activity in spinal halfcenters[END_REF]

H i f (u) = K(dz|u, V i-1 )G(V i |u, V i-1 , z)f (z) = p ∆ (V i , z|V i-1 , u)f (z)dz.
Then π n can be expressed recursively by

π n f = π n-1 H n f π n-1 H n 1 = µ(U 0 ) n i=1 p ∆ (V i , U i |V i-1 , U i-1 )f (U n )dU 0:n µ(U 0 ) n i=1 p ∆ (V i , U i |V i-1 , U i-1 )dU 0:n . (5)
Note that the denominator of ( 5) is

µH 1 • • • H n 1 = p ∆ (V 0:n ),
which is different from 0 since it has support the real line. Thus, the filtering problem is well-posed.

The kernels H i are extensions of the kernels considered by [START_REF] Del Moral | The Monte-Carlo method for filtering with discrete-time observations[END_REF] in the context of two-dimensional SDEs with hidden coordinate U t autonomous (and thus Markovian). We do not extend their particle filter since it is based on simulation of both V i and U i with transition kernel p ∆ (V i , U i |V i-1 , U i-1 ). They avoid the degeneracy of the weights by introducing an instrumental function ψ and the weights are computed as ψ(V

(k) i -V i ).
The choice of this instrumental function may influence the numerical properties of the filter. Therefore, we adopt the general filter proposed by [START_REF] Andrieu | Sequential Monte Carlo methods for optimal filtering[END_REF] for more general dynamic system, that we recall here.

The SMC algorithm provides a set of K particles (U

(k) 0:n ) k=1...K and weights (W (k) 0:n ) k=1...K approximating the conditional smoothing distri- bution p ∆ (U 0:n |V 0:n ) [see 16]. The SMC method relies on proposal distri- butions q(U i |V i , V i-1 , U i-1 )
to sample what we call particles from these distributions. We write V 0:i = (V 0 , . . . , V i ) and likewise for U 0:i . Algorithm 1 (SMC algorithm)

• At time i = 0: ∀ k = 1, . . . , K 1. sample U (k) 0 from p(U 0 |V 0 ) 2.

compute and normalize the weights:

w 0 U (k) 0 = p V 0 , U (k) 0 , W 0 U (k) 0 = w 0 U (k) 0 K k=1 w 0 U (k) 0 • At time i = 1, . . . , n: ∀ k = 1, . . . , K 1. Sample indices A (k) i-1 ∼ r(•|W i-1 (U (1) 0:i-1 ), . . . , W i-1 (U (K) 0:i-1 )) and set U (k) 0:i-1 = U (A (k) i-1 ) 0:i-1 2. sample U (k) i ∼ q •|V i-1:i , U (k) i-1 and set U (k) 0:i = (U (k) 0:i-1 , U (k) i )

compute and normalize the weights

W i (U (k) 0:i ) = w i U (k) 0:i K k=1 w i U (k) 0:i with w i U (k) 0:i = p ∆ V 0:i ,U (k) 0:i p ∆ V 0:i-1 ,U (k) 0:i-1 q U (k) i |V i-1:i ,U (k) 0:i-1
Finally, the SMC algorithm provides an empirical measure

Ψ K n = K k=1 W n (U (k) 0:n )1 U (k) 0:n
which is an approximation to the smoothing distribution p ∆ (U 0:n |V 0:n ). A draw from this distribution can be obtained by sampling an index k from a multinomial distribution with probabilities

W n (U (1) 0:n ), . . . , W n (U (K) 0:n ) and setting the draw U 0:n equal to U 0:n = U (k) 0:n . The variable A (k)
i-1 plays an important role to discard the samples with small weights and multiply those with large weights [START_REF] Gordon | Novel approach to nonlinear/non-gaussian bayesian state estimation[END_REF]. It generates a number of offspring N

( ) i-1 , = 1, . . . , K, such that K =1 N ( ) i-1 = K and E(N ( ) i-1 ) = K W i-1 (U (l) 0:i-1
). Many schemes for r have been presented in the literature, including multinomial sampling [START_REF] Gordon | Novel approach to nonlinear/non-gaussian bayesian state estimation[END_REF], residual sampling [START_REF] Liu | Sequential Monte Carlo methods for dynamic systems[END_REF] or stratified resampling [START_REF] Doucet | On sequential Monte Carlo sampling methods for bayesian filterin[END_REF]. They differ in terms of var(N

( ) i-1 ) [see 15].
The key property that we need in order to prove the deviation inequality is that E(1

{A (k) i-1 = } ) = W i-1 (U (l) 0:i-1 ).
Since our model is not a HMM, the weights w i U (k) 0:i cannot be written in terms of a Markov transition kernel of the hidden path as is usually done. It follows that the proposal q, which is crucial to ensure good convergence properties, has to depend on

V i . The first classical choice of q is q(U i |V i-1:i , U i-1 ) = p ∆ (U i |V i-1 , U i-1
), i.e. the transition density. In this case, the weight reduces to

w i U (k) 0:i = p ∆ (V i |V i-1 , U (k) 0:i ). A second choice for the proposal is q(U i |V i-1:i , U i-1 ) = p ∆ (U i |V i-1:i , U i-1
), i.e. the conditional distribution. In this case, the weight reduces to

w i U (k) 0:i = p ∆ (V i |V i-1 , U (k) 0:i-1
). Transition densities and conditional distributions are detailed in Appendix A. When the two Brownian motions are independent, as we assume, the two choices are equivalent.

This SMC algorithm is plugged into the EM algorithm to estimate the parameters. We thus need non-asymptotic convergence results on the SMC algorithm to ensure the convergence of the EM algorithm. This is discussed in the next section.

Deviation inequality.

In the literature, deviation inequalities for SMC algorithms only appear for HMM. To our knowledge, this is the first non-asymptotic result proposed for a SMC applied to a non-Markovian hidden path. The only result of this type with SDEs has been obtained by [START_REF] Del Moral | The Monte-Carlo method for filtering with discrete-time observations[END_REF], with autonomous second coordinate. Here, we generalize their deviation inequality to a non-autonomous hidden path.

For a bounded Borel function f , denote

Ψ K n f = K k=1 f (U (k) n )W n,θ (U (k)
0:n ), the conditional expectation of f under the empirical measure Ψ K n,θ obtained by the SMC algorithm for a given value of θ. We have: Proposition 1. Under assumption (SMC3), for any ε > 0, and for any bounded Borel function f on R, there exist constants C 1 and C 2 , independent of θ, such that

P Ψ K n,θ f -π n,θ f ≥ ε ≤ C 1 exp -K ε 2 C 2 f 2 (6)
where f is the sup-norm of f . The proof is provided in Appendix D. A similar result can be obtained with respect to the exact smoothing distribution of the exact diffusion model, under assumptions on the number of particles and the step size of the Euler approximation.

Estimation method.

4.1. SAEM algorithm. The EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] is useful in situations where the direct maximization of the marginal likelihood θ → p ∆ (V 0:n ; θ) is more difficult than the maximization of the conditional expectation of the complete likelihood

Q(θ|θ ) = E ∆ [log p ∆ (V 0:n , U 0:n ; θ)|V 0:n ; θ ] ,
where p ∆ (V 0:n , U 0:n ; θ) is the likelihood of the complete data (V 0:n , U 0:n ) of the approximate model ( 2) and the expectation is under the conditional distribution of U 0:n given V 0:n with density p ∆ (U 0:n |V 0:n ; θ ). The EM algorithm is an iterative procedure: at the mth iteration, given the current value θ m-1 , the E-step is the evaluation of Q m (θ) = Q(θ | θ m-1 ), while the M-step updates θ m-1 by maximizing Q m (θ). To fulfill convergence conditions of the algorithm, we consider the particular case of a distribution from an exponential family. More precisely, we assume:

(M1)
The parameter space Θ is an open subset of R p . The complete likelihood p ∆ (V 0:n , U 0:n ; θ) belongs to a curved exponential family, i.e.

log p ∆ (V 0:n , U 0:n ; θ) = -ψ(θ) + S(V 0:n , U 0:n ), ν(θ)
, where ψ and ν are two functions of θ, S(V 0:n , U 0:n ) is known as the minimal sufficient statistic of the complete model, taking its value in a subset S of R d , and •, • is the scalar product on R d .

The approximate Morris-Lecar model (2) satisfies this assumption when the scaling parameters V 1 , V 2 , V 3 and V 4 are known. Details of the sufficient statistic S are given in Appendix B.

Under assumption (M1), the E-step reduces to the computation of

E ∆ S(V 0:n , U 0:n )|V 0:n ; θ m-1 .
When this expectation has no closed form, [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF] propose the Stochastic Approximation EM algorithm (SAEM) replacing the E-step by a stochastic approximation of Q m (θ). The E-step is then divided into a simulation step (S-step) of the non-observed data (U (m) 0:n ) with the conditional density p ∆ (U 0:n |V 0:n ; θ m-1 ) and a stochastic approximation step (SA-step) of E ∆ S(V 0:n , U 0:n )|V 0:n ; θ m-1 with a sequence of positive numbers (a m ) m∈N decreasing to zero. We write s m for the approximation of this expectation. At the S-step, the simulation under the smoothing distribution is done by SMC, as explained in Section 3. We call this algorithm the SAEM-SMC algorithm. Iterations of the SAEM-SMC algorithm are written as follows: SA-Step: update s m-1 using the stochastic approximation:

(7) s m = s m-1 + a m-1 S(V 0:n , U (m) 0:n ) -s m-1 M-Step: update of θ m by θ m = arg max θ∈Θ (-ψ(θ) + s m , ν(θ) ) .
Standard errors of the estimators can be evaluated through the Fisher information matrix. Details are given in Appendix C. An advantage of the SAEM algorithm is the low-level dependence on the initialization θ 0 , due to the stochastic approximation of the E-step. The other advantage over a Monte-Carlo EM algorithm is the computational time. Indeed, only one simulation of the hidden variables U 0:n is needed in the simulation step while an increasing number of simulated hidden variables is required in a Monte-Carlo EM algorithm.

Convergence of the SAEM-SMC algorithm.

The SAEM algorithm we propose in this paper is based on an approximate simulation step performed with an SMC algorithm. We prove that even if this simulation is not exact, the SAEM algorithm still converges towards the maximum of the likelihood of the approximated diffusion (2). This is true because the SMC algorithm has good convergence properties.

Let us be more precise. We introduce a set of convergence assumptions which are the classic ones for the SAEM algorithm [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF].

)|F m ) = f (u)p ∆ (u|v, θ m )du.
Assumptions (M1)-(M5) ensure the convergence of the EM algorithm when the E-step is exact [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF]. Assumptions (M1)-(M5) and (SAEM1)-(SAEM4) together with the additional assumption that (s m ) m≥0 takes its values in a compact subset of S ensure the convergence of the SAEM estimates to a stationary point of the observed likelihood p ∆ (V 0:n ; θ) when the simulation step is exact [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF].

Here the simulation step is not exact and we have three additional assumptions on the SMC algorithm to bound the error induced by this algorithm and prove the convergence of the SAEM-SMC algorithm. Theorem 1 is proved in Appendix D. Note that assumption (SAEM4) is not needed thanks to the conditional independence of the particles generated by the SMC algorithm, as detailed in the proof. Similarly, the additional assumption that (s m ) m≥0 takes its values in a compact subset of S is not needed, as it is directly satisfied under assumption (SMC2).

We deduce that the SAEM algorithm converges to a (local) maximum of the likelihood under standard additional assumptions (LOC1)-(LOC3) proposed by [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF] on the regularity of the log-likelihood ∆ (V 0:n ; θ) that we do not recall here. COROLLARY 1. Under the assumptions of Theorem 1 and additional assumptions (LOC1)-(LOC3), the sequence θ m converges with probability 1 to a (local) maximum of the likelihood p ∆ (V 0:n ; θ).

The classical assumptions (M1)-(M5) are usually satisfied. Assumption (SAEM1) is easily satisfied by choosing properly the sequence (a m ). Assumptions (SAEM2) and (SAEM3) depend on the regularity of the model. They are satisfied for the approximate Morris-Lecar model.

In practice, the SAEM algorithm is implemented with an increasing number equal to the iteration number, which satisfies Assumption (SMC1). Assumption (SMC2) is satisfied for the approximate Morris-Lecar model because the variables U are bounded between 0 and 1 and the variables V are fixed at their observed values. This would not have been the case with the filter of [START_REF] Del Moral | The Monte-Carlo method for filtering with discrete-time observations[END_REF], which resimulates the variables V at each iteration. Assumption (SMC3) is satisfied if we require that γ is strictly bounded away from zero; γ ≥ > 0.

Properties of the approximate diffusion.

The SAEM-SMC algorithm provides a sequence which converges to the set of stationary points of the log-likelihood ∆ (θ) = log p ∆ (V 0:n ; θ). The following result aims at comparing this likelihood, which corresponds to the Euler approximate model ( 2), with the true likelihood p(V 0:n ; θ). The result is based on the bound of the Euler approximation proved by [START_REF] Gobet | Sharp estimates for the convergence of the density of the Euler scheme in small time[END_REF]. Their result holds under the following assumption (H1) Functions f , b, σ are 2 times differentiable with bounded derivatives with respect to u and v of all orders up to 2.

Let us assume we apply the SAEM algorithm on an approximate model obtained with an Euler scheme of step size δ = ∆/L. Then we have THEOREM 2. Under assumption (H1), there exists a constant C, independent of θ, such that for any θ ∈ Θ, and any vector

V 0:n |p(V 0:n ; θ) -p δ (V 0:n ; θ)| ≤ C 1 L n∆
Proof is given in Appendix D. Assumption (H1) is a strong assumption, which is sufficient and not necessary. It does not hold for the Morris-Lecar model. Different sets of weaker assumptions have been proposed to prove the convergence of the Euler scheme in the strong sense (expectation of the absolute error between the exact and approximated process). The proofs are mainly based on localization arguments, see [START_REF] Kloeden | Convergence of numerical methods for stochastic differential equations in mathematical finance[END_REF] for a review paper. The convergence of the densities has been less studied, and it is beyond the scope of this paper.

Intra-cellular recordings from a turtle motoneuron.

The membrane potential from a spinal motoneuron in segment D10 of an adult red-eared turtle (Trachemys scripta elegans) was recorded while a periodic mechanical stimulus was applied to selected regions of the carapace with a sampling step of 0.1 ms (for details see [START_REF] Berg | Balanced inhibition and excitation drive spike activity in spinal halfcenters[END_REF][START_REF] Berg | Intense synaptic activity enhances temporal resolution in spinal motoneurons[END_REF]). The turtle responds to the stimulus with a reflex movement of a limb known as the scratch reflex, causing an intense synaptic input to the recorded neuron. Due to the time varying stimulus, a model for the complete data set needs to incorporate the time-inhomogeneity, as done in [START_REF] Jahn | Motoneuron membrane potentials follow a time inhomogeneous jump diffusion process[END_REF]. However, in [START_REF] Jahn | Motoneuron membrane potentials follow a time inhomogeneous jump diffusion process[END_REF] only one-dimensional diffusions are considered, and spikes are modeled as single points in time by adding a jump term with statedependent intensity function to the SDE, ignoring the detailed dynamics during spikes. In this paper we aim at estimating parameters during spiking activity by explicit modeling of time-varying conductances. Therefore, we only analyze four traces during on-cycles (following [START_REF] Jahn | Motoneuron membrane potentials follow a time inhomogeneous jump diffusion process[END_REF]) where spikes occur. Furthermore, in these time windows, the input is approximately constant, which is required for the Morris-Lecar model with constant parameters. An example of the analyzed data is plotted in Fig. 2, together with a filtered trace of the unobserved coordinate. 2 Parameter estimates obtained from observations of the membrane potential of a spinal motoneuron of an adult red-eared turtle during 600 ms for three different values of σ.

With V1 = -2.4 mV, V2 = 36 mV, V3 = 4 mV, V4 = 60 mV fixed. First trace.

First the model was fitted with the values of the scaling parameters V 1 -V 4 given in [START_REF] Rinzel | Methods in Neural Modeling, chapter Analysis of neural excitability and oscillations[END_REF] and used in Section 6 below, see Table 1 for one of the traces. Most of the estimates are reasonable and in agreement with the expected order of magnitudes for the parameter values, except for the V Ca reversal potential, which in the literature is reported to be around 100-150 mV (estimated to 44.7 mV), and the leak conductance, which is estimated to be negative. Conductances are always non-negative. This is probably due to wrong choices of the scaling constants V 1 -V 4 . For the parameters given in [START_REF] Rinzel | Methods in Neural Modeling, chapter Analysis of neural excitability and oscillations[END_REF], the average of the membrane potential V t between spikes is around -26 mV, whereas the average of the experimental trace between spikes is around -56 mV, a factor two larger. We therefore rerun the estimation procedure fixing V 1 -V 4 to twice the value from before, which provides approximately the same values of the normalized Ca 2+ conductance, m ∞ (•), and the rates of opening and closing of K + ion channels, α(•) and β(•), as in the theoretical model when V t is at its equilibrium value. In this case all parameters are reasonable and in agreement with the expected order of magnitudes. To check the robustness to misspecifications in the diffusion parameter σ of the unobserved coordinate, we fitted the model for three different values of σ, see Table 2. Results are stable and suggest that σ is primarily affecting the subthreshold fluctuations of the channel dynamics, and mainly the spiking dynamics of the unobserved coordinate influences the first coordinate.

Final results for all four traces are presented in Table 3. It is reassuring that the parameter estimates seem so reproducible over different traces; the largest variation was below 10%. This is not due to starting values, e.g. the starting value for g L was 0.1, and all four estimates ended up between 1.3 and 1.4, and the starting value for V K was -55, and all four estimates ended up between -75.4 and -74.4.

Simulation study.

Parameter values of the Morris-Lecar model used in the simulations are the same as those of [START_REF] Rinzel | Methods in Neural Modeling, chapter Analysis of neural excitability and oscillations[END_REF][START_REF] Tateno | Random dynamics of the Morris-Lecar neural model[END_REF] for a class II membrane, except that we set the membrane capacitance constant to C = 1 µF/cm 2 , which is the standard value reported in the literature. Conductances and input current were correspondingly changed, and thus, the two models are the same. The values are:

V K = -84 mV, V L = -60 mV, V Ca = 120 mV, C = 1 µF/cm 2 , g L = 0.1 µS/cm 2 , g Ca = 0.22 µS/cm 2 , g K = 0.4 µS/cm 2 , V 1 = -1.2 mV, V 2 = 18 mV, V 3 = 2 mV, V 4 = 30 mV, φ = 0.04 ms -1 .
Input is chosen to be I = 4.5 µA/cm 2 . Initial conditions of the Morris-Lecar model are V t 0 = -26 mV, U t 0 = 0.2. The volatility parameters are γ = 1 mV ms -1/2 , σ = 0.03 ms -1/2 . Tra- jectories are simulated with time step δ = 0.01 ms and n = 2000 points are subsampled with observations time step ∆ = 10δ. Then θ is estimated on each simulated trajectory. A hundred repetitions are used to evaluate the performance of the estimators. An example of a simulated trajectory (for n = 10000) is given in Figure 1.

6.1. Filtering results. The Particle filter aims at filtering the hidden process (U t ) from the observed process (V t ). We illustrate its performance on a simulated trajectory, with θ fixed at its true value. The SMC Particle filter algorithm is implemented with K = 100 particles and the transition density as proposal, see Figure 3. The true hidden process, the mean filtered signal and its 95% confidence interval are plotted. The filtered process appears satisfactory. The confidence interval includes the true hidden process (U t ).

Estimation results.

The performance of the SAEM-SMC algorithm is illustrated on 100 simulated trajectories. The SAEM algorithm is implemented with m = 200 iterations and a sequence (a m ) equal to 1 during the 100 first iterations and equal to a m = 1/(m -100) 0.8 for m > 100. The SMC algorithm is implemented with K(m) = min(m, 100) particles at each iteration of the SAEM algorithm. The SAEM algorithm is initialized by a random draw of θ 0 not centered around the true value: θ 0 = θ true + 0.1 + θ true /3 N (0, 1).

An example of the convergence of the SAEM algorithm for one of the iterations is presented in Fig. 4. It is seen that the algorithm converges for most of the parameters in few iterations to a neighborhood of the true value, even if the initial values are far from the true ones. Only for φ more iterations are needed, which is expected since this parameter appears in the second, non-observed coordinate.

boxplots of 100 estimates A B C 0.20 0.25 q q q q q q q q q q gCa A B C 100 120 140 160 q q q q q q q q q q q q q q q q VCa A B C 0.5 1.0 q q q q q q q q q q q q q q q q q q gK A B C -100 -80 -60 q q q q q q q q q q q q q qq q q q q q q q q VK A B C 0.00 0.05 0.10 0.15 q q q q q q q q q q q q q q gL A B C 2 4 6 q q q q q q q q q q q I A B C 0.02 0.04 0.06 q q q q q q q q q phi A B C 0.96 0.98 1.00 1.02 q q q q gamma SAEM iteration The SAEM estimator is compared with the pseudo maximum likelihood estimator obtained if both V t and U t were observed. Results are given in Table 4. The parameters are well estimated in this ideal case. The estimation of φ, which is the only parameter in the drift of the hidden coordinate U t , is good and does not deteriorate the estimation of the other parameters. In Fig. 4 we show boxplots of the estimates of the eight parameters for the three estimation settings; both coordinates observed, or only one observed with σ fixed at either the true or a wrong value. All parameters appear As expected, the variance of the estimator of φ hugely increases when only one coordinate is observed, but interestingly, the variance of the parameters of the observed coordinate do not seem much affected by this loss of information.

The SAEM-SMC algorithm provides estimates of the standard errors (SE) of the estimators (see Appendix C). These should be close to the RMSE obtained from the 100 simulated datasets. As an example, the SE for one dataset estimated by SAEM are reported in the last line of Table 4. The estimated SE are satisfactory for most of the parameters, but tends to underestimate.

Discussion.

The main contribution of this paper is an algorithm to handle a more general model than a HMM, and show non-asymptotic convergence results for the method. It out that some of the common problems encountered with particle filters is not present in our case, namely the filter does not degenerate, and we run the algorithm on large data sets of 6000 observations points in reasonable time (35 minutes on a standard portable computer for one of the simulated data sets).

To the authors' knowledge, this is the first time the rate parameter of the unobserved coordinate, φ, is estimated from experimental data.

It is comforting to observe that the estimated value do not seem to be very sensitive to the choice of scaling parameters. Other parameters, like the conductances and the reversal potentials, are more sensitive to this choice, and should be interpreted with care.

The estimation procedure builds on the pseudo likelihood, which approximates the true likelihood by an Euler scheme. This approximation is only valid for small sampling step, i.e. for high frequency data, which is the case for the type of neuronal data considered here. If data were sampled less often, a possibility could be to simulate diffusion bridges between the observed points, and apply the estimation procedure to an augmented data set consisting of the observed data and the imputed values.

There are several issues that deserve further study. First, it is important to understand the influence of the scaling parameters V 1 -V 4 , and how to estimate them for a given data set. The model is not exponential in these parameters (assumption (M1)) and new estimation procedures have to be considered. Secondly, one should be aware of the possible misspecification of the model. More detailed models incorporating further types of ion channels could be explored, but increasing the model complexity might deteriorate the estimates, since the information contained in only observing the membrane potential is limited. Furthermore, the sensitivity on the choice of tuning parameters of the algorithm, like the decreasing sequence of the stochastic approximation, (a m ), and the number of SAEM iterations, needs further investigation. Finally, an automated procedure to find starting values for the procedure is warranted. is the sufficient statistic vector corresponding to the parameters ν 1 (θ) = (g L , g Ca , g K , g K V K , g L V L + I, g Ca V Ca ), where denotes transposition.

The sufficient statistics corresponding to ν 2 (θ) = 1/γ 2 are n i=1

(V i -V i-1 ) U i-1 , n i=1 U 2 i-1 , n i=1 (V i -V i-1 ) V i-1 m ∞ (V i-1 ), n i=1 (V i -V i-1 ) U i-1 V i-1 , n i=1 U 2 i-1 V 2 i-1 .
The sufficient statistics corresponding to φ is also explicit but more complex and not detailed here.

APPENDIX C: FISHER INFORMATION MATRIX

The standard errors (SE) of the parameter estimators can be evaluated from the diagonal elements of the inverse of the Fisher information matrix estimate. Its evaluation is difficult because it has no analytic form. We adapt the estimation of the Fisher information matrix, proposed by [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF] and based on the Louis missing information principle.

The Hessian of the log-likelihood ∆ (θ) can be expressed as:

∂ 2 θ ∆ (θ) = E ∂ 2 θ L(S(V 0:n , U 0:n ), θ)|V 0:n , θ + E ∂ θ L(S(V 0:n , U 0:n ), θ) (∂ θ L(S(V 0:n , U 0:n ), θ)) |V 0:n , θ -E [∂ θ L(S(V 0:n , U 0:n ), θ)|V 0:n , θ] E [∂ θ L(S(V 0:n , U 0:n ), θ)|V 0:n , θ] .
where denotes transposition. The derivatives ∂ θ L(S(V 0:n , U 0:n ), θ) and ∂ 2 θ L(S(V 0:n , U 0:n ), θ) are explicit for the Euler approximation of the Morris-Lecar model. Therefore we implement their estimation using the stochastic approximation procedure of the SAEM algorithm. At the mth iteration of the algorithm, we evaluate the three following quantities:

G m+1 = G m + a m ∂ θ L(S(V 0:n , U (m) 0:n ), θ) -G m H m+1 = H m + a m ∂ 2 θ L(S(V 0:n , U (m) 0:n ), θ) + ∂ θ L(S(V 0:n , U (m) 0:n ), θ) (∂ θ L(S(V 0:n , U (m) 0:n ), θ)) -H m F m+1 = H m+1 -G m+1 (G m+1 ) .
As the sequence ( θm ) m converges to the maximum of the likelihood, the sequence (F m ) m converges to the Fisher information matrix.

APPENDIX D: PROOF OF THE CONVERGENCE RESULTS

D.1. Convergence results of Proposition 1.

We omit θ in the proof for clarity. The conditional expectation π n f is given by (5) and the kernels H i from R into itself are defined by (4). We write ν n = µH 1 • • • H n 1 for the constant conditioned on the observed values V 0:n . Also (4) is bounded, i.e. H i 1(u) ≤ C for all u ∈ R and i = 1, . . . , n, for some constant C. It directly follows that

µH 1 • • • H i-1 1 ≤ C i-1 .
Furthermore, we obtain the bound

µH 1 • • • H i 1 ≥ µH 1 • • • H i+1 1 C ≥ • • • ≥ ν n C n-i .
Using the above bounds and that π i-1 is a transition measure, we obtain

ν n C n-1 ≤ π i-1 H i 1 ≤ C. (10)
Define the two empirical measures obtained at time i:

Ψ K i = 1 K K k=1 1 U (k) 0:i and Ψ K i = K k=1 W i (U (k) 0:i )1 U (k) 0:i
. We also decompose the weights and write

Υ K i f = 1 K K k=1 f (U (k) i )w i (U (k) 0:i ). Then W i (U (k) 0:i ) = w i (U (k) 0:i )/(KΥ K i 1) and Ψ K i f = Υ K i f /Υ K i 1.
Recall the following general result [START_REF] Del Moral | The Monte-Carlo method for filtering with discrete-time observations[END_REF] for ξ 1 , . . . , ξ K random variables, which conditioned on a σ-field G are independent, centered and bounded |ξ k | ≤ a. Then for any ε > 0 we have

P 1 K K k=1 ξ k ≥ ε ≤ 2 exp -K ε 2 2a 2 . ( 11
)
Let f be a bounded function on R. Then under assumption (SMC3)

Ψ K i f -Ψ K i f = 1 K K k=1 f (U (k) i ) -Ψ K i f = 1 K K k=1 ξ k fulfills the conditions for (11) to hold with a = 2 f , since E(f (U (k) i )|G) = Ψ K i f
, where G is the σ-algebra generated by U (k) 0:i . Thus, for any ε > 0,

P Ψ K i f -Ψ K i f ≥ ε ≤ 2 exp -K ε 2 8 f 2 . (12) Define Q i (f )(u) = q(u |V i , V i-1 , u)f (u )du . By definition of the unnor- malized weights in step 3 of the SMC algorithm, w i (u, u ) = p ∆ (V i , V i-1 , u, u )/ p ∆ (V i-1 , u)q(u |V i , V i-1 , u), so that Q i (f w i )(u) = p ∆ (V i , u |V i-1 , u)f (u )du = H i f (u). We therefore have Υ K i f -Ψ K i-1 H i f = 1 K K k=1 f (U (k) i )w i (U (k) 0:i ) -Q i (f w i )(U (k) i-1 ) = 1 K K k=1 ξ k
which fulfills the conditions for [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF] to hold, now with a = 2C f and G is the σ-algebra generated by U (k)

0:i-1 , since U (k) i is drawn from q(•|V i-1:i , U (k)
i-1 ), see step 2 of the SMC algorithm. Hence, for any ε > 0 we obtain

P Υ K i f -Ψ K i-1 H i f ≥ ε ≤ 2 exp -K ε 2 8C 2 f 2 . ( 13
)
We want to show the following two bounds

P Ψ K i f -π i f ≥ ε ≤ 2I i exp -K ε 2 8J i f 2 , i = 1, . . . , n (14)
P Ψ K i f -π i f ≥ ε ≤ 2I i exp -K ε 2 8J i f 2 , i = 0, 1, . . . , n (15)
by induction on i, for some constants I i , I i , J i , J i increasing with i to be computed later. Note first that since π 0 = µ and U (k) 0 are i.i.d. with law µ, then [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF] with ξ k = f (U (k) 0 )-µ(f ) yields [START_REF] Doucet | On sequential Monte Carlo sampling methods for bayesian filterin[END_REF] for i = 0 with I i = J i = 1. Let i ≥ 1 and assume (15) holds for i -1. We can write

Ψ K i f -π i f = 1 π i-1 H i 1 Υ K i f Υ K i 1 (π i-1 H 1 1 -Υ K i 1) + (Υ K i f -π i-1 H i f ) .
Note that Υ K i 1 > 0 because the weights w i are strictly positive. Define L i f = Υ K i f -π i-1 H i f and use that |Υ K i f | ≤ f Υ K i 1 (because f is bounded) and [START_REF] Del Moral | The Monte-Carlo method for filtering with discrete-time observations[END_REF] to see that

|Ψ K i f -π i f | ≤ C n-1 ν n ( f |L i 1| + |L i f |)
and

|L i f | ≤ |Υ K i f -Ψ K i-1 H i f | + |Ψ K i-1 H i f -π i-1 H i f |.
Assuming that (15) holds for i -1 and using [START_REF] Ditlevsen | The Morris-Lecar neuron model embeds a leaky integrate-and-fire model[END_REF] and that H i f ≤ C f yield

P (|L i f | ≥ ε) ≤ 2 exp -K ε 2 32C 2 f 2 + 2I i-1 exp -K ε 2 32J i-1 C 2 f 2 .
We obtain

P Ψ K i f -π i f ≥ ε ≤ P |L i 1| ≥ εν n 2C n-1 f + P |L i f | ≥ εν n 2C n-1 ≤ 4 exp -K ε 2 ν 2 n 128C 2n f 2 + 4I i-1 exp -K ε 2 ν 2 n 128J i-1 C 2n f 2
Hence, ( 14) holds with I i ≥ 2(1 + I i-1 ) and J i ≥ 16C 2n J i-1 /ν 2 n ≥ 16J i-1 since ν n ≤ C n . By [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] and ( 14) we then conclude that (15) also holds for i if I i = 1 + I i and J i = 4J i . These conditions are fulfilled by choosing I i = 3 i+1 -3 and J i = 16 i . Thus, (6) holds with C 1 = 6(3 n -1) and C 2 = 8 • 16 n . This concludes the proof.

D.2. Proof of Theorem 1.

To prove the convergence of the SAEM-SMC algorithm, we study the stochastic approximation scheme used during the SA step. The scheme [START_REF] Borg-Graham | Visual input evokes transient and strong shunting inhibition in visual cortical neurons[END_REF] can be decomposed into: S the expectation of the sufficient statistic S under the empirical measure obtained with the SMC algorithm with K(m) particles and current value of parameters θ(s m ) at iteration m of the SAEM-SMC algorithm. Following Theorem 2 of [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF] on the convergence of the Robbins-Monro scheme, the convergence of the SAEM-SMC algorithm is ensured if we prove the following assertions:

1. The sequence (s m ) m≥0 takes its values in a compact set. 2. The function V (s) = -∆ ( θ(s)) is such that for all s ∈ S, F (s) = ∂ s V (s), h(s) ≤ 0 and such that the set V ({s, F (s) = 0}) is of zero measure.

lim m→∞ m

=1 a e exists and is finite with probability 1. 4. lim m→∞ r m = 0 with probability 1.

Assertion 1 follows from assumption (SMC2) and by construction of s m in formula [START_REF] Borg-Graham | Visual input evokes transient and strong shunting inhibition in visual cortical neurons[END_REF]. Assertion 2 is proved by Lemma 2 of [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF] under assumptions (M1)-(M5) and (SAEM2). Assertion 3 is proved similarly as Theorem 5 of [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF]. By construction of the SMC algorithm, the equivalent of assumption (SAEM3) is checked for the expectation taken under the approximate empirical measure Ψ K(m) n; θm . Indeed, the assumption of independence of the non-observed variables U To verify assertion 4, we use Proposition 1. Under assumptions (SMC2)-(SMC3) and assertion 1, Proposition 1 yields that for any ε > 0, there exist two constants C 1 , C 2 , independent of θ, such that

≤ n i=1 p (V i , U i |V i-1 , U i-1 ; θ) - n i=1 p δ (V i , U i |V i-1 , U i-1 ; θ) dU 0:n ≤ n i=1 |p (V i , U i |V i-1 , U i-1 ; θ) -p δ (V i , U i |V i-1 , U i-1 ; θ)| i-1 j=1 p (V j , U j |V j-1 , U j-1 ; θ) n j=i+1
p δ (V j , U j |V j-1 , U j-1 ; θ) dU 0:n [START_REF] Gobet | Sharp estimates for the convergence of the density of the Euler scheme in small time[END_REF] provide that under assumption (H1), there exist constants C 1 > 0, C 2 > 0, C 3 > 0, C 4 > 0 independent of θ such that

|p δ (V i , U i |V i-1 , U i-1 ; θ) + p(V i , U i |V i-1 , U i-1 ; θ)| ≤ C 1 e -C 2 (V i ,U i )-(V i-1 ,U i-1 ) 2 |p δ (V i , U i |V i-1 , U i-1 ; θ) -p(V i , U i |V i-1 , U i-1 ; θ)| ≤ δ C 3 e -C 4 (V i ,U i )-(V i-1 ,U i-1 ) 2
We deduce that for all i = 1, . . . , n, there exists a constant C > 0 independent of θ such that |p (V i , U i |V i-1 , U i-1 ; θ) -p δ (V i , U i |V i-1 , U i-1 ; θ)| i-1 j=1 p (V j , U j |V j-1 , U j-1 ; θ) × n j=i+1 p δ (V j , U j |V j-1 , U j-1 ; θ) dU 0:n ≤ Cδ Finally, we get |p(V 0:n ; θ) -p δ (V 0:n ; θ)| ≤ Cnδ = C 1 L n∆.

-Lecar model. 2 . 1 .

 21 Exact diffusion model. The stochastic Morris-Lecar model including both current and channel noise is defined as the solution to[START_REF] Aït-Sahalia | Maximum likelihood estimation of discretely sampled diffusions: a closed-form approximation approach[END_REF]

Algorithm 2 (

 2 SAEM-SMC algorithm) • Iteration 0: initialization of θ 0 and set s 0 = 0. • Iteration m ≥ 1: S-Step: simulation of the non-observed data (U (m) 0:n ) with SMC targeting the distribution p ∆ (U 0:n |V 0:n ; θ m-1 ).

( M2 )

 M2 The functions ψ(θ) and ν(θ) are twice continuously differentiable on Θ. (M3) The function s : Θ -→ S defined by s(θ) = S(v, u)p ∆ (u|v; θ)dv du is continuously differentiable on Θ. (M4) The function ∆ (θ) = log p ∆ (v, u, θ) is continuously differentiable on Θ and ∂ θ p ∆ (v, u; θ)dv du = ∂ θ p ∆ (v, u; θ)dv du. (M5) Define L : S × Θ → R by L(s, θ) = -ψ(θ) + s, ν(θ) .There exists a function θ : S → Θ such that ∀θ ∈ Θ, ∀s ∈ S, L(s, θ(s)) ≥ L(s, θ). (SAEM1) The positive decreasing sequence of the stochastic approximation (a m ) m≥1 is such that m a m = ∞ and m a 2 m < ∞. (SAEM2) ∆ : Θ → R and θ : S → Θ are d times differentiable, where d is the dimension of S(v, u). (SAEM3) For all θ ∈ Θ, ||S(v, u)|| 2 p ∆ (u|v; θ)du < ∞ and the function Γ(θ) = Cov θ (S(•, U 0:n )) is continuous, where the covariance is under the conditional distribution p ∆ (U 0:n |V 0:n ; θ). (SAEM4) Let {F m } be the increasing family of σ-algebras generated by the random variables s 0 , U . For any positive Borel function f , E ∆ (f (U

( SMC1 )THEOREM 1 .

 SMC11 The number of particles K used at each iteration of the SAEM algorithm varies along the iteration: there exists a function g(m) → ∞ when m → ∞ such that K(m) ≥ g(m) log(m).(SMC2) The function S is bounded uniformly in u. (SMC3) The functions p ∆ (V i |U i , V i-1 , U i-1 ; θ) are bounded uniformly in θ. Assume that (M1)-(M5), (SAEM1)-(SAEM3), and (SMC1)-(SMC3) hold. Then, with probability 1, lim m→∞ d( θ m , L) = 0 where L = {θ ∈ Θ, ∂ θ ∆ (θ) = 0} is the set of stationary points of the loglikelihood ∆ (θ) = log p ∆ (V 0:n ; θ).

FIG 2 .

 2 FIG 2. Observations of the membrane potential in a spinal motoneuron of an adult redeared turtle during 600 ms (upper panel), and the filtered hidden process of the normalized conductance associated with K + current (lower panel) for the estimated parameters with the scaling parameters fixed at V1 = -2.4 mV, V2 = 36 mV, V3 = 4 mV and V4 = 60 mV.

FIG 3 .

 3 FIG 3. Filtering of (Ut) with the particle filter algorithm (100 particles): hidden simulated trajectory of the Morris-Lecar model (Ut) (black), mean filtered signal (grey full drawn line), 95% confidence interval of filtered signal (grey dashed lines).

FIG 4 .

 4 FIG 4. Estimation results for simulated data. Left panels: Boxplots of 100 estimates from simulated data sets for the 8 parameters. True values used in the simulations are given by the gray lines. A: Both Vt and Ut are observed. B: Only Vt is observed, σ is fixed at the true value 0.03. C: Only Vt is observed, σ is fixed at a wrong value 0.04. Right panels: Convergence of the SAEM algorithm for the 8 estimated parameters on a simulated data set. True values used in the simulation are given by the gray lines.

sS

  m+1 = s m + a m h(s m ) + a m e m + a m r m with h(s m ) = π n, θ(sm) S -s m e m = S(V 0:n , U -π n, θ(sm) Swhere we denote by π n,θ S = E ∆ (S(V 0:n , U 0:n )|V 0:n ; θ) the expectation of the sufficient statistic S under the exact distribution p ∆ (U 0:n |V 0:n ; θ), and by Ψ K(m) n, θ(sm)

  given θ 0 , . . . , θ m is verified. As a consequence, for any positive Borel function f , f . Then m =1 a e is a martingale, bounded in L 2 under assumptions (M5) and (SAEM1)-(SAEM2).

Pε 2 C 2 S 2 . 1 M m=1 1 m

 22211 (|r m | > ε) = M m=1 P Ψ K(m) n, θ(sm) S -π n, θ(sm) S ≥ εFinally, assumptions (SMC1)-(SMC2) imply that there exists a constant C 3 , independent of θ, such thatM m=1 P (|r m | > ε) ≤ C C 3 g(m)ε 2which is finite when M → ∞, proving the a.s. convergence of r m to 0.D.3. Proof of Theorem 2. The Markov property yields|p(V 0:n ; θ) -p δ (V 0:n ; θ)| ≤ |p(V 0:n , U 0:n ; θ) -p δ (V 0:n , U 0:n ; θ)| dU 0:n

  nate and we do not consider measurement noise. We propose to impute the hidden non-Markovian path in the stochastic EM algorithm with a Sequential Monte Carlo (SMC) algorithm. Monte-Carlo methods for non-linear filtering are widely spread, with, among other algorithms, sequential importance sampling, sequential importance sampling with resampling (SISR), auxiliary SISR, and stratified resampling [see 8, for a general presentation]. All SISR algorithms are now called SMC. Most of them are designed for HMM. In the specific setting of multi-dimensional SDEs,[START_REF] Del Moral | The Monte-Carlo method for filtering with discrete-time observations[END_REF] proposes a particle filter for a two-dimensional SDE, where the second equation is autonomous. Although the first coordinate is observed at discrete times, they propose to simulate it at each iteration of the filter.

  The functions α(•) and β(•) model the rates of opening and closing of the K + ion channels. The function m ∞ (•) represents the equilibrium value of the normalized Ca 2+ conductance for a given value of the membrane potential. The parameters V 1 , V 2 , V 3 and V 4 are scaling parameters; g Ca , g K and g L are conductances associated with Ca 2+ , K + and leak currents; V Ca , V K and V L are reversal potentials for Ca 2+ , K + and leak currents; C is the membrane capacitance; φ is a rate scaling parameter for the opening and closing of the K + ion channels; and I is the input current.

		membrane voltage, V(t)		state space
	0.2 0.4 -50 0	normalized conductance, U(t)	U t	0.1 0.2 0.3 0.4 0.5
	0	500	1000	-40 -20 0	20 40
		time		V t
	FIG 1. Simulated trajectory of the stochastic Morris-Lecar model: (Vt) as a function of
	time (left, top), (Ut) as a function of time (left, bottom), and (Ut) against (Vt) (right). Pa-
	rameters are given in Section 6. Time is measured in ms, voltage in mV, the conductance
	is normalized between 0 and 1.

). Processes ( Bt ) t≥t 0 and (B t ) t≥t 0 are independent Brownian motions. The variable V t represents the membrane potential of the neuron at time t, and U t represents the normalized conductance of the K + current. It varies between 0 and 1, and can be interpreted as the probability that a K + ion channel is open at time t. The equation for f (•) describing the dynamics of V t contains four terms, corresponding to Ca 2+ current, K + current, a general leak current, and the input current I.

  Parameter estimates obtained from observations of the membrane potential of a spinal motoneuron of an adult red-eared turtle during 600 ms for two different sets of scaling parameters. With σ = 0.05 fixed. First trace.

	Parameter	gL	gCa	gK	γ	VK	φ	VCa	I
	With V1 = -1.2 mV, V2 = 18 mV, V3 = 2 mV, V4 = 30 mV
	Estimate	-0.296 11.274	6.553 2.801 -124.481 1.989 35.769	-5.024
	SE	0.001	0.028	0.049 0.001	14.563 0.000	0.122	0.052
	With V1 = -2.4 mV, V2 = 36 mV, V3 = 4 mV, V4 = 60 mV
	Estimate	1.046 12.906 20.878 2.466	-67.097 2.153 98.698 -65.403
	SE	0.009	0.008	0.021 0.001	0.227 0.001	0.198	1.204
				TABLE 1
	Parameter	gL	gCa	gK	γ	VK	φ	VCa	I
	σ fixed to 0.02
	Estimate	1.302 12.460 16.550 2.280 -74.301 2.881 99.398 -52.742
	SE	0.003	0.028	0.131 0.000	0.703 0.001	2.384	2.727
	σ fixed to 0.05
	Estimate	1.046 12.906 20.878 2.466 -67.097 2.153 98.698 -65.403
	SE	0.009	0.008	0.021 0.001	0.227 0.001	0.198	1.204
	σ fixed to 0.15
	Estimate	1.308 12.442 16.419 2.301 -74.576 2.911 99.417 -52.341
	SE	0.002	0.007	0.120 0.000	0.021 0.000	4.529	1.558
				TABLE

  Parameter estimates obtained from four different traces of the membrane potential of a spinal motoneuron of an adult red-eared turtle. Each trace has 6000 observations points with a sampling step of 0.1 ms. With V1 = -2.4 mV, V2 = 36 mV, V3 = 4 mV, V4 = 60 mV and σ = 0.05 fixed.

	Parameter	gL	gCa	gK	γ	VK	φ	VCa	I
	First trace
	Estimate	1.046 12.906 20.878 2.466 -67.0.97 2.153	98.698 -65.403
	SE	0.009	0.008	0.021 0.001	0.227 0.001	0.198	1.204
	Second trace
	Estimate	1.430 11.705 15.791 2.253	-73.029 3.138 103.709 -53.183
	SE	0.003	0.008	0.029 0.000	0.641 0.001	1.269	0.651
	Third trace
	Estimate	1.371 11.878 15.379 2.210	-75.024 3.004	99.887 -49.499
	SE	0.002	0.013	0.017 0.000	0.195 0.000	0.464	0.614
	Fourth trace
	Estimate	1.197 11.452 12.521 2.012	-85.982 3.776	99.615 -37.017
	SE	0.002	0.055	0.017 0.000	0.089 0.000	1.466	0.861
				TABLE 3

  Simulation results obtained from 100 simulated Morris-Lecar trajectories (n = 2000, ∆ = 0.1 ms). Two estimators are compared: The pseudo maximum likelihood estimator in the ideal case where both Vt and Ut are observed; and the SAEM estimator when only Vt is observed with the SAEM initialization at a random value not centered around the true value θ. An example of standard errors (SE) estimated with the SAEM-SMC algorithm on one single simulated dataset is also given.

					Parameters
	Estimator	gL	gCa	gK	γ	VK	φ	VCa	I
	true values 0.100 0.220 0.400	1.00	-84.00 0.040	120.00 4.400
	With both Vt and Ut observed (pseudo maximum likelihood estimator)
	mean	0.101 0.219 0.411 0.996	-83.20 0.040	121.97 4.539
	RMSE	0.017 0.019 0.041 0.019	7.61 0.001	8.50 0.560
	With only Vt observed (SAEM estimator)
	mean	0.090 0.225 0.464 1.003 -78.622 0.041 119.677 4.060
	RMSE	0.021 0.024 0.144 0.017	9.459 0.013	10.218 1.028
	estimated
	SE	0.016 0.019 0.042 0.016	4.96 0.001	7.31 0.561
					TABLE 4

Acknowledgments. The authors are grateful to Rune W. Berg for making his experimental data available. We thank E. Gobet for helpful discussions about convergence of the Euler scheme. The work is part of the Dynamical Systems Interdisciplinary Network, University of Copenhagen.

* S. Ditlevsen is supported by grants from the Danish Council for Independent Research | Natural Sciences.

† A. Samson is supported by Grants from the University Paris Descartes PCI.

APPENDIX A: DISTRIBUTIONS OF APPROXIMATE MODEL

Consider the general approximate model (see [START_REF] Andrieu | Sequential Monte Carlo methods for optimal filtering[END_REF])

with ρ the correlation coefficient between the two Brownian motions or perturbations. The distribution of (V i+1 , U i+1 ) conditionally on

The marginal distributions of V i+1 conditionally on

The conditional distributions of V i+1 conditionally on (U i+1 , V i , U i ) and

where

The distributions in ( 8) and ( 9) are equal when the Brownian motions are independent, i.e. when ρ = 0.

APPENDIX B: SUFFICIENT STATISTICS

We here provide the sufficient statistics of the approximate model (2). Consider the n × 6-matrix X = -V 0:(n-1) , -m ∞ (V 0:(n-1) )V 0:(n-1) , -U 0:(n-1) V 0:(n-1) , U 0:(n-1) , 1, m ∞ (V 0:(n-1) ) where 1 is the vector of 1's of size n. Then the vector S 1 (V 0:(n-1) , U 0:(n-1) ) = (X X) -1 X (V 1:n -V 0:(n-1) )