- Email: [email protected]
Phase dynamics in the biological neural networks
Physica A 288 (2000) 380–396
www.elsevier.com/locate/physa
Phase dynamics in the biological neural networks Seunghwan Kim ∗ , Sang Gui Lee Department of Physics, Pohang University of Science and Technology, Kyungpook, Pohang, 790-784, South Korea
Abstract The simpli ed models of neural networks based on biophysical Hodgkin–Huxley neurons are studied with a focus on coherent-phase dynamics. In our approach, each neuron is considered as a nonlinear oscillator, and collective dynamics of a mesoscopic network of neural oscillators are studied using the methods of nonlinear dynamics. We explore the mechanisms for synchrony, clustering and their breakup in the synaptic parameter space and discuss implications c 2000 Elsevier Science B.V. All rights to temporal aspects of neural-information processing. reserved. Keywords: Biological neural networks; Nonlinear oscillations; Synchrony; Phase models
1. Introduction Nonlinear oscillations have been frequently found in physics, chemistry, and biology including neurobiological problems. Examples include motion of a pendulum, chemical oscillations, oscillations in Josephson junction arrays, synchronization in re ies, and so on. Though the subject of synchrony in coupled oscillators has been studied for a long time, it has generated a renewed interest due to its implications in information processing in the brain [1– 4]. Neural oscillations have been observed in various parts of the nervous system, where rhythmic spike generation and synchronized neural oscillations occur in groups of interconnected neurons. The origins of oscillations and their possible functional roles in simple neural network models have been much studied in the context of neuro-information processing including sensory segmentation and binding, and motor control [5 –8]. ∗
Corresponding author. Fax: 82-562-279-3099. E-mail addresses: [email protected] (S. Kim), [email protected] (S.G. Lee).
c 2000 Elsevier Science B.V. All rights reserved. 0378-4371/00/$ - see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 0 0 ) 0 0 4 3 5 - 0
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
381
Recent extensive experiments on primary visual cortex suggest that synchronized oscillations could be crucial in feature linking in di erent parts of visual elds [1– 4]. The observation of 40 –50 Hz synchronous oscillations among the neurons widely separated over the visual cortex suggests that temporal correlation could be an important mechanism for feature binding [5,6]. As a result, it is becoming increasingly important to study the synchronization limit, opposite to the statistical limit much studied with models of arti cial neural networks, where the time (or phase) of individual action potentials plays a crucial role. In this paper, we take a dynamical approach regarding the neuron as a nonlinear oscillator and the neural network as a system of coupled oscillators. The main question here is to understand typical mechanisms for synchrony and determine the conditions for its breakup. When synchrony breaks down, clustering is often observed, for example, in a globally coupled network of nonlinear elements [9 –12] and a system with pulse-coupled oscillators [13]. Our approach is mesoscopic in the sense that we focus on dynamics of a few coupled neurons, which allow us a systematic study of the involved dynamical states and associated bifurcations in the system parameter space. We nd that the study of a small system of neurons, in particular, two and three coupled neurons, yields an important insight into phenomena of synchrony, multistability, and clustering that are often found in much larger networks. Our models also provide a simple neuro-paradigm for exploring the question of how macroscopic, collective, cooperative patterns emerge from microscopic interactions of nonlinear oscillatory elements. Since oscillatory neural networks can be considered as dynamical systems, we can take advantage of techniques of nonlinear dynamics in analyzing neural oscillatory dynamics including a powerful method of phase reduction in dealing with weakly coupled oscillators [14]. The properties of the synaptic coupling between neurons a ect collective nature of oscillations signi cantly. In particular, the synaptic time delay, for example, arising from the axonal propagation delay, is critical in neural systems. The time delay plays a key role in the post-synaptic inhibition in the hippocampus for delays of 2.5 –8.5 ms [15] and it is also critical in the ight system of the locust and the sound localization of an owl [16,17]. The time delay is crucial since it may completely modify the network dynamics in phase oscillators [18–20], providing a mechanism for desynchronization [21]. Recently the network with a time delay has been studied as a possible mechanism for the perception of ambiguous gures [22]. In this paper, we choose a Hodgkin–Huxley neuron model, which represents a typical biophysical model for tonically spiking neurons. First, we study dynamical responses of a single Hodgkin–Huxley neuron under various stimuli with a focus on phase dynamics and mode-locking phenomena. Next, we study a small network of a few Hodgkin–Huxley neurons that are globally coupled through synapses, focusing on the phenomena of synchronization and the role of synaptic time delay in its breakup. The techniques of phase reduction and bifurcation analysis have been extensively used to study the phase dynamics and explore the parameter space of the coupled neurons systematically.
382
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
2. Hodgkin–Huxley neuron models A typical neuron consists of three functionally di erent compartments; dendrites, a soma and an axon. Each neuron receives signals through dendrites and integrates these signals in the soma to produce a nerve impulse at the axon hillock. The nerve impulse, begins with a slight reduction in the negative potential across the membrane of the axon at its junction with the cell body, which opens sodium channels, shifting the voltage further. The in ow of sodium ions accelerates until the inner surface of the membrane is locally positive. The voltage reversal closes the sodium channels and opens the potassium channels. The out ow of potassium ions quickly restores the negative potential. This voltage reversal, called an action potential, has a height of about 100 mV and a duration of about 2 ms and propagates down the axon. After a brief refractory period, a second impulse can follow. When the train of action potentials reaches the terminal, called the synapse, neurotransmitters are released and di use through the gap, producing a potential change in a dendrite of the post-synaptic neuron. Dynamics of the membrane potential has been most thoroughly studied for a giant squid axon by Hodgkin and Huxley in 1952, and the resulting model of the neuron, called the Hodgkin–Huxley neuron model [23], serves as a simple, representative paradigm for tonically spiking neurons based on the voltage-dependent nonlinear membrane conductances. It also provides a good starting point for abstracting and building other kinds of neuron models, for example, the integrate- re model, the FitzHugh– Nagumo model [24,25], the Morris–Lecar model [26] and so on. In the Hodgkin–Huxley neuron model, dynamics of the membrane potential V is given by Kirchho ’s law of current conservation, dV = Iext + Iion + Isyn ; (1) dt where C is the membrane capacitance, Iext the external stimulus current, Iion the current through ion channels, and Isyn the synaptic current. The voltage-dependent ionic current is the sum of currents through the sodium channel (Na), the potassium channel (K), and the leakage channel (l). C
Iion = − gNa m3 h(V − VNa ) − gK n4 (V − VK ) − gl (V − Vl ) ;
(2)
where g’s are the maximal conductances for ion channels and V ’s with subscripts are the corresponding reversal potentials. The neurons are coupled to each other by synapses, the current through which is typically modeled by an -function-type coupling [27] that characterizes a quick rise and a long decay of the post-synaptic potential after the action potential spike arrives. That is, the synaptic current from a synapse is given by X (t − ti − d )(V (t) − Vsyn ) ; (3) Isyn (t) = − gsyn i
where the -function (t) = (t= )e−t= with the characteristic time of the synaptic current, the maximal synaptic conductance gsyn , the activation time ti of the ith potential
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
383
spike in the pre-synaptic neuron, the time delay d generated by the transmission time of the signal through the axon, the synaptic potential Vsyn , and the membrane potential V (t) in the post-synaptic neuron. The summation index i is for summing all the events of action potentials in the pre-synaptic neuron at the synapses. The characteristic time for the synapse, , is chosen to be 2 ms [28]. Typically, the synaptic coupling is called excitatory when Vsyn ¿ Veq and inhibitory when Vsyn ¡ Veq , with equilibrium potential Veq = − 65 mV in Hodgkin–Huxley neuron [29]. The time-varying gate variables, m, h, and n are associated with channel dynamics of sodium activation, sodium inactivation, and potassium activation, respectively, and satisfy the following equations: dx (4) x (V ) = x∞ (V ) − x ; dt where x denotes each of the gate variables, and x and x∞ are the relaxation time and the stationary value, respectively. The detailed values of the parameters can be found elsewhere [28,30]. Note that the Hodgkin–Huxley model is a set of nonlinear coupled ordinary di erential equations with four variables, V , m, n and h, which can be integrated using typical numerical integrators. We use the fourth-order Runge–Kutta method with the time step of 0:02 ms, which is sucient for most of our numerical experiments. 3. Dynamics of a single Hodgkin–Huxley neuron The biological neuron can be considered as a nonlinear oscillator as a rst dynamical approximation. Though neural oscillations can arise from various mechanisms, in the case of the Hodgkin–Huxley neuron they arise through a saddle-node bifurcation of limit cycles. Consider a single Hodgkin–Huxley neuron in natural seawater. It can re periodic action potentials of height 100 mV and duration about 2 ms when it is subjected to DC current stimuli above Idc = 6:2 A=cm2 . The frequency of action potential ring shows a slight dependence on Idc , which, for example, is about 68.5 Hz for Idc = 10 A=cm2 . Under AC current stimuli, a squid giant axon starts to re above a threshold amplitude and its ring rate can lock into the AC driving frequency producing complicated mode-locking structure and chaos [31,32]. Numerical study of the Hodgkin–Huxley model under periodic stimuli yields similar phenomena with ring onset transition and mode-locking ring responses [33], where the periodicity of rings is locked to that of the external stimulus in integer ratios. The dynamical responses of the single neuron can be characterized by the ÿring rate, which is de ned as the average number of spikes per period of the external stimulus, and the largest Lyapunov exponent. Using these order parameters, single neuron dynamics is studied in the parameter space of the stimulus. The phase diagram in the parameter space of the stimulus frequency ! and the stimulus amplitude A from numerical computation of the above order parameters is shown in Fig. 1(a), which contains regions with various mode locked states and chaotic states.
384
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
Fig. 1. (a) The phase diagram in the parameter space of the forcing amplitude A and the frequency ! for a single Hodgkin–Huxley neuron under sinusoidal current stimuli. Solid curves represent phase boundaries for the various mode-locked ring states, dotted curves those for the non ring region, and plus(+) symbols the parameter values for chaotic states. (b) Bifurcation structure of the phase diagram in (a). Solid curves denote inverse- ip bifurcations, dashed curves period-doubling cascades, and dotted curves saddle-node bifurcations, respectively.
The plus(+) symbols in the phase diagram in Fig. 1(a) represent the parameter values for chaotic dynamics determined from the positivity of the largest Lyapunov exponent. The plot of the ring rate as a function of the sinusoidal forcing frequency ! for a xed A gives a staircase-like structure, which is similar to the devil’s staircase in the typical quasi-periodic problem [34]. Each step in the staircase represents the mode-locked region with the same rational ring rate. The mode-locked regions in the
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
385
parameter space are denoted by the corresponding ratios of the period of the stimuli and the period of the ring. In Fig. 1(a), there exists a U -shaped ring onset curve above which the neuron starts to re with a non-zero ring rate. Note that the minimum threshold for the ring onset occurs near the natural frequency of the Hodgkin–Huxley neuron. This structure of the ring onset curve is important in understanding stochastic dynamics of the Hodgkin–Huxley neuron under noise [35]. Detailed numerical investigations reveal the existence of di erent ring onset mechanisms; for example, the non-chaotic ring onset and the chaotic ring onset. In the case of the chaotic ring onset, we nd at least three di erent routes to chaos; the period-doubling route to chaos and two types of the intermittency route to chaos [32]. The details of bifurcations at the phase boundaries can be analyzed through numerical computation of the stability of the periodic orbits. We nd that a major portion of the phase boundaries consist of “inverse ip (IF) bifurcations”, “period-doubling (PD) cascades”, and “saddle-node(SN) bifurcations” as in Fig. 1(b), which is in agreement with the phase diagram in Fig. 1(a). These bifurcation curves explain why di erent types of intermittency occurs between two mode-locked regions and why there exists a chaotic ring onset. 4. Dynamics of two and three coupled Hodgkin–Huxley neurons We take a mesoscopic approach with a few Hodgkin–Huxley neurons with synaptic coupling and ask rst how new dynamical possibilities arise and dynamical complexity grows as the system size increases. To address this, we rst study dynamics of two or three coupled Hodgkin–Huxley neurons under various current stimuli. For simplicity, we consider only the case that all neurons are globally coupled and are under identical suprathreshold DC current stimuli. A mesoscopic network of coupled neurons can display rich dynamical states characterized by phase di erences between spikes of the neurons, the details of which depend on various synaptic parameters including the synaptic time delay. 4.1. Basic phase states When two Hodgkin–Huxley neurons are coupled by a synapse, their ring times (or phases) can develop a certain relationship as in Fig. 2, where ring times of the action potentials are plotted as a sequence of bars (time coding) such as the synchronous state(S) with no phase di erence. Examples include various phase states the out-of-phase(OP) state with a non-zero phase di erence, or the anti-phase (AP) state ◦ with an exact 180 phase di erence. The synchronization between spikes of coupled neurons is the simplest, yet most important collective dynamics which can arise in a system with at least two neural oscillators. The synchronization typically occurs when the synaptic coupling is excitatory, that is, generation of an action potential in one neuron tends to induce a simultaneous ring of the other neuron. The system can be in the
386
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
Fig. 2. Firing time coding for typical phase states of two synaptically coupled Hodgkin–Huxley neurons. (a) The synchronous state (S), (b) the anti-phase state (AP), (c) the period doubled anti-phase state (PDA), and (d) the period four anti-phase state (PFA).
anti-phase state when the synaptic coupling is inhibitory. For strong enough coupling, synchrony may break up through di erent bifurcation processes and there occur other types of dynamical responses including the period doubled anti-phase state (PDA) in Fig. 2(c) and the period-four anti-phase state in Fig. 2(d), where the period of ring times is doubled through the period-doubling bifurcations, respectively. When three neurons are coupled globally, phase states can be even more complex with a hierarchy of clustering, where two or three neurons can be clustered into di erent groups of identical phases or ring rates [36,37]. This clustering phenomenon may play an important role in representation of neuro-information by a group of neurons, in particular, in temporal segmentation of overlapping patterns [38– 40]. Typical phase states of three globally coupled Hodgkin–Huxley neurons are shown in Fig. 3. As in the case of two coupled neurons, there exists a synchronous state (S) (Fig. 3(a)) when the coupling is excitatory. For inhibitory coupling, the clustered anti-phase (CAP) state may occur, where neurons are grouped into two synchronous neurons and the other in anti-phase (Fig. 3 (b)). The CAP state in three coupled neurons plays the role of the anti-phase state of two coupled neurons. In addition to these phase states, there exists a symmetric three-cluster state (SC) (Fig. 3(c)), where the phase of each neuron is shifted exactly by 2=3 [11]. Since the reduced phase space for the phase model of three weakly coupled Hodgkin–Huxley neurons is a torus, a state with quasi-periodic dynamics (QP) is possible for a non-zero time delay as in Fig. 3(d). In the quasi-periodic state, a slow switching between three CAP states may occur, which is also found in
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
387
Fig. 3. Firing time coding for typical phase states of three globally coupled Hodgkin–Huxley neurons. (a) The synchronous state (S), (b) the clustered anti-phase state (CAP), (c) the symmetric three-cluster state (SC), and (d) the quasi-periodic state (QP).
larger networks. For stronger coupling, there occurs a variety of nontrivial complex states including chaotic states that cannot be described by phase dynamics only. 4.2. Phase diagrams In the study of coupled neurons, an important order parameter used to label dynamics is the average phase di erence. The phase diagram in the parameter space of the time delay and the synaptic coupling strength for two Hodgkin–Huxley neurons coupled by an excitatory synapse with Vsyn = −30 mV is shown in Fig. 4(a), where a given symbol labels a dynamical state x and the curve SNx , a corresponding saddle-node bifurcation curve. As expected, for small synaptic coupling strength gsyn and small time delay d we nd a large region of synchrony. However, as the time delay is increased, the synchrony breaks down at the critical time delay of 2.3 ms, where a stable AP state is created by the saddle-node bifurcation denoted as SNAP in Fig. 4(a). Due to the coexistence of the synchronous state and the AP state, there exists a relatively large region of hysteresis involving these two states over a time delay of 2– 4 ms. When we increase the time delay further, the synchronous state breaks down at d ∼ 4 ms by the saddle-node bifurcation beyond which only the AP state is stable. Note that the critical value of the time delay for the creation of the AP state does not depend on gsyn much, but, the critical time delay for the breakup of the synchronous state depends strongly on gsyn . For large coupling, there exists another mechanism for breaking up synchrony through the period doubling bifurcations, where period-doubled anti-phase (PDA) state
388
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
Fig. 4. (a)The phase diagram for two coupled Hodgkin–Huxley neurons in the parameter space of the time delay d and the synaptic coupling strength gsyn . (b) The phase diagram for three globally coupled Hodgkin–Huxley neurons in the parameter space of the time delay and the synaptic coupling. The notations SNx denote saddle-node bifurcations associated with the state x, respectively.
and other complex states such as the period-four anti-phase state (PFA) or the chaotic state (C) can be found. The phase diagram for three coupled Hodgkin–Huxley neurons in the parameter space of gsyn and d is shown with Vsyn = − 30 mV in Fig 4(b). The global bifurcation structure in this case is similar to that of two coupled Hodgkin–Huxley neurons for small gsyn . For weak coupling, we nd the synchronous state for a small time delay and the CAP state for a large time delay. The saddle-node bifurcation curve of the synchronous state here is exactly the same as one for two coupled neurons. However, the saddle-node bifurcation curve of the CAP state is shifted to the right by about 1 ms. Note that the QP state is present in the region of d ∼ 2 − 3 ms, which is created by a Hopf bifurcation around the SC state. For a time delay of 2– 4 ms, there appear large regions of hysteresis between the synchronous state and the QP state and between the synchronous state and the CAP state. Note that the critical value of gsyn is about 0.3 above which a period doubling cascade(PDC) occurs and predictions of phase models fail qualitatively. Above this critical value, we can nd a variety of nontrivial complex states including chaotic states. 4.3. Phase model analysis A powerful phase reduction method [14], which allows a reduction in degrees of freedom for a system of weakly coupled oscillators, can be applied to our problem of coupled Hodgkin–Huxley neurons, predicting the existence of various phase states, and determining their stability properties. When the synaptic coupling is weak enough, the amplitude degree of freedom in a limit cycle can be neglected and the system can be reduced to phase oscillators with an e ective interaction given by the -function [14]. Using the phase reduction method, we get equations of motion for the asymptotic phases of the neuron, i ; i = 1; 2; : : : ; N , N d i X = dt j6=i
j i( i
−
j)
;
(5)
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
-function is given by Z 1 T dt j j ( t + i )gsyn dt {Vi ( t + i ( i − j) = − T 0 dV X j j − tspk t+ − d : ×
j
389
where the
i)
j − Vsyn }
(6)
tspk
j denote spiking times of action potentials in the pre-synaptic neuron The times tspk and the angular frequency of periodic oscillation in the HH neuron when they are uncoupled. Note that e ective phase dynamics is solely described by a -function, which is a 2-periodic function of phase di erences only. For two synaptically coupled Hodgkin–Huxley neurons with no time delay, we get evolution equations for the phase di erence = 1 − 2 from phase equations
d = dt
2 1 ( )
1 2 (− )
−
:
In the case of symmetric coupling, d =2 dt
− ( )
;
(7) 2 1
=
1 2
= , we get (8)
where − ( ) = 12 { ( ) − (− )}, the anti-symmetric part of . The phase states correspond to the stable xed points ∗ of Eq. (8), that is, the zeros of − ( ) with a negative slope. For symmetric coupling, the synchronized state ( ∗ = 0) and the anti-phase state ( ∗ = ± ) are always solutions, though their stability depends on synaptic parameters. The synchronized solution becomes unstable if the coupling or the dc input is not symmetric. In our numerical calculations, we choose Iext = Idc = 10 A=cm2 . As a phase origin, we selected an instant when the membrane potential crosses Veq = − 65 mV. The -function is calculated from a numerical convolution integral in Eq. (6) over one period, T . Next the -function is transformed into a Fourier series and truncated up to the rst ve modes, which produces an approximate -function within a numerical error of the order of 10−3 . The shape of the 2-periodic function depends on a number of parameters such as Vsyn , d , gsyn , and Iext . Note that for the symmetric coupling, gsyn produces no qualitative changes in analysis of the phase model since it a ects only the overall magnitude of . The time delay d simply causes a “phase shift” of the -functions in the negative direction with an amount of ! d , which may change the shape of − ( ) and, therefore, the stability of xed points. Fig. 5 shows the phase diagrams for asymptotically stable phase states in the parameter space of Vsyn and d . This gure is constructed by using the bisection method nding all stable xed points of − ( ) varying d from 0 to ! d = for a set of xed values of Vsyn . Each phase state is depicted in a di erent shade or pattern. The white denotes the region in the synchronized state, the black the anti-phase state, the gray the out-of-phase state, and the hatched regions the various bistable states with hysteresis.
390
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
Fig. 5. The phase diagram for two Hodgkin–Huxley neurons with synaptic coupling and identical dc input in the parameter space of Vsyn and d from the phase model analysis.
The phase diagram consists mostly of two basic phase states: the synchronized state, (S) and the anti-phase state, (AP). When d is suciently small, two neurons are synchronized under an excitatory coupling, while they are in the anti-phase state under an inhibitory coupling. As d increases, the phase state in the excitatory case exhibits a transition from the synchronous state to the AP state through bistability. This means that the nature of the coupling is changed from phase attraction to phase repulsion. For the inhibitory case, the anti-phase state makes a transition to the synchronous state via the out-of-phase state as d increases. The transition for the excitatory synapse is mainly induced by the sign change of the coecients of the rst and the second sine terms in the Fourier representation of − ( ), and one for the inhibitory synapse is due to that of higher modes. The phase model predicts that the typical critical time delay above which the nature of coupling changes is about 2– 4 ms depending on the direction of the parameter scan. Note that the critical time delay becomes smaller for larger synaptic coupling strength. The results of phase model analysis are consistent with the results of direct numerical simulations of two weakly coupled Hodgkin–Huxley neurons for gsyn = 0:1 [41]. In the case of three globally coupled Hodgkin–Huxley neurons, the equations of motion for phase di erences, 1 = 1 − 2 and 2 = 1 − 3 , for zero time delay
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
391
Fig. 6. The phase diagram for three globally coupled Hodgkin–Huxley neurons with identical dc stimuli in the parameter space of Vsyn and d from the phase model analysis.
d = 0 becomes 1 d 1 = ( ( 1 ) − (− 1 ) + ( 2 ) − ( 2 − 1 )) ; dt 2 d 2 1 = ( ( 2 ) − (− 2 ) + ( 1 ) − ( 1 − 2 )) ; dt 2
(9)
where the function is the same as one used for two coupled Hodgkin–Huxley neurons. If d 6= 0, the arguments of all -functions on the right-hand side of Eq. (9) are uniformly displaced by the amount d = ! d . Typical phase states for three coupled neurons are given by; ( 1 ; 2 ) = (0; 0) for the synchronous state, ( 1 ; 2 ) = (±; 0) or ( 1 ; 2 ) = (0; ±) for the clustered anti-phase (CAP) state, and ( 1 ; 2 ) = (±2=3; ∓2=3) for the symmetric three cluster(SC) state. Since the e ective phase space is two dimensional, a periodic orbit is possible and may produce quasi-periodic dynamics with a slow switching between the CAP states as in Fig. 3(d) [28,42]. Fig. 6 shows the phase diagram in the parameter space of Vsyn and d for the phase model of three globally coupled Hodgkin–Huxley neurons under dc stimuli of Idc = 10 A=cm2 . The basic phase states corresponding to stable xed points are found using the Newton–Raphson method, and the QP state is identi ed by direct integration. The visualization of phase portraits using the saddle-manifold method [42] has been ecient in locating the QP state and other aspects of the phase portraits. Note that
392
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
many of the phase states are located on the invariant manifolds, that is, under various symmetry operations. The regions of stability for basic phase states in Fig. 6 are denoted in di erent patterns and shades. Note that the overall structure of the phase diagram is somewhat similar to that of two coupled neurons. For excitatory coupling and small d , the typical phase state is the synchronous state. Though the SC state is also stable, its basin of attraction bounded by three saddle manifolds is much smaller than that of the synchronous state. The stability boundaries of the synchronous state is the same as those for two coupled neurons. For large enough d , the synchronous state becomes unstable, and turns into the clustered anti-phase state (CAP). However, the region in the parameter space for the CAP state is much smaller than its counterpart for two coupled neurons. The SC state may reappear for large d . For inhibitory coupling, the CAP state is a typical phase state. As d increases, a sequence of states including the quasi-periodic states, the symmetric 3-cluster state appears, which turns into the synchrony for large d . 5. Dynamics of many coupled Hodgkin–Huxley neurons In this section, we study the existence and the stability of basic phase states including the symmetric N -clustered state and the clustered anti-phase state in a large network of globally coupled phase oscillators. Suppose that N neural oscillators are grouped into n clusters with N=n oscillators per cluster. The equations of motion for each neural oscillator in the cluster k are given by 1 ˙ k = n
n−1 X
( l − k ) ;
(10)
l=0
where k = 0; 1; : : : ; n − 1. Then each equation has a solution with the phase k = !t + Pn−1 (2k=n). The stability of the symmetric N -cluster state can be 2k=n and ! = 1n determined by computing the following eigenvalues [11]: n−1 1 X 0 2k (N -n multiplicity) ; (11) = n n n−1
p =
1X n
0
2k n
(1 − e−2ikp=n ) (p = 0; 1; : : : ; n − 1) :
(12)
The eigenvalues in the rst equation are for the intra-cluster phase perturbations and the eigenvalue in the second equation for inter-cluster phase perturbations. Numerical stability analysis yields the stable regions for the symmetric N -cluster state with N = 3; 4; 5; 6 in the parameter space of the time delay and the synaptic reversal potential Vsyn in Fig. 7. Note that the symmetric 3-cluster state can be stable for no time delay. However, the symmetric 4-cluster state is stable for excitatory coupling with a large time delay or inhibitory coupling with an intermediate time delay. The symmetric 5 and 6 cluster states are more dicult to obtain since they can be stable only near the
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
393
Fig. 7. The regions in the parameter space of Vsyn and d for stable symmetric N -cluster states with (a) N = 3, (b) N = 4, (c) N = 5, and (d) N = 6 are shown for the phase model of N globally coupled Hodgkin–Huxley neurons.
rest potential and a large time delay. The size of the stability regions is sharply reduced as the number of the clusters increases. We have not found the symmetric N -cluster state with N ¿ 7 in the parameter region explored in Fig. 7. This may be due to the existence of the nite refractory period which makes the uniform distribution of the phases during one period of the ring more unstable. For example, for the symmetric 6-cluster state the time (phase) interval is about T=n ∼ 2:43 ms, which is comparable to the synaptic time constant 2 ms. The symmetric N -cluster state persists even under small symmetry-breaking perturbations and is expected to occur generically in large neural networks. Next we explore the existence and the stability of the clustered anti-phase (CAP) state in a large neural network through the phase model of the globally coupled Hodgkin–Huxley neurons. In this case, equations of motion for two phases 1 and 2 of two groups of phase-locked neural oscillators are given by N − N1 N1 − 1 (0) + ( 2 − 1 ) ; ˙ 1 = N N
(13)
N1 N − N1 − 1 ( 1 − 2 ) + (0) ; (14) ˙ 2 = N N where N1 and N2 are the numbers of oscillators in each cluster, respectively. If the CAP state is a solution, ˙ 1 = ˙ 2 . Therefore, we get the condition for the phase di erence
394
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
Fig. 8. The dependence of the phase di erence between two clusters on the time delay for a set of N1 =N = 0:1; 0:2; 0:3; 0:4, and 0:5 for (a) and (c) the excitatory synapse Vsyn = − 30 mV, and (b) and (d) the inhibitory synapse Vsyn = − 85 mV. The shade regions in (a), (b), (c), and (d) denote ones for stable clustered anti-phase (CAP) states.
= 1 − 2 , N − N1 N1 N − 2N1 (0) + (− ) + ( ) = 0 : N N N The solution of this equation for a set of N1 is calculated numerically using a root- nding method and its stability is identi ed from the eigenvalues 1 =
N1 N
0
(0) +
N − N1 N
0
(− ) ;
(15)
2 =
N1 N
0
( ) +
N − N1 N
0
(16)
(0) ;
N1 0 N − N1 0 ( ) + (− ) ; (17) N N where 1 and 2 are the stability multipliers for intra-cluster perturbations of the clusters 1 and 2, respectively, and 3 is the stability multiplier for inter-cluster phase perturbation. In Fig. 8, we show the dependence of the phase di erence between two clusters as a function of the time delay for a set of N1 for the excitatory synapse, Vsyn = − 30 mV (Fig. 8 (a)), and the inhibitory synapse Vsyn = − 85 mV (Fig. 8(b)). The shaded regions in Fig. 8 denote the regions for the stable CAP state. For the 3 =
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
395
excitatory synapse, the CAP state is not stable for small time delay. However, for the inhibitory synapse, the CAP state can be stable for zero or small time delay. Note that for both the excitatory and inhibitory synapse, the size of the shaded region for the stable CAP state gradually decreases to zero as N1 =N decreases, which is best seen in Fig. 8(c) and (d). 6. Summary and discussions The ring responses and collective phase dynamics of a Hodgkin–Huxley neuron under current stimuli and synaptically coupled Hodgkin–Huxley neurons have been studied through numerical simulations, phase reduction method, and nonlinear dynamical techniques. It would be interesting to study how our results on synchrony and clustering can be extended to a network with various connection and input con gurations, and how the mesoscopic results are generalized for larger networks. It would be also interesting to explore dynamical responses of the networks under various current stimuli including asymmetric, random, noisy, spatio-temporally correlated currents [28,43]. The vast majority of synapses in the mammalian cortex originate from within cortex. However, models of cortical development and processing do not re ect this basic fact of cortical anatomy and remain (almost) totally dominated by feedforward networks. On the basis of anatomy, a new paradigm for neuro-information processing with biological neural networks needs to be established, in which a small sensory input can be ampli ed via massive feedback and lateral connections. In particular, the conditions for reduction of complex realistic cell models into simple mathematical ones should be systematically studied, for example, through a comparative study of complex neuron models and their simpli ed versions. The phenomena of chaotic synchronization need to be explored further in connection with feature binding in cat visual cortex, etc. It would be also interesting to implement models of biological neural networks in analog VLSI and explore possibilities for their applications. Acknowledgements The work was supported by the Korea Research Foundation through the Basic Science Research Institute program and the Ministry of Science and Technology. References [1] C. Gray, W. Singer, Stimulus-Speci c Neuronal Oscillations in the Cat Visual Cortex, A Cortical Functional Unit, Soc. Neurosci. Abstracts 13 (1987) 404. [2] C. Gray, W. Singer, Proc. Natl. Acad. Sci. USA 86 (1989) 1698. [3] A.K. Engel, P. Konig, C.M. Gray, W. Singer, European J. Neurosci. 2 (1990) 588. [4] C.M. Gray, P. Koenig, A.K. Engel, W. Singer, Nature 338 (1989) 334.
396 [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43]
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396 C. von der Marlsburg, W. Schneider, Biol. Cybernet. 54 (1986) 29. H.G. Schuster, Nonlinear Dynamics and Neuronal Networks, VCH, New York, 1991. H.G. Schuster, P. Wagner, Biol. Cybernet. 64 (1990) 83. H. Sompolinsky, D. Golomb, D. Kleinfeld, Phys. Rev. A 43 (1991) 6990. K. Kaneko, Physica D 75 (1994) 74. V. Hakim, W.J. Rappel, Phys. Rev. A 46 (1992) R7347. K. Okuda, Physica D 63 (1993) 424. N. Nakagawa, Y. Kuramoto, Physica D 75 (1994) 74. U. Ernst, K. Pawelzik, T. Geisel, Phys. Rev. Lett. 74 (1995) 1570. Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence, Springer, New York, 1984. R. Miles, R.D. Traub, R.K.S. Wong, J. Neurophysiol. 60 (1988) 1471. M. Robertson, K.G. Pearson, J. Neurophysiol. 53 (1985) 110. C.E. Carr, M. Konishi, Proc. Natl. Acad. Sci. 85 (1988) 8311. E. Niebur, H.G. Schuster, D.M. Kammen, Phys. Rev. Lett. 67 (1991) 2753. H.G. Schuster, P. Wagner, Prog. Theor. Phys. 81 (1989) 939. S.H. Park, S. Kim, H.B. Pyo, S. Lee, Phys. Rev. E 60 (1999) 4962. C.V. Vrreswijk, L.F. Abbott, J. Comput. Neurosci. 1 (1994) 313. S. Kim, S.H. Park, C.S. Ryu, Phys. Rev. Lett. 79 (1997) 2911. A.L. Hodgkin, A.F. Huxley, J. Physiol. (London) 117 (1952) 500. R. Fitz Hugh, Biophys. J. 1 (1961) 445. J.S. Nagumo, S. Arimoto, S. Yoshizawa, Proc. IRE 50 (1961) 2061. C. Morris, H. Lecar, Biophys. J. 35 (1981) 193. E.R. Kandel, J.H. Schwartz, T.M. Jessel, Principles of Neural Science, 3rd Edition, Appleton & Lange, Norwalk, 1991. D. Hansel, G. Mato, C. Meunier, Phys. Rev. E 48 (1993) 3470. C. Koch, I. Segev, Methods in Neuronal Modeling, 2nd Edition, MIT Press, Cambridge, 1998. S.G. Lee, S. Kim, H. Kook, Int. J. Bifurc. Chaos 7 (1997) 889. K. Aihara, T. Numajiri, G. Matsumoto, M. Kotani, Phys. Lett. A 116 (1986) 313. N. Takahashi, Y. Hanyu, T. Musha, R. Kubo, G. Matsumoto, Physica D 43 (1990) 318. S.G. Lee, S. Kim, unpublished. P. Bak, Phys. Today 39 (1986) 38. S.G. Lee, S. Kim, Phys. Rev. E 60 (1999) 826. K. Kaneko, Physica D 41 (1990) 137. Y. Kuramoto, Prog. Theo. Phys. 87 (1992) 1119. D.L. Wang, J.M. Buhmann, C. von der Malsburg, Neural Comput. 2 (1990) 94. D. Horn, I. Opher, Neural Comput. 8 (1996) 373. S.K. Han, W.S. Kim, H. Kook, Phys. Rev. E 58 (1998) 2325. M.H. Park, S. Kim, J. Kor. Phys. Soc. 29 (1996) 9. C. Baesens, J. Guckenheimer, S. Kim, R.S. MacKay, Physica D 49 (1991) 387. E. Grannan, D. Kleifeld, H. Sompolinsky, Neural Comput. 4 (1993) 550.
www.elsevier.com/locate/physa
Phase dynamics in the biological neural networks Seunghwan Kim ∗ , Sang Gui Lee Department of Physics, Pohang University of Science and Technology, Kyungpook, Pohang, 790-784, South Korea
Abstract The simpli ed models of neural networks based on biophysical Hodgkin–Huxley neurons are studied with a focus on coherent-phase dynamics. In our approach, each neuron is considered as a nonlinear oscillator, and collective dynamics of a mesoscopic network of neural oscillators are studied using the methods of nonlinear dynamics. We explore the mechanisms for synchrony, clustering and their breakup in the synaptic parameter space and discuss implications c 2000 Elsevier Science B.V. All rights to temporal aspects of neural-information processing. reserved. Keywords: Biological neural networks; Nonlinear oscillations; Synchrony; Phase models
1. Introduction Nonlinear oscillations have been frequently found in physics, chemistry, and biology including neurobiological problems. Examples include motion of a pendulum, chemical oscillations, oscillations in Josephson junction arrays, synchronization in re ies, and so on. Though the subject of synchrony in coupled oscillators has been studied for a long time, it has generated a renewed interest due to its implications in information processing in the brain [1– 4]. Neural oscillations have been observed in various parts of the nervous system, where rhythmic spike generation and synchronized neural oscillations occur in groups of interconnected neurons. The origins of oscillations and their possible functional roles in simple neural network models have been much studied in the context of neuro-information processing including sensory segmentation and binding, and motor control [5 –8]. ∗
Corresponding author. Fax: 82-562-279-3099. E-mail addresses: [email protected] (S. Kim), [email protected] (S.G. Lee).
c 2000 Elsevier Science B.V. All rights reserved. 0378-4371/00/$ - see front matter PII: S 0 3 7 8 - 4 3 7 1 ( 0 0 ) 0 0 4 3 5 - 0
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
381
Recent extensive experiments on primary visual cortex suggest that synchronized oscillations could be crucial in feature linking in di erent parts of visual elds [1– 4]. The observation of 40 –50 Hz synchronous oscillations among the neurons widely separated over the visual cortex suggests that temporal correlation could be an important mechanism for feature binding [5,6]. As a result, it is becoming increasingly important to study the synchronization limit, opposite to the statistical limit much studied with models of arti cial neural networks, where the time (or phase) of individual action potentials plays a crucial role. In this paper, we take a dynamical approach regarding the neuron as a nonlinear oscillator and the neural network as a system of coupled oscillators. The main question here is to understand typical mechanisms for synchrony and determine the conditions for its breakup. When synchrony breaks down, clustering is often observed, for example, in a globally coupled network of nonlinear elements [9 –12] and a system with pulse-coupled oscillators [13]. Our approach is mesoscopic in the sense that we focus on dynamics of a few coupled neurons, which allow us a systematic study of the involved dynamical states and associated bifurcations in the system parameter space. We nd that the study of a small system of neurons, in particular, two and three coupled neurons, yields an important insight into phenomena of synchrony, multistability, and clustering that are often found in much larger networks. Our models also provide a simple neuro-paradigm for exploring the question of how macroscopic, collective, cooperative patterns emerge from microscopic interactions of nonlinear oscillatory elements. Since oscillatory neural networks can be considered as dynamical systems, we can take advantage of techniques of nonlinear dynamics in analyzing neural oscillatory dynamics including a powerful method of phase reduction in dealing with weakly coupled oscillators [14]. The properties of the synaptic coupling between neurons a ect collective nature of oscillations signi cantly. In particular, the synaptic time delay, for example, arising from the axonal propagation delay, is critical in neural systems. The time delay plays a key role in the post-synaptic inhibition in the hippocampus for delays of 2.5 –8.5 ms [15] and it is also critical in the ight system of the locust and the sound localization of an owl [16,17]. The time delay is crucial since it may completely modify the network dynamics in phase oscillators [18–20], providing a mechanism for desynchronization [21]. Recently the network with a time delay has been studied as a possible mechanism for the perception of ambiguous gures [22]. In this paper, we choose a Hodgkin–Huxley neuron model, which represents a typical biophysical model for tonically spiking neurons. First, we study dynamical responses of a single Hodgkin–Huxley neuron under various stimuli with a focus on phase dynamics and mode-locking phenomena. Next, we study a small network of a few Hodgkin–Huxley neurons that are globally coupled through synapses, focusing on the phenomena of synchronization and the role of synaptic time delay in its breakup. The techniques of phase reduction and bifurcation analysis have been extensively used to study the phase dynamics and explore the parameter space of the coupled neurons systematically.
382
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
2. Hodgkin–Huxley neuron models A typical neuron consists of three functionally di erent compartments; dendrites, a soma and an axon. Each neuron receives signals through dendrites and integrates these signals in the soma to produce a nerve impulse at the axon hillock. The nerve impulse, begins with a slight reduction in the negative potential across the membrane of the axon at its junction with the cell body, which opens sodium channels, shifting the voltage further. The in ow of sodium ions accelerates until the inner surface of the membrane is locally positive. The voltage reversal closes the sodium channels and opens the potassium channels. The out ow of potassium ions quickly restores the negative potential. This voltage reversal, called an action potential, has a height of about 100 mV and a duration of about 2 ms and propagates down the axon. After a brief refractory period, a second impulse can follow. When the train of action potentials reaches the terminal, called the synapse, neurotransmitters are released and di use through the gap, producing a potential change in a dendrite of the post-synaptic neuron. Dynamics of the membrane potential has been most thoroughly studied for a giant squid axon by Hodgkin and Huxley in 1952, and the resulting model of the neuron, called the Hodgkin–Huxley neuron model [23], serves as a simple, representative paradigm for tonically spiking neurons based on the voltage-dependent nonlinear membrane conductances. It also provides a good starting point for abstracting and building other kinds of neuron models, for example, the integrate- re model, the FitzHugh– Nagumo model [24,25], the Morris–Lecar model [26] and so on. In the Hodgkin–Huxley neuron model, dynamics of the membrane potential V is given by Kirchho ’s law of current conservation, dV = Iext + Iion + Isyn ; (1) dt where C is the membrane capacitance, Iext the external stimulus current, Iion the current through ion channels, and Isyn the synaptic current. The voltage-dependent ionic current is the sum of currents through the sodium channel (Na), the potassium channel (K), and the leakage channel (l). C
Iion = − gNa m3 h(V − VNa ) − gK n4 (V − VK ) − gl (V − Vl ) ;
(2)
where g’s are the maximal conductances for ion channels and V ’s with subscripts are the corresponding reversal potentials. The neurons are coupled to each other by synapses, the current through which is typically modeled by an -function-type coupling [27] that characterizes a quick rise and a long decay of the post-synaptic potential after the action potential spike arrives. That is, the synaptic current from a synapse is given by X (t − ti − d )(V (t) − Vsyn ) ; (3) Isyn (t) = − gsyn i
where the -function (t) = (t= )e−t= with the characteristic time of the synaptic current, the maximal synaptic conductance gsyn , the activation time ti of the ith potential
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
383
spike in the pre-synaptic neuron, the time delay d generated by the transmission time of the signal through the axon, the synaptic potential Vsyn , and the membrane potential V (t) in the post-synaptic neuron. The summation index i is for summing all the events of action potentials in the pre-synaptic neuron at the synapses. The characteristic time for the synapse, , is chosen to be 2 ms [28]. Typically, the synaptic coupling is called excitatory when Vsyn ¿ Veq and inhibitory when Vsyn ¡ Veq , with equilibrium potential Veq = − 65 mV in Hodgkin–Huxley neuron [29]. The time-varying gate variables, m, h, and n are associated with channel dynamics of sodium activation, sodium inactivation, and potassium activation, respectively, and satisfy the following equations: dx (4) x (V ) = x∞ (V ) − x ; dt where x denotes each of the gate variables, and x and x∞ are the relaxation time and the stationary value, respectively. The detailed values of the parameters can be found elsewhere [28,30]. Note that the Hodgkin–Huxley model is a set of nonlinear coupled ordinary di erential equations with four variables, V , m, n and h, which can be integrated using typical numerical integrators. We use the fourth-order Runge–Kutta method with the time step of 0:02 ms, which is sucient for most of our numerical experiments. 3. Dynamics of a single Hodgkin–Huxley neuron The biological neuron can be considered as a nonlinear oscillator as a rst dynamical approximation. Though neural oscillations can arise from various mechanisms, in the case of the Hodgkin–Huxley neuron they arise through a saddle-node bifurcation of limit cycles. Consider a single Hodgkin–Huxley neuron in natural seawater. It can re periodic action potentials of height 100 mV and duration about 2 ms when it is subjected to DC current stimuli above Idc = 6:2 A=cm2 . The frequency of action potential ring shows a slight dependence on Idc , which, for example, is about 68.5 Hz for Idc = 10 A=cm2 . Under AC current stimuli, a squid giant axon starts to re above a threshold amplitude and its ring rate can lock into the AC driving frequency producing complicated mode-locking structure and chaos [31,32]. Numerical study of the Hodgkin–Huxley model under periodic stimuli yields similar phenomena with ring onset transition and mode-locking ring responses [33], where the periodicity of rings is locked to that of the external stimulus in integer ratios. The dynamical responses of the single neuron can be characterized by the ÿring rate, which is de ned as the average number of spikes per period of the external stimulus, and the largest Lyapunov exponent. Using these order parameters, single neuron dynamics is studied in the parameter space of the stimulus. The phase diagram in the parameter space of the stimulus frequency ! and the stimulus amplitude A from numerical computation of the above order parameters is shown in Fig. 1(a), which contains regions with various mode locked states and chaotic states.
384
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
Fig. 1. (a) The phase diagram in the parameter space of the forcing amplitude A and the frequency ! for a single Hodgkin–Huxley neuron under sinusoidal current stimuli. Solid curves represent phase boundaries for the various mode-locked ring states, dotted curves those for the non ring region, and plus(+) symbols the parameter values for chaotic states. (b) Bifurcation structure of the phase diagram in (a). Solid curves denote inverse- ip bifurcations, dashed curves period-doubling cascades, and dotted curves saddle-node bifurcations, respectively.
The plus(+) symbols in the phase diagram in Fig. 1(a) represent the parameter values for chaotic dynamics determined from the positivity of the largest Lyapunov exponent. The plot of the ring rate as a function of the sinusoidal forcing frequency ! for a xed A gives a staircase-like structure, which is similar to the devil’s staircase in the typical quasi-periodic problem [34]. Each step in the staircase represents the mode-locked region with the same rational ring rate. The mode-locked regions in the
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
385
parameter space are denoted by the corresponding ratios of the period of the stimuli and the period of the ring. In Fig. 1(a), there exists a U -shaped ring onset curve above which the neuron starts to re with a non-zero ring rate. Note that the minimum threshold for the ring onset occurs near the natural frequency of the Hodgkin–Huxley neuron. This structure of the ring onset curve is important in understanding stochastic dynamics of the Hodgkin–Huxley neuron under noise [35]. Detailed numerical investigations reveal the existence of di erent ring onset mechanisms; for example, the non-chaotic ring onset and the chaotic ring onset. In the case of the chaotic ring onset, we nd at least three di erent routes to chaos; the period-doubling route to chaos and two types of the intermittency route to chaos [32]. The details of bifurcations at the phase boundaries can be analyzed through numerical computation of the stability of the periodic orbits. We nd that a major portion of the phase boundaries consist of “inverse ip (IF) bifurcations”, “period-doubling (PD) cascades”, and “saddle-node(SN) bifurcations” as in Fig. 1(b), which is in agreement with the phase diagram in Fig. 1(a). These bifurcation curves explain why di erent types of intermittency occurs between two mode-locked regions and why there exists a chaotic ring onset. 4. Dynamics of two and three coupled Hodgkin–Huxley neurons We take a mesoscopic approach with a few Hodgkin–Huxley neurons with synaptic coupling and ask rst how new dynamical possibilities arise and dynamical complexity grows as the system size increases. To address this, we rst study dynamics of two or three coupled Hodgkin–Huxley neurons under various current stimuli. For simplicity, we consider only the case that all neurons are globally coupled and are under identical suprathreshold DC current stimuli. A mesoscopic network of coupled neurons can display rich dynamical states characterized by phase di erences between spikes of the neurons, the details of which depend on various synaptic parameters including the synaptic time delay. 4.1. Basic phase states When two Hodgkin–Huxley neurons are coupled by a synapse, their ring times (or phases) can develop a certain relationship as in Fig. 2, where ring times of the action potentials are plotted as a sequence of bars (time coding) such as the synchronous state(S) with no phase di erence. Examples include various phase states the out-of-phase(OP) state with a non-zero phase di erence, or the anti-phase (AP) state ◦ with an exact 180 phase di erence. The synchronization between spikes of coupled neurons is the simplest, yet most important collective dynamics which can arise in a system with at least two neural oscillators. The synchronization typically occurs when the synaptic coupling is excitatory, that is, generation of an action potential in one neuron tends to induce a simultaneous ring of the other neuron. The system can be in the
386
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
Fig. 2. Firing time coding for typical phase states of two synaptically coupled Hodgkin–Huxley neurons. (a) The synchronous state (S), (b) the anti-phase state (AP), (c) the period doubled anti-phase state (PDA), and (d) the period four anti-phase state (PFA).
anti-phase state when the synaptic coupling is inhibitory. For strong enough coupling, synchrony may break up through di erent bifurcation processes and there occur other types of dynamical responses including the period doubled anti-phase state (PDA) in Fig. 2(c) and the period-four anti-phase state in Fig. 2(d), where the period of ring times is doubled through the period-doubling bifurcations, respectively. When three neurons are coupled globally, phase states can be even more complex with a hierarchy of clustering, where two or three neurons can be clustered into di erent groups of identical phases or ring rates [36,37]. This clustering phenomenon may play an important role in representation of neuro-information by a group of neurons, in particular, in temporal segmentation of overlapping patterns [38– 40]. Typical phase states of three globally coupled Hodgkin–Huxley neurons are shown in Fig. 3. As in the case of two coupled neurons, there exists a synchronous state (S) (Fig. 3(a)) when the coupling is excitatory. For inhibitory coupling, the clustered anti-phase (CAP) state may occur, where neurons are grouped into two synchronous neurons and the other in anti-phase (Fig. 3 (b)). The CAP state in three coupled neurons plays the role of the anti-phase state of two coupled neurons. In addition to these phase states, there exists a symmetric three-cluster state (SC) (Fig. 3(c)), where the phase of each neuron is shifted exactly by 2=3 [11]. Since the reduced phase space for the phase model of three weakly coupled Hodgkin–Huxley neurons is a torus, a state with quasi-periodic dynamics (QP) is possible for a non-zero time delay as in Fig. 3(d). In the quasi-periodic state, a slow switching between three CAP states may occur, which is also found in
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
387
Fig. 3. Firing time coding for typical phase states of three globally coupled Hodgkin–Huxley neurons. (a) The synchronous state (S), (b) the clustered anti-phase state (CAP), (c) the symmetric three-cluster state (SC), and (d) the quasi-periodic state (QP).
larger networks. For stronger coupling, there occurs a variety of nontrivial complex states including chaotic states that cannot be described by phase dynamics only. 4.2. Phase diagrams In the study of coupled neurons, an important order parameter used to label dynamics is the average phase di erence. The phase diagram in the parameter space of the time delay and the synaptic coupling strength for two Hodgkin–Huxley neurons coupled by an excitatory synapse with Vsyn = −30 mV is shown in Fig. 4(a), where a given symbol labels a dynamical state x and the curve SNx , a corresponding saddle-node bifurcation curve. As expected, for small synaptic coupling strength gsyn and small time delay d we nd a large region of synchrony. However, as the time delay is increased, the synchrony breaks down at the critical time delay of 2.3 ms, where a stable AP state is created by the saddle-node bifurcation denoted as SNAP in Fig. 4(a). Due to the coexistence of the synchronous state and the AP state, there exists a relatively large region of hysteresis involving these two states over a time delay of 2– 4 ms. When we increase the time delay further, the synchronous state breaks down at d ∼ 4 ms by the saddle-node bifurcation beyond which only the AP state is stable. Note that the critical value of the time delay for the creation of the AP state does not depend on gsyn much, but, the critical time delay for the breakup of the synchronous state depends strongly on gsyn . For large coupling, there exists another mechanism for breaking up synchrony through the period doubling bifurcations, where period-doubled anti-phase (PDA) state
388
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
Fig. 4. (a)The phase diagram for two coupled Hodgkin–Huxley neurons in the parameter space of the time delay d and the synaptic coupling strength gsyn . (b) The phase diagram for three globally coupled Hodgkin–Huxley neurons in the parameter space of the time delay and the synaptic coupling. The notations SNx denote saddle-node bifurcations associated with the state x, respectively.
and other complex states such as the period-four anti-phase state (PFA) or the chaotic state (C) can be found. The phase diagram for three coupled Hodgkin–Huxley neurons in the parameter space of gsyn and d is shown with Vsyn = − 30 mV in Fig 4(b). The global bifurcation structure in this case is similar to that of two coupled Hodgkin–Huxley neurons for small gsyn . For weak coupling, we nd the synchronous state for a small time delay and the CAP state for a large time delay. The saddle-node bifurcation curve of the synchronous state here is exactly the same as one for two coupled neurons. However, the saddle-node bifurcation curve of the CAP state is shifted to the right by about 1 ms. Note that the QP state is present in the region of d ∼ 2 − 3 ms, which is created by a Hopf bifurcation around the SC state. For a time delay of 2– 4 ms, there appear large regions of hysteresis between the synchronous state and the QP state and between the synchronous state and the CAP state. Note that the critical value of gsyn is about 0.3 above which a period doubling cascade(PDC) occurs and predictions of phase models fail qualitatively. Above this critical value, we can nd a variety of nontrivial complex states including chaotic states. 4.3. Phase model analysis A powerful phase reduction method [14], which allows a reduction in degrees of freedom for a system of weakly coupled oscillators, can be applied to our problem of coupled Hodgkin–Huxley neurons, predicting the existence of various phase states, and determining their stability properties. When the synaptic coupling is weak enough, the amplitude degree of freedom in a limit cycle can be neglected and the system can be reduced to phase oscillators with an e ective interaction given by the -function [14]. Using the phase reduction method, we get equations of motion for the asymptotic phases of the neuron, i ; i = 1; 2; : : : ; N , N d i X = dt j6=i
j i( i
−
j)
;
(5)
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
-function is given by Z 1 T dt j j ( t + i )gsyn dt {Vi ( t + i ( i − j) = − T 0 dV X j j − tspk t+ − d : ×
j
389
where the
i)
j − Vsyn }
(6)
tspk
j denote spiking times of action potentials in the pre-synaptic neuron The times tspk and the angular frequency of periodic oscillation in the HH neuron when they are uncoupled. Note that e ective phase dynamics is solely described by a -function, which is a 2-periodic function of phase di erences only. For two synaptically coupled Hodgkin–Huxley neurons with no time delay, we get evolution equations for the phase di erence = 1 − 2 from phase equations
d = dt
2 1 ( )
1 2 (− )
−
:
In the case of symmetric coupling, d =2 dt
− ( )
;
(7) 2 1
=
1 2
= , we get (8)
where − ( ) = 12 { ( ) − (− )}, the anti-symmetric part of . The phase states correspond to the stable xed points ∗ of Eq. (8), that is, the zeros of − ( ) with a negative slope. For symmetric coupling, the synchronized state ( ∗ = 0) and the anti-phase state ( ∗ = ± ) are always solutions, though their stability depends on synaptic parameters. The synchronized solution becomes unstable if the coupling or the dc input is not symmetric. In our numerical calculations, we choose Iext = Idc = 10 A=cm2 . As a phase origin, we selected an instant when the membrane potential crosses Veq = − 65 mV. The -function is calculated from a numerical convolution integral in Eq. (6) over one period, T . Next the -function is transformed into a Fourier series and truncated up to the rst ve modes, which produces an approximate -function within a numerical error of the order of 10−3 . The shape of the 2-periodic function depends on a number of parameters such as Vsyn , d , gsyn , and Iext . Note that for the symmetric coupling, gsyn produces no qualitative changes in analysis of the phase model since it a ects only the overall magnitude of . The time delay d simply causes a “phase shift” of the -functions in the negative direction with an amount of ! d , which may change the shape of − ( ) and, therefore, the stability of xed points. Fig. 5 shows the phase diagrams for asymptotically stable phase states in the parameter space of Vsyn and d . This gure is constructed by using the bisection method nding all stable xed points of − ( ) varying d from 0 to ! d = for a set of xed values of Vsyn . Each phase state is depicted in a di erent shade or pattern. The white denotes the region in the synchronized state, the black the anti-phase state, the gray the out-of-phase state, and the hatched regions the various bistable states with hysteresis.
390
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
Fig. 5. The phase diagram for two Hodgkin–Huxley neurons with synaptic coupling and identical dc input in the parameter space of Vsyn and d from the phase model analysis.
The phase diagram consists mostly of two basic phase states: the synchronized state, (S) and the anti-phase state, (AP). When d is suciently small, two neurons are synchronized under an excitatory coupling, while they are in the anti-phase state under an inhibitory coupling. As d increases, the phase state in the excitatory case exhibits a transition from the synchronous state to the AP state through bistability. This means that the nature of the coupling is changed from phase attraction to phase repulsion. For the inhibitory case, the anti-phase state makes a transition to the synchronous state via the out-of-phase state as d increases. The transition for the excitatory synapse is mainly induced by the sign change of the coecients of the rst and the second sine terms in the Fourier representation of − ( ), and one for the inhibitory synapse is due to that of higher modes. The phase model predicts that the typical critical time delay above which the nature of coupling changes is about 2– 4 ms depending on the direction of the parameter scan. Note that the critical time delay becomes smaller for larger synaptic coupling strength. The results of phase model analysis are consistent with the results of direct numerical simulations of two weakly coupled Hodgkin–Huxley neurons for gsyn = 0:1 [41]. In the case of three globally coupled Hodgkin–Huxley neurons, the equations of motion for phase di erences, 1 = 1 − 2 and 2 = 1 − 3 , for zero time delay
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
391
Fig. 6. The phase diagram for three globally coupled Hodgkin–Huxley neurons with identical dc stimuli in the parameter space of Vsyn and d from the phase model analysis.
d = 0 becomes 1 d 1 = ( ( 1 ) − (− 1 ) + ( 2 ) − ( 2 − 1 )) ; dt 2 d 2 1 = ( ( 2 ) − (− 2 ) + ( 1 ) − ( 1 − 2 )) ; dt 2
(9)
where the function is the same as one used for two coupled Hodgkin–Huxley neurons. If d 6= 0, the arguments of all -functions on the right-hand side of Eq. (9) are uniformly displaced by the amount d = ! d . Typical phase states for three coupled neurons are given by; ( 1 ; 2 ) = (0; 0) for the synchronous state, ( 1 ; 2 ) = (±; 0) or ( 1 ; 2 ) = (0; ±) for the clustered anti-phase (CAP) state, and ( 1 ; 2 ) = (±2=3; ∓2=3) for the symmetric three cluster(SC) state. Since the e ective phase space is two dimensional, a periodic orbit is possible and may produce quasi-periodic dynamics with a slow switching between the CAP states as in Fig. 3(d) [28,42]. Fig. 6 shows the phase diagram in the parameter space of Vsyn and d for the phase model of three globally coupled Hodgkin–Huxley neurons under dc stimuli of Idc = 10 A=cm2 . The basic phase states corresponding to stable xed points are found using the Newton–Raphson method, and the QP state is identi ed by direct integration. The visualization of phase portraits using the saddle-manifold method [42] has been ecient in locating the QP state and other aspects of the phase portraits. Note that
392
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
many of the phase states are located on the invariant manifolds, that is, under various symmetry operations. The regions of stability for basic phase states in Fig. 6 are denoted in di erent patterns and shades. Note that the overall structure of the phase diagram is somewhat similar to that of two coupled neurons. For excitatory coupling and small d , the typical phase state is the synchronous state. Though the SC state is also stable, its basin of attraction bounded by three saddle manifolds is much smaller than that of the synchronous state. The stability boundaries of the synchronous state is the same as those for two coupled neurons. For large enough d , the synchronous state becomes unstable, and turns into the clustered anti-phase state (CAP). However, the region in the parameter space for the CAP state is much smaller than its counterpart for two coupled neurons. The SC state may reappear for large d . For inhibitory coupling, the CAP state is a typical phase state. As d increases, a sequence of states including the quasi-periodic states, the symmetric 3-cluster state appears, which turns into the synchrony for large d . 5. Dynamics of many coupled Hodgkin–Huxley neurons In this section, we study the existence and the stability of basic phase states including the symmetric N -clustered state and the clustered anti-phase state in a large network of globally coupled phase oscillators. Suppose that N neural oscillators are grouped into n clusters with N=n oscillators per cluster. The equations of motion for each neural oscillator in the cluster k are given by 1 ˙ k = n
n−1 X
( l − k ) ;
(10)
l=0
where k = 0; 1; : : : ; n − 1. Then each equation has a solution with the phase k = !t + Pn−1 (2k=n). The stability of the symmetric N -cluster state can be 2k=n and ! = 1n determined by computing the following eigenvalues [11]: n−1 1 X 0 2k (N -n multiplicity) ; (11) = n n n−1
p =
1X n
0
2k n
(1 − e−2ikp=n ) (p = 0; 1; : : : ; n − 1) :
(12)
The eigenvalues in the rst equation are for the intra-cluster phase perturbations and the eigenvalue in the second equation for inter-cluster phase perturbations. Numerical stability analysis yields the stable regions for the symmetric N -cluster state with N = 3; 4; 5; 6 in the parameter space of the time delay and the synaptic reversal potential Vsyn in Fig. 7. Note that the symmetric 3-cluster state can be stable for no time delay. However, the symmetric 4-cluster state is stable for excitatory coupling with a large time delay or inhibitory coupling with an intermediate time delay. The symmetric 5 and 6 cluster states are more dicult to obtain since they can be stable only near the
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
393
Fig. 7. The regions in the parameter space of Vsyn and d for stable symmetric N -cluster states with (a) N = 3, (b) N = 4, (c) N = 5, and (d) N = 6 are shown for the phase model of N globally coupled Hodgkin–Huxley neurons.
rest potential and a large time delay. The size of the stability regions is sharply reduced as the number of the clusters increases. We have not found the symmetric N -cluster state with N ¿ 7 in the parameter region explored in Fig. 7. This may be due to the existence of the nite refractory period which makes the uniform distribution of the phases during one period of the ring more unstable. For example, for the symmetric 6-cluster state the time (phase) interval is about T=n ∼ 2:43 ms, which is comparable to the synaptic time constant 2 ms. The symmetric N -cluster state persists even under small symmetry-breaking perturbations and is expected to occur generically in large neural networks. Next we explore the existence and the stability of the clustered anti-phase (CAP) state in a large neural network through the phase model of the globally coupled Hodgkin–Huxley neurons. In this case, equations of motion for two phases 1 and 2 of two groups of phase-locked neural oscillators are given by N − N1 N1 − 1 (0) + ( 2 − 1 ) ; ˙ 1 = N N
(13)
N1 N − N1 − 1 ( 1 − 2 ) + (0) ; (14) ˙ 2 = N N where N1 and N2 are the numbers of oscillators in each cluster, respectively. If the CAP state is a solution, ˙ 1 = ˙ 2 . Therefore, we get the condition for the phase di erence
394
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
Fig. 8. The dependence of the phase di erence between two clusters on the time delay for a set of N1 =N = 0:1; 0:2; 0:3; 0:4, and 0:5 for (a) and (c) the excitatory synapse Vsyn = − 30 mV, and (b) and (d) the inhibitory synapse Vsyn = − 85 mV. The shade regions in (a), (b), (c), and (d) denote ones for stable clustered anti-phase (CAP) states.
= 1 − 2 , N − N1 N1 N − 2N1 (0) + (− ) + ( ) = 0 : N N N The solution of this equation for a set of N1 is calculated numerically using a root- nding method and its stability is identi ed from the eigenvalues 1 =
N1 N
0
(0) +
N − N1 N
0
(− ) ;
(15)
2 =
N1 N
0
( ) +
N − N1 N
0
(16)
(0) ;
N1 0 N − N1 0 ( ) + (− ) ; (17) N N where 1 and 2 are the stability multipliers for intra-cluster perturbations of the clusters 1 and 2, respectively, and 3 is the stability multiplier for inter-cluster phase perturbation. In Fig. 8, we show the dependence of the phase di erence between two clusters as a function of the time delay for a set of N1 for the excitatory synapse, Vsyn = − 30 mV (Fig. 8 (a)), and the inhibitory synapse Vsyn = − 85 mV (Fig. 8(b)). The shaded regions in Fig. 8 denote the regions for the stable CAP state. For the 3 =
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396
395
excitatory synapse, the CAP state is not stable for small time delay. However, for the inhibitory synapse, the CAP state can be stable for zero or small time delay. Note that for both the excitatory and inhibitory synapse, the size of the shaded region for the stable CAP state gradually decreases to zero as N1 =N decreases, which is best seen in Fig. 8(c) and (d). 6. Summary and discussions The ring responses and collective phase dynamics of a Hodgkin–Huxley neuron under current stimuli and synaptically coupled Hodgkin–Huxley neurons have been studied through numerical simulations, phase reduction method, and nonlinear dynamical techniques. It would be interesting to study how our results on synchrony and clustering can be extended to a network with various connection and input con gurations, and how the mesoscopic results are generalized for larger networks. It would be also interesting to explore dynamical responses of the networks under various current stimuli including asymmetric, random, noisy, spatio-temporally correlated currents [28,43]. The vast majority of synapses in the mammalian cortex originate from within cortex. However, models of cortical development and processing do not re ect this basic fact of cortical anatomy and remain (almost) totally dominated by feedforward networks. On the basis of anatomy, a new paradigm for neuro-information processing with biological neural networks needs to be established, in which a small sensory input can be ampli ed via massive feedback and lateral connections. In particular, the conditions for reduction of complex realistic cell models into simple mathematical ones should be systematically studied, for example, through a comparative study of complex neuron models and their simpli ed versions. The phenomena of chaotic synchronization need to be explored further in connection with feature binding in cat visual cortex, etc. It would be also interesting to implement models of biological neural networks in analog VLSI and explore possibilities for their applications. Acknowledgements The work was supported by the Korea Research Foundation through the Basic Science Research Institute program and the Ministry of Science and Technology. References [1] C. Gray, W. Singer, Stimulus-Speci c Neuronal Oscillations in the Cat Visual Cortex, A Cortical Functional Unit, Soc. Neurosci. Abstracts 13 (1987) 404. [2] C. Gray, W. Singer, Proc. Natl. Acad. Sci. USA 86 (1989) 1698. [3] A.K. Engel, P. Konig, C.M. Gray, W. Singer, European J. Neurosci. 2 (1990) 588. [4] C.M. Gray, P. Koenig, A.K. Engel, W. Singer, Nature 338 (1989) 334.
396 [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43]
S. Kim, S.G. Lee / Physica A 288 (2000) 380–396 C. von der Marlsburg, W. Schneider, Biol. Cybernet. 54 (1986) 29. H.G. Schuster, Nonlinear Dynamics and Neuronal Networks, VCH, New York, 1991. H.G. Schuster, P. Wagner, Biol. Cybernet. 64 (1990) 83. H. Sompolinsky, D. Golomb, D. Kleinfeld, Phys. Rev. A 43 (1991) 6990. K. Kaneko, Physica D 75 (1994) 74. V. Hakim, W.J. Rappel, Phys. Rev. A 46 (1992) R7347. K. Okuda, Physica D 63 (1993) 424. N. Nakagawa, Y. Kuramoto, Physica D 75 (1994) 74. U. Ernst, K. Pawelzik, T. Geisel, Phys. Rev. Lett. 74 (1995) 1570. Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence, Springer, New York, 1984. R. Miles, R.D. Traub, R.K.S. Wong, J. Neurophysiol. 60 (1988) 1471. M. Robertson, K.G. Pearson, J. Neurophysiol. 53 (1985) 110. C.E. Carr, M. Konishi, Proc. Natl. Acad. Sci. 85 (1988) 8311. E. Niebur, H.G. Schuster, D.M. Kammen, Phys. Rev. Lett. 67 (1991) 2753. H.G. Schuster, P. Wagner, Prog. Theor. Phys. 81 (1989) 939. S.H. Park, S. Kim, H.B. Pyo, S. Lee, Phys. Rev. E 60 (1999) 4962. C.V. Vrreswijk, L.F. Abbott, J. Comput. Neurosci. 1 (1994) 313. S. Kim, S.H. Park, C.S. Ryu, Phys. Rev. Lett. 79 (1997) 2911. A.L. Hodgkin, A.F. Huxley, J. Physiol. (London) 117 (1952) 500. R. Fitz Hugh, Biophys. J. 1 (1961) 445. J.S. Nagumo, S. Arimoto, S. Yoshizawa, Proc. IRE 50 (1961) 2061. C. Morris, H. Lecar, Biophys. J. 35 (1981) 193. E.R. Kandel, J.H. Schwartz, T.M. Jessel, Principles of Neural Science, 3rd Edition, Appleton & Lange, Norwalk, 1991. D. Hansel, G. Mato, C. Meunier, Phys. Rev. E 48 (1993) 3470. C. Koch, I. Segev, Methods in Neuronal Modeling, 2nd Edition, MIT Press, Cambridge, 1998. S.G. Lee, S. Kim, H. Kook, Int. J. Bifurc. Chaos 7 (1997) 889. K. Aihara, T. Numajiri, G. Matsumoto, M. Kotani, Phys. Lett. A 116 (1986) 313. N. Takahashi, Y. Hanyu, T. Musha, R. Kubo, G. Matsumoto, Physica D 43 (1990) 318. S.G. Lee, S. Kim, unpublished. P. Bak, Phys. Today 39 (1986) 38. S.G. Lee, S. Kim, Phys. Rev. E 60 (1999) 826. K. Kaneko, Physica D 41 (1990) 137. Y. Kuramoto, Prog. Theo. Phys. 87 (1992) 1119. D.L. Wang, J.M. Buhmann, C. von der Malsburg, Neural Comput. 2 (1990) 94. D. Horn, I. Opher, Neural Comput. 8 (1996) 373. S.K. Han, W.S. Kim, H. Kook, Phys. Rev. E 58 (1998) 2325. M.H. Park, S. Kim, J. Kor. Phys. Soc. 29 (1996) 9. C. Baesens, J. Guckenheimer, S. Kim, R.S. MacKay, Physica D 49 (1991) 387. E. Grannan, D. Kleifeld, H. Sompolinsky, Neural Comput. 4 (1993) 550.