Lecture Notes in Computer Science

3

− y − z + I

˙

y = (a + α)x

2

− y

˙z = μ(bx + c

− z)

where x represents the membrane potential y and z are associated with fast

and slow currents, respectively. I is an applied current, and a, α, μ, b and c are

constant parameters. We rewrite the single HR neuron as a vectorized form:

(S

0

) : ˙

w = h(w) + Ξ(x, z)θ

where

w =

⎡

⎣

x

y

z

⎤

⎦

T

, h(w) =

⎡

⎣

−(x

3

+ y + z)

−y

0

⎤

⎦, Ξ(x, z) =

⎡

⎣

x

2

1 0 0 0 0

0 0 x

2

0 0 0

0 0 0 x 1

−z

⎤

⎦,

θ = θ

1

, θ

2

, θ

3

, θ

4

, θ

5

, θ

6

T

= a, I, a + α, μb, μc, μ

T

.

3.2

Numerical Examples

The HR model shows a large variety of behaviors with respect to the parameter

values in the differential equations[12]. Thus, we can characterize the dynamic

behaviors with respect to different values of the parameters. We focus on the

parameters a and I. The parameter a is an internal parameter in the single

neuron and I is an external depolarizing current. For the fixed I = 0.05, the

HR model shows a tonic bursting with a

∈ [1.8, 2.85] and a tonic spiking with

a ≥ 2.9. On the other hand, for the fixed a = 2.8, the HR model shows a tonic

bursting with I

∈ [0, 0.18] and a tonic spiking with a ∈ [0.2, 5].

86

K. Mitsunaga, Y. Totoki, and T. Matsuo

0

200

400

600

800

1000

−1.5

−1

−0.5

0

0.5

1

1.5

2

time

x

Fig. 1. The response of x in the tonic

bursting

−2

−1

0

1

2

0

2

4

6

8

−1

−0.9

−0.8

−0.7

−0.6

−0.5

−0.4

−0.3

x

y

z

Fig. 2. 3-D surface of x, y, z in the

tonic bursting

0

200

400

600

800

1000

−1.5

−1

−0.5

0

0.5

1

1.5

2

time

x

Fig. 3. The response of x in the tonic

spiking

−2

−1

0

1

2

0

2

4

6

8

−1

−0.9

−0.8

−0.7

−0.6

−0.5

x

y

z

Fig. 4. 3-D surface of x, y, z in the

tonic spiking

0

200

400

600

800

1000

−1

−0.5

0

0.5

1

1.5

2

time

x

Fig. 5. The response of x

1

in the in-

trinsic bursting neuron

−1

0

1

2

0

2

4

6

8

−1

−0.8

−0.6

−0.4

−0.2

x

y

z

Fig. 6. 3-D surface of x, y, z in the in-

trinsic bursting neuron

The parameters of the HR model in the tonic bursting (TB) are given by

a = 2.8, α = 1.6, c = 5, b = 9, μ = 0.001, I = 0.05.

Figure 1 shows the response of x. Figure 2 shows the 3 dimensional surface of

x, y, z. We call this neuron the intrinsic bursting neuron(IBN).

Firing Pattern Estimation of Biological Neuron Models

87

The parameters of the HR model in the tonic spiking (TS) are given by

a = 3.0, α = 1.6, c = 5, b = 9, μ = 0.001, I = 0.05.

The difference between the tonic bursting and the tonic spiking is only the

value of the parameter a. Figure 3 shows the responses of x. Figure 4 shows the

3 dimensional surface of x, y, z. We also call this neuron the intrinsic spiking

neuron(ISN).

When the external current change from I = 0.05 to I = 0.2, the IBN shows a

tonic spiking. Figure 5 shows the responses of x. Figure 6 shows the 3 dimensional

surface of x, y, z.

4

Synaptically Coupled Model of HR Neuron

4.1

Dynamical Equations

Consider the following synaptically coupled two HR neurons[1]:

˙

x

1

= a

1

x

2

1

− x

3

1

− y

1

− z

1

− g

s

(x

1

− V

s1

)Γ (x

2

)

˙

y

1

= (a

1

+ α

1

)x

2

1

− y

1

,

˙

z

1

= μ

1

(b

1

x

1

+ c

1

− z

1

)

˙

x

2

= a

2

x

2

2

− x

3

2

− y

2

− z

2

− g

s

(x

2

− V

s2

)Γ (x

1

)

˙

y

2

= (a

2

+ α

2

)x

2

2

− y

2

,

˙

z

2

= μ

2

(b

2

x

2

+ c

2

− z

2

)

where Γ (x) is the sigmoid function given by

