Sequential estimation of intrinsic activity and synaptic input in single neurons by particle filtering with optimal importance density

This paper deals with the problem of inferring the signals and parameters that cause neural activity to occur. The ultimate challenge being to unveil brain's connectivity, here we focus on a microscopic vision of the problem, where single neurons (potentially connected to a network of peers) are at the core of our study. The sole observation available are noisy, sampled voltage traces obtained from intracellular recordings. We design algorithms and inference methods using the tools provided by stochastic filtering, that allow a probabilistic interpretation and treatment of the problem. Using particle filtering we are able to reconstruct traces of voltages and estimate the time course of auxiliary variables. By extending the algorithm, through PMCMC methodology, we are able to estimate hidden physiological parameters as well, like intrinsic conductances or reversal potentials. Last, but not least, the method is applied to estimate synaptic conductances arriving at a target cell, thus reconstructing the synaptic excitatory/inhibitory input traces. Notably, these estimations have a bound-achieving performance even in spiking regimes.


I. INTRODUCTION
N EUROSCIENCE is the science that delves into the understanding of the nervous system. It is one of the most interdisciplinary sciences, gathering together experts from a vast variety of fields of knowledge. Neuroscience is a rather broad discipline and encompasses many aspects related to the Central Nervous System (CNS). The different topics in neuroscience can be studied from various perspectives depending on the prism used to focus the problem. This ranges from understanding the internal mechanisms that cause spiking of a single cell (a neuron), to explaining the dynamics occurring in populations of neurons that are interconnected. Even more macroscopically, one could consider the analysis of functional parts in the brain which are ultimately composed of individual neurons grouped together. In this work, we are interested in the microscopic view of the problem.
Measurements of membrane potential traces constitute the main observable quantities to derive a biophysical neuron model. In particular, the dynamics of auxiliary variables and the model parameters are inferred from voltage traces, in a costly process that typically entails a variety of channel blocks and clamping techniques (see for instance [1]), as well as some uncertainty in the parameter values due to noise in the signal. Recent works in the literature deal with the problem of inferring hidden parameters of the model, see for instance [2]- [4] and, for an exhaustive review, [5]. In the same line, it is worth to mention attempts to extract connectivity in networks of neurons from calcium imaging [6].
Apart from inferring intrinsic parameters of the model, voltage traces are also useful to obtain valuable information about synaptic input, an inverse problem with some satisfactory (see for instance [7]- [13]) but no complete solutions yet. The main shortcomings are the requirement of multiple (supposedly identical) trials for some methods to be applied, and the need of avoiding signals obtained when ionic currents are active. The latter constraint arises from the fact that many methods rely on the linearity of the signal and this is not possible to achieve under quite general situations, like spiking regimes (see [14]) or subthreshold regimes when specific currents (e.g., AHP, LTS, etc.) are active (see [15]). An ideal method should be able to sequentially infer the time-course of the membrane potential and its intrinsic/extrinsic activity from noisy observations of a voltage trace. The main features of the envisaged algorithm are fivefold i) Single-trial: the method should be able to estimate the desired signals and parameters from a single voltage trace, thus avoiding the experimental variability among trials; ii) Sequential: the algorithm should provide estimates each time a new observation is recorded, thus avoiding re-processing of all data stream each time; iii) Spike regime: contrary to most solutions operating only under the subthreshold assumption, the method should be able to operate in the presence of spikes as well; iv) Robust: if the method is model-dependent, thus implying knowledge of the model parameters, then the algorithm should be provided with enhancements to adaptively learn these parameters; and v) Efficient: the performance of the method should be close to the theoretical lower bounds, meaning that the estimation error cannot be substantially reduced. Notice that the focus here is not on reducing the computational cost of the inference method and thus we allow ourselves to use resource-consuming algorithms. Indeed, the target application does not demand (at least as a main requirement) real-time operation, and thus we prioritized performance (i.e., estimation accuracy and the rest of features described earlier) in our developments.
According to the above desired features, in this work, that substantially extends our previous contributions [16], [17], we are interested in methods that can provide on-line estimation and avoid the need of repetitions that could be contaminated by neuronal variability. Particularly, we concentrate on methods to extract intrinsic activity of ionic channels, namely the probabilities of opening and closing ionic channels, and the contribution of synaptic conductances. We built a method based on Bayesian theory to sequentially infer these quantities from single-trace, noisy membrane potentials. The material therein includes a discussion of the discrete state-space representation of the problem and the model inaccuracies due to mismodeling effects. We present two sequential inference algorithms: i) a method based on particle filtering (PF) to estimate the time-evolving states of a neuron under the assumption of perfect model knowledge; and ii) an enhanced version where model parameters are jointly estimated, and thus the rather strong assumption of perfect model knowledge is relaxed. We provide exhaustive computer simulation results to validate the algorithms and observe that they are attaining the theoretical lower bounds of accuracy, which are derived in the paper as well.
The results show the validity of the approach and its statistical efficiency. Although we used for convenience a specific neuron model (i.e., the Morris-Lecar) in the computer simulations, the proposed procedure can be applied to any neuron model without loss of generality.
In this paper we use the powerful tools of PF to make inferences in general state-space models. PF are a set of methods able to sample from the marginal filtering distribution in situations where analytical solutions are hard to work out, or simply impossible. In the recent years PFs played an important role in many research areas such as signal detection and demodulation, target tracking, positioning, Bayesian inference, audio processing, financial modeling, computer vision, robotics, control or biology [18]- [24]. At a glance, PF approximate the filtering distribution of states given measurements by a set of random points, properly weighted according to the Bayes' rule. The generation of the random particles can be done through a variety of distributions, known as importance density. Particularly, we formulate the problem at hand and observe that it is possible to use the optimal importance density [25]. This distribution generates particles close to the target distribution and thus it can be shown to reduce the variance of the particles. As a consequence, for a fix number of particles, usage of this approach (not always possible) leads to better accuracy results than other choices. To the author's knowledge, the utilization of such sampling distribution is novel in the context of neural model filtering. Similar works have used PF to track neural dynamics, but with no optimal importance density (see [2], [3]) or to estimate synaptic input from subthreshold recordings [11], as opposite to our proposed approach where we aim at providing estimates during the, highly nonlinear, spike regime. These references use the expectation-maximization algorithm to estimate the model parameters. Lighter filtering methods based on the Gaussian assumption were considered in the literature (see [12], [13], [26] for instance), but the assumption might not hold in general. For instance, due to outliers in the membrane measurements or if more sophisticated models for the synaptic conductances are considered. In these situations, a PF approach seems more appropriate. As mentioned earlier, the focus here is on highly efficient and reliable filtering methods, rather than on computationally light inference methods.
The remainder of the article is organized as follows. Section II provides a brief overview of the main biological concepts regarding neural activity, which are relevant to this work. This section has been included for the sake of completeness and the reader familiar with neuroscience can skip it. In Section III we expose the problem and present the statistical model, essentially a discretization of the well-known Morris-Lecar model, and we analyze the model inaccuracies as well. Next, in Section IV we present the different inference algorithms we apply depending on the knowledge of the system. Results are given in Section V, where we tackle three inference problems: i) when the parameters defining the model are known; ii) when the parameters of the model are unknown, and thus they need to be estimated; and, iii) estimation of synaptic conductances from voltage traces assuming unknown model parameters. Finally, Section VI concludes the paper with final remarks.

