state and parameter estimation for canonic models ... - Semantic Scholar
June 5, 2010 9:54 00235
International Journal of Neural Systems, Vol. 20, No. 3 (2010) 193–207 c World Scientific Publishing Company DOI: 10.1142/S0129065710002358
STATE AND PARAMETER ESTIMATION FOR CANONIC MODELS OF NEURAL OSCILLATORS IVAN TYUKIN Department of Mathematics, University of Leicester University Road Leicester, LE1 7RH United Kingdom [email protected] ERIK STEUR∗ and HENK NIJMEIJER† Department of Mechanical Engineering Eindhoven University of Technology Eindhoven, P.O. Box 513 5600 MP The Netherlands ∗ [email protected] † [email protected] DAVID FAIRHURST Department of Mathematics, University of Leicester University Road, Leicester, LE1 7RH United Kingdom [email protected] INSEON SONG‡ and ALEXEY SEMYANOV§ Neuronal Circuit Mechanisms Research Group RIKEN Brain Science Institute Wako-shi, 351-0198, Japan ‡ [email protected] § [email protected] CEES VAN LEEUWEN Laboratory for Perceptual Dynamics RIKEN Brain Science Institute Wako-shi, 351-0198, Japan [email protected] We consider the problem of how to recover the state and parameter values of typical model neurons, such as Hindmarsh-Rose, FitzHugh-Nagumo, Morris-Lecar, from in-vitro measurements of membrane potentials. In control theory, in terms of observer design, model neurons qualify as locally observable. However, unlike most models traditionally addressed in control theory, no parameter-independent diffeomorphism exists, such that the original model equations can be transformed into adaptive canonic observer form. For a large class of model neurons, however, state and parameter reconstruction is possible nevertheless. We propose a method which, subject to mild conditions on the richness of the measured signal, allows model parameters and state variables to be reconstructed up to an equivalence class. Keywords: State and parameter estimation; universal observers; neural oscillators; nonlinear parameterization; unstable convergence; non-uniform small-gain theorem.
193
June 5, 2010 9:54 00235
194
1.
I. Tyukin et al.
Notations and Nomenclature
• The Euclidian norm of x ∈ Rn is denoted by x. • Let h : Rn → R be a differentiable function and f : Rn → Rn . Then Lf h(x), or simply Lf h, is the Lie derivative of h with respect to f : Lf h(x) = ∂h ∂x f (x) • Let f, g : Rn → Rn be differentiable vector-fields. Then the symbol [f, g] stands for the Lie bracket: ∂g [f, g] = ∂f ∂x g − ∂x f . • The adjoint representation of the Lie bracket is defined as ad0f g = g, adkf g = [f, adfk−1 g] • Let A be a subset of Rn , then for all x ∈ Rn , we define dist(A, x) = inf q∈A x − q. • Finally, let ∈ R>0 , then x stands for the following: x − , x > , x = 0, x ≤ . 2.
Introduction
2.1. Motivation and state-of-the art The brain consists of billions of neurons, each of which contains thousands of ion channels. This sheer complexity poses insurmountable difficulties towards full-scale in-vivo experiments and analysis. To overcome these issues, mathematical modeling combined with large-scale computer simulation of brain function is an inexorable tool.14 Successful examples include but are not limited to the simulation of memory functions,5 mechanisms of phase-resetting in olivo-cerebellar networks,15 and demonstrating the role of refractory period for synchronization in neural populations.23 Many model neurons to date share the following two features: they are systems of differential equations describing the response of a biological neuron to stimulation; their parameters characterize variables such as time constants, conductances, response thresholds, etc. To apply these models for describing the specific behavior of a biological neuron, their parameter values need to be determined from available data. Fitting parameters of systems of nonlinear differential equations, however, is a non-trivial operation. a
Even in the case of a system of ordinary differential equations it is recognized as a hard computational problem6 that “has not been yet treated in full generality”.17 Hand-tuning and automatic but heuristic fitting of model neurons to measured data conventionally involves trial-and-error or exhaustive search methods26 and thus suffers from the curse of dimensionality; only a few model parameters can be varied simultaneously, and the variation takes place over a coarse grid. Such coarseness leads to non-uniqueness of signal representation and the inability to detect subtle changes in the cells. Rigorous methods have been developed for computationally efficient state and parameter estimation in various classes of ordinary differential equations in system identification and control theory.3, 16 Exponentially fast on-line reconstruction is achievable for linearly parameterized systems that can be transformed into canonic adaptive observer form.a,4 State-of-the-art methods for parameter estimation of model neurons, e.g. Refs. 6 and 11, require that the model can be reduced to this canonic representation with linear parametrization. Typical models of spike generation, however, do not belong to the class of linearly parameterized systems. They belong to a more general class of equations13 : ϕj (v, t)pj (r)θj + I(t); v˙ = j
r˙i = −ai (v, θ, t)ri + bi (v, θ, t);
(1)
θ = (θ1 , θ2 , · · ·). Variable v in (1) is the membrane potential, and ri are the currents the values of which are not available for direct observation. Functions ai (·), bi (·), pj (·) are known, yet ai (·), bi (·) may depend on the unknown parameter vector θ. In light of the afore-mentioned computational difficulties in genuine search-based approaches, we wish to pursue observer-based solutions further, because of the possibility of exponential convergence of observer-based estimation even though equations (1) are not in the adaptive observer canonic form. Attempts have previously been made to develop observer-based methods for inferring parameters of
A canonic adaptive observer form is the following class of equations: x˙ = Ax + ϕ(y)θ + g(y)u, y = Cx, where x ∈ Rn is the state variable, y is the measured output, θ is the vector of unknown parameters, and u is the input (usually defined as a piece-wise continuous function). The function ϕ is known and continuous, matrices A, C are known and satisfy the following observability rank condition: rank [C T , AT C T , . . . , (AT )n−1 C T ] = n.
June 5, 2010 9:54 00235
State and Parameter Estimation for Canonic Models of Neural Oscillators
(1) from measured data such as membrane potential, v, and input current, I. For example, in Ref. 29, heuristic algorithms were proposed to estimate parameters of Hindmarsh-Rose model neurons. Yet, as the authors concluded, there is no systematic way to date for designing observers for parameter estimation of even the simplest model neurons. Developing techniques for state and parameter estimation with guaranteed convergence properties is the main goal of this study. 2.2. Contribution and structure of the paper We consider the following class of equations: x˙ 0 = θ0T φ0 (x0 , t) +
n
xi + ξ0 (t)
i=1
x˙ i = −τi xi + θiT φi (x0 , pi , t) + ξi (t), y = x0 ,
(2)
195
universal estimator that is capable of reconstructing unknown parameters of system (2) from the values of x0 (t) when t → ∞. Our method requires neither differentiation of measured evoked potentials, x0 (t), nor availability of internal currents for direct observation is required. It ensures not only asymptotic reconstruction of unknown model parameters, but also allows estimating the accuracy of fit as a function of the “richness” of the available data. Here richness is defined as a simple model-dependent condition on the non-singularity of a matrix that can easily be constructed from the values of x0 (t) and knowledge of the functions φi (·) in the right-hand side of (2). The paper is organized as follows. In Section 3 we present relevant concepts used to tackle the problem of parameter estimation, Section 4 contains main results, in Section 5 we provide an illustrative example, and Section 6 concludes the paper.
xi (t0 ) = xi,0
where φi : R × R≥0 → Rdi , di ∈ N/{0}, i = {0, . . . , n} are continuous functions modelling dynamic conductances of ion channels. Variable x0 in system (2) models the dynamics of the cell’s membrane potential, and variables xi , i ≥ 1 correspond to the ion currents in the cell. Parameters θi ∈ Rdi , pi ∈ R characterize membrane conductances. Parameters τi ∈ R>0 are the time constants of the currents. Variables ξi (t) constitute the unmodeled dynamics, or noise. The values of variable x0 (t) are assumed to be available for direct observation, and functions φ0 (x0 , t), φi (x0 , pi , t) are known. The variables xi (t), i = {1, . . . , n} are assumed to be not available for direct observation. Furthermore, true values of parameters θ0 , . . . , θn , τ1 , . . . , τn , p1 , . . . , pn are allowed to be unknown. We assume, however, that domains of admissible values of θi , pi , τi are known. In particular, we consider the case when θi,j ∈ [θi,min , θi,max ], τi ∈ [τi,min , τi,max ], pi ∈ [pi,min , pi,max ] and the values of θi,min , θi,max , τi,min , τi,max , pi,min , pi,max are available. Let us denote θ = (θ0 , θ1 , . . . , θn ),
3.
General Concepts of Parameter Inference for Ordinary Differential Equations
We will first specify the evaluation criteria needed for fitting the model parameters to data. Model neurons are systems of differential equations. Solutions of these systems are functions of space and time, and hence can be characterized by (1) their global asymptotic properties (e.g. by sets to which these solutions converge with time), and (2) by their local behavior in space and time. Generally, neither global nor local criteria alone give satisfactory results, especially if measurements are corrupted by noise (see Fig. 1). Plane (x0 , x1 ) in Fig. 1 represents the system’s state space. Black curves correspond to trajectories generated by the model, and gray curves depict
x0
x1
λ = (τ1 , p1 , . . . , τn , pn ),
and domains of admissible values for θ, λ are denoted by symbols Ωθ and Ωλ respectively. In this article we provide a method for systematic design of algorithms for asymptotic parameter reconstruction of model neurons that belong to a broad subclass of models (1). In particular, we present a
o
t1 t2 Fig. 1. time.
t3
t
Diagram of model and data trajectories in space-
June 5, 2010 9:54 00235
196
I. Tyukin et al.
a sample of actual data. Plane (x0 , t) is the space of observables, and x1 is a variable of which the values are not available for direct observation. In both planes there are points (marked by filled circles) at which distinct orbits in the space-time have identical or close projections on the state space, and viceversa. Hence neither the (x0 , x1 ) nor the (x0 , t) plane alone characterize data or model satisfactorily. 3.1. Observers Design of observers is the principal concept for asymptotic reconstruction of state and parameters of dynamical systems from knowledge of the history of trajectories and signals. In our case these trajectories and signals are the values of membrane potential, v(t), and input current I(t). In adaptive observer design, the original reconstruction problem is stated as follows. Let us suppose that, in compact notation, model (2) of evoked responses of a cell to stimulation I(t) is given by x˙ = g(x, I(t), θ, λ), θ˙ = 0, λ˙ = 0, v(t) = h(x(t)), and v(t), θ, λ are the measured potential and parameter vectors of which the values are unknown. Then the problem of inferring inaccessible variables of the model is viewed as a matter of finding a computational procedure in the form of a system: z˙ = f (z, t, I(t), v(t)),
ζ0 = z(t0 ) ∈ Ωz ,
(3)
such that inaccessible state variables and parameters are tracked by known functions of the internal state z of the observer. In the case of a state and parameter reconstruction problem, we will seek for (3) and funcˆ ˆ tions x ˆ(z), λ(z), θ(z) such that for sufficiently small given constant εx , ελ , εθ > 0, the following properties hold: lim sup ˆ x(z(t, ζ0 )) − x(t) ≤ εx , t→∞
ˆ lim sup λ(z(t, ζ0 )) − λ ≤ ελ , t→∞
ˆ ζ0 )) − θ ≤ εθ lim sup θ(z(t,
(4)
t→∞
∀ ξ0 ∈ Ωz . along the trajectories of a coupled model (data) observer system: x˙ = g(x, I(t), θ, λ) (5) z˙ = f (z, t, I(t), v(t)).
Traditionally xˆ(z) is viewed as an estimate of the ˆ ˆ are the estimates model state vector, x, and λ(z), θ(z) of parameter vectors λ, θ. The concept of observer-based estimation allows a general perspective on the problem of finding algorithms for fitting parameters of (1) to data. This perspective enables us to use the tools of analysis familiar from dynamical systems theory. In Section 4 we outline main ideas leading to the construction of observer-based parameter estimators for models (2).
3.2. Weakly attracting sets Current methods to determine algorithms (3), (4) (i.e. observers) are based on the notion of Lyapunov stability of solutions. If a globally Lyapunov-stable solution exists in the state space of (5) such that requirements (4) are satisfied, then (3) is the algorithm needed. Hence finding such algorithms may be achieved by employing the convenience offered by the method of Lyapunov functions. However, according to the results of simulation studies in Ref. 26, parameter values of general conductance models (1) corresponding to nearby trajectories, are not unique and are distributed across the whole parameter space. Therefore, if all we had available were the concepts of Lyapunov stability and the method of Lyapunov functions, we would have to deal with multiple separated attracting invariant sets. In such systems unstable sets inevitably occur at the boundaries of the attractor basins. Hence only local stability is generally guaranteed. It is known, however, that the method of Lyapunov functions is conservative when it comes to the local analysis in the problem of observer design for (1): its success and accuracy depend on how good the choice of the Lyapunov candidate is. To overcome locality and non-uniqueness issues, it is advantageous to replace the concept of Lyapunov stability with the more relaxed concepts of convergence to weakly attracting sets 2, 22 and relaxation times.8 A diagram illustrating the benefits of this less restrictive notion is provided in Fig. 2. In Fig. 2 Panel a depicts a diagram of standard stable attractors, panel b illustrates weak attractors, and panel c shows observer concept for (1). Domains of stable attractors are neighborhoods containing A1 , A2 . Estimates of sizes of these domains depend on parameters S0,1 , S0,2 , S0,3 . These estimates are
June 5, 2010 9:54 00235
State and Parameter Estimation for Canonic Models of Neural Oscillators
S 0,1 S 0,2 S 0,3 A1
A2
(a)
S 0,1 S 0,2 S 0,3 A1 A2
eo ac sp eters b Su ram pa
o pl Ex
ra
tio
n
Con
tr a c
Subspa ce of linea r pa r a meter s
ar
ne
nli
o fn
Subspace of linear p a r a meter s
(b)
tion
A
Con
tr a c
tion
t1
A
t2
(c)
Fig. 2.
Attractor and observer concepts.
depicted as closed curves around A1 , A2 . Domains of attraction for weak (Milnor) attracting sets are not neighborhoods. Hence, even in case of an unfortunate initial guess trajectories will simply converge to another (alternative) attractor via unstable paths. As part of our theoretical framework, we exploit the notion of weakly attracting sets in the problem of observer design and advocate the use of associated analysis results31 to solve this problem, as a useful complementary tool to the method of Lyapunov functions. 3.3. A-posteriori bifurcation analysis As follows from Refs. 1, 31, 28 and 7, parameter reconstruction of model neurons from data is b
197
achievable in principle. When the model is capable of reproducing data without errors, our algorithm ensures that model parameters can be reconstructed such that, over a given finite time interval, mismatches between trajectories of the model and data could be made negligibly small.b This case is equivalent to the assumption that the unmodeled dynamics, ξi (t), in (2) are negligible or non-existent. In practice, such a case will hardly ever be encountered; in general, non-zero mismatches between assumed model and data are to be expected. Model solutions are generally sensitive to parameter values; small perturbations of these values will result in models exhibiting different regimes with distinct invariant sets and asymptotic properties. Because of this, the presence of unmodeled dynamics may result in estimates that fail to generate plausible trajectories. Consider Class-I excitable model neurons, a class to which all the following model neurons belong: Morris-Lecar, Hindmarsh-Rose, FizhughNagumo, and modified Hodgkin-Huxley models.13 In these models the spiking mechanism is described by the homoclinic or saddle-on-invariant circle bifurcation. Figure 3 illustrates this mechanism for the homoclinic bifurcation. A model close to the bifurcation point is generating low-frequency oscillations (Fig. 3(c)). As shown in Fig. 3(b), a small perturbation of the excitation parameter may lead to abrupt changes of the solutions: instead of being periodic, almost all solutions will asymptotically converge to the stable equilibrium. The situation becomes even more complicated in practice for, according to most recent studies,27 ,25 popular models such as Morris-Lecar or HindmarshRose inherit more than one type of excitability. This means that not only the topology of the invariant sets changes with parameters (e.g. resting vs spiking) but also the spiking mechanisms become different. Hence in general, in order to account for model feasibility regions a-posteriori tailoring of the results of estimation is needed. The most straightforward solution to this problem is to map the results of estimation onto those domains in which the estimates are feasible. For planar models (or those reducible to this class) the regions of feasibility can be located by examining
This is a consequence of the continuity property of solutions of ordinary differential equations with respect to parameters and initial conditions.
June 5, 2010 9:54 00235
198
I. Tyukin et al.
3
3
3
2
2
1 1 (a)
(b)
(c)
Fig. 3. Typical transition of invariant sets in class-I model neurons. When the excitation parameter value is sufficiently small (Panel a), the system has three distinct equilibria: a stable one labeled 1 and two unstable ones labeled 2 and 3, respectively. Increasing the excitation parameter value (Panel b) leads to the collision of stable and unstable Equilibria 1 and 2, such that their stable and unstable manifolds coincide. Further increase (Panel c) gives rise to stable periodic orbits with variable oscillation periods.
local stability of model equilibria.10, 13 Projection of estimated parameters onto the feasibility regions corresponding to typical bifurcations (such as e.g. b in Fig. 3) can be achieved by searching for the feasible values of model parameters in the vicinity of the ˆ λ. ˆ obtained estimates θ, In order to do so, it is possible to identify the values of model parameters corresponding to the specific bifurcation that leads to unfolding of the desired dynamic regime of the system; in our case, the homoclinic loop (fold) bifurcation. Here we show that, at least for the Hindmarsh-Rose model neurons, MonteCarlo search procedures can be used successfully to locate the values of model parameters corresponding to the bifurcation. 4.
Observer-based Estimation of Unknown Model Parameters
Let us consider the following class of dynamical systems x˙ = f (x, θ) + g(x, θ)u(t), y = h(x),
x(t0 ) ∈ Ωx ⊂ Rn
x ∈ Rn , θ ∈ Rd , y ∈ R
(6)
where f, g : Rn → Rn , h : Rn → R, u : R → R are continuous and differentiable functions. Variable x is the state vector, u is the known input (representing e.g. the stimulation current I(t) in (1)), θ is the vector of unknown parameters, and y is the output of (6). Recall that y(t, 0, x0 , u) = h(x(t, 0, x0 , u)) denotes the output of (6) for u(·) and initial state x(0) = x0 . System (6) includes equations (1) as a subclass and in this respect can be considered
as a plausible generalization. Obviously, conclusions about (6) should be valid for systems (1) as well. Definition 1 (Observability24 ). Two states x1 , x2 ∈ Rn are said to be indistinguishable (denoted by x1 Ix2 ) for (6) if for every admissible input function u the output function t → y(t, 0, x1 , u), t ≥ 0 of the system for initial state x(0) = x1 , and the output function t → y(t, 0, x2 , u), t ≥ 0 of the system for initial state x(0) = x2 , are identical on their common domain of definition. The system is called observable if x1 Ix2 implies x1 = x2 . According to Definition 1, observability of a dynamical system implies that the values of its state, x(t) ∈ Rn , t ∈ [t1 , t2 ] are completely determined by inputs and outputs u(t), y(t) over [t1 , t2 ], although this definition of observability does not imply that every input function distinguishes points of Rn . Remark 1. Although Definition 1 does not account for the observability of any unknown parameter vectors, one can easily see that the very same definition can be used for parameterized systems as well. Indeed, extending the original equations (6) by including the parameter vector θ as a component of the extended state vector x ˜ = (xT , θT )T results in x ˜˙ = f˜(˜ x) + g˜(˜ x)u(t), ˜ y = h(˜ x), x˜(t0 ) = (x(t0 ), θ) ∈ Ωx˜ ⊂ R T
(7) n+d
where f˜(˜ x) = (f (x)T , 0)T , g˜(˜ x) = (g(x)T , 0)T , and T T ˜ x) = (h(x) , 0) . All uncertainties in (6), includh(˜ ing the parameter vector θ, are now combined into
June 5, 2010 9:54 00235
State and Parameter Estimation for Canonic Models of Neural Oscillators
the state vector of (7). Hence the problem of state and parameter reconstruction of (6) can be viewed as that of recovering the values of state for (7). Throughout this paper we shall consider input u as a constant (this corresponds I(t) = const in (1)), and combining the unknown parameter vector with the state vector as described in remark 1 enables us to write the system in the following form x˜˙ = f˜(˜ x), ˜ x), y = h(˜
x ˜ ∈ Rn˜ , y ∈ R
(8)
24
Proposotion 1 (Local observability test ). The system y = h(x), x ∈ Rn , y ∈ R
(9)
with smooth functions f, h is observable locally to the point x∗ if the following rank condition holds at that point rank
∂ (h(x∗ ) Lf h(x∗ ) · · · Ln−1 h(x∗ ))T = n f ∂x (10)
In7 it is shown that typical models neurons (Hindmarsh-Rose, Morris-Lecar) are locally observable systems.
The system with output measurements to be observed, state estimator (11), and parameter adaptation mechanism (12) are given below x(t) ˙ = A1 x(t) + φ0 (y(t), u(t)) + b
βi (y(t), u(t))θi
i=1
y(t) = C1 x(t) 0 1 0 ··· 0 0 1 ··· A1 = · · · ··· 0 0 0 ··· 0 0 0 ···
C1 = 1 0 0 · · · 0 x(t) ∈ R , n
(12)
K=
1 (A1 b + λb) = (k1 , . . . , kn )T b1
˙ ˆ) sign (b1 ) θˆ = Γβ(t)(y − C1 x
(13)
with Γ an arbitrary symmetric positive definite matrix and λ an arbitrary positive real. The n × 1 vector, b, is Hurwitz with b1 = 0 (see also21 where an observer is presented for systems in which the vector b is an arbitrary vector from Rn ). Theorem 1.18 There exists a local change of coordinates, z = Φ(x), transforming x˙ = f (x) +
p
θi (t)qi (x),
i=1 n
y = x1
x ∈ R , y ∈ R, θi ∈ R, qi : Rn → Rn , n ≥ 2
z˙ = A1 z + ψ(y) +
p
θi (t)ψi (y),
i=1
y = C1 z
z ∈ Rn , ψi : R → Rn
with (A1 , C1 ) in canonical observer form (11), if and only if (i) [adif g, adjf g] = 0, 0 ≤ i, j ≤ n − 1 (ii) [qi , adjf g] = 0, 0 ≤ j ≤ n − 2, 1 ≤ i ≤ p
0 0 · , 1 0
y(t) ∈ R,
x(t) + φ0 (y(t), u(t)) x ˆ˙ (t) = (A1 − KC1 )ˆ p βi (y(t), u(t))θˆi + Ky(t) +b
with h(xo ) = 0 and (f, h) an observable pair, into the system
4.1. Canonical adaptive observers
p
known, bounded and piecewise continuous functions of y(t), u(t). The column vector b ∈ Rn is assumed to be Hurwitz with b1 = 0 (e.g. all roots of the polynomial b1 pn−1 + · · · + bn−1 p + bn should have negative real part). An adaptive observer as presented in Ref. 19 is given as
i=1
Systems (8), in turn, admit the following simple observability test:
x˙ = f (x),
199
where the vector field, g(x), is uniquely defined by ∂ n−1 T (h(x), Lf h(x), . . . , Lf h(x)) , g(x) ∂x T
= (0, 0, . . . , 1) βi (·, ·) : R × R → R (11)
In (11) x(t) ∈ Rn is the state vector, and y = x1 (t) is the output; θ ∈ Rp = (θ1 , . . . , θp )T is the vector of unknown parameters. The functions βi are
(14)
Consider, for example, the problem of fitting parameters of the conventional HindmarshRose oscillator to measured data. It is shown in Refs. 28 and 7 that the Hindmarsh-Rose model, after translation and stretching of the state space and time
June 5, 2010 9:54 00235
I. Tyukin et al.
200
axis, can be written in vector-matrix notation as follows:
x1 x˙ 1 (15) = A(λ) + ϕ(x1 )θ x˙ 2 x2
cannot be transformed by diffeomorphic change of coordinates, z = φ(x), into z˙ = A1 z + ψ0 (y) +
A(λ)
=
ϕ(x1 ) θ
=
0 1 0 −λ x31
x21
0
0
0
1
0
0
0
0
x31
x21
x1
(16)
= (θ1,3 , θ1,2 , θ1,1 , θ1,0 , θ2,3 , θ2,2 , θ2,1 )
One of the main issues is that the original equations of neural dynamics are not written in a canonical form for which the reconstruction algorithms are available. The question, therefore, is if there exists an invertible coordinate transformation such that the model equations can be rendered canonic. Below we demonstrate that this is generally not the case if the transformation is parameter-independent. However, if we allow our transformation to be both parameter and time-dependent, a relevant class of models with polynomial right-hand sides can be transformed into one of the canonic forms. Let us consider a class of systems that can be described by (15). Clearly this system is not in a canonical adaptive observer form because A(λ) depends on the unknown parameter explicitly (see (16)). The question, however, is if there exists a differentiable coordinate transformation z = Φ(x), z1 = x1 such that in the new coordinates the equations of system (15) satisfy the canonic descriptions. We show that the answer to this question is negative, and it follows from the following slightly more general statement Theorem 2.7 The system x˙ = f (x) +
p
z ∈ Rn ,
y ∈ R,
ψi : R → Rn
(17)
with A1 , C1 in canonical observer form (11), if either (i) n > 2 or (ii) there exists i ∈ {1, . . . , p}, j ∈ {2, . . . , n} such that ∂qi /∂xj = 0.
, x1
θi ψi (y),
i=1
y = C1 z
where
p
Let us now consider the case in which the transformation z = Φ(x, λ, θ, t) is allowed to depend on unknown parameters and time. As we show below, this class of transformations is much more flexible. In principle it allows us to solve the problem of partial state and parameter reconstruction for an important class of oscillators with polynomial right-hand side and time-invariant time constants. We start by searching for a transformation Φ : q = T (λ)x, |T (λ)| = 0 such that
1 (18) T (λ)A(λ)T −1 (λ) =
0 where the matrix A(λ) is defined as in (16). It is easy to see that the transformation satisfying this constraint exists, and it is determined by
1 0 T (λ) = . (19) λ 1 According to (18), (19) and (16) equations of (13) in the coordinates q can be written as q˙ = A1 q + ψ(q1 )η(θ, λ), where
A1 =
ψ(q1 ) =
0 1 0 0 q13 0
(20)
q12 0
,
(21) q1 0
1 0 0 q13
0 q12
0 q1
0 1
,
η = (θ1,3 , θ1,2 , θ1,1 − λ, θ1,0 , λθ1,3 + θ2,3 , λθ1,2 + θ2,2 , λθ1,1 + θ2,1 , λθ1,0 )T
θi qi (x),
Remark 2. Notice that
i=1
y ∈ R,
y = h(x) x ∈ R , n
qi : R → R , n
with
0 1 1 ··· 1 0 0 0 ··· 0 f (x) = · · · · · · · x, 0 0 0 ··· 0
h(x) = 1 0 0 · · · 0 x
n
n≥2
• availability of the parameter vector η in (20), (21), expressed as a function of θ, λ, implies the availability of θ, λ if θ1,0 = 0. Indeed, in this case the value of λ = η8 /η4 and the values of all θi,j are uniquely defined by ηi ; • condition θ1,0 = 0 is sufficient for reconstructing the values of x2 provided that q and η are available; indeed in this case x = T −1 (λ)q.
June 5, 2010 9:54 00235
State and Parameter Estimation for Canonic Models of Neural Oscillators
As follows from Remark 2 the problem of state and parameter reconstruction of (13) from measured data x1 (t) amounts to solving the problem of state and parameter reconstruction of (20). In order to solve this problem we shall employ yet another coordinate transform: z1 = q1 , z2 = q2 + ζ T (t)η
(22)
in which the functions ζ(t) are some differentiable functions of time. Coordinate transformation (22) is clearly time-dependent. The role of this additional transformation is to transform the equations of system (20) into the form for which a solution already exists. Define the vector-function ζ(t) in (22) as follows: ζ˙i = −kζi + kψ1,i (q1 ) − ψ2,i (q1 ),
k ∈ R>0
In this case we have z˙2 =
8
ψ2,i (q1 )ηi + k −
i=1
−
8
ψ2,i (q1 )ηi = k
i=1
8
ζi ηi + ψ1,i (q1 )ηi
i=1 8
(−ζi + ψ1,i (q1 ))ηi
i=1
Taking (22) into account and expressing q2 as q2 = z2 − ζ T (t)η(θ, λ) we obtain
z˙1 = z2 +
8
(−ζi + ψ1,i (q1 ))ηi
i=1
z˙2 = k
8
(23)
(−ζi + ψ1,i (q1 ))ηi
i=1
Notice that ψ1,4 (q1 ) = ψ2,8 (q1 ) = 1, hence −ζ4 + ψ1,4 (q1 ) and −ζ8 + ψ1,8 (q1 ) converge to constants in R exponentially fast as t → ∞. Moreover, the sum −ζ4 + ψ1,4 (q1 ) is converging to zero, and the sum −ζ8 + ψ1,8 (q1 ) is converging to −1/k as t → ∞. Taking these facts into account we can conclude that system (23) can be rewritten in the following (reduced) form z˙ = A1 z + bφT (z1 , t)υ(θ, λ) + bε(t),
0 1 1 A1 = , b= , 0 0 k
201
−ζ1 + ψ1,1 (z1 ) −ζ + ψ (z ) 1,2 1 2 −ζ + ψ (z ) 3 1,3 1 φ(z1 , t) = −ζ8 + ψ1,8 (z1 ) −ζ5 + ψ1,5 (z1 ) −ζ6 + ψ1,6 (z1 ) −ζ7 + ψ1,7 (z1 ) υ(θ, λ) = (θ1,3 , θ1,2 , θ1,1 , θ¯1,0 , λθ1,3 + θ2,3 ,
λθ1,2 + θ2,2 , λθ1,1 + θ2,1 )T ,
(24)
where ε(t) is an exponentially decaying term. System (24) is clearly in adaptive canonic observer form. Hence it admits the following adaptive observer z − z) + bφT (z1 , t)ˆ υ zˆ˙ = A1 zˆ + L(ˆ
−l1 0 (25) L = , l1 , l2 ∈ R > 0 −l2 0 z1 − z1 )φi (z1 , t), γ ∈ R>0 υˆ˙ i = −γ(ˆ of which the asymptotic properties are specified in.7 In particular, the state and parameters are shown to converge to their true value provided that ∃ L, δ ∈ R>0 : t+L φ(z1 (τ ), τ )φ(z1 (τ ), τ )T dτ ≥ δI,
(26)
t
The function φ(z1 (t), t) is determined by the measured data, and condition (26) can easily be tested. A detailed explanation on how to get the values of θ and λ from the values of υˆ can be found in Ref. 7. So far we have considered special cases of (1) in which the time constants of unmeasured variables were unknown and parametrization of the right-hand side was linear. In the next section we discuss an approach that for a large subclass of (1) allows us to produce an observer solving the problem of state and parameter reconstruction from membrane potential measurement. The structure of this observer does not depend significantly on specific equations describing dynamics of the observed system. For this reason, and similarly to Ref. 12, we refer to this class of observers as universal adaptive observers. Detailed definitions of this class of observers and formal statements regarding their asymptotic properties are presented in Ref. 30. 4.2. Universal observers The main idea of our observer-based estimation method is to design uniformly stable estimators for
June 5, 2010 9:54 00235
202
I. Tyukin et al.
the parameters θi in (2) if time-constants and nonlinear parameters are known, followed by an automated search for the remaining parameters.c Thus the observer comprises of two compartments: stable, contracting and exploring, weakly attracting (interplay between these contracting-exploring compartments is illustrated in Fig. 2, Panel c). Contracting compartment. Consider (2), if parameters τi and pi would be known then the model would be in the adaptive canonic observer form. For these systems, estimation procedures are available that allow exponentially fast convergence of estimates θˆi of θi to small neighborhoods of θi cf.4, 20, 21 Let us combine variables xi and vectors θi into a single vector as follows: z = (x0 , . . . , xn , θ0T , . . . , θnT ), and let zˆ be its estimate. Then, according to4 and under mild non-degeneracy assumptions, similar to the non-singularity requirement for a matrix composed of filtered measured data in Ref. 11, there will exist a computational scheme x ˆ˙ = f (ˆ x, x0 (t), I(t)) (27) zˆ = (χ0 (ˆ x), . . . , χn (ˆ x), θˆ0 (ˆ x), . . . , θˆn (ˆ x)), where f (ˆ x, x0 (t), I(t)), χi (ˆ x), θˆi (ˆ x), are known functions, such that the dynamics of z − zˆ will obey the following estimated : z(t) − zˆ(t) ≤ De−ρ(t−t0 ) z(t0 ) − zˆ(t0 ) + ∆. (28) The presence of ∆ in (28) is due to measurement noise, ξi (t), and D ∈ R>0 is a constant of which the value can be computed a-priori. Observe that property (28) implies that the unperturbed system, i.e. when ξi (t) = 0, has a unique asymptotically stable solution in the sense of Lyapunov. Hence, we call this computational scheme the stable compartment of our observer. When the values of pi , τi are unknown, i.e. the vector λ is not measured, we replace them with their ˆ = (ˆ estimates, λ τ1 , pˆ1 , . . . , τˆn , pˆn ). Then it is possible to show that the dynamics of z(t) − zˆ(t) will satisfy the following inequality: z(t) − zˆ(t) ≤ De−ρ(t−t0 ) z(t0 ) − zˆ(t0 ) ˆ ) + C∆, + B sup λ − λ(τ τ ∈[t0 ,t]
c
(29)
where B, C are positive constants (see Refs. 30 and 32). Exploring compartment. In order to reconstruct the values of λ we extend the estimation procedure (27) by adding a component which searches for the values of λ by exploring all points in the domain Ωλ to which the values of unknown λ are supposed to belong. This explorative search is realized by a dynamical system of which the solutions form a (dense) grid in Ωλ (see Ref. 30 for details). Briefly, its working can be explained as follows. Let, for instance, Ωλ be a unit interval, and let λ(h) be a dense trajectory in Ωλ parameterized by h.e A possible candidate ˆ for the dense trajectory in [0, 1] is λ(h) = (sin(h) + 1)/2 (see Fig. 4). This trajectory can be modeled (up to a constant phase shift) by the following system: ˆ 1 ) = (1 + h1 )/2; λ(h h˙ 1 = h2 h˙ 2 = −h1 , h21 (t0 ) + h22 (t2 ) = 1. ˆ 1) If no further adjustments are made then λ(h will periodically visit every point in [0, 1]. In order to be able to stop at the point that corresponds to x0 (t)) = 0, those values of λ for which limt→∞ (x0 (t)−ˆ where xˆ0 (t) is the estimate of x0 (t) provided by the observer, we make the speed of exploration along the chosen dense trajectory proportional to the absoˆ0 (t)|. lute value of the trajectory error |x0 (t) − x Thus, we impose the following regulatory feedback mechanism for the explorative part of our observer: ˆ0 (t)|, γ ∈ R>0 . Then the exploh˙ ∼ γ|x0 (t) − x rative compartment of the estimator can be modeled by the following system of ordinary differential equations: ˆ 1 ) = (1 + h1 )/2; λ(h
(30)
ˆ0 (t)|h2 h˙ 1 = γ|x0 (t) − x
(31)
h˙ 2 = −γ|x0 (t) − xˆ0 (t)|h1 ,
(32)
h21 (t0 ) + h22 (t2 ) = 1.
General structure, formal definition and asymptotic properties of universal observers are provided in Refs. 30 and 7. Because of space limitations we do not reproduce these results here. d We refer the reader to Ref. 30 for details of derivations needed to obtain expressions in this section. e Specifically, we require that for all λ ∈ Ωλ and any δ ∈ R>0 there is a diverging sequence s1 , s2 ,. . . such that λ − λ(si ) ≤ δ.
June 5, 2010 9:54 00235
State and Parameter Estimation for Canonic Models of Neural Oscillators
λ
>
λ
203
e
1 ∧ λ=
h2
λ
h’2 0
λ
h
h’1
h1
ˆ as a function of h. Trajectory λ(h) ˆ Fig. 4. Explorative part of the observer. On the left - λ is dense in [0, 1]. In fact, ˆ i ) = λ. On the right - a diagram of for any λ ∈ [0, 1] there will be an infinite set of points si , i = 1, 2, . . . such that λ(s ˆ0 (t). Two invariant sets are present in the phase space of (31), (32). Variable e denotes instantaneous difference x0 (t) − x ˆ 1 ) = λ holds. system, both corresponding to the values of h1 for which the target equivalence λ(h
The values of γ are to be chosen such that suitable estimates form (weakly) attracting sets in the system state space (see Fig. 2). The evidence that such a choice can actually be made and the estimates of domains to which the values of γ should belong are provided in Ref. 31. In the next section we present universal observers designed specifically for dealing with state and parameter estimation of typical model neurons considered in this article. In comparison with other search-based methods our estimation technique has two distinct features that make its application advantageous: (1) Explorative search is conducted in a lowerdimensional subspace (one - dimensional for 2D Hindmarsh-Rose models (33)). This property allows us to combine efficiently the arsenal of conventional Lyapunov-based tools for designing stable compartments of the observer with the flexibility of the more general concept of weakly attracting sets. The result is – a substantial improvement in computational performance from trial-and-error approaches. (2) Exploration in Ωλ occurs simultaneously with exponentially fast estimation of θ. The explorative search does not need to wait until the exponentially fast compartment converges. Hence the amount of time needed to explore the domain of model parameters is reduced even further.
5.
Example: Fitting the Hindmarsh-Rose Model Neuron to Data
In this section we illustrate how the observer-based approach for state and parameter estimation of neural oscillators can be implemented for fitting models of neural oscillators to actual electrophysiological data. The data we used for this demonstration is shown in Fig. 5 (lower panel). Notice that the measured spikes are pretty sharp and narrow, which makes using the Hindmarsh-Rose model a preferred option over other low-dimensional models such as the Fitzhugh-Nagumo or Morris-Lecar equations.9 Consider a two dimensional Hindmarsh-Rose model written in the following form: x˙ 0 = θ0,3 x30 + θ0,2 x20 + θ0,1 x0 + x1 + θ0,0 + I(t) x˙ 1 = −λ1 x1 + θ1,2 x20 + θ1,1 x0 + θ1,0 .
(33)
To reconstruct the dynamics of a single spike we use those trials in which the input function I(t) is a square impulse. Within the time span of a singe burst this input function can be replaced with a constant. The value of θ1,0 can be aggregated into the parameter θ0,0 . Thus, instead of (33), we obtain the following equations: x˙ 0 = θ0,3 x30 + θ0,2 x20 + θ0,1 x0 + x1 + θ0,0 (34) x˙ 1 = −λ1 x1 + θ1,2 x20 + θ1,1 x0 .
June 5, 2010 9:54 00235
204
I. Tyukin et al.
International Journal of Neural Systems, Vol. 20, No. 3 (2010) 193–207 c World Scientific Publishing Company DOI: 10.1142/S0129065710002358
STATE AND PARAMETER ESTIMATION FOR CANONIC MODELS OF NEURAL OSCILLATORS IVAN TYUKIN Department of Mathematics, University of Leicester University Road Leicester, LE1 7RH United Kingdom [email protected] ERIK STEUR∗ and HENK NIJMEIJER† Department of Mechanical Engineering Eindhoven University of Technology Eindhoven, P.O. Box 513 5600 MP The Netherlands ∗ [email protected] † [email protected] DAVID FAIRHURST Department of Mathematics, University of Leicester University Road, Leicester, LE1 7RH United Kingdom [email protected] INSEON SONG‡ and ALEXEY SEMYANOV§ Neuronal Circuit Mechanisms Research Group RIKEN Brain Science Institute Wako-shi, 351-0198, Japan ‡ [email protected] § [email protected] CEES VAN LEEUWEN Laboratory for Perceptual Dynamics RIKEN Brain Science Institute Wako-shi, 351-0198, Japan [email protected] We consider the problem of how to recover the state and parameter values of typical model neurons, such as Hindmarsh-Rose, FitzHugh-Nagumo, Morris-Lecar, from in-vitro measurements of membrane potentials. In control theory, in terms of observer design, model neurons qualify as locally observable. However, unlike most models traditionally addressed in control theory, no parameter-independent diffeomorphism exists, such that the original model equations can be transformed into adaptive canonic observer form. For a large class of model neurons, however, state and parameter reconstruction is possible nevertheless. We propose a method which, subject to mild conditions on the richness of the measured signal, allows model parameters and state variables to be reconstructed up to an equivalence class. Keywords: State and parameter estimation; universal observers; neural oscillators; nonlinear parameterization; unstable convergence; non-uniform small-gain theorem.
193
June 5, 2010 9:54 00235
194
1.
I. Tyukin et al.
Notations and Nomenclature
• The Euclidian norm of x ∈ Rn is denoted by x. • Let h : Rn → R be a differentiable function and f : Rn → Rn . Then Lf h(x), or simply Lf h, is the Lie derivative of h with respect to f : Lf h(x) = ∂h ∂x f (x) • Let f, g : Rn → Rn be differentiable vector-fields. Then the symbol [f, g] stands for the Lie bracket: ∂g [f, g] = ∂f ∂x g − ∂x f . • The adjoint representation of the Lie bracket is defined as ad0f g = g, adkf g = [f, adfk−1 g] • Let A be a subset of Rn , then for all x ∈ Rn , we define dist(A, x) = inf q∈A x − q. • Finally, let ∈ R>0 , then x stands for the following: x − , x > , x = 0, x ≤ . 2.
Introduction
2.1. Motivation and state-of-the art The brain consists of billions of neurons, each of which contains thousands of ion channels. This sheer complexity poses insurmountable difficulties towards full-scale in-vivo experiments and analysis. To overcome these issues, mathematical modeling combined with large-scale computer simulation of brain function is an inexorable tool.14 Successful examples include but are not limited to the simulation of memory functions,5 mechanisms of phase-resetting in olivo-cerebellar networks,15 and demonstrating the role of refractory period for synchronization in neural populations.23 Many model neurons to date share the following two features: they are systems of differential equations describing the response of a biological neuron to stimulation; their parameters characterize variables such as time constants, conductances, response thresholds, etc. To apply these models for describing the specific behavior of a biological neuron, their parameter values need to be determined from available data. Fitting parameters of systems of nonlinear differential equations, however, is a non-trivial operation. a
Even in the case of a system of ordinary differential equations it is recognized as a hard computational problem6 that “has not been yet treated in full generality”.17 Hand-tuning and automatic but heuristic fitting of model neurons to measured data conventionally involves trial-and-error or exhaustive search methods26 and thus suffers from the curse of dimensionality; only a few model parameters can be varied simultaneously, and the variation takes place over a coarse grid. Such coarseness leads to non-uniqueness of signal representation and the inability to detect subtle changes in the cells. Rigorous methods have been developed for computationally efficient state and parameter estimation in various classes of ordinary differential equations in system identification and control theory.3, 16 Exponentially fast on-line reconstruction is achievable for linearly parameterized systems that can be transformed into canonic adaptive observer form.a,4 State-of-the-art methods for parameter estimation of model neurons, e.g. Refs. 6 and 11, require that the model can be reduced to this canonic representation with linear parametrization. Typical models of spike generation, however, do not belong to the class of linearly parameterized systems. They belong to a more general class of equations13 : ϕj (v, t)pj (r)θj + I(t); v˙ = j
r˙i = −ai (v, θ, t)ri + bi (v, θ, t);
(1)
θ = (θ1 , θ2 , · · ·). Variable v in (1) is the membrane potential, and ri are the currents the values of which are not available for direct observation. Functions ai (·), bi (·), pj (·) are known, yet ai (·), bi (·) may depend on the unknown parameter vector θ. In light of the afore-mentioned computational difficulties in genuine search-based approaches, we wish to pursue observer-based solutions further, because of the possibility of exponential convergence of observer-based estimation even though equations (1) are not in the adaptive observer canonic form. Attempts have previously been made to develop observer-based methods for inferring parameters of
A canonic adaptive observer form is the following class of equations: x˙ = Ax + ϕ(y)θ + g(y)u, y = Cx, where x ∈ Rn is the state variable, y is the measured output, θ is the vector of unknown parameters, and u is the input (usually defined as a piece-wise continuous function). The function ϕ is known and continuous, matrices A, C are known and satisfy the following observability rank condition: rank [C T , AT C T , . . . , (AT )n−1 C T ] = n.
June 5, 2010 9:54 00235
State and Parameter Estimation for Canonic Models of Neural Oscillators
(1) from measured data such as membrane potential, v, and input current, I. For example, in Ref. 29, heuristic algorithms were proposed to estimate parameters of Hindmarsh-Rose model neurons. Yet, as the authors concluded, there is no systematic way to date for designing observers for parameter estimation of even the simplest model neurons. Developing techniques for state and parameter estimation with guaranteed convergence properties is the main goal of this study. 2.2. Contribution and structure of the paper We consider the following class of equations: x˙ 0 = θ0T φ0 (x0 , t) +
n
xi + ξ0 (t)
i=1
x˙ i = −τi xi + θiT φi (x0 , pi , t) + ξi (t), y = x0 ,
(2)
195
universal estimator that is capable of reconstructing unknown parameters of system (2) from the values of x0 (t) when t → ∞. Our method requires neither differentiation of measured evoked potentials, x0 (t), nor availability of internal currents for direct observation is required. It ensures not only asymptotic reconstruction of unknown model parameters, but also allows estimating the accuracy of fit as a function of the “richness” of the available data. Here richness is defined as a simple model-dependent condition on the non-singularity of a matrix that can easily be constructed from the values of x0 (t) and knowledge of the functions φi (·) in the right-hand side of (2). The paper is organized as follows. In Section 3 we present relevant concepts used to tackle the problem of parameter estimation, Section 4 contains main results, in Section 5 we provide an illustrative example, and Section 6 concludes the paper.
xi (t0 ) = xi,0
where φi : R × R≥0 → Rdi , di ∈ N/{0}, i = {0, . . . , n} are continuous functions modelling dynamic conductances of ion channels. Variable x0 in system (2) models the dynamics of the cell’s membrane potential, and variables xi , i ≥ 1 correspond to the ion currents in the cell. Parameters θi ∈ Rdi , pi ∈ R characterize membrane conductances. Parameters τi ∈ R>0 are the time constants of the currents. Variables ξi (t) constitute the unmodeled dynamics, or noise. The values of variable x0 (t) are assumed to be available for direct observation, and functions φ0 (x0 , t), φi (x0 , pi , t) are known. The variables xi (t), i = {1, . . . , n} are assumed to be not available for direct observation. Furthermore, true values of parameters θ0 , . . . , θn , τ1 , . . . , τn , p1 , . . . , pn are allowed to be unknown. We assume, however, that domains of admissible values of θi , pi , τi are known. In particular, we consider the case when θi,j ∈ [θi,min , θi,max ], τi ∈ [τi,min , τi,max ], pi ∈ [pi,min , pi,max ] and the values of θi,min , θi,max , τi,min , τi,max , pi,min , pi,max are available. Let us denote θ = (θ0 , θ1 , . . . , θn ),
3.
General Concepts of Parameter Inference for Ordinary Differential Equations
We will first specify the evaluation criteria needed for fitting the model parameters to data. Model neurons are systems of differential equations. Solutions of these systems are functions of space and time, and hence can be characterized by (1) their global asymptotic properties (e.g. by sets to which these solutions converge with time), and (2) by their local behavior in space and time. Generally, neither global nor local criteria alone give satisfactory results, especially if measurements are corrupted by noise (see Fig. 1). Plane (x0 , x1 ) in Fig. 1 represents the system’s state space. Black curves correspond to trajectories generated by the model, and gray curves depict
x0
x1
λ = (τ1 , p1 , . . . , τn , pn ),
and domains of admissible values for θ, λ are denoted by symbols Ωθ and Ωλ respectively. In this article we provide a method for systematic design of algorithms for asymptotic parameter reconstruction of model neurons that belong to a broad subclass of models (1). In particular, we present a
o
t1 t2 Fig. 1. time.
t3
t
Diagram of model and data trajectories in space-
June 5, 2010 9:54 00235
196
I. Tyukin et al.
a sample of actual data. Plane (x0 , t) is the space of observables, and x1 is a variable of which the values are not available for direct observation. In both planes there are points (marked by filled circles) at which distinct orbits in the space-time have identical or close projections on the state space, and viceversa. Hence neither the (x0 , x1 ) nor the (x0 , t) plane alone characterize data or model satisfactorily. 3.1. Observers Design of observers is the principal concept for asymptotic reconstruction of state and parameters of dynamical systems from knowledge of the history of trajectories and signals. In our case these trajectories and signals are the values of membrane potential, v(t), and input current I(t). In adaptive observer design, the original reconstruction problem is stated as follows. Let us suppose that, in compact notation, model (2) of evoked responses of a cell to stimulation I(t) is given by x˙ = g(x, I(t), θ, λ), θ˙ = 0, λ˙ = 0, v(t) = h(x(t)), and v(t), θ, λ are the measured potential and parameter vectors of which the values are unknown. Then the problem of inferring inaccessible variables of the model is viewed as a matter of finding a computational procedure in the form of a system: z˙ = f (z, t, I(t), v(t)),
ζ0 = z(t0 ) ∈ Ωz ,
(3)
such that inaccessible state variables and parameters are tracked by known functions of the internal state z of the observer. In the case of a state and parameter reconstruction problem, we will seek for (3) and funcˆ ˆ tions x ˆ(z), λ(z), θ(z) such that for sufficiently small given constant εx , ελ , εθ > 0, the following properties hold: lim sup ˆ x(z(t, ζ0 )) − x(t) ≤ εx , t→∞
ˆ lim sup λ(z(t, ζ0 )) − λ ≤ ελ , t→∞
ˆ ζ0 )) − θ ≤ εθ lim sup θ(z(t,
(4)
t→∞
∀ ξ0 ∈ Ωz . along the trajectories of a coupled model (data) observer system: x˙ = g(x, I(t), θ, λ) (5) z˙ = f (z, t, I(t), v(t)).
Traditionally xˆ(z) is viewed as an estimate of the ˆ ˆ are the estimates model state vector, x, and λ(z), θ(z) of parameter vectors λ, θ. The concept of observer-based estimation allows a general perspective on the problem of finding algorithms for fitting parameters of (1) to data. This perspective enables us to use the tools of analysis familiar from dynamical systems theory. In Section 4 we outline main ideas leading to the construction of observer-based parameter estimators for models (2).
3.2. Weakly attracting sets Current methods to determine algorithms (3), (4) (i.e. observers) are based on the notion of Lyapunov stability of solutions. If a globally Lyapunov-stable solution exists in the state space of (5) such that requirements (4) are satisfied, then (3) is the algorithm needed. Hence finding such algorithms may be achieved by employing the convenience offered by the method of Lyapunov functions. However, according to the results of simulation studies in Ref. 26, parameter values of general conductance models (1) corresponding to nearby trajectories, are not unique and are distributed across the whole parameter space. Therefore, if all we had available were the concepts of Lyapunov stability and the method of Lyapunov functions, we would have to deal with multiple separated attracting invariant sets. In such systems unstable sets inevitably occur at the boundaries of the attractor basins. Hence only local stability is generally guaranteed. It is known, however, that the method of Lyapunov functions is conservative when it comes to the local analysis in the problem of observer design for (1): its success and accuracy depend on how good the choice of the Lyapunov candidate is. To overcome locality and non-uniqueness issues, it is advantageous to replace the concept of Lyapunov stability with the more relaxed concepts of convergence to weakly attracting sets 2, 22 and relaxation times.8 A diagram illustrating the benefits of this less restrictive notion is provided in Fig. 2. In Fig. 2 Panel a depicts a diagram of standard stable attractors, panel b illustrates weak attractors, and panel c shows observer concept for (1). Domains of stable attractors are neighborhoods containing A1 , A2 . Estimates of sizes of these domains depend on parameters S0,1 , S0,2 , S0,3 . These estimates are
June 5, 2010 9:54 00235
State and Parameter Estimation for Canonic Models of Neural Oscillators
S 0,1 S 0,2 S 0,3 A1
A2
(a)
S 0,1 S 0,2 S 0,3 A1 A2
eo ac sp eters b Su ram pa
o pl Ex
ra
tio
n
Con
tr a c
Subspa ce of linea r pa r a meter s
ar
ne
nli
o fn
Subspace of linear p a r a meter s
(b)
tion
A
Con
tr a c
tion
t1
A
t2
(c)
Fig. 2.
Attractor and observer concepts.
depicted as closed curves around A1 , A2 . Domains of attraction for weak (Milnor) attracting sets are not neighborhoods. Hence, even in case of an unfortunate initial guess trajectories will simply converge to another (alternative) attractor via unstable paths. As part of our theoretical framework, we exploit the notion of weakly attracting sets in the problem of observer design and advocate the use of associated analysis results31 to solve this problem, as a useful complementary tool to the method of Lyapunov functions. 3.3. A-posteriori bifurcation analysis As follows from Refs. 1, 31, 28 and 7, parameter reconstruction of model neurons from data is b
197
achievable in principle. When the model is capable of reproducing data without errors, our algorithm ensures that model parameters can be reconstructed such that, over a given finite time interval, mismatches between trajectories of the model and data could be made negligibly small.b This case is equivalent to the assumption that the unmodeled dynamics, ξi (t), in (2) are negligible or non-existent. In practice, such a case will hardly ever be encountered; in general, non-zero mismatches between assumed model and data are to be expected. Model solutions are generally sensitive to parameter values; small perturbations of these values will result in models exhibiting different regimes with distinct invariant sets and asymptotic properties. Because of this, the presence of unmodeled dynamics may result in estimates that fail to generate plausible trajectories. Consider Class-I excitable model neurons, a class to which all the following model neurons belong: Morris-Lecar, Hindmarsh-Rose, FizhughNagumo, and modified Hodgkin-Huxley models.13 In these models the spiking mechanism is described by the homoclinic or saddle-on-invariant circle bifurcation. Figure 3 illustrates this mechanism for the homoclinic bifurcation. A model close to the bifurcation point is generating low-frequency oscillations (Fig. 3(c)). As shown in Fig. 3(b), a small perturbation of the excitation parameter may lead to abrupt changes of the solutions: instead of being periodic, almost all solutions will asymptotically converge to the stable equilibrium. The situation becomes even more complicated in practice for, according to most recent studies,27 ,25 popular models such as Morris-Lecar or HindmarshRose inherit more than one type of excitability. This means that not only the topology of the invariant sets changes with parameters (e.g. resting vs spiking) but also the spiking mechanisms become different. Hence in general, in order to account for model feasibility regions a-posteriori tailoring of the results of estimation is needed. The most straightforward solution to this problem is to map the results of estimation onto those domains in which the estimates are feasible. For planar models (or those reducible to this class) the regions of feasibility can be located by examining
This is a consequence of the continuity property of solutions of ordinary differential equations with respect to parameters and initial conditions.
June 5, 2010 9:54 00235
198
I. Tyukin et al.
3
3
3
2
2
1 1 (a)
(b)
(c)
Fig. 3. Typical transition of invariant sets in class-I model neurons. When the excitation parameter value is sufficiently small (Panel a), the system has three distinct equilibria: a stable one labeled 1 and two unstable ones labeled 2 and 3, respectively. Increasing the excitation parameter value (Panel b) leads to the collision of stable and unstable Equilibria 1 and 2, such that their stable and unstable manifolds coincide. Further increase (Panel c) gives rise to stable periodic orbits with variable oscillation periods.
local stability of model equilibria.10, 13 Projection of estimated parameters onto the feasibility regions corresponding to typical bifurcations (such as e.g. b in Fig. 3) can be achieved by searching for the feasible values of model parameters in the vicinity of the ˆ λ. ˆ obtained estimates θ, In order to do so, it is possible to identify the values of model parameters corresponding to the specific bifurcation that leads to unfolding of the desired dynamic regime of the system; in our case, the homoclinic loop (fold) bifurcation. Here we show that, at least for the Hindmarsh-Rose model neurons, MonteCarlo search procedures can be used successfully to locate the values of model parameters corresponding to the bifurcation. 4.
Observer-based Estimation of Unknown Model Parameters
Let us consider the following class of dynamical systems x˙ = f (x, θ) + g(x, θ)u(t), y = h(x),
x(t0 ) ∈ Ωx ⊂ Rn
x ∈ Rn , θ ∈ Rd , y ∈ R
(6)
where f, g : Rn → Rn , h : Rn → R, u : R → R are continuous and differentiable functions. Variable x is the state vector, u is the known input (representing e.g. the stimulation current I(t) in (1)), θ is the vector of unknown parameters, and y is the output of (6). Recall that y(t, 0, x0 , u) = h(x(t, 0, x0 , u)) denotes the output of (6) for u(·) and initial state x(0) = x0 . System (6) includes equations (1) as a subclass and in this respect can be considered
as a plausible generalization. Obviously, conclusions about (6) should be valid for systems (1) as well. Definition 1 (Observability24 ). Two states x1 , x2 ∈ Rn are said to be indistinguishable (denoted by x1 Ix2 ) for (6) if for every admissible input function u the output function t → y(t, 0, x1 , u), t ≥ 0 of the system for initial state x(0) = x1 , and the output function t → y(t, 0, x2 , u), t ≥ 0 of the system for initial state x(0) = x2 , are identical on their common domain of definition. The system is called observable if x1 Ix2 implies x1 = x2 . According to Definition 1, observability of a dynamical system implies that the values of its state, x(t) ∈ Rn , t ∈ [t1 , t2 ] are completely determined by inputs and outputs u(t), y(t) over [t1 , t2 ], although this definition of observability does not imply that every input function distinguishes points of Rn . Remark 1. Although Definition 1 does not account for the observability of any unknown parameter vectors, one can easily see that the very same definition can be used for parameterized systems as well. Indeed, extending the original equations (6) by including the parameter vector θ as a component of the extended state vector x ˜ = (xT , θT )T results in x ˜˙ = f˜(˜ x) + g˜(˜ x)u(t), ˜ y = h(˜ x), x˜(t0 ) = (x(t0 ), θ) ∈ Ωx˜ ⊂ R T
(7) n+d
where f˜(˜ x) = (f (x)T , 0)T , g˜(˜ x) = (g(x)T , 0)T , and T T ˜ x) = (h(x) , 0) . All uncertainties in (6), includh(˜ ing the parameter vector θ, are now combined into
June 5, 2010 9:54 00235
State and Parameter Estimation for Canonic Models of Neural Oscillators
the state vector of (7). Hence the problem of state and parameter reconstruction of (6) can be viewed as that of recovering the values of state for (7). Throughout this paper we shall consider input u as a constant (this corresponds I(t) = const in (1)), and combining the unknown parameter vector with the state vector as described in remark 1 enables us to write the system in the following form x˜˙ = f˜(˜ x), ˜ x), y = h(˜
x ˜ ∈ Rn˜ , y ∈ R
(8)
24
Proposotion 1 (Local observability test ). The system y = h(x), x ∈ Rn , y ∈ R
(9)
with smooth functions f, h is observable locally to the point x∗ if the following rank condition holds at that point rank
∂ (h(x∗ ) Lf h(x∗ ) · · · Ln−1 h(x∗ ))T = n f ∂x (10)
In7 it is shown that typical models neurons (Hindmarsh-Rose, Morris-Lecar) are locally observable systems.
The system with output measurements to be observed, state estimator (11), and parameter adaptation mechanism (12) are given below x(t) ˙ = A1 x(t) + φ0 (y(t), u(t)) + b
βi (y(t), u(t))θi
i=1
y(t) = C1 x(t) 0 1 0 ··· 0 0 1 ··· A1 = · · · ··· 0 0 0 ··· 0 0 0 ···
C1 = 1 0 0 · · · 0 x(t) ∈ R , n
(12)
K=
1 (A1 b + λb) = (k1 , . . . , kn )T b1
˙ ˆ) sign (b1 ) θˆ = Γβ(t)(y − C1 x
(13)
with Γ an arbitrary symmetric positive definite matrix and λ an arbitrary positive real. The n × 1 vector, b, is Hurwitz with b1 = 0 (see also21 where an observer is presented for systems in which the vector b is an arbitrary vector from Rn ). Theorem 1.18 There exists a local change of coordinates, z = Φ(x), transforming x˙ = f (x) +
p
θi (t)qi (x),
i=1 n
y = x1
x ∈ R , y ∈ R, θi ∈ R, qi : Rn → Rn , n ≥ 2
z˙ = A1 z + ψ(y) +
p
θi (t)ψi (y),
i=1
y = C1 z
z ∈ Rn , ψi : R → Rn
with (A1 , C1 ) in canonical observer form (11), if and only if (i) [adif g, adjf g] = 0, 0 ≤ i, j ≤ n − 1 (ii) [qi , adjf g] = 0, 0 ≤ j ≤ n − 2, 1 ≤ i ≤ p
0 0 · , 1 0
y(t) ∈ R,
x(t) + φ0 (y(t), u(t)) x ˆ˙ (t) = (A1 − KC1 )ˆ p βi (y(t), u(t))θˆi + Ky(t) +b
with h(xo ) = 0 and (f, h) an observable pair, into the system
4.1. Canonical adaptive observers
p
known, bounded and piecewise continuous functions of y(t), u(t). The column vector b ∈ Rn is assumed to be Hurwitz with b1 = 0 (e.g. all roots of the polynomial b1 pn−1 + · · · + bn−1 p + bn should have negative real part). An adaptive observer as presented in Ref. 19 is given as
i=1
Systems (8), in turn, admit the following simple observability test:
x˙ = f (x),
199
where the vector field, g(x), is uniquely defined by ∂ n−1 T (h(x), Lf h(x), . . . , Lf h(x)) , g(x) ∂x T
= (0, 0, . . . , 1) βi (·, ·) : R × R → R (11)
In (11) x(t) ∈ Rn is the state vector, and y = x1 (t) is the output; θ ∈ Rp = (θ1 , . . . , θp )T is the vector of unknown parameters. The functions βi are
(14)
Consider, for example, the problem of fitting parameters of the conventional HindmarshRose oscillator to measured data. It is shown in Refs. 28 and 7 that the Hindmarsh-Rose model, after translation and stretching of the state space and time
June 5, 2010 9:54 00235
I. Tyukin et al.
200
axis, can be written in vector-matrix notation as follows:
x1 x˙ 1 (15) = A(λ) + ϕ(x1 )θ x˙ 2 x2
cannot be transformed by diffeomorphic change of coordinates, z = φ(x), into z˙ = A1 z + ψ0 (y) +
A(λ)
=
ϕ(x1 ) θ
=
0 1 0 −λ x31
x21
0
0
0
1
0
0
0
0
x31
x21
x1
(16)
= (θ1,3 , θ1,2 , θ1,1 , θ1,0 , θ2,3 , θ2,2 , θ2,1 )
One of the main issues is that the original equations of neural dynamics are not written in a canonical form for which the reconstruction algorithms are available. The question, therefore, is if there exists an invertible coordinate transformation such that the model equations can be rendered canonic. Below we demonstrate that this is generally not the case if the transformation is parameter-independent. However, if we allow our transformation to be both parameter and time-dependent, a relevant class of models with polynomial right-hand sides can be transformed into one of the canonic forms. Let us consider a class of systems that can be described by (15). Clearly this system is not in a canonical adaptive observer form because A(λ) depends on the unknown parameter explicitly (see (16)). The question, however, is if there exists a differentiable coordinate transformation z = Φ(x), z1 = x1 such that in the new coordinates the equations of system (15) satisfy the canonic descriptions. We show that the answer to this question is negative, and it follows from the following slightly more general statement Theorem 2.7 The system x˙ = f (x) +
p
z ∈ Rn ,
y ∈ R,
ψi : R → Rn
(17)
with A1 , C1 in canonical observer form (11), if either (i) n > 2 or (ii) there exists i ∈ {1, . . . , p}, j ∈ {2, . . . , n} such that ∂qi /∂xj = 0.
, x1
θi ψi (y),
i=1
y = C1 z
where
p
Let us now consider the case in which the transformation z = Φ(x, λ, θ, t) is allowed to depend on unknown parameters and time. As we show below, this class of transformations is much more flexible. In principle it allows us to solve the problem of partial state and parameter reconstruction for an important class of oscillators with polynomial right-hand side and time-invariant time constants. We start by searching for a transformation Φ : q = T (λ)x, |T (λ)| = 0 such that
1 (18) T (λ)A(λ)T −1 (λ) =
0 where the matrix A(λ) is defined as in (16). It is easy to see that the transformation satisfying this constraint exists, and it is determined by
1 0 T (λ) = . (19) λ 1 According to (18), (19) and (16) equations of (13) in the coordinates q can be written as q˙ = A1 q + ψ(q1 )η(θ, λ), where
A1 =
ψ(q1 ) =
0 1 0 0 q13 0
(20)
q12 0
,
(21) q1 0
1 0 0 q13
0 q12
0 q1
0 1
,
η = (θ1,3 , θ1,2 , θ1,1 − λ, θ1,0 , λθ1,3 + θ2,3 , λθ1,2 + θ2,2 , λθ1,1 + θ2,1 , λθ1,0 )T
θi qi (x),
Remark 2. Notice that
i=1
y ∈ R,
y = h(x) x ∈ R , n
qi : R → R , n
with
0 1 1 ··· 1 0 0 0 ··· 0 f (x) = · · · · · · · x, 0 0 0 ··· 0
h(x) = 1 0 0 · · · 0 x
n
n≥2
• availability of the parameter vector η in (20), (21), expressed as a function of θ, λ, implies the availability of θ, λ if θ1,0 = 0. Indeed, in this case the value of λ = η8 /η4 and the values of all θi,j are uniquely defined by ηi ; • condition θ1,0 = 0 is sufficient for reconstructing the values of x2 provided that q and η are available; indeed in this case x = T −1 (λ)q.
June 5, 2010 9:54 00235
State and Parameter Estimation for Canonic Models of Neural Oscillators
As follows from Remark 2 the problem of state and parameter reconstruction of (13) from measured data x1 (t) amounts to solving the problem of state and parameter reconstruction of (20). In order to solve this problem we shall employ yet another coordinate transform: z1 = q1 , z2 = q2 + ζ T (t)η
(22)
in which the functions ζ(t) are some differentiable functions of time. Coordinate transformation (22) is clearly time-dependent. The role of this additional transformation is to transform the equations of system (20) into the form for which a solution already exists. Define the vector-function ζ(t) in (22) as follows: ζ˙i = −kζi + kψ1,i (q1 ) − ψ2,i (q1 ),
k ∈ R>0
In this case we have z˙2 =
8
ψ2,i (q1 )ηi + k −
i=1
−
8
ψ2,i (q1 )ηi = k
i=1
8
ζi ηi + ψ1,i (q1 )ηi
i=1 8
(−ζi + ψ1,i (q1 ))ηi
i=1
Taking (22) into account and expressing q2 as q2 = z2 − ζ T (t)η(θ, λ) we obtain
z˙1 = z2 +
8
(−ζi + ψ1,i (q1 ))ηi
i=1
z˙2 = k
8
(23)
(−ζi + ψ1,i (q1 ))ηi
i=1
Notice that ψ1,4 (q1 ) = ψ2,8 (q1 ) = 1, hence −ζ4 + ψ1,4 (q1 ) and −ζ8 + ψ1,8 (q1 ) converge to constants in R exponentially fast as t → ∞. Moreover, the sum −ζ4 + ψ1,4 (q1 ) is converging to zero, and the sum −ζ8 + ψ1,8 (q1 ) is converging to −1/k as t → ∞. Taking these facts into account we can conclude that system (23) can be rewritten in the following (reduced) form z˙ = A1 z + bφT (z1 , t)υ(θ, λ) + bε(t),
0 1 1 A1 = , b= , 0 0 k
201
−ζ1 + ψ1,1 (z1 ) −ζ + ψ (z ) 1,2 1 2 −ζ + ψ (z ) 3 1,3 1 φ(z1 , t) = −ζ8 + ψ1,8 (z1 ) −ζ5 + ψ1,5 (z1 ) −ζ6 + ψ1,6 (z1 ) −ζ7 + ψ1,7 (z1 ) υ(θ, λ) = (θ1,3 , θ1,2 , θ1,1 , θ¯1,0 , λθ1,3 + θ2,3 ,
λθ1,2 + θ2,2 , λθ1,1 + θ2,1 )T ,
(24)
where ε(t) is an exponentially decaying term. System (24) is clearly in adaptive canonic observer form. Hence it admits the following adaptive observer z − z) + bφT (z1 , t)ˆ υ zˆ˙ = A1 zˆ + L(ˆ
−l1 0 (25) L = , l1 , l2 ∈ R > 0 −l2 0 z1 − z1 )φi (z1 , t), γ ∈ R>0 υˆ˙ i = −γ(ˆ of which the asymptotic properties are specified in.7 In particular, the state and parameters are shown to converge to their true value provided that ∃ L, δ ∈ R>0 : t+L φ(z1 (τ ), τ )φ(z1 (τ ), τ )T dτ ≥ δI,
(26)
t
The function φ(z1 (t), t) is determined by the measured data, and condition (26) can easily be tested. A detailed explanation on how to get the values of θ and λ from the values of υˆ can be found in Ref. 7. So far we have considered special cases of (1) in which the time constants of unmeasured variables were unknown and parametrization of the right-hand side was linear. In the next section we discuss an approach that for a large subclass of (1) allows us to produce an observer solving the problem of state and parameter reconstruction from membrane potential measurement. The structure of this observer does not depend significantly on specific equations describing dynamics of the observed system. For this reason, and similarly to Ref. 12, we refer to this class of observers as universal adaptive observers. Detailed definitions of this class of observers and formal statements regarding their asymptotic properties are presented in Ref. 30. 4.2. Universal observers The main idea of our observer-based estimation method is to design uniformly stable estimators for
June 5, 2010 9:54 00235
202
I. Tyukin et al.
the parameters θi in (2) if time-constants and nonlinear parameters are known, followed by an automated search for the remaining parameters.c Thus the observer comprises of two compartments: stable, contracting and exploring, weakly attracting (interplay between these contracting-exploring compartments is illustrated in Fig. 2, Panel c). Contracting compartment. Consider (2), if parameters τi and pi would be known then the model would be in the adaptive canonic observer form. For these systems, estimation procedures are available that allow exponentially fast convergence of estimates θˆi of θi to small neighborhoods of θi cf.4, 20, 21 Let us combine variables xi and vectors θi into a single vector as follows: z = (x0 , . . . , xn , θ0T , . . . , θnT ), and let zˆ be its estimate. Then, according to4 and under mild non-degeneracy assumptions, similar to the non-singularity requirement for a matrix composed of filtered measured data in Ref. 11, there will exist a computational scheme x ˆ˙ = f (ˆ x, x0 (t), I(t)) (27) zˆ = (χ0 (ˆ x), . . . , χn (ˆ x), θˆ0 (ˆ x), . . . , θˆn (ˆ x)), where f (ˆ x, x0 (t), I(t)), χi (ˆ x), θˆi (ˆ x), are known functions, such that the dynamics of z − zˆ will obey the following estimated : z(t) − zˆ(t) ≤ De−ρ(t−t0 ) z(t0 ) − zˆ(t0 ) + ∆. (28) The presence of ∆ in (28) is due to measurement noise, ξi (t), and D ∈ R>0 is a constant of which the value can be computed a-priori. Observe that property (28) implies that the unperturbed system, i.e. when ξi (t) = 0, has a unique asymptotically stable solution in the sense of Lyapunov. Hence, we call this computational scheme the stable compartment of our observer. When the values of pi , τi are unknown, i.e. the vector λ is not measured, we replace them with their ˆ = (ˆ estimates, λ τ1 , pˆ1 , . . . , τˆn , pˆn ). Then it is possible to show that the dynamics of z(t) − zˆ(t) will satisfy the following inequality: z(t) − zˆ(t) ≤ De−ρ(t−t0 ) z(t0 ) − zˆ(t0 ) ˆ ) + C∆, + B sup λ − λ(τ τ ∈[t0 ,t]
c
(29)
where B, C are positive constants (see Refs. 30 and 32). Exploring compartment. In order to reconstruct the values of λ we extend the estimation procedure (27) by adding a component which searches for the values of λ by exploring all points in the domain Ωλ to which the values of unknown λ are supposed to belong. This explorative search is realized by a dynamical system of which the solutions form a (dense) grid in Ωλ (see Ref. 30 for details). Briefly, its working can be explained as follows. Let, for instance, Ωλ be a unit interval, and let λ(h) be a dense trajectory in Ωλ parameterized by h.e A possible candidate ˆ for the dense trajectory in [0, 1] is λ(h) = (sin(h) + 1)/2 (see Fig. 4). This trajectory can be modeled (up to a constant phase shift) by the following system: ˆ 1 ) = (1 + h1 )/2; λ(h h˙ 1 = h2 h˙ 2 = −h1 , h21 (t0 ) + h22 (t2 ) = 1. ˆ 1) If no further adjustments are made then λ(h will periodically visit every point in [0, 1]. In order to be able to stop at the point that corresponds to x0 (t)) = 0, those values of λ for which limt→∞ (x0 (t)−ˆ where xˆ0 (t) is the estimate of x0 (t) provided by the observer, we make the speed of exploration along the chosen dense trajectory proportional to the absoˆ0 (t)|. lute value of the trajectory error |x0 (t) − x Thus, we impose the following regulatory feedback mechanism for the explorative part of our observer: ˆ0 (t)|, γ ∈ R>0 . Then the exploh˙ ∼ γ|x0 (t) − x rative compartment of the estimator can be modeled by the following system of ordinary differential equations: ˆ 1 ) = (1 + h1 )/2; λ(h
(30)
ˆ0 (t)|h2 h˙ 1 = γ|x0 (t) − x
(31)
h˙ 2 = −γ|x0 (t) − xˆ0 (t)|h1 ,
(32)
h21 (t0 ) + h22 (t2 ) = 1.
General structure, formal definition and asymptotic properties of universal observers are provided in Refs. 30 and 7. Because of space limitations we do not reproduce these results here. d We refer the reader to Ref. 30 for details of derivations needed to obtain expressions in this section. e Specifically, we require that for all λ ∈ Ωλ and any δ ∈ R>0 there is a diverging sequence s1 , s2 ,. . . such that λ − λ(si ) ≤ δ.
June 5, 2010 9:54 00235
State and Parameter Estimation for Canonic Models of Neural Oscillators
λ
>
λ
203
e
1 ∧ λ=
h2
λ
h’2 0
λ
h
h’1
h1
ˆ as a function of h. Trajectory λ(h) ˆ Fig. 4. Explorative part of the observer. On the left - λ is dense in [0, 1]. In fact, ˆ i ) = λ. On the right - a diagram of for any λ ∈ [0, 1] there will be an infinite set of points si , i = 1, 2, . . . such that λ(s ˆ0 (t). Two invariant sets are present in the phase space of (31), (32). Variable e denotes instantaneous difference x0 (t) − x ˆ 1 ) = λ holds. system, both corresponding to the values of h1 for which the target equivalence λ(h
The values of γ are to be chosen such that suitable estimates form (weakly) attracting sets in the system state space (see Fig. 2). The evidence that such a choice can actually be made and the estimates of domains to which the values of γ should belong are provided in Ref. 31. In the next section we present universal observers designed specifically for dealing with state and parameter estimation of typical model neurons considered in this article. In comparison with other search-based methods our estimation technique has two distinct features that make its application advantageous: (1) Explorative search is conducted in a lowerdimensional subspace (one - dimensional for 2D Hindmarsh-Rose models (33)). This property allows us to combine efficiently the arsenal of conventional Lyapunov-based tools for designing stable compartments of the observer with the flexibility of the more general concept of weakly attracting sets. The result is – a substantial improvement in computational performance from trial-and-error approaches. (2) Exploration in Ωλ occurs simultaneously with exponentially fast estimation of θ. The explorative search does not need to wait until the exponentially fast compartment converges. Hence the amount of time needed to explore the domain of model parameters is reduced even further.
5.
Example: Fitting the Hindmarsh-Rose Model Neuron to Data
In this section we illustrate how the observer-based approach for state and parameter estimation of neural oscillators can be implemented for fitting models of neural oscillators to actual electrophysiological data. The data we used for this demonstration is shown in Fig. 5 (lower panel). Notice that the measured spikes are pretty sharp and narrow, which makes using the Hindmarsh-Rose model a preferred option over other low-dimensional models such as the Fitzhugh-Nagumo or Morris-Lecar equations.9 Consider a two dimensional Hindmarsh-Rose model written in the following form: x˙ 0 = θ0,3 x30 + θ0,2 x20 + θ0,1 x0 + x1 + θ0,0 + I(t) x˙ 1 = −λ1 x1 + θ1,2 x20 + θ1,1 x0 + θ1,0 .
(33)
To reconstruct the dynamics of a single spike we use those trials in which the input function I(t) is a square impulse. Within the time span of a singe burst this input function can be replaced with a constant. The value of θ1,0 can be aggregated into the parameter θ0,0 . Thus, instead of (33), we obtain the following equations: x˙ 0 = θ0,3 x30 + θ0,2 x20 + θ0,1 x0 + x1 + θ0,0 (34) x˙ 1 = −λ1 x1 + θ1,2 x20 + θ1,1 x0 .
June 5, 2010 9:54 00235
204
I. Tyukin et al.
