Bounded-Error Parameter Estimation Using Integro-Differential Equations for Hindmarsh–Rose Model

: A numerical parameter estimation method, based on input-output integro-differential polynomials in a bounded-error framework is investigated in this paper. More precisely, the measurement noise and parameters belong to connected sets (in the proposed work, intervals). First, this method, based on the Rosenfeld–Groebner elimination algorithm, is presented. The latter provides differential equations containing derivatives, sometimes of high order. In order to improve the numerical results, a pretreatment of the differential relations is done and consists in integration. The new relations contain, essentially, integrals depending only on the outputs. In comparison with the initial relations, they are less sensitive to measurement noise. Finally, the impact of the size of the measurement noise domain on the estimated intervals is studied.


Introduction
This paper proposes a parameter estimation procedure, based on integro-differential (ID) relations in the set-membership framework. These relations are based on differential relations obtained, owing to the Rosenfeld-Groebner algorithm implemented in the package DifferentialAlgebra of Maple (see [1], for more details). The main numerical difficulty in using differential relations provided by elimination algorithms, comes from the presence of derivatives, sometimes of high order, which must be estimated from noisy measurements. Several methods have been used in the literature to obtain new relations less sensitive to the noise (see [2][3][4][5]).
In this paper, the main idea is the use of iterative integrals leading to ID relations, sometimes with no derivative. These are used in a bounded-error framework, in order to estimate the model unknown parameters. A bounded-error framework means that all uncertainties (measurement noise, parameters) are considered unknown but belonging to bounded connected sets, and, in the proposed work, intervals. Notice that the least squares method, adapted to the set-membership framework, could be used. However, with intervals, it requires the inversion of interval matrices, and an interval matrix is invertible if all punctual matrices in this interval matrix are invertible, which is a very hard constraint.
The dynamic systems considered in this work are given by the following form: ẋ(t, θ) = f (x(t, θ), θ) + u(t)g(x(t, θ), θ), y(t, θ) = h(x(t, θ), θ). (1) x(t, θ) ∈ R n represents the state variables, y(t, θ) ∈ R m represents the model outputs, and u(t) ∈ R l represents the input variables. u can be considered equal to 0, in the case of uncontrolled systems. The set of model parameters to be estimated are given by θ ∈ U p (U p is an open subset in R s ). The functions f (., θ), g(., θ) and h(., θ) are real, rational and analytic, for every θ ∈ U p on M (a connected open subset of R n , such that x(t, θ) ∈ M for every θ ∈ U p and every t ∈ [0, T]). In the following work, we let x 0 = (x 0,i ) i=1,...,n , the vector of initial conditions for x(t, θ); note that some components can depend on the parameters to be estimated.
is the smallest connected box belonging to U p . A box is an interval vector (a vector with intervals components) and may, equivalently, be seen as a Cartesian product of intervals. In this framework, a real interval [u] = [u, u] is a closed and connected subset of R, where u (respectively, u) represents the lower (respectively the upper) bound of [u]. Similarly, an interval matrix [A] is a matrix with interval components (see [6], for more details on interval analysis).
Considering the order eliminating, first, the unobserved state variables, then the model outputs and the model parameters, the Rosenfeld-Groebner algorithm applied on system (1) provides input-output polynomials (IO polynomials). Then, using iterative integrals, integro-differential input-output polynomials (ID-IO) can be obtained and used to propose a parameter estimation procedure less sensitive to noise, compared to the initial IO polynomials [7]. Indeed, contrary to the latter, ID-IO polynomials may contain derivatives of lower order.
To illustrate the method, the Hindmarsh-Rose (HR) model, resulting from a simplification and a generalization of the Hodgkin-Huxley model, is considered (see [8,9]). The HR model was developed to better understand neuron activity, from a simpler model that can be studied mathematically (for example [10][11][12]). Its particularity is to reproduce different dynamics of neurons. For example, it presents a bifurcation with respect to its slow-fast parameter [13]. In order to recover the behaviors of a neuron from the HR model, several parameter-estimation procedures were proposed in the literature. In [14], nonlinear optimization is used and exploits the particular structure of the relevant cost function. In [12], the authors propose two approaches. The first one deals with a synchronization-based parameter estimation and a least squares problem, subject to constraints. The second one is based on adaptive observers as in [15]. This method aims to find a dynamical system, so that it synchronizes with the measured voltages from a real neuron. However, none of the cited papers considered the noise on the data and the impact of the size of the measurement noise domain on the estimated intervals. The first work evoking this question can be found in paper [5]. It deals with a method based on ID relations, to estimate, first, the HR model parameters, then the probability that the results permit to reproduce the correct behavior of the model output near an equilibrium point. The probability is calculated from the evaluation, M times, of the model and a classical floating-point method, to create stochastically disturbed measurements. To complete this study, the bounded-error framework is considered in the proposed paper. In contrast to [5], this article does not use a probabilistic interpretation of measurement noise. Instead, the system outputs are assumed to be disturbed by bounded uncertainty, with unknown probability distributions in their interior. To handle this task, the concept of interval arithmetic is used. It is employed to analyze the accuracy of the parameter identification algorithm, relying on the same integration relations as in [5]. The results of the proposed new approach are interval bounds for all parameters to be identified, in which the true parameters are located with 100% certainty, if the measurement noise is assumed to be bounded. This kind of result is impossible to be obtained, with the algorithm published in [5]. Moreover, these results allow for detecting the structural changes of the system dynamics with certainty, by checking whether the parameter values corresponding to the bifurcation point are included or not in the estimated parameter ranges.
This paper is organized as follows. In Section 2, the parameter estimation procedure based on the ID-IO polynomials in the bounded-error framework is explained, and Section 3 presents both the HR Model and the numerical results. Section 4 concludes the paper.