II. NEURON ELECTROPHYSIOLOGY IN A NUTSHELL
Neurons are the basic information processing structures in the CNS [27]- [29]. The main function of a neuron is to receive input information from other neurons, to process that information, and to send output information to other neurons. Synapses are connections between neurons, through which they communicate this information. It is controversial how this information is encoded, but it is quite accepted that information produces changes in the electrical activity of the neurons, seen as voltage changes in the membrane potential (i.e., the difference in electrical potential between the interior and the exterior of a biological cell).
The basic constituents of a neuron are the soma, which contains the nucleus of the cell, it is the body of the neuron where most of the information processing is carried; the dendrites, that are extensions of the soma which connect the neuron to neighboring neurons and capture their stimuli; the axon, which is the largest part of a neuron where the information is transmitted in form of an electrical current. A cell might have only one axon or more. The physiological meaning for the propagation of the voltage through the axon can be understood in terms of voltage-gated ionic channels located in the axon membrane; and the synapses, located at the axon terminal are in charge of the electrochemical reactions that cause neuron communications to happen. More precisely, the membrane potential (an electrical phenomenon) traveling through the axon, when reaching the synapse, activates the emission of neurotransmitters (a chemical phenomenon) from the neuron to the receptors of the target neurons. This chemical reaction is transformed again into electrical impulses in the dendrites of the receiving neurons.
We are specially interested in understanding the phenomena through which an electrical voltage travels the axon from the soma to the synapse. The basic idea is that the membrane covering the axon is essentially impermeable to most charged molecules. This makes the axon to act as a capacitor (in terms of electrical circuits) that separates the inner and outer parts of the neuron's axon. This is combined with the so-called ionicchannels, that allow the exchange of intracellular/extracellular ions through electrochemical gradients. This exchange of ions is responsible for the generation of an electrical pulse called action potential, that travels along the neuron's axon. Ionicchannels are found throughout the axon and are typically voltage-dependent, which is primarily how the action potential propagates.
The most common ionic species involved in the generation of the action potential are sodium (Na + ), potassium (K + ), chloride (Cl − ), and calcium (Ca 2+ ). For each ionic species, the corresponding ionic-channel aims at balancing the concentration and electrical potential gradients, which are opposite forces regulating the exchange of ions through the gate. The point at which both forces counterbalance is known as Nernst equilibrium potential.
For the sake of simplicity and without loss of generality, we consider the evolution of the membrane potential at a specific site of the axon. Therefore, v v(t) denotes the continuoustime membrane potential at a point in the axon. Accounting that the membrane potential is seen as a capacitor, the currentvoltage relation allows us to express the total current flowing in the membrane as proportional to the time derivative of the voltage. Then, we obtain the so-called conductance-based mathematical model for the evolution of membrane potentials, where C m is the membrane capacitance and I app represents the externally applied currents, for instance injected via an electrode and used to perform a controlled experiment. The time-varying ionic current for the i-th ionic species, is a constant called the maximal conductance, which is fixed for each ionic species. The variable p i p i (v) is the average proportion of channels of the type i ∈ I in the open state. Notice that the proportion is in general voltage-dependent (i.e., sensitive to the membrane potential) and thus they are said to be voltage-gated. p i can be further classified into gates that activate (i.e., gates that open the i-th ionic channel) and those that inactivate (i.e., gates that close the i-th ionic channel). Mathematically, omitting the dependence on i, p = m a h b , where a is the number of activation gates and 0 < m m i (v) < 1 is the probability of activating gate being in the open state. Similarly, b is the number of inactivation gates and 0 < h h i (v) < 1 the probability of inactivating gate being in the open state. We refer to m and h as the gating variables of the ionic channel. The dynamics of these gating variables are responsible for the membrane potential generation The leakage term in (1),ḡ L (v−E L ), is mathematically used to gather all ionic channels that are not explicitly modeled. The maximal conductance of the leakage,ḡ L , is considered constant and it is adjusted to match the membrane potential at resting state. Similarly, the equilibrium potential E L has to be estimated at rest.
I syn gathers the contribution of neighboring neurons and it is referred to as the synaptic current. The most general model for I syn considers decomposition in 2 independent components: corresponding to excitatory (the most common being α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic (AMPA) neuroreceptors) and inhibitory (the most common being γaminobutyric acid (GABA) neuroreceptors) terms, respectively. Roughly speaking, whereas the excitatory synaptic term makes the postsynaptic neuron more likely to generate a spike, the inhibitory term makes the postsynaptic neuron less likely to generate an action potential. E E and E I are the corresponding reverse potentials, typically close to 0 mV and −80 mV, respectively. A longstanding problem is to characterize the time-varying global excitatory and inhibitory conductances g E (t) and g I (t). There are various mathematical models for the synaptic current that could be considered [30], here we use the so-called effective point-conductance model of [7], [31]. In this model, the excitatory/inhibitory global conductances are treated as Ornstein-Uhlenbeck (OU) processeṡ where u = {E, I}. χ(t) is a zero-mean, white noise, Gaussian process with unit variance. Then, the OU process has mean g u,0 , standard deviation σ u , and time constant τ u . This simple model was shown in [31] to yield a valid description of the synaptic noise, capturing the properties of more complex models.
Conductance-based models like (1) are widely used, mostly varying on the number and type of gating variables and the activation/inactivation functions defining the dynamics of the gating variables. The pioneer work by Hodgkin and Huxley in [32] has been followed by a plethora of alternative models such as [33]- [37] to name a few. Fig. 1 shows the basic setup we are dealing with in this article. The neuron under observation has its own dynamics, producing electrical voltage patterns. The generation of action potentials is regulated by internal drivers (e.g., the active gating variables of the neuron) as well as exogenous factors like excitatory and inhibitory synaptic conductances produced by pools of connected neurons. This system is unobservable, in the sense that we cannot measure it directly. The sole observation from this system are the noisy membrane potentials y k . In this experimental scenario, we aim at applying sequential inference methods to extract the following quantities:  1) The time-evolving states characterizing the neuron dynamics, including a filtered membrane potential and the dynamics of the gating variables.
2) The parameters defining the neuron model. It was seen that conductance-based models require the knowledge (or adjustment) of a number of static parameters. It is desirable to have autonomous inference algorithms that estimate these parameters as well, on the top of the timeevolving states.
3) The dynamics of synaptic conductances and its parameters, the final goal being to discern the temporal contributions of global excitation from those of global inhibition, g E (t) and g I (t) respectively.