Γ (x) =

1

1 + exp(

−λ(x − θ

s

))

.

4.2

Numerical Examples

Consider the IBN neuron with a = 2.8 and the ISN neuron with a = 10.8 whose

other parameters are as follows:

α

i

= 1.6, c

i

= 5, b

i

= 9, μ

i

= 0.001, V

si

= 2, θ

s

=

−0.25, λ = 10.

Figures 7,8 show the responses of the membrane potentials in the coupling of

IBN neuron and ISN neuron with the coupling strength g

s

= 0.05, respectively.

Each neuron behaves as an intrinsic single neuron. As increasing the coupling

strength, however, the IBN neuron shows a chaotic behavior. Figures 9,10 show

the responses of the membrane potentials in the coupling of IBN neuron and

ISN neuron with the coupling strength g

s

= 1, respectively. Figure 11 shows

the response of the membrane potentials in the coupling of two same IBNs with

88

K. Mitsunaga, Y. Totoki, and T. Matsuo

0

200

400

600

800

1000

−1.5

−1

−0.5

0

0.5

1

1.5

2

time

x

Fig. 7. The response of x

1

of the IBN

0

200

400

600

800

1000

−4

−2

0

2

4

6

8

10

12

time

x

Fig. 8. The response of x

2

of the ISN

0

200

400

600

800

1000

−1.5

−1

−0.5

0

0.5

1

1.5

2

time

x

Fig. 9. The response of x

1

of the IBN

0

200

400

600

800

1000

−4

−2

0

2

4

6

8

10

12

time

x

Fig. 10. The response of x

2

of the ISN

0

200

400

600

800

1000

−1.5

−1

−0.5

0

0.5

1

1.5

2

time

x

Fig. 11. The response of x

1

of the IBN-

IBN coupling with g

s

= 0.05

0

200

400

600

800

1000

−1

−0.5

0

0.5

1

1.5

2

time

x

Fig. 12. The response of x

1

of the IBN-

IBN coupling with g

s

= 1

the coupling strength g

s

= 0.05. In this case, two IBNs synchronize as bursting

neurons. Figure 12 shows the response of the membrane potentials in the coupling

of two same IBNs with the coupling strength g

s

= 1. Two IBNs synchronize as

spiking neurons.

Firing Pattern Estimation of Biological Neuron Models

89

5

Adaptive Observer with Full States

We present the parameter estimation problem to distinguish the firing patterns

by using early-time dynamic behaviors. In this section, assuming that the full

states are measurable, we present an adaptive observer to estimate all parameters

in the single HR neuron.

5.1

Construction of Adaptive Observer

We present an adaptive observer as

(O

0

) : ˙ˆ

w = W ( ˆ

w − w) + h(w) + Ξˆθ

where ˆ

w = [ˆx,

ˆ

y,

ˆ

z] is an estimate of the states, ˆθ is an estimate of the

unknown parameters and W is selected as a stable matrix. Using the standard

adaptive control theory[11], the parameter update law is given by

˙ˆθ = Γ Ξ

T

P (w − ˆ

w).

where P is a positive definite solution of the following Lyapunov equation for a

positive definite matrix Q:

W

T

P + P W = −Q.

5.2

Numerical Examples

We will show the simulation results of single IBN case. The parameters in the

tonic bursting are given by

a = 2.8, α = 1.6, c = 5, b = 9, μ = 0.001, I = 0.05.

The parameters of the adaptive observers are selected as

W = −10I

3

, Γ = diag{100, 50, 300}.

Figure 13 shows the estimation behavior of a (solid line) and I (dotted line).

The estimates ˆ

a and ˆ

I converge to the true values of a and I.

6

Adaptive Observer with a Partial State

We assume that the membrane potential x is available, but the others are im-

measurable. In this case, we consider following problems:

– Estimate y and z using the available signal x;

– Estimate the parameter a or I to distinguish the firing patterns by using

early-time dynamic behaviors.

90

K. Mitsunaga, Y. Totoki, and T. Matsuo

6.1

Construction of Adaptive Observer

The parameters a and I are key parameters that determine the firing pattern.

The HR model can be rewritten by the following three forms[9]:

(S

1

) : ˙

w = Aw + h

1

(x) + b

1

(x

2

a)

(1)

(S

2

) : ˙

w = Aw + h

2

(x) + b

2

(I)

(2)

(S

3

) : ˙

w = Aw + h

3

(x) + b

2

(θ

T

ξ)

(3)

where

A =

⎡

⎣

0

−1 −1

0

−1 0

μb 0 −μ

⎤

⎦ , h

1

=

⎡

⎣

−x

3

+ I