Differential Algebra
Parameter estimation is done through relations linking the inputs, outputs and parameters of the system (1). These relations are obtained using the Rosenfeld-Groebner algorithm, implemented in the package DifferentialAlgebra of Maple [1]. This algorithm provides relations, called IO polynomials, from an elimination order, consisting of eliminating unobservable variables. The IO polynomials have the following form, for i from 1 to m: and (m l,i ) 1≤l≤q i are differential polynomials, with respect to y and u. m i 0 = 0 and i from 1 to m. According to [16], the number of the relations is the number of observations. Afterwards, only one output is considered, and the index i is omitted to lighten the notations.

Estimation Procedure
A numerical method, deduced from (2), was proposed to estimate the unknown constant parameters in [17] and is first recalled.
In the numerical applications, the measurement y is supposed to be described by , and θ * represents the "true" parameter vector value. Denote {y k = y(t k ), k = 1, . . . , N}, the set of measurements at (t k ) 1≤k≤N and u k = u(t k ) the associated inputs.
The parameter vector θ belongs to Θ, where Θ is an interval vector. Consider Γ k (Θ), the associated expression of γ k (θ) defined in the polynomial (2), where θ is substituted by Θ. Then, the following system, whose interval vector Θ is unknown, can be deduced: Notice that (3) is linear, with respect to Γ 1 (Θ), . . . , Γ q (Θ). Solving the previous system comes back to solving 0 ∈ However, some derivatives of high order can be involved in the linear system in using elimination algorithms, since the IO polynomials are deduced from the model equations in using addition and multiplication by any polynomials in x, u, y and θ as well as differentiation in time. Integrating these relations will permit not only to decrease the derivative order but also to attenuate the structured noise, whose amplitude is unknown (see [18], for more details).
Afterwards, we present an improvement of the method proposed in [17], which is based on the integrated relations obtained from (2).
Let f , a real-valued function, and I ν ( f ), ν ≥ 0, the integrated function, such that Using the linearity of the integral, a new relation is obtained from P and can be rewritten: I ν (P) = I ν (m 0 (y, u)) + q ∑ l=1γ l (θ)I ν (m l (y, u)).
I ν (P) is called the ID-IO polynomial. The approximated value of I ν (P), by a numerical procedure at the measurement points t k (k = 1, . . . , N), will be denoted I k ν (P).
In the same way, as previously, evaluating the expression I ν (P) at each t k leads to the linear systems.
where In the numerical applications, System (4) will be solved using the algorithm SIVIA, as presented in the following section.