III. PROBLEM STATEMENT AND MODEL
The membrane potential, obtained from intracellular recordings, is one of the most valuable signals of neurons' activity. Most of the neuron models have been derived from fine measurements and allow the progress of "in silico" experiments. The recording of the membrane potential is a physical process, including some approximations/inaccuracies involving: 1) Voltage observations are noisy. This is due, in part, to the thermal noise at the sensing device, non-ideal conditions in experimental setups, etc. 2) Recorded observations are discrete. All sensing devices record data by sampling at regular time intervals the continuous-time natural phenomena. This is the task of an Analog-to-Digital Converter (ADC). Moreover, these samples are typically digitized, i.e. expressed by a finite number of bits. This latter issue is not tackled in the work as we assume that modern computer capabilities allow us to sample with relatively large number of bits per sample.
The problem investigated in this paper considers recordings of noisy voltage traces to infer the hidden gating variables of the neuron model, as well as filtered voltage estimates, model parameters and input synaptic conductances. Data is recorded at discrete time-instants k at a sampling frequency f s = 1/T s . The problem can thus be posed in the form of a discrete-time, state-space model. The observations are with v k representing the nominal membrane potential and σ 2 y,k modeling the noise variance due to the sensor or the instrumentation inaccuracies when performing the experiment. To provide comparable results, we define the signal-to-noise ratio (SNR) as SNR = P s /P n , with P s being the average signal power and P n = σ 2 y,k the noise power. The methods presented in this paper rely on models for the evolution of the voltage-traces and the hidden variables of a neuron. For instance, the conductance-based neuron models and the OU processes for the synaptic conductances introduced in Section II.
Concerning neuron models, in this work, we focus on the Morris-Lecar model [35] for the sake of clarity. This model is very popular in computational neuroscience as it is able to model a wide variety of neural dynamics. Details of the Morris-Lecar model can be consulted in Appendix A. The unknown state vector in this case is composed of the membrane potential and the K + gating variable, Notice that the Morris-Lecar neuron model is defined by a system of continuous-time, ordinary differential equations (ODEs) of the formẋ = f (x). In general, mathematical models for neurons are of this type. However, due to the sampled recording of measurements, we are interested in expressing the model in the form of a discrete state-space, where ν k ∼ N (0, Σ x,k ) is the process noise which includes the model inaccuracies. The covariance matrix Σ x,k is used to quantify our confidence in the model that maps In general, obtaining a closedform analytical expression for f k without approximations is not possible. In such a case, we could use the Euler method: where ∆t = T s is the sampling period. Thus, we can write (6) as which is of the Markovian type. If we focus on the Morris-Lecar model, the resulting discrete version of the ODE system in (39)-(40) is: with m ∞ (v k ), n ∞ (v k ), and τ n (v k ) defined as in (44)- (46), see Appendix A. Then, (9) and (10) can be interpreted as The goal is to express the inference problem in state-space formulation and apply stochastic filtering tools learned from signal processing. The final ingredient to do so is to introduce the so-called process noise in the state equation where the noise terms ν v,k and ν n,k are assumed jointly Gaussian with covariance matrix Σ x,k . Further details of this matrix are discussed in Section III-A. This general form of Markovian type is preserved when the model is extended with a couple of OU processes associated to the excitatory and inhibitory synaptic conductances. In this case, the resulting state-space model is composed by the Morris-Lecar model used so far (with the peculiarity that the term −I syn is added to (39)), plus the OU stochastic process in (3) describing I syn . Therefore, the continuous-time state is x = (v, n, g E , g I ) . The discrete version of the state-space is used again.