αx

2

μc

⎤

⎦ , b

1

=

⎡

⎣

1

1

0

⎤

⎦ , b

2

=

⎡

⎣

1

0

0

⎤

⎦ ,

h

2

=

⎡

⎣

−x

3

+ ax

2

(a + α)x

2

μc

⎤

⎦ , h

3

=

⎡

⎣

−x

3

δx

2

μc

⎤

⎦ , θ = a

I

, ξ =

x

2

1

.

In (S

1

) and (S

2

), the unknown parameters are assumed to be a and I, respec-

tively. In (S

3

), we assume that the parameter δ = a + α is known, and a and I

are unknown.

Since the measurable signal is x, the output equation is given by

x = cw = 1 0 0 w

We present adaptive observers that estimate the parameter for each system

(S

i

), i = 1, 2, 3, as follows:

(O

1

) : ˙ˆ

w

1

= A ˆ

w

1

+ h

1

(x) + b

1

(x

2

ˆ

a) + g(x − ˆx)

(4)

(O

2

) : ˙ˆ

w

2

= A ˆ

w

2

+ h

2

(x) + b

2

( ˆ

I) + g(x − ˆx)

(5)

(O

3

) : ˙ˆ

w

2

= A ˆ

w

2

+ h

3

(x) + b

2

(ˆ

θ

T

ξ) + g(x − ˆx)

(6)

where g is selected such that A

− gc is a stable. Since (A, b

1

, c) and (A, b

2

, c)

are strictly positive real, the parameter estimation laws are given as

˙ˆa = γ

1

x

2

(x

− ˆx),

˙ˆ

I = γ

2

(x

− ˆx).

(7)

Using the Kalman-Yakubovich (KY) lemma, we can show the asymptotic sta-

bility of the error system based on the standard adaptive control theory[11].

6.2

Numerical Examples

We will show the simulation results of single IBN case. The parameters in the

tonic spiking are same as in the previous simulation. Figures 14 and 15 show

the estimated parameters by the adaptive observers (O

1

) and (O

2

), respectively.

Figures 16 and 17 show the responses of y (solid line) and its estimate ˆ

y (dot-

ted line) the adaptive observers (O

1

) for t

≤ 500 and for t ≤ 20, respectively.

Figure 18 shows the responses of z (solid line) and its estimate ˆ

z (dotted line) by

the adaptive observers (O

1

). The simulation results of other cases are omitted.

The states and parameters can be asymptotically estimated.

Firing Pattern Estimation of Biological Neuron Models

91

0

50

100

150

200

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

time

ˆ

a

ˆ

I

ˆ

a

ˆ

I

Fig. 13. ˆa(solid line) and ˆ

I(dotted line) in

the adaptive observer (O

0

) with full states

0

20

40

60

80

100

1

1.5

2

2.5

3

3.5

4

4.5

5

time

ˆ

a

Fig. 14. ˆa(solid line) in the adaptive ob-

server (O

1

) with x

0

20

40

60

80

100

−0.5

−0.4

−0.3

−0.2

−0.1

0

0.1

time

ˆ

I

Fig. 15. ˆ

I(solid line) in the adaptive ob-

server (O

2

) with x

0

100

200

300

400

500

0

1

2

3

4

5

6

7

time

y

ˆ

y

y

ˆ

y

Fig. 16. y(solid line) and ˆ

y(dottedline) in

the adaptive observer (O

1

) (t

≤ 500)

0

5

10

15

20

0

1

2

3

4

5

6

7

time

y

ˆ

y

y

ˆ

y

Fig. 17. y(solid line) and ˆ

y(dottedline) in

the adaptive observer (O

1

) (t

≤ 20)

0

5

10

15

20

−1.5

−1

−0.5

0

time

z

ˆ

z

z

ˆ

z

Fig. 18. z(solid line) and ˆz(dottedline) in

the adaptive observer (O

1

) (t

≤ 20)

92

K. Mitsunaga, Y. Totoki, and T. Matsuo

7

Conclusion

We presented estimators of the parameters of the HR model using the adaptive

observer technique with the output measurement data such as the membrane

potential. The proposed observers allow us to distinguish the firing pattern in

early time and to recover the immeasurable internal states.

References

1. Belykh, I., de Lange, E., Hasler, M.: Synchronization of Bursting Neurons: What

Matters in the Network Topology. Phys. Rev. Lett. 94, 101–188 (2005)

2. Izhikevich, E.M.: Simple Model of Spiking Neurons. IEEE Trans on Neural Net-

works 14(6), 1569–1572 (2003)

3. Izhikevich, E.M.: Which model to use for cortical spiking neurons? IEEE Trans on

