Parameter Estimation of a Class of Neural Systems with Limit Cycles

: This work addresses parameter estimation of a class of neural systems with limit cycles. An identification model is formulated based on the discretized neural model. To estimate the parameter vector in the identification model, the recursive least-squares and stochastic gradient algorithms including their multi-innovation versions by introducing an innovation vector are proposed. The simulation results of the FitzHugh–Nagumo model indicate that the proposed algorithms perform according to the expected effectiveness.


Background
The nervous system is a complex and huge system, composed of a large number of nerve cells.The performance of neurons is closely related to the function of the nervous system.Oscillatory activity is ubiquitous in the nervous system and its synchronous behavior has become an increasing topic in neuroinformatics [1,2].Although great progress has been made in information coding and understanding nonlinear dynamics of neural oscillator systems in neuroscience and from the system control perspective, there is relatively little research on the modeling and control of neural oscillator systems.In a typical neural system, the measurable state is limited by using voltage clamp technique, while the internal parameters of the system can only be estimated by the measurable state output sequence such as membrane voltage.Sometimes, it is even impossible for biologists to get the modeling information due to various technical difficulties.For example, in the FitzHugh-Nagumo model, the membrane voltage of the neuron is measurable and the model control parameters are unknown, which relies on the parameter identification methods in control theory.In particular, parameter estimation for neural models are used to estimate the uncertain parameters from the available input and output data.

Parameter Estimation in Neural Model
System parameter estimation has received substantial attention in several decades [3][4][5][6][7].One reason is its broad application to mechanical systems [8,9], neuron systems [10], communication systems [11], structural systems [12], power systems [13], etc. Existing black-box parameter identification methods for nonlinear systems are simple and universal, but the identification accuracy and convergence of the identification algorithm are not of universal significance.Therefore, accurately identifying the system parameters of the spiking neuron models under finite measurement presents a far more challenging modeling problem.Examples of parameter identification methods for neural models include the time/frequency-domain method [14], maximum likelihood method [15], self-organizing state-space-model approach [16], heuristic optimization method [10], and so on.In [14], a time-domain method and a frequency-domain method are applied to estimate parameters of a neuronal model consisting of a soma coupled to a uniform dendritic cylinder, where the time-domain method is more prone to estimation errors in the cable parameters.Vavoulis et al. [16] presented an adaptive sampling algorithm using the self-organizing state-space model to estimate parameters in a Hodgkin-Huxley-type model of single neurons and to achieve reduced variance of parameter estimates.Mullowney et al. [15] proposed a maximum likelihood methodology for parameter estimation in a leaky integrate-and-fire neuronal model with the only available data being the interspike intervals or the times between firings.In [10], a global heuristic search method using in-vitro and in-vivo electrophysiological data was explored for identification of an Izhikevich-type neuron model.
In spite of these advances, little attention on parameter estimation of the FitzHugh-Nagumo (FHN) model has been received, including Bayesian statistical approaches [17,18], and least-squares algorithms [19,20].In [17], a Bayesian framework was proposed for drift parameter estimation of the stochastic FHN model.Arnold and Lloyd [18] employed nonlinear filtering for periodic, time-varying parameter estimation, which was illustrated in estimation of the external voltage parameter in the FHN model.In comparison, from the recursive estimation point of view, we perform parameter inference for the FHN model with data generated from the model in this paper.Concha and Garrido [19] presented two methodologies based on the least-squares algorithm to estimate the parameters of the FHN model.The two methods were only suitable for the case with noncontinuous input current stimulus, at the cost of a linear integral filter.Che et al. [20] employed the recursive least-squares algorithm for parameter estimation of the FHN model, which requires the first and second time derivatives of the membrane potential and the input current stimulus being continuously differentiable.It is well known that the stochastic gradient algorithm is a class of important stochastic approximation methods, which have received much attention and have been widely used in different systems, such as Hammerstein systems [21], Wiener systems [22] and sampled systems [23].Though the stochastic gradient algorithm has a slower convergence rate than the least-squares algorithms, but it requires less computational effort.In this paper, to improve the convergence rate of the stochastic gradient algorithm, we extend the innovation concept in [24] and explore the multiinnovation stochastic gradient algorithm for parameter estimation of the FHN neuron system.The parameterized FHN model in [20] can also be performed by using the proposed multiinnovation recursive least-squares algorithm in this paper and a better accuracy of the parameter estimation will be obtained due to the use of past innovations and repeated available data.However, from the aspect of the computing style, identification algorithms with qualified identification accuracy and convergence are desired for parameter estimation of the neural models.In particular, it seems that little effort has been made toward performing parameter estimation of the FHN model using the stochastic gradient methods.

Contributions
Inspired by these works, we propose four parameter estimation algorithms for the FHN neuron system with limit cycle and external disturbance.The contributions of this paper include the following:

•
We formulate the FHN neuron system as an identification model based on the explicit forward Euler method.

•
We propose a recursive least-squares algorithm and a stochastic gradient algorithm to estimate the unknown parameters of the model.

•
We extend the innovation concept in [24], and explore the multiinnovation recursive least-squares algorithm and multiinnovation stochastic gradient algorithm for parameter estimation of the FHN neuron system.

•
We show that a faster convergence rate and better accuracy can be achieved using the innovation and repeated available data.

Organization
The organization of the paper is as follows.Section 2 describes the FHN neural model.Section 3 formulates the identification model and presents the parameter estimation problem.The proposed algorithms are given in Section 4, where we estimate the unknown parameters using four different algorithms and summarize corresponding algorithm procedures.Section 5 provides the computational simulations.Finally, we draw some conclusions in Section 6.

The Spiking Neuron Model
In this section, we introduce the general spiking neuron models.A conductance-based spiking neuron model is described by where x(t) ∈ R n is the state which usually consists of the membrane potential, gating variable, recovery variable and adaptation variable; ξ(t) ∈ R n denotes the disturbance in the neuron system; u(t) ∈ R is the external input injected to the neuron; B ∈ R n is a constant vector, which is typically taken as ] when the external input acts on the membrane potential.
The system described by ( 1) is quite general as it includes many spiking neuron models such as the Hodgkin-Huxley model, the Hindmarsh-Rose (HR) model, the FHN model and so on.In this paper, we shall focus on parameter estimation of the FHN model using different identification algorithms, but the proposed approach is applicable to other neuron models in a similar manner.The FHN model in a dimensionless form [25,26] can be written in the form of ( 1) where v denotes the voltage potential of the neuron membrane and w denotes the inactivation of the sodium channels; ξ ∈ R 2 denotes the disturbance in the neuron system and is the white noise with zero mean.The parameters a, b, c 1 , c 2 , µ are unknown to be estimated later in detail.In this paper, we consider a constant external input u(t) ≡ J resulting in periodic spiking dynamics and the external input J will also be estimated.Note that when the parameters are specified as a = 0.1, b = 1, c 1 = 1, c 2 = 0.5, µ = 100 and the input current value is taken as J = 0.5, the neural system from (v(0), w(0)) = (−0.3,0.6) without disturbance exhibits periodic spiking dynamics and converge to a limit cycle (see Figure 1).In the real nervous system, the parameters of the spiking neurons are difficult to be measured or determined, the identification of system parameters has always been an important topic in the field of neural computing and system control.Motivated by this, the present work is directed towards developing efficient algorithms to identify the parameters of spiking neuron models.

The Identification Model of Spiking Neurons
Usually, the neuron membrane and the inactivation of the sodium channels are discretely sampled in experiment.Hence, let us consider the discretized system of the FHN model using the explicit forward Euler method with step size T as follows where T is the sampling period which is sufficiently small, k represents the kT time instant and Define Therefore, Equation (4) can be rewritten as which is called as the identification model.Apparently, to estimate the parameters a, b, µ, J, c 1 , c 2 , it is equivalent to estimate the parameter vector θ.

Least-Squares Estimation Algorithms
To estimate the parameters θ, let us consider the cost function where L d denotes the data length.Using the least-squares principle to solve the optimization problem, we explore the following recursive least-squares (RLS) parameter estimation formulas [3,27] to estimate the parameter vector θ(k).
To initialize the RLS algorithm, the initial value θ(0) is generally taken to be a zero vector or a small real vector, e.g., θ(0) = 10 −6 1 ñ with 1 ñ being an ñ-dimensional column vector whose elements are 1.Now, the RLS parameter estimation algorithm for the spiking neuron model is summarized in Algorithm 1.

Algorithm 1 RLS algorithm
(1) Discretize the FHN model based on the explicit forward Euler method with step size T.
(3) Collect the measurement state data and determine a data length L d .Then, form the output data y(k) and the information vector φ(k).(4) Compute L(k) and P(k) by 6) If k < L d , increase k by 1 and go to step (3); otherwise, stop the procedure and output the estimate θ(L d ) of the parameter vector θ.
To obtain better estimation accuracy, we modify the algorithm (6) and apply a multiinnovation recursive least-squares (MIRLS) parameter estimation algorithm [7] for (5) by introducing an innovation length p.By defining the information matrix Φ(p, k) and stacked output vector Y(p, k) as the innovation vector E (p, k) can be expressed as From here, we present the following MIRLS iterative formulas The MIRLS parameter estimation algorithm for the FHN model is summarized in Algorithm 2.