A. Model inaccuracies
The proposed estimation method relies on the fact that the neuron model is known. This is true to some extent, but most of the parameters in the Morris-Lecar model discussed are to be estimated. Typically this model calibration is done beforehand, but as we will see later in Section IV-B this could be done in parallel to the filtering process. Therefore, the robustness of the method to possible inaccuracies should be assessed. In this section, we point out possible causes of mismodeling. Computer simulations are later used to characterize the performance of the proposed methods under these impairments.
In the single-neuron model considered, three major sources of inaccuracies can be identified. First, the applied current I app can be itself noisy, with a variance depending on the quality of the instrumentation used and the experiment itself. We model the actual applied current as a random variable where I o is the nominal current applied and σ 2 I the variance around this value. Plugging (12) into (9) we obtain that the contribution of I app to the noise term is T s Secondly, when the conductance of the leakage term is estimated beforehand, some inaccuracies might be taken into account. In general, this term is considered constant in the models although it gathers relatively distinct phenomena that can potentially be time-varying. The maximal conductance of the leakage term is therefore inaccurate and modeled as whereḡ o L is the nominal, estimated conductance and σ 2 g the variance of this estimate. Similarly, inserting (13) into (9) we see that the contribution ofḡ L to the noise term is T s , and τ n (v k ) are to be estimated. In general, these parameters can be properly obtained off-line by standard methods, see [28]. However, as they are estimates, a residual error typically remains. To account for these inaccuracies, we consider that the equation governing the evolution of gating variables is corrupted by a zero-mean additive white Gaussian process with variance σ 2 n . This analysis allows us to construct a realistic Σ x,k , as the contribution of the aforementioned inaccuracies. In a practical setup, in order to compute the noise variance due to leakage, we need to use the approximationv k−1 ≈ v k−1 , wherev k−1 is estimated by the filtering method. We construct the covariance matrix of the model as where we used that the overall noise in the voltage model is and as an estimate of σ 2 v , provided accurate knowledge of σ 2 I and σ 2 g . Otherwise, the covariance matrix of the process could be estimated by other means, as the ones presented in Section IV-B for mixed state-parameter estimation in nonlinear filtering problems.

IV. METHODS
Two filtering methods are proposed here, depending on the knowledge regarding the underlying dynamical model. Section IV-A presents an algorithm able to estimate the states in x k by a PF methodology, the particularity being that an optimal distribution is used to draw the random samples characterizing the joint filtering distribution of interest. This method assumes knowledge of the parameter values of the system model, although we account for some inaccuracies as detailed in Section III-A. An enhanced version of this method is presented in Section IV-B, where we relax the assumption of knowing the parameter values. Leveraging on a Particle Markov-Chain Monte-Carlo (PMCMC) algorithm and the use of the optimal importance density as in the first method, we present a method that is able to filter x k while estimating the values describing the neuron model.

