Parameter Estimation for Hindmarsh–Rose Neurons

Abstract

In the paper, a new adaptive model of a neuron based on the Hindmarsh–Rose third-order model of a single neuron is proposed. The learning algorithm for adaptive identification of the neuron parameters is proposed and analyzed both theoretically and by computer simulation. The proposed algorithm is based on the Lyapunov functions approach and reduced adaptive observer. It allows one to estimate parameters of the population of the neurons if they are synchronized. The rigorous stability conditions for synchronization and identification are presented.

1. Introduction

The human brain is one of the most complex systems existing on Earth. Many researchers have been working hard for many years in order to understand how the human brain is functioning, to create its model, and to redesign its construction [1,2]. Recently, the study of the whole brain dynamics models has started [3,4]. However, the overwhelming complexity demands a variety of simplified models of different complexity degrees. The simplest models are represented as the networks of biological neuron models connected by couplings. Typical models of this class are based on the celebrated Hodgkin–Huxley model and its simplifications: FitzHugh–Nagumo (FHN), Morris–Lecar (ML), Hindmarsh–Rose (HR), etc. New models still appear in the literature, e.g., [5]. The HR model is the simplest one that can exhibit most of the biological neuron’s behavior, such as spiking or bursting.

Chong et al. [6] studied the activity of neurons at the macroscopic level, known in the literature as “neural mass models”. To estimate the neurons’ activity, a robust circle criterion observer was proposed, and its application to estimate the average membrane potential of neuron populations in the single cortical column model was demonstrated. Further steps were based on using network models. However, well-understood linear network models [7,8] are not suitable for modeling networks of biological neurons. A number of approaches to state estimation for neural populations are based on synchronization; see [9,10] and the references therein.

An important part of brain network modeling is the identification of the model parameters. Its first step is the identification of the building blocks—the models of the single neurons or their subpopulations. Ideally, it should be performed based on the neuron state measurements. However, there exists the problem of uncertainty and incomplete measurements. A number of such approaches based on adaptation and learning are described in the literature. Dong et al. [11] dealt with an identification of the FHN model dynamics, employing the deterministic learning and interpolation method. For global identification, the FHN model was transformed in [11] into a set of ordinary differential equations (ODEs), and the dynamics of the approximation system were then identified by employing deterministic learning. In [12], an approach to solving the problem of identifying topology and parameters in HR neural networks was proposed. For this purpose, the so-called generalized extremal optimization (GEO) was introduced, and a heuristic identification algorithm was employed. Identifying HR neural networks’ topology was also considered by Zhao et al. [13], who employed the sinusoidal disturbance to identify the topology at the stage when the complex network achieves synchronization. It was demonstrated by the simulations that, compared with the disturbance of all the nodes, the disturbance of the key nodes alone can achieve a very good effect. In [14], an adaptive observer for asymptotical estimation of the parameters and states of a model of interconnected cortical columns was presented. The model adopted is capable of realistically reproducing the patterns seen on (intracranial) electroencephalograms (EEGs). The assessment of the parameters and status allows a better understanding of the mechanisms underlying neurological phenomena and can be used to predict the onset of epileptic seizures. Tang et al. [15] studied the effect of electromagnetic induction on the electrical activity of neurons, and the variable for magnetic flow was used to improve the HR neuron model. By the simulations, it was demonstrated that the neuron model, proposed in [15], can show multiple modes of electrical activity, which is dependent on the time delay and external forcing current. In [16], an approach based on adaptive observers was developed for partial identification HR model parameters. Malik and Mir [17] studied the synchronization of HR neurons, demonstrating that the coupled system shows several behaviors depending on the parameters of the HR model and coupling function. Recently, Xu [18] proposed to use an impulse response identification experiment with dynamical observations with an increasing data length for capturing the real-time information of systems and serving for online identification. In Xu [18], a separable Newton recursive parameter estimation approach was developed, and its efficacy was demonstrated by the Monte Carlo tests.

The models of biological neural networks and single neurons have many applications, e.g., in brain–computer interfaces [19,20,21]. A practical application of coupled HR neural networks was also discussed in [17]. In particular, it was demonstrated that the spiking network successfully encodes and decodes a time-varying input. Synchronization of the artificial HR neurons was also studied in [22,23]. An adaptive controller that provides synchronization of two connected HR neurons using only the output signal of the reference neuron was suggested in [22]. Andreev and Maksimenko [24] considered synchronization in a coupled neural network with inhibitory coupling. It was shown in [24] that in the case of a discrete neuron model, the periodic dynamics are manifested in the alternate excitation of various neural ensembles, whereas periodic modulation of the synchronization index of neural ensembles was observed in the continuous-time model.