Neural Networks 15(5), 1063–1070 (2004)

4. Watts, L.: A Tour of NeuraLOG and Spike - Tools for Simulating Networks of

Spiking Neurons (1993), http://www.lloydwatts.com/SpikeBrochure.pdf

5. Hindmarsh, J.L., Rose, R.M.: A model of the nerve impulse using two first order

differential equations. Nature 296, 162–164 (1982)

6. Hindmarsh, J.L., Rose, R.M.: A model of neuronal bursting using three coupled

first order differential equations. Proc. R. Soc. Lond. B. 221, 87–102 (1984)

7. Carroll, T.L.: Chaotic systems that are robust to added noise, CHAOS, 15, 013901

(2005)

8. Meunier, N., Narion-Poll, R., Lansky, P., Rospars, J.O.: Estimation of the Individ-

ual Firing Frequencies of Two Neurons Recorded with a Single Electrode. Chem.

Senses 28, 671–679 (2003)

9. Yu, H., Liu, Y.: Chaotic synchronization based on stability criterion of linear sys-

tems. Physics Letters A 314, 292–298 (2003)

10. Marino, R.: Adaptive Observers for Single Output Nonlinear Systems. IEEE Trans.

on Automatic Control 35(9), 1054–1058 (1990)

11. Narendra, K.S., Annaswamy, A.M.: Stable Adaptive Systems. Prentice Hall Inc.,

Englewood Cliffs (1989)

12. Arena, P., Fortuna, L., Frasca, M., Rosa, M.L.: Locally active Hindmarsh-Rose

neurons. Chaos, Soliton and Fractals 27, 405–412 (2006)

13. Tokuda, I., Parlitz, U., Illing, L., Kennel, M., Abarbanel, H.: Parameter estimation

for neuron models. In: Proc. of the 7th Experimental Chaos Conference (2002),

http://www.physik3.gwdg.de/

∼

ulli/pdf/TPIKA02 pre.pdf

14. Steur, E.: Parameter Estimation in Hindmarsh-Rose Neurons (2006),

http://alexandria.tue.nl/repository/books/626834.pdf

15. Fujikawa, H., Mitsunaga, K., Suemitsu, H., Matsuo, T.: Parameter Estimation of

Biological Neuron Models with Bursting and Spiking. In: Proc. of SICE-ICASE

International Joint Conference 2006 CD-ROM, pp. 4487–4492 (2006)

Thouless-Anderson-Palmer Equation for

Associative Memory Neural Network Models

with Fluctuating Couplings

Akihisa Ichiki and Masatoshi Shiino

Department of Applied Physics, Faculty of Science, Tokyo Institute of Technology,

2-12-2 Ohokayama Meguro-ku Tokyo, Japan

Abstract. We derive Thouless-Anderson-Palmer (TAP) equations and

order parameter equations for stochastic analog neural network models

with fluctuating synaptic couplings. Such systems with finite number

of neurons originally have no energy concept. Thus they defy the use

of the replica method or the cavity method, which require the energy

concept. However for some realizations of synaptic noise, the systems

have the effective Hamiltonian and the cavity method becomes applicable

to derive the TAP equations.

1

Introduction

The replica method [1] for random spin systems has been successfully employed

in neural network models of associative memory to have the order parameters

and the storage capacity [2] and the cavity method [3] has been employed to de-

rive the Thouless-Anderson-Palmer (TAP) equations [4,5]. However these tech-

niques require the energy concept. On the other hand, various types of neural

network models which have no energy concept, such as networks with tempo-

rally fluctuating synaptic couplings, may exist. The alternative approach to the

replica method to derive the order parameter equations, called the self-consistent

signal-to-noise analysis (SCSNA), is closely related to the cavity concept in the

case where networks have free energy [6,7]. An advantage to apply the SCSNA

to neural networks is that the energy concept is not required to derive the order

parameter equations once the TAP equations are obtained. The SCSNA, which

was originally proposed for deriving a set of order parameter equations for deter-

ministic analog neural networks, becomes applicable to stochastic networks by

noting that the TAP equations define the deterministic networks. Furthermore,

the coefficients of the Onsager reaction terms characteristic to the TAP equations

which determine the form of the transfer functions in analog networks are self-

consistently obtained through the concept shared by the cavity method and the

SCSNA. Thus the TAP equations as well as the order parameter equations are

derived self-consistently by the hybrid use of the cavity method and the SCSNA

in the case where the energy concept exists. However the networks with synaptic

noise, which have no energy concept, defy the use of the cavity method to have

the TAP equations. On the other hand, as in [8], the network with a specific