A. Sequential estimation of voltage traces and gating variables
The type of problems we are interested in involve the estimation of time-evolving signals that can be expressed through a state-space formulation. Particularly, estimation of the states in a single-neuron model from noisy voltage traces can be readily seen as a filtering problem. Bayesian theory provides the mathematical tools to deal with such problems in a systematic manner. The focus is on sequential methods that can incorporate new available measurements as they are recorded without the need for reprocessing all past data.
Bayesian filtering involves the recursive estimation of states x k ∈ R nx given measurements y k ∈ R at time t based on all available measurements, y 1:k = {y 1 , . . . , y k }. To that aim, we are interested in the filtering distribution p(x k |y 1:k ), which can be recursively expressed as with p(y k |x k ) and p(x k |x k−1 ) referred to as the likelihood and the prior distributions, respectively. Unfortunately, (16) can only be obtained in closed-form in some special cases and in more general setups we should resort to more sophisticated methods. In this paper we consider PF to cope with the nonlinearity of the model. Although other lighter approaches might be possible as well [22], we seek the maximum accuracy regardless the involved computational cost. Theoretically, for sufficiently large number of particles, particle filters offer such performances.
Particle filters, see [18], [20], [21], [23], approximate the filtering distribution by a set of N weighted random samples, forming the random measure x . These random samples are drawn from the importance density distribution, and weighted according to the general formulation The importance density from which particles are drawn is a key issue in designing efficient PFs. It is well-known that the optimal importance density is π(x k |x (i) 0:k−1 , y 1:k ) = p(x k |x (i) k−1 , y k ), in the sense that it minimizes the variance of importance weights. Weights are then computed using (18) as w . This choice requires the ability to draw from p(x k |x (i) k−1 , y k ) and to evaluate p(y k |x (i) k−1 ). In general, the two requirements cannot be met and one needs to resort to suboptimal choices. However, we are able to use the optimal importance density since the state-space model assumed here is Gaussian, with nonlinear process equations but related linearly to observations [25]. The equations are: with and the importance weights can be updated using (1, 0) . The PF provides a discrete approximation of the filtering distribution of the form p(x k |y 1:k , which gather all information from x k that the measurements up to time k provide. The minimum mean square error estimator can be obtained aŝ wherex k = (v k ,n k ) . Recall that the method discussed in this section could be easily adapted to other neuron models by simply substituting the corresponding transition function f k and constructing the state vector x k accordingly. As a final step, PFs incorporate a resampling strategy to avoid collapse of particles into a single state point. Resampling consists in eliminating particles with low importance weights and replicating those in high-probability regions [38]. The overall algorithm can be consulted in Algorithm 1 at instance k. Notice that this version of the algorithm requires knowledge of noise statistics and all the model parameters, which for the Morris-Lecar model are When we add the dynamics of the synaptic conductances, the vector Θ of model parameters also includes τ E , τ I , g E,0 , g I,0 , σ E and σ I .

B. Joint estimation of states and model parameters
In practice the parameters in (24) might not be known. It is reasonable to assume that Θ, or a subset of these parameters θ ⊆ Θ, are unknown and need to be estimated at the same time the filtering method in Algorithm 1 is executed. Therefore, the ultimate goal in this case is to estimate jointly the time evolving states and the unknown parameters of the model, x 1:k and θ respectively. Joint estimation of states and parameters is a longstanding problem in Bayesian filtering, and specially hard to handle in the context of PFs. Refer to [39]- [41] and their references for a complete survey. Here, we follow the approach in [42] and make us of the so-called PMCMC to enhance the presented PF algorithm with parameter estimation capabilities. In the remainder of this section we provide the basic ideas to use the algorithm, following a similar approach as in [24].
Following the Bayesian philosophy we adopt here, the problem fundamentally reduces to assigning an a priori distribution for the unknown parameter θ ∈ R n θ and extending the statespace model (here we adopt its probabilistic representation) and initial state distribution x 0 ∼ p(x 0 |θ). Applying Bayes' rule, the full posterior distribution can be expressed as p(x 0:T , θ|y 1:T ) = p(y 1:T |x 0:T , θ)p(x 0:T |θ)p(θ) p(y 1:T ) with p(y 1: p(x 0: Notice here that we are dealing with a finite horizon of observations T . Then, from a Bayesian perspective, the estimation of θ is equivalent to obtaining its marginal posterior distribution p(θ|y 1:T ) = p(x 0:T , θ|y 1:T )dx 0:T . However, this is in general extremely hard to compute analytically and one needs to find workarounds. Evaluation of the full posterior turns to be not only computationally prohibitive, but useless if states cannot be marginalized out analytically. Alternative methods resort on the factorization of the parameter marginal distribution as p(θ|y 1:T ) = p(y 1:T |θ)p(θ) and how Bayesian filters can be transformed to provide characterizations of the marginal likelihood distribution p(y 1:T |θ) and related quantities. The marginal likelihood distribution can be recursively factorized in terms of the predictive distributions of the observations: p(y 1:T |θ) = T k=1 p(y k |y 1:k−1 , θ), with p(y k |y 1:k−1 , θ) = p(y k |x k , θ)p(x k |y 1:k−1 , θ)dx k obtained straightforwardly as a byproduct of any of the Bayesian filtering methods, see [22].
A useful transformation of the marginal likelihood is the so-called energy function, which is sometimes more convenient to deal with. The energy function is defined as ϕ T (θ) = − ln p(y 1:T |θ) − ln p(θ) or, equivalently, p(θ|y 1:T ) ∝ exp(−ϕ T (θ)). The energy function can then be recursively computed as a function of the predictive distribution Then, the basic problem is to obtain an estimate of the predictive distribution p(y k |y 1:k−1 , θ) from the PF we have designed in Section IV-A and use it in conjunction with p(θ) to infer the marginal distribution p(θ|y 1:T ) of interest. This latter step can be performed in several ways, from which we choose to use the Markov-Chain Monte-Carlo (MCMC) methodology to continue with a fully Bayesian solution. Besides, it is known to be the solution that provides best results when used in a PF. Next, we detail how ϕ k (θ) can be obtained from a PF algorithm, we present the MCMC method for parameter estimation, and finally we sketch the overall algorithm. 1) Computing the energy function from particle filters: The modification needed is rather small. Actually, it is non-invasive in the sense that the algorithm remains the same and the energy function can be computed adding some extra formulae. Recall that the predictive distribution p(y k |y 1:k−1 , θ) is composed of two distributions and that the PF provides characterizations of these two distributions. Then, one could use a particle approximation p(y k |y 1:k−1 , θ) ≈p(y k |y 1: k−1 as in the original algorithm and Then, it is straightforward to identify the energy function approximation as which can be computed recursively in the PF algorithm. 2) The Particle Markov-Chain Monte-Carlo algorithm: Once an approximation of the energy function is available, we can apply MCMC to infer the marginal distribution of the parameters. MCMC methods constitute a general methodology to generate samples recursively from a given distribution by randomly simulating from a Markov chain [43]- [47]. There are many algorithms implementing the MCMC concept, being one of the most popular the Metropolis-Hastings (MH) algorithm. At the j-th iteration, the MH algorithm samples a candidate point θ * from a proposal distribution q(θ * |θ (j−1) ) based on the previous sample θ (j−1) . Starting from an arbitrary value θ (0) , the MH algorithm accepts the new candidate point (meaning that it was generated from the target distribution, p(θ|y 1:T )) using the rule where u is drawn randomly from a uniform distribution, u ∼ U(0, 1), and is referred to as the acceptance probability. It is critical for the performance of the algorithm the choice of the proposal density. A common choice is q(θ|θ (j−1) ) = N (θ; θ (j−1) , Σ (j−1) ) with the selection of the transitional covariance remaining as the tuning Σ (j−1) parameter. This covariance can be adapted as iterations of the MCMC method progress. In this work we have adopted the Robust Adaptive Metropolis (RAM) algorithm [48]. The RAM algorithm can be consulted in Algorithm 2. We use the notation that S = Chol (A) denotes the Cholesky factorization of an Hermitian positive-definite matrix A such that A = SS , and S is a lower triangular matrix [49]. , which can be used to approximate (after neglecting the first samples corresponding to the transient phase of the algorithm) it as and one can obtain the desired statistics from the characterization of the marginal distribution. For instance, point estimates of the parameter likê The main assumption in Algorithm 2 is the ability of evaluating the energy function, ϕ T (·). We have seen earlier how this can be done in a PF. Roughly speaking, the PMCMC algorithm consists on putting together these two algorithms [42]. The resulting PMCMC method can be consulted in Algorithm 3. Compute θ * = θ (j−1) + S j−1 a

