Bifurcations in a Silicon Neuron

Viewer
Transcript

Bifurcations in a Silicon Neuron Arindam Basu∗, Csaba Petre† and Paul Hasler‡ School of Electrical and Computer Engineering Georgia Institute of Technology, Atlanta, Georgia 30332–0250 Email: ∗ [email protected];†[email protected]; ‡ [email protected] Abstract— In this paper, we describe the bifurcations occurring in a silicon neuron with one sodium and one potassium channel. The channels are designed to model the physics of ion flow in actual biological channels instead of modeling a particular set of equations. We show a pair of subcritical Hopf-bifurcation with increase in current stimulus which is characteristic of class 2 excitability in Hodgkin-Huxley neurons. Theoretical analysis of the bifurcations lead to conditions for designing and biasing the circuit. The circuit is very compact, comprising six transistors and three capacitors, lending itself to easy integration. The parameters are set using floating-gate transistors and can be programmed as desired. We hope to study more complicated dynamics of large networks of these neurons, a task which might be beyond a typical digital computer.

I. M ODELING B IOLOGICAL N EURONS

IN

M1

ENa MNa

CZ

I τh

Vfg M3

Vout M2

Iin

CNa

Vτh

Iamp Vmem

Cmem

Vτn

CK

S ILICON

Historically there has been a trend of studying neuronal models on an analog platform, the earliest instance of which dates back to Richard Fitzhugh in the sixties. Also with the progress in semiconductor technology, power supply levels have become similar to biology opening a new possibility of the same silicon model to interface with live neurons [1]. Almost all the implementations use the Hodgkin-Huxley (HH) formalism, some of them using a set of equations simpler than the original H-H equations [2] [3] while others have modelled in detail the full set of equations [4]. However, it is believed that computational properties of neurons are based on the bifurcations exhibited by these dynamical systems in response to some changing stimulus [5] leading one to believe that all models which present the same set of bifurcations should be equally good in analyzing and modeling neurons. In other words, it is not necessary to make a silicon model to reproduce the exact H-H equations; any model which generates a ﬂow that is topologically equivalent to the ﬂow generated by the equations governing the physical system is good enough. This is evidently attractive as it might be simpler to use the nonlinearities in transistors to make an equivalent nonlinear ﬂow than to implement the original H-H equations. We exploit the similarity between the biological and transistor channels by using a mosfet to represent a channel [6]. The gates of these “channel” transistors are controlled by ampliﬁers so that their voltage clamp responses resemble biology. A combination of one sodium and one potassium channel makes a basic neuron as shown in Fig. 1. Of course, this does not rigorously prove that the equations corresponding to this circuit produce topologically equivalent ﬂows as the H-H equations and in general, it is difﬁcult to prove such a proposition. We rather concentrate on the bifurcations ex-

978-1-4244-1684-4/08/$25.00 ©2008 IEEE

Vamp

MK

VK

M4

I τn

Vgk

EK Fig. 1: A Neuron with two channels: A simple neuron with one sodium channel (MNa ) and one potassium channel (MK ). The sodium channel transistor is controlled by an ampliﬁer with bandpass characteristics giving it the fast activation and slow inactivation dynamics. The potassium channel transistor is controlled by a lowpass ampliﬁer providing it with slow activation.

hibited by this circuit and use it as a metric to validate this model. We hope that this analysis not only helps in the design of the silicon neuron but also results in an alternate model for computational neuroscientists that can be run on a digital computer. In section II, we describe the operation of the neuron and derive the differential equations governing its dynamics. We also ﬁnd ﬁxed points and conditions for bifurcation. Section III deals with two particular Hopf-bifurcations observed in this model. In section IV, we present some results from a reconﬁgurable IC and describe other chips with similar structures that have been fabricated. Finally we present the conclusions in the last section. II. C IRCUIT O PERATION AND M ODEL We start by discussing the main parameters used in the model to study bifurcations and their relation to circuit properties. As shown in ﬁg. 1, the channels consist of a channel transistor and a control ampliﬁer. The high frequency corner in

428

Authorized licensed use limited to: Georgia Institute of Technology. Downloaded on July 08,2010 at 08:38:47 UTC from IEEE Xplore. Restrictions apply.

800 700

100

100

80

80

60

60

40

40

Increasing Iamp

500

300

Γ

Γ

400

Γ

600

20

200 100

20

0

0

0

Increasing Iτ h

−100 −200

0

10

20

30 iin

40

Increasing Iτ n

−20

50

60

−40

0

5

10

15 iin

(a)

20

−20

25

30

(b)

−40

0

10