Algorithm 2 MIRLS algorithm
(1) Discretize the FHN model based on the explicit forward Euler method with step size T.
(3) Collect the measurement state data and determine a data length L d .Then, form the output data y(k) and the information vector φ(k).(4) Given an innovation length p, form Y(p, k) by ( 7) and Φ(p, k) by ( 8).
(5) Compute L(k) and P(k) by increase k by 1 and go to step (3); otherwise, stop the procedure and output the estimate θ(L d ) of the parameter vector θ.

Stochastic Gradient Estimation Algorithms
Let X 2 := tr[XX ] denote the norm of the matrix X.For the model ( 5), we present the following stochastic gradient (SG) parameter estimation formulas to estimate the parameter vector θ(k) where α is the forgotten factor.Note that a smaller α results in a faster convergence but the price we paid is a large parameter fluctuation.In this paper, since the SG algorithm converges slowly in parameter estimation of the FHN neuron, we choose a moderate α = 0.8.The SG parameter estimation algorithm for the spiking neuron model is summarized in Algorithm 3.

Algorithm 3 SG algorithm
(1) Discretize the FHN model based on the explicit forward Euler method with step size T.
, reset a large forgotten factor α ∈ (0, 1] and go to step (7); otherwise, go to step (7).(7) If k < L d , increase k by 1 and go to step (3); otherwise, stop the procedure and output the estimate θ(L d ) of the parameter vector θ.
Similarly, to obtain better estimation accuracy, we can derive multiinnovation stochastic gradient (MISG) parameter estimation formulas [28] using the Y(p, k) and Φ(p, k) in Section 4.1 as follows Comparing with the SG algorithm in (10) using only the current data, the innovation Φ in (11) uses not only the current innovation, but also the past innovations, which thus can improve the convergence rates compared with the SG algorithm.Moreover, the available data are repeatedly used in the MISG algorithm, and such a treatment can enhance accuracy of the parameter estimation [24].The MISG parameter estimation algorithm for the FHN model is summarized in Algorithm 4.

Simulations
For purpose of illustration, we consider the FHN model [25,26] described in the dimensionless form (2). In simulation, the parameters are taken as a = 0.1, b = 1, c 1 = 1, c 2 = 0.5, µ = 100 and J = 0.5.The FHN model is firstly discretized based on the explicit forward Euler method with step size T = 10ms.Then, consider the initial value v(0) = −0.3,w(0) = 0.6, and perform the simulations based on the identification algorithms in Section 4. Taking the data length L d = 20, 000, the forgotten factor λ is chosen as 0.99 for RLS and MIRLS algorithms.For both SG and MISG algorithms, we take the forgotten factor α = 0.8.To compare the estimation performance of the RLS, SG, MIRLS and MISG algorithms, apply these four algorithms to estimate the parameter vector θ of the neuron system, respectively.Note that, when p = 1, the MIRLS algorithm reduces to the RLS algorithm and the MISG algorithm reduces to the SG algorithm.For the MIRLS and MISG algorithms, two cases with different innovation lengths p = 3 and p = 5 are considered, respectively.To quantify the estimation accuracy and clearly compare the performance of the RLS, SG, MIRLS and MISG algorithms, we consider two different noise levels, i.e., σ 2 = 0.2 2 and σ 2 = 0.5 2 .The parameter estimates and their errors are shown in Tables 1-4 with different noise variances.The estimation errors δ versus k are shown in Figures 2  and 3, where δ := θ − θ / θ .From Tables 1 and 2, it is seen that both the RLS and MIRLS algorithms have a fast convergence rate and a high estimation accuracy, which can also be observed from Figures 2 and 3. On the other hand, by using the SG and MISG algorithms, it is shown in Tables 3 and 4 that it requires a data length around L d = 20, 000 to achieve an acceptable estimation accuracy.Though not shown in table, the MIRLS algorithm with p = 5 enjoys the most accurate parameter estimate and the fastest convergence rate.From Figures 2 and 3, we can see that for the same batch of data, compared with the RLS algorithm, the SG algorithm extracts less information from the measured data and uses information inefficiently, which results in a much slower convergence.Hence, to show the effect clearly, the results by the RLS and SG algorithms are shown in two figures, respectively.Meanwhile, from Figure 2, it can be found that due to the high efficiency of the RLS algorithm, the MIRLS algorithm has limited improvement potential for the accuracy of parameter estimation.However, the computation efficiency of the MIRLS algorithm will be revealed in the case of data missing, which remains our future work.From Figure 3, it is clear that the MISG algorithm has a faster convergence rate than the SG algorithm or the MISG algorithm with p = 1 and the MISG estimates with p = 3 and p = 5 have higher accuracy than the SG estimates.In fact, the parameter estimation errors by the MISG algorithm become smaller and smaller as the innovation length p increases.

Figure 1 .
Figure 1.A limit cycle of the FHN model.

Table 1 .
The RLS estimates and errors.

Table 3 .
The SG estimates and errors.