5:
Run the PF in Algorithm 1 with model parameters set to θ * . Required outputs: State filteringx * 1:T as in (23) Energy functionφ T (θ * ) as in (35) 6:  Compute η (j) = j −γ 15: Three sets of simulations are discussed. First, we validated the filtering method considering perfect knowledge of the model. In this case, the method in Algorithm 1 was used. Secondly, the model assumptions were relaxed in the sense that the parameters of the model were not known by the method. We analyzed the capabilities of the proposed method to infer both the time-evolving states of the system and some of the parameters defining the model. In this case, the method in Algorithm 3 was used. Finally, we validated the performance of the proposed methods in inferring the synaptic conductances. We tested both PF and PMCMC methods, that is with and without full knowledge of the model respectively.

A. Model parameters are known
In the simulations we considered the model inaccuracies described in Section III-A. To excite the neuron into spiking activity a nominal applied current was injected with I o = 110 µA/cm 2 and two values for σ I were considered, namely 1% and 10% of I o . The nominal conductance used in the model wasḡ L = 2 mS/cm 2 , whereas the underlying neuron had a zero-mean Gaussian error with standard deviation σḡ L . Two variance values were considered as well, 1% and 10% ofḡ L . Finally, we considered σ n = 10 −3 in the dynamics of the gating variable.
In order to evaluate the efficiency of the proposed estimation method, we computed the Posterior Cramér-Rao Bound (PCRB) [50] in Appendix B. We plot the PCRB as a benchmark for the root mean square error (RMSE) curves obtained by computer simulations, obtained after averaging 200 independent Monte Carlo trials. For a generic time series w k , the RMSE of an estimatorŵ k is defined as whereŵ j,k denotes the estimate of w k at the j-th realization and M the number of independent Monte Carlo trials used to approximate the mathematical expectation. Figures 2 and 3 show the time course of the RMSE using N = {500, 1000} particles. We see that in both scenarios, our method efficiently attains the PCRB. We measure the efficiency (η ≥ 1) of the method as the quotient between the RMSE and the PCRB, averaged over the entire simulation time. The worse efficiency on estimating v k was 1.43 corresponding to 500 particles and 10% of inaccuracies (see Fig. 3), the best was 1.11 for 1000 particles and 1% of errors (see Fig. 2). In estimating n k the discrepancy was even lower, 1.06 and 1.03 for maximum and minimum η. As a conclusion, the PF tends to the PCRB with the number of particles. Also, the performance (both theoretical and empirical) could be improved if model inaccuracies are reduced, i.e., if the model parameters are better estimated at a previous stage. For the sake of completeness, we summarize the results in Table I, where the average RMSE and PCRB along the 500 ms simulation can be consulted. It is apparent that increasing the number of particles from N = 500 to N = 1000 does not improve significantly the performance of the method.  To give some intuition on the operation and performance of the PF method in Algorithm 1, we show the results for a single realization in Fig. 4. The results are for 500 particles and two different values of σ 2 y,k , corresponding to 0 and 32 dB respectively. Even in very low SNR regimes, the method is able to operate and provide reliable filtering results.