In this paper, the problem of HR neuron model parameters’ identification is considered. To solve this, the reduced-order adaptation algorithm based on the speed gradient (SG) method and feedback Kalman–Yakubovich lemma (FKYL) [25,26] is proposed. A rigorous statement about the convergence of the parameter estimates to the true values was formulated, and the proof is given. The performance of the identification procedure was analyzed by computer simulation.

The remainder of the paper is organized as follows. The problem statement is presented in Section 2. Section 3 presents the design of the adaptation algorithm. The main results are given in Section 4. Computer simulation results are described in Section 5. Concluding remarks and the future work intentions of Section 6 finalize the paper.

2. Problem Statement

Consider the classical Hindmarsh–Rose (HR) neuron model [27]:

$$\begin{array}{c}{\dot{x}}_1\left(t\right) = x_2\left(t\right)-a{x}_1^3\left(t\right) + b{x}_1^2\left(t\right)-x_3\left(t\right) + I\hfill \\ {\dot{x}}_2\left(t\right) = c-d{x}_1^2\left(t\right)-x_2\left(t\right)\hfill \\ {\dot{x}}_3\left(t\right) = \epsilon s{x}_1\left(t\right)-\epsilon sr-\epsilon x_3\left(t\right)\hfill \end{array}$$

(1)

The variable $x_1\left(t\right)$ is the membrane potential, while $x_2\left(t\right), x_3\left(t\right)$ are the fast and slow ionic currents, respectively. Variables $x_1\left(t\right), x_2\left(t\right), x_3\left(t\right)$ represent the state vector $x\left(t\right)$. The values $a, b, c, d, \epsilon , s, r, I$ are model parameters. We considered all quantities to be dimensionless.

Assume that $x_1\left(t\right)$, $x_2\left(t\right)$, $x_3\left(t\right)$ can be measured. In some cases (e.g., in in vitro experiments), this assumption is quite realistic. In order to estimate the state and parameters of the HR neuron model (1), introduce an auxiliary system—adaptive model—as follows:

$$\begin{array}{c}\hfill \begin{array}{c}{\dot{z}}_1\left(t\right) = z_2\left(t\right)-\widehat{a}{x}_1^3\left(t\right) + \widehat{b}{x}_1^2\left(t\right)-z_3\left(t\right) + \widehat{I}-k(x_1\left(t\right)-z_1\left(t\right))\hfill \\ {\dot{z}}_2\left(t\right) = \widehat{c}-\widehat{d}{x}_1^2\left(t\right)-z_2\left(t\right)\hfill \\ {\dot{z}}_3\left(t\right) = \epsilon s{x}_1\left(t\right)-\epsilon sr-\epsilon z_3\left(t\right)\hfill \end{array}\end{array}$$

(2)

The variables $z_1\left(t\right), z_2\left(t\right), z_3\left(t\right)$ represent the state vector $z\left(t\right)$ of the adaptive model (2). The term $k_1(x_1\left(t\right)-z_1\left(t\right))$ is the stabilizing term, typical for observers.

In this study, we assumed for simplicity that the value $\epsilon$ characterizing the relative rate of fast and slow currents is known. Then, one may denote the vector of tunable parameters $\theta = {[{\theta }_1,\dots ,{\theta }_7]}^{\mathrm{T}}$, where ${\theta }_1 = \widehat{a}$, ${\theta }_2 = \widehat{b}$, ${\theta }_3 = \widehat{I}$, ${\theta }_4 = \widehat{c}$, ${\theta }_5 = \widehat{d}$, ${\theta }_6 = \widehat{\epsilon s}$, ${\theta }_7 = \widehat{\epsilon sr}$. The problem is to design the adaptation/learning algorithm for $\theta$ ensuring the goals:

$$\underset{t\to \infty }{lim}\parallel z\left(t\right)-x\left(t\right)\parallel = 0,$$

(3)

$$\underset{t\to \infty }{lim}\parallel \theta \left(t\right)-{\theta }^{\ast }\parallel = 0,$$

(4)

where ${\theta }^{\ast }={[a,b,I,c,d,\epsilon s,\epsilon sr]}^{\mathrm{T}}$.

3. Adaptation/Learning Algorithm Design

The equations for the state estimation error $e\left(t\right) = x\left(t\right)-z\left(t\right)$ are as follows:

$$\begin{array}{c}\hfill \dot{e} = A_{k}e + B_1({\theta }_1(-x_1^3) + {\theta }_2{x}_1^2 + {\theta }_3·1) + B_2({\theta }_4 + {\theta }_5(-x_1^2)) + B_3({\theta }_6{x}_1 + {\theta }_7(-1)),\end{array}$$

(5)

where $A_k = \left[\begin{array}{ccc}-k& 1& -1\\ 0& -1& 0\\ 0& 0& -\epsilon \end{array}\right]$, $B_1 = \left[\begin{array}{c}1\\ 0\\ 0\end{array}\right]$, $B_2 = \left[\begin{array}{c}0\\ 1\\ 0\end{array}\right]$, $B_3 = \left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]$

For the design of the adaptation/learning algorithm, we used the speed gradient method [28], suggesting to change the tunable vector $\theta \left(t\right)$ in the direction of the gradient in $\theta$ of the speed of changing the goal function $Q\left(e\right)$, where $e = x-z$ along the trajectories of the system (1), (2). It is suggested to choose quadratic goal function $Q\left(e\right) = e^{\mathrm{T}}Pe$, where $P = P^{\mathrm{T}}>0$ is a positive-definite matrix to be determined later. Since $Q\left(e\right)$ depends only on e, the error model (5) can be used instead of (1), (2) at this stage. Applying the SG methodology [25,28] yields the following algorithm:

$$\begin{array}{c}{\dot{\widehat{\theta }}}_1 = -\gamma e^{\mathrm{T}}P{B}_1·(-x_1^3)\hfill \\ {\dot{\widehat{\theta }}}_2 = -\gamma e^{\mathrm{T}}P{B}_1·x_1^2\hfill \\ {\dot{\widehat{\theta }}}_3 = -\gamma e^{\mathrm{T}}P{B}_1·1\hfill \\ {\dot{\widehat{\theta }}}_4 = -\gamma e^{\mathrm{T}}P{B}_2·1\hfill \\ {\dot{\widehat{\theta }}}_5=-\gamma e^{\mathrm{T}}P{B}_2·(-x_1^2)\hfill \\ {\dot{\widehat{\theta }}}_6=-\gamma e^{\mathrm{T}}P{B}_3·x_1\hfill \\ {\dot{\widehat{\theta }}}_7=-\gamma e^{\mathrm{T}}P{B}_3·(-1)\hfill \end{array}$$

(6)

4. Main Results

The following proposition concerning the convergence of the adaptive model to the true one holds.

Theorem 1.

Let the parameter ε of the HR neuron model be known. Let $k_1>0$ and $\gamma >0$. Then, the goals (3), (4) are achieved for any initial conditions of the system (1), (2), (6).

The proof is based on exploiting the Lyapunov function:

$$\begin{array}{c}\hfill V\left(e\right) = e^{\mathrm{T}}Pe + {\parallel \theta \parallel }^2/\left(\gamma \right)\end{array}$$

(7)

Evaluating the derivative of (7) along the system (1), (2), (6), one obtains:

$$\begin{array}{c}\hfill \dot{V}\left(e\right) = e^{\mathrm{T}}(A_k^{\mathrm{T}}P + P{A}_k)e + 2e^{\mathrm{T}}P{B}_1\left(-{\theta }_1{x}_1^3 + {\theta }_2{x}_1^2 + {\theta }_3\right)\\ \hfill +2e^{\mathrm{T}}P{B}_2\left({\theta }_4-{\theta }_5{x}_1^2\right)+2e^{\mathrm{T}}P{B}_3\left({\theta }_6{x}_1-{\theta }_7\right)+{\displaystyle \sum _{i=1}^7}{\theta }_i{\dot{\theta }}_i/\gamma \end{array}$$

(8)

Applying (6) to (8) yields $\dot{V}\left(e\right)=e^{\mathrm{T}}(A_k^{\mathrm{T}}P+P{A}_k)e$. Condition $\dot{V}\left(e\right)<0$ is fulfilled $\forall e:e\ne 0$ when:

$$A_k^{\mathrm{T}}P+P{A}_k=-Q,$$

(9)

where $Q=Q^{\mathrm{T}}>0$. Since the eigenvalues of $A_k$ have a negative real part, P, introduced before, is the solution of Lyapunov Equation (9). The fact that $V\left(e\right(t\left)\right)$ is bounded is obtained from [29]. It is shown that the solutions of the system (1) are bounded, and the synchronization error is also bounded.

Definition 1.

Time-varying vector function $f:[0,+\infty ]\to {\mathbb{R}}^m$ is persistent excitation (PE) if it is bounded and there exists $T>0$ and $\alpha >0$ such that $\forall t>0$:

$$\underset{t}{\overset{t+T}{\int }}f\left(s\right)f{\left(s\right)}^{\mathrm{T}}ds⩾\alpha E.$$

(10)

Theorem 5.1 in Fradkov [25], the conditions for achieving goal (4) are derived. His theorem is reformulated below for the problem of interest.

Theorem 2.

If trajectories of the systems (1), (2) are bounded as far as $e\left(t\right)$, $\theta \left(t\right)$ are bounded, all roots of $det(\lambda E-A_k)$ have a negative real part, and functions $f_1\left(t\right)={[-x_1^3{x}_1^{2}1]}^{\mathrm{T}}$, $f_2\left(t\right)={[1-x_1^2]}^{\mathrm{T}}$, $f_3\left(t\right) = {[x_1-1]}^{\mathrm{T}}$ are PE, then the tunable parameters converge to the regulator, achieving the goals (3), (4).

Persistent excitation of vector $\mathsf{\Phi }\left(t\right)$ is equivalent to the existence of $\alpha >0,T>0,t_0>0$ such as $\forall h\in \mathbb{R}\left(\right|h| = 1),t>t_0$ satisfying $\underset{[t, t + T]}{max}\left|{\mathsf{\Phi }}^{\mathrm{T}}h\right|>\alpha$ [28]. The proof of PE for the signals mentioned before is quite simple. Let at some moment t, $[-x_1{\left(t\right)}^3{x}_1{\left(t\right)}^{2}1]t = 0\iff -x_1{\left(t\right)}^3{h}_1 + x_1{\left(t\right)}^2{h}_2 + h_3 = 0$. This equation has no more then three separate solutions. Then, for an arbitrary close moment, $t + \Delta t$ equation does not hold. The PE of the rest of the equations is proven similarly.

Remark 1.

The identification algorithm (6) requires the measurement of all state system variables. It is aimed at studying real neural cells “in vitro”, i.e., during real experiments with single cells where the measurement of all auxiliary variables can be implemented.

5. Computer Simulation

5.1. Neuron Modeling

For the systems (1), (2) and algorithm (6), we conducted mathematical modeling to ensure that the approach adequately meets the formulated problem. We built the neuron model (1) and observed different dynamical behaviors, which are presented in Figure 1, Figure 2 and Figure 3.

We chose numerical values of quantities $a, b, c, d, e, s, r, I$ from [30] where the bifurcation analysis of the HR neuron was performed. $a = 1, c = 1, d = 5, s = 4, r = -1$, and parameters $b, r, \epsilon$, and I are in the figure captions.

5.2. Regular Neuron Modeling

Secondly, we modeled the difference of the systems (1) and (2) in order to perform the identification process by Algorithm (6). The parameters of the observer were calculated just in time by the adaptation algorithm. One can see the identification process in Figure 4. It is seen that the observer parameters converge to the reference neuron parameters. However, synchronization error converges to zero, which is shown in Figure 5

Further simulations showed that the convergence depends on the small parameter $\epsilon$, and too small values may result in the divergence of the identification process. The reason is the liner part of error Equation (5). One may subtract regularization term $k{e}_3$ from the third equation of (2). The linear part of $\dot{e}$ will take the form:

$$A_k^{\prime }=\left[\begin{array}{ccc}-k& -1& 1\\ 0& -1& 0\\ 0& 0& -\epsilon -k\end{array}\right]$$

(11)

After this modification, convergence becomes faster. You may see the accelerated identification process in Figure 6 and Figure 7. The synchronization error for irregular bursting is shown on Figure 8.

5.3. Robustness of Identification with Respect to Noise

Suppose time-varying ultimately bounded disturbance $v\left(t\right)$. Let us propagate $v\left(t\right)$ in the first equation on Model (1). If $v_{max}$ is small enough [31], then the adaptive system will have globally bounded solutions. The identification process for the system with disturbance is presented in Figure 9.

6. Conclusions

A new algorithm for Hindmarsh–Rose neuron model parameter estimation based on the Lyapunov functions approach and reduced adaptive observer was proposed. It allows one to estimate the parameters of the population of the neurons if they are synchronized. The rigorous stability conditions for synchronization and identification were presented. Future work will be aimed at taking into account disturbances and measurement noises. It would be also interesting to use the Hindmarsh–Rose model instead of the FitzHugh–Nagumo model for studying gamma oscillation activity in the brain [32]. In future research, it is also planned to apply the results of the works [33], where the state filter for a time delay state-space system with unknown parameters from noisy observation information was proposed, and [34], where a gradient approach for the adaptive filter design based on the fractional-order derivative and a linear filter was developed.