Interval Set Inversion
This section recalls the algorithm SIVIA (Set Inverter Via Interval Analysis), well known in the interval-analysis community. This algorithm, proposed in [19], leads to characterize the solution set of a system of nonlinear real constraints, by enclosing it between internal and external unions of interval boxes (pavings).
Consider the problem of determining a solution set S for the unknown quantities u, belonging to an a priori search set U, defined by: where [y] is a priori known, and f is a nonlinear function, not necessarily invertible in the classical sense. (5) involves computing the reciprocal image of f and is known as a set inversion problem, which can be solved using the algorithm. SIVIA is a recursive algorithm, which explores all the search space without losing any solution. This algorithm makes it possible to derive a guaranteed enclosure of the solution set S, as follows: The inner enclosure S is composed of the boxes that have been proven feasible. To prove that a box [u] is feasible, it is sufficient to prove that f ([u]) ⊆ [y]. Reversely, if it can be proven that f ([u]) ∩ [y] = ∅, then the box [u] is unfeasible. Otherwise, no conclusion can be reached, and the box [u] is said to be undetermined. The latter is, then, bisected and tested, again, until its size reaches a user-specified precision threshold ε > 0. Such a termination criterion ensures that SIVIA terminates after a finite number of iterations.
Thus, the algorithm SIVIA allows to obtain these two subpavings, with a required precision ε, based on an inclusion test. The relation between the two subpavings can be characterized as: where ∆S is called the inclusion test uncertainty, in which no decision can be made during the test. The properties of solutions are: • if S = ∅ the problem (5) has no solution; • if S = ∅, there exists at least one verified solution for (5). A further alternative to the parameter identification proposed above consists of using interval methods, with a subsequent subdivision of parameter domains, in order to reliably identify the implausible parameter subintervals. This method, however, is much more computationally demanding and may, significantly, be affected by the wrapping effect of interval analysis, if specific properties, such as cooperativity, are not satisfied by [20] or [21].

Hindmarsh-Rose Model
The model of Hindmarsh-Rose (HR) results from a simplification and a generalization of the Hodgkin-Huxley model [8,9]. From this slow-fast model, rich dynamics of a neuron can be reproduced, such as spiking, bursting and chaotic behaviors. The HR model [9] reads as follows ( θ = (a, b, d, ε) T ): where • x 1 describes the membrane potential; • x 2 is the recovery variable, associated with the fast current, due to the passage of the Na + or K + ions; • x 3 is the adaptation current, associated with the slow current, due to the passage of the Ca + ions.
The variable supposed to be observed is the membrane potential, and we denote y = x 1 . u corresponds to the applied current (in amperes), and, afterwards, it is supposed constant. u generates the opening or closing of ion channels at one point in the membrane, which produces a local change in the membrane potential. Notice that the experimental data can be obtained in vivo, by using the current stimulus to generate a potential difference (see [22,23], for more details). Parameters a, b and d are, experimentally, determined from measurements of membrane potentials, while c x 1 is the x 1 -coordinate of the leftmost equilibrium of the two-dimensional system, given by the first two equations of (8), when u(t) = 0 and x 3 (t, θ) = 0, thus, θ = [a, b, d, ] T . Finally, parameter ε represents the ratio of time scales between fast and slow fluxes, across the membrane of a neuron, and controls the speed of variation of the slow variable x 3 .
It has been proven in [5] that the proposed model is globally identifiable (in a stochastic framework), thus, the identifiability property and, consequently, uniqueness of the parameters is maintained in a set-membership framework [24]. Relying on global identifiability, we know that this interval can be as small as possible (due to the threshold chosen in the bisection algorithm: SIVIA).
This model permits to obtain different dynamics, with respect to the parameter values. For example, for a = 3, b = 4, d = 5, u(t) = 3.25A, the authors of [13] prove that the parameter ε presents a Hopf bifurcation. When a Hopf bifurcation occurs, a local periodic solution near an equilibrium point appears or disappears, with the change of one parameter value. The time series of system (8), for the two parameter values ε = 0.12 and ε = 0.13, are represented in Figure 1. The chosen case is the second one, which is ε = 0.20.