B. Model parameters are unknown
In this section we validate the algorithm presented in Section IV-B. According to the previous analysis, we deem that 500 particles are enough for the filter to provide reliable results. The parameters of the PMCMC algorithm were set to γ = 0.9 andᾱ * = 0.234. Fig. 5 shows the results for a single realization when a number of parameters in the nominal model are unknown. We considered 1, 2, and 4 unknown parameters. Each of the plots In these plots we omitted the results for the gating variable for the sake of clarity. The true and initial values used in the experiments, as well as the initial covariances assumed, can be consulted in Table II. From the plots we can observe that the method performs reasonably well even in the case of estimating the model parameters at the same time it is filtering out the noise in the membrane voltage traces. A biologically meaningful signal is the leakage current. In general, the leakage gathers those ionic channels that are not explicitly modeled and other non-modeled sources of activity. The parameters driving the leak current areḡ L and E L . We tested and validated the proposed PMCMC in an experiment where the leak parameters were estimated at the same time the filtering solution was computed. Moreover, the statistics of the process noise were estimated as well, Σ x,k . In this case, we iterated the PMCMC method 1000 times and average the  results over 100 Monte Carlo independent trials. The results can be consulted in Fig. 6, where the RMSE performance of the PMCMC method is compared to the performance of the original PF with perfect knowledge of the model.
It can be observed that the filtering performances with perfect knowledge of the model and with estimation of parameters by PMCMC are similar. Moreover, both approaches attain the theoretical lower bound of accuracy given by the PCRB.
In Fig. 7, validation results for the parameter estimation capabilities of the PMCMC are shown. Particularly, we plotted in Fig. 7(a)  . We observe that all of them converge to the true values of the parameter. Recall that these true values were θ = (ḡ L , E L ) = (2, −60) . In Fig. 7(b) and 7(d), the average of these realizations can be consulted, where the aforementioned convergence to the true parameter is highlighted.

C. Estimation of synaptic conductances
Finally, once the methods to estimate state variables and unknown parameters were consolidated, we proceeded to  1  10  0.5   TABLE II  TRUE VALUE, INITIAL  test the methods to our ultimate goal: estimating jointly the intrinsic states of the neuron and the extrinsic inputs (i.e., the synaptic conductances). First, the method with perfect knowledge of the model was validated in Fig. 8. It can be observed in Fig. 8(a) that the intrinsic signals can be effectively recovered as before where synaptic inputs were not accounted for. The estimation of g E (t) and g I (t) is seen in Fig. 8(b). We see that the estimation of the the excitatory and inhibitory terms is quite accurate, and that the presence of spikes does not degrade the estimation capabilities of the method.
The PMCMC algorithm was tested similarly. In this case, we assumed that the model parameters related to v k and n k were accurately estimated, for instance using an off-line procedure or that analyzed in Section V-B. Therefore, we focused on the estimation of those parameters that describe the OU process of each of the synaptic terms. Particularly, we considered the values in Table III. The results can be consulted in Fig. 9 and compared to those in Fig. 8. We observe that  g L and E L were unknown.