20 iin

30

40

(c)

Fig. 2: Effect of varying parameters: Condition for a Hopf-bifurcation is a pair of imaginary eigen-values with zero real part. The zero crossings of Γ in eq. 6 provide necessary conditions for such a case when (a) Iτ h is varied (b) Iτ n is varied and (c) Iamp is varied. We want to bias the neuron such that we get at least two such intersections which could lead to two Hopf-bifurcations.

the sodium ampliﬁer is controlled by the bias current Iamp set by M2 and the low corner is set by the current Iτ h through M3. M3 acts as a nonlinear resistor in feedback across the ampliﬁer. Different nonlinear resistors like an NMOS or a diode connected MOS can also be used and are similar in the small-signal sense but give rise to different interesting large-signal dynamics. The gain of this ampliﬁer is set to approximately 8 by the ratio of the capacitors CN a and CK . Iamp is responsible for the activation of the sodium channel and Iτ h provides inactivation. The potassium channel is controlled by a lowpass ampliﬁer with its corner being set by Iτ n , the current through M4. Again different nonlinear resistors could be used here. The currents from MN a and MK are added on the membrane capacitor, Cmem on the node Vmem . The voltage Vmem varies between EN a and EK which have a difference of 200mV that is close to biology. The power supply for the sodium ampliﬁer, Vamp ﬁxes a bias current, IN a0 through MN a . Also, an external stimulus current, Iin is considered as a bifurcation parameter. We now develop the dynamical equations for the system by using KCL at the Vmem , Vk , Vf g and Vout nodes respectively. The variables x, n, h and m are deﬁned as follows: δVmem δVk x= ;n = UT UT δVf g δVout h= ;m = , (1) UT UT where the δV variables represent changes in the voltages from steady state and UT is thermal voltage. Notice that there are four variables as in the H-H equations which is natural as both of them model the same system. UT x˙ =

IN a0 −m Iin x−4 Iin e {1 − (1 + )e }+ − Cmem IN a0 Cmem

IK0 x−n e Cmem Iτ UT n˙ = UT x˙ − n f (n) CK

Iamp −h {e − em/50 } 8CZ 1.125Iamp −h Iτ UT m ˙ = UT x˙ + {e − em/50 } + h g(h, m) CZ CZ (2) UT h˙ = UT x˙ +

where f and g represent the currents through M4 and M3 respectively. This generalization allows for analyzing a general nonlinear resistor in that place. For this case the functions are Iin x−4 f = en −1 and g = eh −em . Also, note the term 1− IN e a0 in the ﬁrst equation which models the sodium transistor entering ohmic region which is the essential term for the bifurcations at large currents. Ohmic effects of MK , M 1, M 2 are ignored for this particular case but are needed in a more general setting. To reduce the number of parameters, we do another renormalization as follows: tIN a0 Cmem UT Ci Ci = Cmem Ij Ij = IN a0 n = n − x, h = h − x, m = m − x, τ=

(3)

and then dropping the primes for economy of notation we get: x˙ = e−(w+x) {1 − (1 + Iin )ex−4 } + Iin − IK0 e−y Iτ y˙ = − n {ey+x − 1} CK Iamp −(z+x) z˙ = {e − e(w+x)/50 } 8CZ 1.125Iamp −(z+x) Iτ w˙ = {e − e(w+x)/50 } + h {ez+x − ew+x } CZ CZ (4) where IK0 = (1 + Iin )(1 − e−4 ) includes the effect of input stimulus.

429

Authorized licensed use limited to: Georgia Institute of Technology. Downloaded on July 08,2010 at 08:38:47 UTC from IEEE Xplore. Restrictions apply.

8

8

2

6

6

1.95

4

4

2

2

0

Vmem (V)

Im(λ)

Im(λ)

1.9

0

1.85

−2

−2

1.8

−4

−4

1.75

−6

−6

−8 −50

−8 −1

1.7

0 Real(λ)

50

0 Real(λ)

1

Fig. 3: Spectrum of the stability matrix: The eigen-values of the stability matrix are plotted with increasing Iin for Iτh = Iτ n = 0.8, Iamp = 15. The ﬁgure on the right shows a closeup of the complex eigen-values demonstrating two Hopf-bifurcations.