ID-IO Polynomial
The package DifferentialAlgebra of Maple is used, in order to obtain the IO polynomial of the HR model. u is supposed to be an input of the system. To eliminate the variables x 1 , x 2 and x 3 as well as to acquire relations between y, u and θ, the elimination order [θ] ≺ [y, u] ≺ [x 1 , x 2 , x 3 ] is chosen. The Rosenfeld-Groebner algorithm provides a polynomial with a derivative of order 3, with 24 expressions. To obtain a simpler polynomial, let In integrating the second Equation of (8), one gets x 2 (t, θ) = (x 2 (0, θ) − 1)w 1 (t) + 1 − dv 1 (t). The following system is, now, considered: System (8) completed with the initial conditions (x 1 (0, θ), x 2 (0, θ), x 3 (0, θ)) is equivalent to System (9) completed with (x 1 (0, θ), Since the estimation of derivatives from noisy measurements is an ill-posed problem, some technicals were proposed to decrease the derivative orders of these polynomials. The most natural one is the integration of this IO polynomial, and the use of integration by parts.
The following relation, which does not contain any derivative, is, also, obtained: Evaluating this ID-IO polynomial at each t k , k = 1, . . . , N gives the following linear system.
In [5], the estimated parameters are obtained in solving the linear systemÃx =b with the QR factorization, whereas, we propose in this paper to estimate the enclosure of the parameters in solving (10) with the SIVIA algorithm. The important differences, with respect to the previous paper, are summarized in the two structure diagrams 2. In the left diagram, the stochastic procedure is presented. M is the number of iterations, to calculate the probability that the system reproduces the expected behavior of the model, given a numerical procedure and a noise. The right diagram summarizes the method presented in this paper, to estimate an enclosure of the parameters.
Let a parameter vector be given, for which the model solution is near an equilibrium point, such that one of the parameters is a bifurcation parameter. In our case, ε is the bifurcation parameter (see Section 3). The two works aim, also, to detect if the parameter value obtained by an optimization procedure leads to the expected behavior of the model output or not. The first method is based on the use of probabilistic tools. In contrast to this, the proposed method, based on interval arithmetic, is used to analyze the accuracy of the parameter identification algorithm. This new approach provides interval bounds, for all parameters to be identified, and in which the true parameters are located with 100% certainty, since the measurement noise is assumed to be bound. Consequently, this method allows for detecting, with certainty, the structural changes of the dynamics, by checking whether the parameter values corresponding to the bifurcation point are included or not in the estimated parameter ranges.
(a) (b) Figure 2. These two diagrams summarize the differences between the procedure presented in [5] and the one presented in this paper. (a) Diagram summarizing the stochastic procedure, presented in [5].
(b) Diagram summarizing the interval procedure, presented in this paper.

Parameter Estimation
In this section, enclosures of different integrals in the set-membership framework are obtained, by the interval extension of the trapezoidal classic method.
The integrals are evaluated, by using 29 points. In the output equation, [η(t)] is given by three successive intervals: [0.00008, 0.00012], [0.0004, 0.0006] and [0.0005, 0.0015], then System (10) is solved by using SIVIA's algorithm (see Section 2.3). For each instance of this value, the estimated intervals ofã,b,d,ε and the widths of the estimated intervals are given in Tables 1 and 2. [η(t)]εãbd All the intervals in Tables 1 and 2 contain the true values. Tables 1 and 2 highlight the impact of the width of the measurement noise interval on the width of the estimated; the larger the size of this interval is, the larger the size of the estimates.

Conclusions
In this paper, a new approach to estimate unknown parameters of the HR model is considered in the set-membership framework, to detect, with certainty, a behavior change in the dynamic of the system. Indeed, unlike probabilistic methods, which only give the probability that the estimated parameter leads to the expected behavior, the method we propose certifies whether the parameter values corresponding to the bifurcation parameter are included or not in the estimated parameter ranges. It takes the benefit of the differential algebra-based method, integration and the SIVIA algorithm. Numerical simulations highlight the interest of the proposed approach, in terms of estimated interval widths, compared to the previous approaches that are based on IO polynomials. A sensitivity analysis of estimates to the noisy data will be the subject of future works, as this could, indeed, improve the estimation of the parameters, based on the works of [25].