VI. CONCLUSIONS
In this paper we propose a filtering method that is able to sequentially infer the time-course of the membrane potential, the intrinsic activity of ionic channels, and the input synaptic conductances from noisy observations of voltage traces. The method works both for subthreshold and spiking regimes. It is based on the PF methodology and features an optimal im- portance density, providing enhanced use of the particles that characterize the filtering distribution. In addition, we tackle the problem of joint parameter estimation and state filtering by extending the designed PF with an MCMC procedure in an iterative method known as PMCMC. Another distinctive contribution with respect to the other works in the literature is that here we provide accuracy bounds for the problem at hand, given by the PCRB. The RMSEs of our methods are then compared to the bound and, therefore, we can assess the efficiency of the proposed inference methods.
Filtering methods of different types (e.g., PF or sigma-point Kalman filtering) have been used in other recent contributions to similar problems, see [7]- [13]. From a methodological perspective, the novelty of this paper is in the use of an optimal importance density to generate particles, fact that increases the estimation accuracy for a given budget of particles. This technical detail only applies to PF methods. Although Gaussian methods (e.g. the family of sigma-point Kalman filters) have a lower computational cost in general, they require Gaussianity of the measures, whereas PFs do not. This is an advantage that we think can be crucial in estimating synaptic conductances, since the assumption of Gaussianity is generally assumed in the literature ( [4], [31]) but there are no conclusive evidences to assert this assumption. In this paper, we have still applied the PF to an Ornstein-Uhlenbeck process in order to check that basic results can be attained, but further research will go in the direction of assuming other types of distributions for the synaptic conductances. The use of more complex distributions nicely fits within the framework of our PF-based method. Another advantage of PFs versus Gaussian filters is their enhanced robustness to outliers [51], for instance due to recording artifacts, future applications shall also incorporate this feature.
We have found excellent estimations of the synaptic conductances, even in spiking regimes. Estimating synaptic conductances in spiking regimes is a challenge which is far to be solved. It is well-known that linear estimations of synaptic conductances are not trustable in this situation when data is extracted intracellularly from spiking activity of neurons, see [14]. In experiments, thus, caution has to be taken in eliminating part of the voltage traces, thus loosing also part of the temporal information of both excitatory and inhibitory conductances. Our method is able to perform well in this regime. This information is highly valuable in problems (epilepsy, schizophrenia, Alzheimer's disease,. . . ) where a debate on the balance of excitation and inhibition is open, see the introduction of [12] for a rather complete overview on this feature.
The results show the validity of the approach when applied to Morris-Lecar type of neuron. However, the procedure is general and could be applied to any neuron model, exhibiting more complex dynamics like bursting, mixed-mode oscillations,. . . Forthcoming applications would be validating the method using real data recordings, both for inferring parameters of the model and synaptic conductances. The latter problem is a challenging hot topic in the neuroscience literature, which is recently focusing on methods to extract the conductances from single-trace measurements. We think that our PF method would give useful and interesting results to physiologists that aim at inferring brain's activation rules from neurons' activities. Actually, knowing the excitatoryinhibitory time course separation can help in getting important conclusions about brain's functional connectivity (see [52]- [54]).
We have not tried to obtain estimations when subthreshold ionic currents are active, where the presence of nonlinearities could also contaminate the estimations, see [15]. According to the excellent performance in spiking regimes, where nonlinearities are stronger, we expect also a good agreement between the estimated data and the prescribed synaptic conductances.
Other extensions of the model can be devoted to incorporate the dendro-somatic interaction (see for instance [55]), by considering multi-compartmental neuron models, thus taking into account the morphology and the functional properties of the cell. This is another big challenge for which we think that our method can account for.

APPENDIX A MORRIS-LECAR NEURON MODEL
From the myriad of existing single-neuron models, we consider without loss of generality the Morris-Lecar model proposed in [35]. The model can be related (see [28]) to the I Na,p + I K -model (pronounced persistent sodium plus potassium model). The dynamics of the neuron is modeled by a continuous-time dynamical system composed of the currentbalance equation for the membrane potential, v = v(t), and the K + gating variable 0 ≤ n = n(t) ≤ 1, which represents the probability of the K + ionic channel to be active. Then, the system of differential equations is where C m is the membrane capacitance and φ a nondimensional constant. I app represents the (externally) applied current. For the time being, we have neglected I syn in (39). The leakage, calcium, and potassium currents are of the form respectively.ḡ L ,ḡ Ca , andḡ K are the maximal conductances of each current. E L , E Ca , and E K denote the Nernst equilibrium potentials, for which the corresponding current is zero, a.k.a. reverse potentials. The dynamics of the activation variable m is considered at the steady state, and thus we write m = m ∞ (v). On the other hand, the time constant τ n (v) for the gating variable n cannot be considered that fast and the corresponding differential equation needs to be considered. The formulae for these functions is which parameters V 1 , V 2 , V 3 , and V 4 can be measured experimentally [28]. The knowledgeable reader would have noticed that the Morris-Lecar model is a Hodgin-Huxley type-model with the usual considerations, where the following two extra assumptions were made: the depolarizing current is generated by Ca 2+ ionic channels (or Na + depending on the type of neuron modeled), whereas hyperpolarization is carried by K + ions; and that m = m ∞ (v). The Morris-Lecar model is very popular in computational neuroscience as it models a large variety of neural dynamics while its phase-plane analysis is more manageable as it involves only two states [56].
The Morris-Lecar, although simple to formulate, results in a very interesting model as it can produce a number of different dynamics. For instance, for given values of its parameters, we encounter a subcritic Hopf bifurcation for I app = 93.86 µA/cm 2 . On the other hand, for another set of parameter values, the system of equations has a Saddle-Node on an Invariant Circle (SNIC) bifurcation at I app = 39.96 µA/cm 2 .

APPENDIX B PCRB IN MORRIS-LECAR MODELS
This appendix is devoted to the derivation of the PCRB estimation bound for the Morris-Lecar model used to benchmark the proposed methods in the simulations. We follow the sequential procedure given in [57], accounting that we have nonlinear functions in the state evolution and linear measurements, both with additive Gaussian noise. We are interested in an estimation error bound of the type of wherex k (y 1:k ) represents an estimator of x k given y 1:k .
Recall that the state-space we are dealing with is of the form where h = (1, 0), x k = (v k , n k ) , and f k−1 (x k−1 ) defined by (9) and (10). The noise terms are of the form In this case, the PCRB can be computed recursively by virtue of the result in [57] by computing the following terms x,kF k (51) and plugging them into for some initial J 0 . Notice that, in our case, D 22 k becomes deterministic, but the rest of terms involving expectations should be computed by Monte Carlo integration over independent state trajectories.
Since the state function is nonlinear, we use the Jacobian evaluated at the true value of x k instead, that is where functions f 1 and f 2 are as in (9) and (10), respectively. Therefore, to evaluate the bound we need to compute the derivatives in the Jacobian. These are,