A. Finding Fixed Points The ﬁxed points of the set of equations is found by setting the RHS to zero which gives the origin as an unique ﬁxed point. This is intuitively good as you expect the deviations from steady state to be zero. It should be noted that the increase of resting membrane potential with increasing stimulus current does not show up because of the way we deﬁne IK0 . An easy practical way of ﬁnding the equilibria is to sweep Vmem and ﬁnd the steady state current through the voltage source. The number of zero crossings of this I-V curve give the number of equilibria. B. Conditions for Hopf-Bifurcations As mentioned earlier, in this paper we concentrate on two Hopf-bifurcations exhibited by neurons. In fact the deﬁning property of Class 2 excitability is the existence of a stable equilibrium for small Iin but at a critical value of Iin , the neuron exhibits spontaneous ﬁring with large amplitude and non-zero frequency. This corresponds to a subcritical Hopfbifurcation. At a very large value of the input, the equilibrium regains stability by a second Hopf bifurcation which may be subcritical or supercritical. So we want to ﬁnd the region of parameter space where our circuit can produce two Hopfbifurcations with increasing Iin . The relevant conditions for Hopf-bifurcation in our case may be stated as follows: (1) Two of the eigen-values of the stability matrix of the equilibrium must be complex conjugates with zero real part and non-zero imaginary part at the bifurcation. (2) The derivative of the real part of the complex eigenvalues with respect to the parameter must be non-zero at the bifurcation [7]. In our case the two real eigen-values will be negative as we want the resting state to be stable before the bifurcation. To use the above conditions, we equate the characteristic polynomial of the stability matrix of the equilibrium to the desired form with two real and two complex eigen-values with

0

0.2

0.4

0.6 0.8 Iin (A)

1

1.2

1.4 −7

x 10

Fig. 4: Bifurcation in SPICE simulation: Simulation of the neuron with Iin being a slowly increasing ramp. We can clearly see the loss of stability of the equilibrium and spontaneous oscillations at around Iin = 2nA. At around Iin = 120nA the limit cycle disappears.

zero real part at equilibrium: x4 +a3 x3 +a2 x2 +a1 x+a0 = (x2 +ω02 )(x2 +(λ1 +λ2 )x+λ1 λ2 ), (5) where −λ1 and −λ2 are the real eigen-values. Then the conditions for Hopf-bifurcation become: a3 = λ1 + λ2 > 0 a1 = ω02 > 0 a3 a1 a0 a3 Γ = a2 − + =0 a3 a1

(6)

These conditions are necessary but not sufﬁcient as the condition for non-zero derivative of the real part of the complex eigen-values needs to be appended. Fig. 2 shows plots of Γ with varying Iin for different Iτ h , Iτ n and Iamp . The zerocrossings of these curves along with the ﬁrst two conditions in 6 give possible bifurcation points. Fig. 3 shows a plot of the eigen-values with changing stimulus current and Iτh = Iτ n = 0.8, Iamp = 15 clearly showing two Hopf-bifurcations. III. B IFURCATIONS

IN THE

N EURON

Having established the regions in parameter space where we need to bias the circuit, we can directly look at the response of the system when the input stimulus is slowly increased. Fig. 4 shows a SPICE simulation of the neuron with Iτh = Iτ n = 0.8, Iamp = 15. The input was a slow ramp of current. The circuit exhibits spontaneous large oscillations when Iin is large enough, a classic case of Class 2 excitability. At a much larger value of the stimulus the equilibrium becomes stable again. This happens as the resting potential becomes close to EN a which makes MN a ohmic. The correctness of the mathematical model was veriﬁed by numerically integrating the equations to generate a similar bifurcation diagram as shown in ﬁg. 5. The critical value of Iin is different for the second bifurcation as the value of EN a − EK was kept at 300mV in the circuit as compared to 200mV in the numerical analysis.

430

Authorized licensed use limited to: Georgia Institute of Technology. Downloaded on July 08,2010 at 08:38:47 UTC from IEEE Xplore. Restrictions apply.

4

1.5

0.5 Input Output

3

1

1

-0.5

-1

0 −1

-1.5

−2 −3

Vmem(V)

2 V(V)

Vmem(normalised)

0

0

5

10

15 20 Iin(normalised)

25

30

35

Fig. 5: Bifurcation in numerical integration: The bifurcation of the theoretical model is observed by integrating the differential equations numerically. The equilibrium was perturbed to check its stability and ﬁnd bifurcations. We can see two subcritical Hopf-bifurcations.

It should be noted that the two bifurcation diagrams are different. Fig. 5 is a true bifurcation plot of the equilibrium showing the hard loss of stability. For the ramp experiment though, the second bifurcation we see is actually a fold bifurcation of the limit-cycle. The Hopf-bifurcation of the equilibrium occurred earlier but was not visible as the solution was in the basin of attraction of the limit-cycle. The unstable limit cycle acts as the threshold or basin boundary between their basins of attraction. Thus the co-existence of stable attractors leads to a kind of hysteresis in real experiments. The existence of the unstable limit cycle causing the Hopfbifurcations can be proved by making a Poincar´e map of the system of equations, or by using a continuation method or by integrating the equations in reverse time. To rigorously prove that the bifurcation is indeed a subcritical Hopf, we project the system on a suspended center manifold at equilibrium to get the normal form coefﬁcients [7]: r˙ = r(0.113μ − 0.03μ2 + (0.61 + 0.061μ − 0.012μ2)r2 ) θ˙ = 2.42 + 0.5μ − 0.02μ2 + (−0.015 − 0.034μ + 0.014μ2), (7) where μ = Iin − Iin,b is the deviation of Iin from its critical value. The cubic coefﬁcient in the equation for amplitude is positive at bifurcation which indicates that the r = 0 solution is unstable at bifurcation and there exists an unstable limit cycle prior to bifurcation. This conclusively shows that the bifurcation was indeed a subcritical Hopf. The details of the normal-form calculation are omitted here for lack of space. A similar center-manifold expansion at the second bifurcation conﬁrms it to be subcritical too. However, there are parameter values which result in a supercritical Hopf, details of which are being studied currently. IV. H ARDWARE P LATFORM The neuron was compiled on a Field Programmable Analog Array (FPAA) [8] and initial experiments showing the bandpass response of a sodium control ampliﬁer and spiking

0.5

0 -0.5

0

0.2

0.4 0.6 0.8 t(sec) (a)

1

-1

0

0.1

0.2

0.3

(b) t(sec)

0.4

0.5

Fig. 6: (a) Measured step response of the sodium ampliﬁer showing two time constants. (b) Measured spiking waveforms generated by the neuron.

behavior is demonstrated in ﬁg. 6. A test chip having one sodium and one potassium channels is in fabrication and expected to be back soon which will enable us to conduct a detailed study of the bifurcations. V. C ONCLUSION In this paper, we modelled the dynamics of a silicon neuron and demonstrated the bifurcations exhibited by it. The principle of the design is based on similarities between transport of charge carriers in transistor channels and biological ion channels. The gates of these channel transistors are driven by ampliﬁers designed to give voltage clamp responses resembling biology. A neuron with one such sodium and potassium channel show dynamics very similar to biology and those exhibited by Hodgkin-Huxley equations. In particular, we show subcritical Hopf-bifurcations which is the trademark of Class 2 neural excitability. The compactness of the circuit allows a large density of channels on a chip and opens up the possibility exploring the dynamics of complicated networks. R EFERENCES [1] T. Bal S. Le Masson, A. Laﬂaquiere and G. Le Masson, “Analog circuits for modeling biological neural networks: Design and applications,” IEEE Transactions on Biomedical Engineering, vol. 46, no. 6, pp. 638–645, June 1999. [2] A. Rodriguez-Vazquez B.Linares-Barranco, E. Sanchez-Sinencio and J.L. Huertas, “A cmos implementation of ﬁtzhugh-nagumo neuron model,” IEEE Journal of Solid-State Circuits, vol. 26, no. 7, pp. 956–965, July 1991. [3] R.L. Calabrese G.N. Patel, G.S. Cymbalyuk and S.P. DeWeerth, “Bifurcation analysis of a silicon neuron,” in Proceedings of the Neural Information Processing Systems, 1999, pp. 731–737. [4] M.F. Simoni, G.S. Cymbalyuk, M.E. Sorensen, R.L. Calabrese, and S.P. DeWeerth, “A multiconductance silicon neuron with biologically matched dynamics,” IEEE Transactions on Biomedical Engineering, vol. 51, no. 2, pp. 342–354, Feb 2004. [5] E.M. Izhikevich, “Neural excitability, spiking and bursting,” International Journal of Bifurcations and Chaos, vol. 10, pp. 1171–1266, 2000. [6] E. Farquhar and P. Hasler, “A bio-physically inspired silicon neuron,” IEEE Transactions on Circuits and Systems I, vol. 52, no. 3, pp. 477– 488, March 2005. [7] J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields, Springer-Verlag, Applied mathematical Sciences 42, New York, 1983. [8] C.M. Twigg and P.E. Hasler, “A large-scale reconﬁgurable analog signal processor ic,” in Proceedings of the IEEE Custom Integrated Circuits Conference, Sept 2006, pp. 5–8.

431

Authorized licensed use limited to: Georgia Institute of Technology. Downloaded on July 08,2010 at 08:38:47 UTC from IEEE Xplore. Restrictions apply.