9512.net
甜梦文库
当前位置:首页 >> >>

Abstract Physica D 118 (1998) 343–370 Phase-locking in weakly heterogeneous neuronal netwo


Physica D 118 (1998) 343–370

Phase-locking in weakly heterogeneous neuronal networks
Carson C. Chow 1
Department of Mathematics and Center for BioDynamics, Boston University, Boston, MA 02215, USA Received 22 September 1997; received in revised form 5 February 1998; accepted 17 February 1998 Communicated by J.D. Meiss

Abstract We examine analytically the existence and stability of phase-locked states in a weakly heterogeneous neuronal network. We consider a model of N neurons with all-to-all synaptic coupling where the heterogeneity is in the intrinsic ?ring frequency of the individual neurons. We consider both inhibitory and excitatory coupling. We derive the conditions under which stable phase-locking is possible. In homogeneous networks, many different periodic phase-locked states are possible. Their stability depends on the dynamics of the neuron and the coupling. For weak heterogeneity, the phase-locked states are perturbed from the homogeneous states and can remain stable if their homogeneous counterparts are stable. For enough heterogeneity, phase-locked solutions either lose stability or are destroyed completely. We analyze the possible states the network can take when phase-locking is broken. ? 1998 Elsevier Science B.V.
PACS: 87.10.+e; 87.22.Jb; 87.22.As; 05.45.+b Keywords: Neural networks; Synchronization; Phase-locking; Oscillations

1. Introduction Phase-locked and synchronized neuronal ?ring is found throughout the brain and this coherent activity may play an important role in behavior and cognition [1–5]. Given their importance, it is of interest to understand the mechanisms responsible for their generation. Here, we consider a network of heterogeneous neurons that are coupled synaptically and determine analytically the conditions under which their ?ring patterns will be phase-locked. Our results could have implications for synchronous phenomena in other areas of biology and physics that are based on pulse-coupled oscillators [6–8]. Much of the previous analytical work on phase-locking and synchronization in pulse-coupled neuronal oscillators have focused on homogeneous networks [9–15]. However, any biophysically relevant neuronal network will likely include heterogeneous effects. It is known that heterogeneity can greatly reduce synchrony in networks of phasecoupled oscillators [6,16,17]. Coexistence of synchronous and asynchronous states has been shown to exist in large networks of heterogeneous pulse-coupled oscillators [18,19]. Recently it has been shown numerically that
1 After

September 1998: Department of Mathematics, University of Pittsburgh, Pittsburgh, PA 15260, USA.

0167-2789/98/$19.00 Copyright ? 1998 Elsevier Science B.V. All rights reserved PII S 0 1 6 7 - 2 7 8 9 ( 9 8 ) 0 0 0 8 2 - 7

344

C.C. Chow / Physica D 118 (1998) 343–370

heterogeneous assemblies of neurons can ?re coherently. However, the coherent activity can be greatly reduced even when the heterogeneity is very small [20,21]. Here, we examine analytically the existence and stability of phase-locked solutions of weakly heterogeneous neuronal networks. We consider an all-to-all coupled network of size N where the heterogeneity is in the intrinsic ?ring frequency of the component neurons. We examine both inhibitory and excitatory connections. Phase-locking in homogeneous synaptically coupled neuronal networks has been studied analytically with different formalisms [10–12,22]. Here, we generalize the formalism developed by Gerstner et al. [12] to examine the existence and stability of phase-locked solutions. We ?nd the conditions for stability of homogeneous phase-locked solutions and examine what happens as we allow the heterogeneity to increase. Our formalism can analyze all periodic phase-locked states but we mainly focus on the synchronized or in-phase state, the anti-synchronized state, and the splay-phase or anti-phase state. We ?nd that stable phase-locking depends crucially on the time course of the synaptic coupling and on the response properties of the neuron to synaptic inputs. Neurons have recently been classi?ed in accordance with their phase response properties [11,23]. Type I neurons have the property that their phase response curves are always positive, i.e. the next ?ring time of the neuron is always advanced in response to a positive (depolarizing) input. In contrast, for Type II neurons the next ?ring time is either delayed or advanced depending on when the input is received. Our results con?rm previous observations that for Type I neurons, excitatory coupling is generally desynchronizing while inhibitory coupling is synchronizing [10–12]. For Type II neurons, the opposite can be true although not necessarily. In previous work, the requirements on the rise and decay time of the synaptic coupling for stable phase-locking have been separately emphasized [10,12,22]. We ?nd that depending on the nature of the synaptic coupling conditions on both are required. An interesting result for homogeneous neurons is that the splay-phase or anti-phase state can be stable for inhibitory coupling for any ?nite sized network. For in?nite N , the splay-phase state is known to be unstable [26,29]. With weak heterogeneity, phase-locked solutions can exist but with small phase dispersion around their homogeneous counterparts. Stability is possible if the original homogeneous solutions are stable. The addition of heterogeneity cannot stabilize an unstable homogeneous phase-locked state. As the heterogeneity is increased, phase-locked states can either lose stability or cease to exist. Heterogeneity can also foster clustering. When stable phase-locking is lost, the network can enter asynchronous, harmonically locked or suppressed states. Fig. 1 shows examples of various states for two coupled inhibitory neurons obtained from numerical simulations of a biophysically based neuron model given in Appendix A. In Section 2 we introduce the model and formalism we use for our analysis. In Section 3 we give the conditions which must be satis?ed for a periodic phase-locked solution. In Section 4 we give general conditions for stability of locked solutions. In Section 5 we apply our conditions for the existence and stability of phase-locked states to two coupled homogeneous and heterogeneous neurons. We also examine suppression and harmonic locking. In Section 6 we examine phase-locking and other dynamics in a network of N neurons. We discuss our results and conclude in Section 7.

2. Model Our network consists of N all-to-all coupled neurons. Each time a given neuron ?res, a synaptic pulse is transmitted to all the other neurons. The effect of the synapse remains over a ?nite amount of time. The dynamics of neurons are generally described by biophysically based (Hodgkin–Huxley type) equations which are a set of differential equations for the current balance of the membrane potential and for the dynamics of the component ion channels [24]. In Appendix A, we give an example set of biophysically based equations [21]. We use these equations in numerical

C.C. Chow / Physica D 118 (1998) 343–370

345

Fig. 1. Examples of possible states for two heterogeneous neurons from numerical simulations of the model given in Appendix A. (a) Asynchrony with I1 = 9.0, I2 = 9.9 ?A/cm2 ; gs = 0.25 mS/cm2 ; τs = 10 ms. (b) Near-synchrony with I1 = 1.6, I2 = 1.78 ?A/cm2 ; gs = 0.25 mS/cm2 ; τs = 10 ms. (c) Harmonic locking with I1 = 9.0, I2 = 9.9 ?A/cm2 ; gs = 0.5 mS/cm2 ; τs = 10 ms. (d) Suppression with I1 = 1.6, I2 = 1.78 ?A/cm2 ; gs = 0.5 mS/cm2 ; τs = 10 ms.

simulations to con?rm our analytical results. These neurons are Type I in contrast to classic Hodgkin–Huxley equations which are Type II [11,23]. Generally, the differential equations for neuron models are complicated and cumbersome for theoretical analysis. Instead, we use the spike response method for our neuronal dynamics [12,25,26]. This formalism is fairly general and can be applied to biophysical neuron models [27]. In this method, the neuron is considered to be a device

346

C.C. Chow / Physica D 118 (1998) 343–370

which ?res when a single scalar variable (membrane potential) crosses a threshold from below. The dynamics of the membrane potential consists of a sum of two sets of kernels akin to an integral expansion: vi (t) =
l

ηi (t ? tli ) +

ij (t jl

? tl ),

j

(2.1)

where tli marks the ?ring times of neuron i (i.e. the l th time neuron i reaches the threshold θi ) and j is summed from 1 to N excluding i . The kernels are zero for t ? tlk ≤ 0. The kernel ηi (t) describes the intrinsic dynamics of the neuron after it has been triggered to ?re. It represents the ensuing spike (action potential) and recovery behavior. The kernel ij (t) represents the response due to synaptic inputs. We do not include the effects of self-coupling in (2.1) but can do so by modifying ηi (t). Fig. 3 gives example sketches of possible kernels. To linear order we can write ij (t) in the form
t ij (t)

=
0

G(s)Iij (t ? s) ds,

(2.2)

where Iij (t) is the synaptic input to neuron i from neuron j , and G(t) is the linear response function. It has been shown that the linear response is often an adequate approximation for the full nonlinear response [12,26,27]. In biophysically based neurons, the synaptic input current is usually given by Iij (t) = gs Sij (t)(vi ? vs ) where vs is a constant reversal potential, gs is the maximal synaptic conductance, and Sij (t) is the synaptic gating variable obtained from a differential equation (see Eq. (A.4) in Appendix A). Models often use the simpli?cation Iij (t) gSij (t) since (vi ? vs ) is approximately constant away from a spike [11]. The dynamics of the synaptic gating variable has previously been approximated by a recursive scheme where Sij (t) is updated each time the neuron receives an input with Sij (t) → Sij (t) + Sf (t ? t j ), where t j is the time of the pre-synaptic spike, and Sf (t) is a stereotypical synaptic response [10,11,29]. We call this type of synaptic model nonsaturating since Sij (t) increases without bound as the ?ring rate increases. On the other hand, the synaptic model given by (A.4) is saturating. Here Sij (t) saturates to a ?nite value as the ?ring rate of the pre-synaptic neuron increases. In the limit of very fast synaptic rise times we can use the approximation S(t) → Sf (t ? t j ) to model the dynamics of (A.4) [28]. In this case the synapse has no memory of previous activity. In our network, we allow for heterogeneity in the thresholds of the neurons which implies heterogeneity in the intrinsic ?ring frequencies. We maintain homogeneity in the time courses of the kernels. We can make the threshold heterogeneity explicit with the transformation vi → vi ? Ii and θi → θ ? Ii , where Ii can be thought of as a normalized applied driving current, and θ is a constant threshold. We make the synaptic strengths of the neurons explicit by setting ij (t) = J (t), where J represents the coupling strength. Heterogeneous coupling could be accounted for by setting J = Jij . When J > 0 we call this excitatory coupling and for J < 0 we call this inhibitory coupling. Putting this in Eq. (2.1) yields vi (t) = Ii +
l

η(t ? tli ) +

j,l

J (t ? tl ).

j

(2.3)

In this form, the network heterogeneity is explicitly expressed in the applied current Ii . We can set the threshold for all the oscillators to 1 by a simple rescaling without any loss of generality. In this formulation, Iij (t) in Eq. (2.2) can be thought of as the post-synaptic current while ij (t) can be thought of as the post-synaptic potential. Often we will refer to the rise time of the synaptic current which will mean the rise time of Iij (t). We note that the synaptic kernel ij (t) will itself have a rise time which can depend on both the rise and decay time of the synaptic current.

C.C. Chow / Physica D 118 (1998) 343–370

347

Remark 1. We have assumed that the synaptic kernel does not depend on the time of arrival or phase of the synaptic input. For the integrate-and-?re model, as will be shown below, this is true. However, in general, there j will be a dependence which we can express as ij = (t ? tli , t ? tm ). For the ensuing analysis, we examine phase-locked solutions where the time since the last spike is ?xed. Our assumption is then valid if |?u (u, v)| |?v (u, v)|. For neurons where the response changes drastically with phase, our analysis will need to be modi?ed. 2.1. Integrate-and-?re model To compare to previous analytical results and to provide a simple concrete example we consider a leaky integrateand-?re model [10–12,30] with dynamics given by dvi = Ii ? vi ? dt δ(t ? tli ) + J Sf (t ? tl ),
j

(2.4)

l

j,l

where the ?ring times tli indicate when vi reaches the threshold value which we have set arbitrarily to 1. Each time the potential reaches threshold, a delta function pulse is subtracted to reset the potential. As an example for a saturating synaptic current, consider a double exponential function Sf (t) = exp(?βt) ? exp(?αt), 0 < t < tf , 0, tf < t, (2.5)

where α > β and tf is the time elapsed to the next synaptic input. For nonsaturating synapses, we can set tf = ∞ [10,11]. Fig. 2 shows numerical simulations of two coupled integrate-and-?re neurons. We see that the integrateand-?re model exhibits ?ring states that are similar to those of the biophysically based neuron model shown in Fig. 1. Eq. (2.4) can be integrated directly [12,26] to obtain the kernels: η(t) = ? exp(?t), t > 0, (2.6)

? 1 1 ? ? [exp(?βt) ? exp(?t)] ? [exp(?αt) ? exp(?t)], 0 < t < tf , ? ? 1?α ? ? 1?β (t) = 1 1 ? [exp(?βtf ) ? exp(?tf )] ? [exp(?αtf ) ? exp(?tf )] ? ? 1?β 1 ? α ? ? ? × exp(?(t ? tf )), tf < t.

(2.7)

We show examples of the kernels for the integrate-and-?re model in Figs. 3(a) and (b). Note that for α > β , (t) > 0 so that the time of the next ?ring is always advanced in response to a positive synaptic input and thus the integrate-and-?re model can be classi?ed as Type I. In Fig. 3(c) we show sketches of possible synaptic kernels for the Type II neurons where the synaptic in?uence can be both positive and negative depending on the time after the spike.

3. Existence of phase-locked solutions We ?rst examine the existence of periodic phase-locked solutions. We derive a condition for which all periodic phase-locked solutions must satisfy. Our condition generalizes that of Gerstner et al. [12] to include heterogeneity

348

C.C. Chow / Physica D 118 (1998) 343–370

Fig. 2. Examples of possible states for two heterogeneous neurons from numerical simulations of the integrate-and-?re model from Section 2.1. The threshold for ?ring is set to 1. The simulations include self-inhibition with strength equal to the coupling strength J . The time axis has been scaled up by a factor of 10. (a) Asynchrony with parameters I1 = 2.5, I2 = 2.4; J = ?0.4; α = 1 β = 0.5. (b) Near-synchrony with I1 = 1.51, I2 = 1.5; J = ?0.4; α = 1 β = 0.5. (c) Harmonic locking with I1 = 2.6, I2 = 2.5; J = ?2; α = 1 β = 0.5. (d) Suppression with I1 = 1.6, I2 = 1.5; J = ?2; α = 1 β = 0.5.

in the applied current and locking at arbitrary phases. Van Vreeswijk et al. [10] derived a locking condition for the speci?c case of two homogeneous integrate-and-?re neurons locked at arbitrary phases. Our strategy is to assume that there is a periodic phase-locked solution and then check for self-consistency, i.e. the membrane potential must be 1 (the threshold value) at the time of the next ?ring. Consider each neuron labeled

C.C. Chow / Physica D 118 (1998) 343–370

349

Fig. 3. Example kernels. (a) An η kernel with a spike. The function used has a spiking part given by η(t) = (exp[50t ] ? exp[?t ])/(1 + 50), t < 0.1, and a recovery part given by η(t) = ? exp[?t + 0.1], t > 0.1. The recovery part has the same form as the integrate-and-?re η kernel of Eq. (2.6). (b) The kernel from Eq. (2.7) for the saturating case (solid line) and the nonsaturating case (dashed line) for parameters α = 1.8, β = 0.6 and tf = 3. (c) Two sketches of possible kernels for Type II neurons. Note that they have both positive and negative parts. These two examples are not derived from any particular model.

by i to be ?ring periodically in the past at tli = (φi ? l)T , where l = 1, 2, . . ., and 0 ≤ φi < 1 is a phase shift. Inserting these ?ring times into Eq. (2.3) gives vi (t) = Ii +
l ≥1

η(t ? (φi ? l)T ) +
j =i,l ≥0

J (t ? (φj ? l )T ),

t ≤ φi T .

(3.1)

This assumes that each neuron will ?re once within a period. We will discuss the situation where this does not occur later. Neuron i will ?re next at t = φi T . Inserting this into Eq. (3.1) and imposing self-consistency yields vi (φi T ) = 1 = Ii +
l ≥1

η(lT ) +
j =i,l ≥0

J (l T + (φi ? φj )T ).

(3.2)

Note that for (φi ? φj ) < 0, the ?rst term of the kernel sum will be equivalently zero by de?nition of the kernels. The interpretation is that within the given period (l = 0), neuron j spikes after neuron i and the synaptic effect does not appear until the next period, l = 1. If Eq. (3.2) is satis?ed then a phase-locked solution exists. The solution will be speci?ed by the set of phases φi and the period T . We note that the phase-locking condition given by Eq. (3.2) consists of N equations for N phases and the period. However, the system is translationally invariant in time, so one of the phases can be arbitrarily ?xed. This then leaves us with a well-de?ned system of N equations and N unknowns (N ? 1 phases and the period T ). For a homogeneous system, the equations are symmetric with respect to the neuron index, implying a large number of possible solutions for the phases. The most notable solutions are the in-phase [6,9–12,16] or synchronized state

350

C.C. Chow / Physica D 118 (1998) 343–370

where all the phases are the same and the splay-phase state [31–34] where the phases are spaced evenly apart over one period. There are also clustered solutions where equal numbers of neurons will ?re at evenly spaced phases. An example is the anti-synchronous [10,35–39] state where half the oscillators ?re a half a period apart. The inclusion of heterogeneity breaks the permutation symmetry. However, for small enough heterogeneity, near synchronous, anti-synchronous and splay-phase solutions that are perturbed from their homogeneous counterparts can exist. 3.1. Phase-locked solutions for two neurons As a speci?c example, we consider the existence of periodic phase-locked solutions for the case of two coupled neurons. We consider neurons with heterogeneous applied current I1 ≥ I2 and homogeneous coupling. We suppose the neurons have ?red at times tli = (φi ? l)T , where l = 1, 2, 3, . . ., i = 1, 2. We choose the origin of time so that φ2 ≥ φ1 , i.e. neuron 2 always ?res after (or with) neuron 1. Applying condition (3.2) for N = 2 yields v1 (φ1 T ) = 1 = I1 +
l ≥1

[η(lT ) + J (lT ? φT )], η(lT ) +
l ≥1 l ≥0

(3.3) (3.4)

v2 (φ2 T ) = 1 = I2 +

J (lT + φT ),

where φ = φ2 ? φ1 ≥ 0 and we have taken into account that neuron 2 ?res after neuron 1. We must solve Eqs. (3.3) and (3.4) for the phase φ and period T . A relation for the phase can be derived by subtracting Eq. (3.4) from (3.3) to give 0 = I1 ? I2 ? G(φ), where G(φ) = J
l ≥1

(3.5)

[ (lT ? T + φT ) ? (lT ? φT )].

(3.6)

For homogeneous neurons, (I1 = I2 ), condition (3.5) for the phase is equivalent to that derived by Van Vreeswijk et al. [10]. For (0) = 0 (which is generally true), we ?nd that G(0) = G(0.5) = G(1) = 0. This implies that for homogeneous neurons φ = 0 and φ = 0.5 are phase-locked solutions which we refer to as the synchronized and anti-synchronized solutions, respectively. In addition, G(φ) is an odd function around G(0.5). Thus, if other phaselocked solutions exist, they must always appear symmetrically in pairs on either side of φ = 0.5. From Eq. (3.5), it is apparent that neurons with heterogeneous Ii cannot lock in a purely synchronous or anti-synchronous state. However, if the heterogeneity is not too large, there can still be phase-locked solutions near to synchrony and near to anti-synchrony. For heterogeneous neurons it is also possible that phase-locked solutions do not exist at all. As a speci?c example, consider the integrate-and-?re model. Using (t) without saturation from (2.7) in condition (3.6), we ?nd G(φ) = eφT ?T ? e?φT eαφT ?αT ? e?αφT J ? 1?α 1 ? e?αT 1 ? e?T ? J 1?β eβφT ?βT ? e?βφT eφT ?T ? e?φT ? 1 ? e?βT 1 ? e?T . (3.7)

Fig. 4 shows some examples of G(φ). The parameters T , α and β control the shape of G(φ). Generally, for moderate rise rate, α , and slow decay rate, β , as compared to the period T , there are only the synchronous and anti-synchronous

C.C. Chow / Physica D 118 (1998) 343–370

351

Fig. 4. Two examples of G(φ) for inhibitory coupling. The top ?gure has parameters: J = ?1, T = 10, α = 0.1, β = 0.05. The bottom ?gure has the same parameters except α = 0.5. The intersection of G(φ) and the horizontal line given by I1 ? I2 gives a phase-locked solution. If |I1 ? I2 | is large enough then no phase-locked solution exists.

solutions as seen in Fig. 4(a). However, if α is increased, a third phase-locked solution can arise as seen in Fig. 4(b). With heterogeneity, we see that the phase-locked solutions move away from their homogeneous counterparts. To specify the period for a given phase, we add Eq. (3.4) to (3.3) and divide by 2 to obtain ?+ 1=I
l ≥1

η(lT ) +

J [ (lT ? φT ) + (lT ? T + φT )] , 2

(3.8)

352

C.C. Chow / Physica D 118 (1998) 343–370

? = (I1 + I2 )/2. Using Eq. (3.8) and assuming the rise time of the synaptic current is very short compared to where I the decay time (α β ), we ?nd that the period for the synchronized state (φ = 0) for the integrate-and-?re model with saturating synapses satis?es [28] ?(1 ? e?T ) + 1=I J (e?βT ? e?T ) . (1 ? β) (3.9)

The period for nonsaturating synapses is given by Eq. (3.9) with an extra factor of (1 ? e?βT ) in the denominator of the second term. The period does not change much for near synchrony where the phase is near zero. It has been found that (3.9) can be used to estimate the period of a synchronized inhibitory network of biophysically based Type I neuron models [28]. We should note that a periodic phase-locked state exists only if the phase condition (3.5) and the period condition (3.8) can be solved simultaneously.

4. Stability of phase-locking We examine the local stability of phase-locked solutions for a network of N all-to-all coupled neurons. Our formalism generalizes that of Ref. [12] to include heterogeneity in the applied current Ii and for locking at arbitrary phases φi . In this strategy, perturbations around the ?ring times of the phase-locked solutions are checked for growth or decay. Based on the perturbation dynamics we obtain suf?cient conditions for stability. The resulting condition is similar to that derived for phase models [40]. For neuron i we consider ?ring in the past at times tli = (φi ? l)T + δi (n ? l), l = 1, 2, 3, . . . , (4.1)

where δi (n ? l) is a small perturbation around each ?ring time and n is the index for the current perturbation. Since the unperturbed neurons ?re periodically, we can translate n to 0 in the periodic part of tli . We want to examine whether δi (n) grows or decays as n → ∞. We will show that δi (n) depends on the perturbations of the previous ?ring times of all the neurons in the network. To derive a mapping for δi (n), we insert the perturbed ?ring times (4.1) into Eq. (2.3) to obtain vi (t) = Ii +
l ≥1

η(t ? (φi ? l)T ? δi (n ? l)) +
j =i,l ≥0

J (t ? (φj ? l )T ? δj (n ? l )).

(4.2)

The neuron will next ?re at t = φi T + δi (n). Hence we can set vi (φi T + δi (n)) = 1, which implies after linearizing in δi (n) that vi (φi T ) 1?v ˙i (φi T )δi (n), (4.3)

where the dot denotes the ?rst derivative with respect to time. We now set t = φi T in Eq. (4.2). Using Eq. (4.3) and linearizing in δk (l) we obtain 1?v ˙i (φi T )δi (n) = Ii +
l ≥1

[η(lT ) ? η(lT ˙ )δi (n ? l)] J [ (l T + (φi ? φj )T ) ? ˙ (l T + (φi ? φj )T )δj (n ? l )]. (4.4)

+
j =i,l ≥0

C.C. Chow / Physica D 118 (1998) 343–370

353

For an unperturbed spiking history Eq. (3.2) holds. This allows us to simplify Eq. (4.4). Solving for δi (n) yields ? ? 1 ? δi (n) = η(lT ˙ )δi (n ? l) + J ˙ (l T + (φi ? φj )T )δj (n ? l )? , (4.5) v ˙i (φi T )
l ≥1 j =i,l ≥0

where v ˙i (φi T ) =
l ≥1

η(lT ˙ )+
j =i,l ≥0

J ˙ (l T + (φi ? φj )T ).

(4.6)

Stability is determined from Eq. (4.5). If the perturbations, δi (n), decay to zero as n approaches in?nity then the system is considered to be stable. The outcome is not immediately apparent since a given perturbation of a neuron depends on previous perturbations including those within the same cycle (i.e have the same index n). One point of concern might be that if perturbations decay in the future then should they not amplify in the past and violate the smallness assumption required for the stability analysis. We show that this does not pose a problem and the perturbed dynamics can be treated as an initial value problem. Although a single perturbation depends on the entire history of previous perturbations, we can choose as an initial condition a state of perfectly phase-locked periodic ?ring (i.e. zero perturbations) in?nitely into the past followed by a single perturbation at one ?ring time. All the ensuing perturbations will then respond according to Eq. (4.5). If the perturbations eventually decay away then the system is stable. By turning the perturbed dynamics into a high-dimensional ?rst return map we obtain a suf?cient condition for stability which we state in the following Theorem 1. Theorem 1. If the series (4.5) can be truncated then a suf?cient condition for stability of a phase-locked state is that all the coef?cients the series in (4.5) are positive. Proof. We ?rst make δi (n) depend only on previous ?rings by substituting for δj (n) recursively on the right-hand side of the mapping (4.5) whenever φi > φj . This substitution can be done systematically if we index the neurons so that they are ordered by phase with φ1 < φ2 < · · · < φN , i.e. neuron 1 ?res ?rst, neuron 2 second, and so forth. With this phase ordering the mapping for δ1 (n) is given by Eq. (4.5) with i = 1 and l ≥ 1. Since neuron 1 ?res ?rst it only depends on the previous ?rings of the other neurons. To obtain the mapping for δ2 (n) we must substitute for δ1 (n) on the right-hand side of mapping (4.5) with i = 2. We continue this recursive substitution scheme for all N neurons where for i = k we must substitute for δj (n) on the right-hand side of Eq. (4.5) for j < k . We write the resulting mapping as
N ∞

δi (n) =
j =1 l =1

Mij (l)δj (n ? l),

(4.7)

where Mij (l) is obtained from applying the prescription described above and the sum over j now includes i . The elements of Mij are composed of products and sums of the coef?cients in the series (4.5). Thus, if the coef?cients are positive then the elements of Mij are positive. Theorem 1 then directly follows from Lemma 1. Lemma 1. If the series (4.5) can be truncated then a suf?cient condition for stability is given by Mij (l) > 0. Proof. By assumption the kernels η(t) and (t) vanish suf?ciently quickly so that we can truncate the series in Eq. (4.5). This allows us to truncate the inner sum in Eq. (4.7) at some l = lM . We can then express the mapping (4.7) in terms of a ?rst return mapping by embedding in a higher dimension with δ (n) = M · δ (n ? 1), (4.8)

354

C.C. Chow / Physica D 118 (1998) 343–370

where in block form, ? δ(n) ? δ(n ? 1) ? ? δ (n) = ? δ(n ? 2) ? . . ? . δ(n ? lM + 1)

? ? ? ? ?, ? ?

M(1) ? I ? ? M=? 0 ? . ? . . 0

?

? M(2) · · · M(lM ? 1) M(lM ) 0 ··· 0 0 ? ? I ··· 0 0 ? ?, ? . . . . . . . . . ? 0 ··· I 0

(4.9)

I denotes the N × N identity matrix, δ(k) is the column vector of length N with elements δi (k) and M(l) is the N × N matrix with elements Mij (l) from Eq. (4.7). The dynamics of the perturbations are determined by repeated matrix multiplications of M. The phase-locked solution is stable if all the eigenvalues of M have absolute value less than unity. For simple enough systems we can calculate the eigenvalues explicitly (see Appendix B). This is usually not the case, so we need to determine bounds on the modulus of the maximum eigenvalue ρ(M) (spectral radius). The spectral radius is less than or equal to any given matrix norm [41]. In particular ρ(M) ≤ M ∞ = max j |Mij |. From Eqs. (4.5)–(4.7) we ?nd that the rows of l M(l) sum to unity, which implies that all the rows M also sum to unity. Thus we immediately ?nd that 1 is an eigenvalue with eigenvector δ = (1, 1, 1, . . . , 1). This corresponds to a global uniform shift in time and is due to the time translational invariance of the dynamics. It can be disregarded for stability considerations [12]. If the elements of M are non-negative then ρ(M) = M ∞ = 1. The spectral radius is the unique maximum modulus eigenvalue if M is a primitive matrix which is satis?ed if Mm > 0 for some m [41]. By inspection of M from (4.9), we see that Mm > 0 if Mij (l) > 0 proving Lemma 1. It is important to note that Mij (l) > 0 is not a necessary condition for stability. We make this explicit in the following corollary. Corollary 1. Stability is possible even if Mij has negative elements. Proof. If M is primitive, we know that one eigenvalue is always 1, it is nondegenerate and all the other eigenvalues are within the unit disc [41]. Stability is lost if an interior eigenvalue crosses the unit boundary, which can occur if we allow some of the elements of Mij to become negative. Because the eigenvalues are continuous functions of the elements, we can make some of the elements negative and still keep the eigenvalues within the unit disc by the intermediate value theorem. Remark 2. Our perturbation scheme generalizes that of Gerstner et al. [12] who arrived at the stability equation (4.8) for the purely synchronized solution (all phases equal). They considered dynamics where ˙ (0) = 0, so only previous perturbations affect the current one. For a large homogeneous network, they proved the necessary and suf?cient condition for stability is J l ˙ (lT ) > 0. In our case with heterogeneity, pure synchrony is not possible. We are thus forced to consider locking at arbitrary phases and use the recursive substitution scheme above to calculate the mapping for the perturbations. 4.1. Stability for two coupled neurons As an example we apply the formalism for N = 2. Eq. (4.5) is δ1 (n) = 1 v ˙1 [η(lT ˙ )δ1 (n ? l) + J ˙ (lT ? φT )δ2 (n ? l)],
l ≥1

(4.10)

C.C. Chow / Physica D 118 (1998) 343–370

355

δ2 (n) =

1 v ˙2

η(lT ˙ )δ2 (n ? l) +
l ≥1

J v ˙2

˙ (l T + φT )δ1 (n ? l ),
l ≥0

(4.11)

where φ = φ2 ? φ1 . Since neuron 2 ?res after neuron 1, we substitute δ1 (n) from Eq. (4.10) into Eq. (4.11) to yield δ2 (n) =
l ≥1

η(lT ˙ ) J 2 ˙ (φT ) ˙ (lT ? φT ) + δ2 (n ? l) v ˙2 v ˙1 v ˙2 ˙ ) J ˙ (lT + φT ) J ˙ (φT )η(lT + δ1 (n ? l) . v ˙2 v ˙1 v ˙2 ? J ˙ (lT ? φT ) ? v ˙1 ?, η(lT ˙ ) J 2 ˙ (φT ) ˙ (lT ? φT ) ? + v ˙2 v ˙1 v ˙2 (4.12)

+

Combining (4.10) with (4.12) we obtain ? η(lT ˙ ) ? v ˙1 M(l) = ? ? J ˙ (lT + φT ) J ˙ (φT )η(lT ˙ ) + v ˙2 v ˙1 v ˙2

(4.13)

which is to be used in Eqs. (4.8) and (4.9). The suf?cient condition for stability of a phase-locked solution is given by Lemma 1 which we denote by M(l) > 0. For the case of two neurons we can also obtain the necessary condition for stability by examining whether or not one neuron can remain locked to the other one independently of whether they are ?ring periodically. This type of stability was considered by van Vreeswijk et al. [10]. We derive their condition using the spike response formalism which we state in the following theorem. Theorem 2. A necessary condition for stability of a periodic phase-locked solution for two neurons is G (φ) > 0, where G(φ) is given by Eq. (3.6) and prime denotes the ?rst derivative with respect to φ . Proof. Suppose that neuron 1 ?res at times tl1 = t1?l and neuron 2 ?res at times tl2 = t1?l + θ + δ(n ? l), where l ≥ 0. We consider the situation where neuron 2 has ?red after neuron 1. Neuron 1 will satisfy the condition v1 (t1 ) = 1 while neuron 2 satis?es the condition v2 (t1 + θ + δ(n)) = 1, from which we derive the linearization ˙2 (t1 + θ)δ(n). Inserting this into Eq. (2.3) and linearizing in δ(n ? l) yields v2 (t1 + θ ) = 1 ? v 1 = I1 +
l ≥1

[η(t1 ? t1?l ) + J (t1 ? t1?l ? θ ) ? J ˙ (t1 ? t1?l ? θ )δ(n ? l)], [η(t1 ? t1?l ) ? η(t ˙ 1 ? t1?l )δ(n ? l) + J (t1 ? t2?l + θ )].
l ≥1

(4.14) (4.15)

1?v ˙2 (t1 + θ)δ(n) = I2 +

Subtracting Eq. (4.15) from Eq. (4.14) and canceling terms, we obtain the condition for synchrony: I 1 ? I2 = J
l ≥1

[ (t1 ? t2?l + θ) ? (t1 ? t1?l ? θ )],

(4.16)

and the perturbation mapping: δ(n) = where v ˙2 (t1 + θ) =
l ≥1

1 v ˙2 (t1 + θ)

[η(t ˙ 1 ? t1?l ) ? J ˙ (t1 ? t1?l ? θ )]δ(n ? l),
l ≥1

(4.17)

[η(t ˙ 1 ? t1?l ) + J ˙ (t1 ? t2?l + θ )].

(4.18)

356

C.C. Chow / Physica D 118 (1998) 343–370

Eq. (4.17) can be rewritten as δ(n) =
l ≥1

A(l)δ(n ? l),

(4.19)

where A(l) = 1 [η(t ˙ 1 ? t1?l ) ? J ˙ (t1 ? t1?l ? θ )]. v ˙2 A(l) < 1. (4.20)

It will be shown in Lemma 2 that a necessary condition for stability is Applying this condition leads to ˙ 1 l ≥1 η(t which implies J
l ≥1

? t1?l ) ? J ˙ (t1 ? t1?l ? θ)

˙ 1 ? t1?l ) + J ˙ (t1 ? t2?l + θ) l ≥1 η(t

< 1,

(4.21)

[ ˙ (t1 ? t1?l ? θ) + ˙ (t1 ? t2?l + θ)] > 0.

(4.22)

Condition (4.22) is the general necessary condition for stable locking. For the speci?c case of periodic ?ring we can set tl = (1 ? l)T and θ = φT in (4.22) which is equivalent to G (φ)/T > 0. This leads to G (φ) > 0 as stated in Theorem 2. Lemma 2. For perturbed dynamics of the form
lM

δ(n) =
l =1

A(l)δ(n ? l),

(4.23) A(l) ≤ 1. If A(l) > 0, for all l , then the condition is suf?cient as well.

a necessary condition for stability is that

Proof. The dynamics can be expressed in matrix notation as in Eq. (4.8) where M(l) is replaced by the scalar A(l) and the identity matrices are replaced by 1. The corresponding matrix M is a companion matrix and the characteristic polynomial is given by [12,41]
lM

1=
l =1

A(l)λ?l .

(4.24)

Let f (λ) = l A(l)λ?l , so that eigenvalues of M satisfy f (λ) = 1. Suppose l A(l) > 1. Then f (1) > 1. However, f (λ) → 0 as λ → ∞. Since f is a polynomial and hence continuous, f (λ) = 1 at some λ > 1 by the intermediate value theorem meaning at least one eigenvalue will be greater than 1. Thus l A(l) ≤ 1 is necessary for all eigenvalues to be less than 1 and hence for stability. If A(l) > 0, then l A(l) < 1 is necessary and suf?cient for stability since f (λ) < 1, for all λ ≥ 1 implying there are no eigenvalues greater than or equal to 1. Remark 3. The condition G (φ) > 0 is not a suf?cient condition for periodic phase-locking. However, if the dynamics are constrained so that A(l) as given in Eq. (4.20) is positive for all l , then l ≥1 A(l) < 1, which can be re-expressed as condition (4.22), is the necessary and suf?cient condition for arbitrary phase-locking. However, this constraint does not ensure periodicity. For instance, it may be possible that two neurons could be locked together in a higher order periodic or even an aperiodic orbit.

C.C. Chow / Physica D 118 (1998) 343–370

357

5. Two neuron network We now apply our formalism to examine the existence and stability of speci?c phase-locked solutions for two coupled neurons. We will revisit the case of two homogeneous neurons and then consider neurons with heterogeneous applied current but with homogeneous inhibitory and excitatory coupling. In Appendix B we consider a simple example for two neurons where only the kernels from the last ?ring are kept. We compute the eigenvalues for stability explicitly and compare with our results for the general case. Applying the stability conditions for phase-locked solutions, we make the following proposition for the stability of two coupled neurons. Although we treat the spiking period T as a parameter in the proposition, it is important to note that it is an emergent timescale that depends on parameters of the network. Proposition 1. A periodic phase-locked solution for two coupled neurons is stable if in (4.13) (i) η(lT ˙ ) > 0 for l ≥ 1, (ii) J ˙ (φT ) ≥ 0 and (iii) J ˙ (lT ± φT ) > 0, for all l ≥ 1. These are suf?cient but not necessary conditions. ˙2 > 0. Proposition 1 is established by analyzing the matrix elements of M(l) in (4.13). By de?nition v ˙1 > 0 and v Conditions (i)–(iii) then imply that M(l) > 0 proving the proposition by Lemma 1. Condition (i) is generally true for most neuron models. Condition (ii) is strongly dependent on the neuron type. Condition (iii) is contingent on the neuron type as well as the neuron parameters, i.e. strength of the applied current, the synaptic strength and the synaptic decay time. Gerstner et al. [12] proved Proposition 1 for the homogeneous synchronous case (φ = 0) for N coupled neurons where N is large. We discuss this case in Section 6.1. Condition (iii) corresponds to a special case of their locking theorem. This puts a condition on how slow the onset of the synapse can be. As will be discussed below, condition (ii) was implicitly assumed in their theorem. Van Vreeswijk et al. [10] emphasized condition (ii) in their criterion for stability. 5.1. Homogeneous neurons The existence and stability of phase-locking for two homogeneous neurons has been analyzed in detail by van Vreeswijk et al. [10] and by Hansel et al. [11] for excitatory connections. These analytical results considered integrate-and-?re dynamics or used weak coupling phase reductions. Our formalism is applicable to a more general class of neuron models. In Section 5.2, we consider heterogeneous neurons. Gerstner et al. [12] considered the stability of the synchronous solution for a large homogeneous network. We will relate the stability conditions of Ref. [12] to those of Ref. [10]. Periodic phase-locked solutions are given by Eqs. (3.5) and (3.8) for I1 = I2 . As discussed in Section 3.1, these include the synchronous (φ = 0) and anti-synchronous (φ = 0.5) solutions as well as a possible third locked solution. 5.1.1. Inhibitory coupling Here we relate the results and formalism of van Vreeswijk et al. [10] to our formalism. First consider synchronous locking. As mentioned, condition (i) in Proposition 1 is satis?ed for most neuron models. For Type I models, (t) is always positive, so condition (ii) can only be satis?ed if ˙ (0) = 0. This condition is equivalent to a nonzero rise time of the synaptic current which was emphasized by Van Vreeswijk et al. [10] as the condition required for synchrony with inhibition. In order for synchrony to occur with inhibition, the effect of the ?rst neuron on the second one cannot be felt immediately. They gave an empirical criterion that the synaptic rise time be slower than the width of the spike. For certain Type II models where (t) is initially negative and then becomes positive J ˙ (0) could be greater than zero (see Fig. 3). Condition (iii) is easily satis?ed for Type I neurons (which includes the integrate-and-?re model). Recall that inhibitory coupling is given by J < 0. Thus, if (t) is positive with a maximum at t = tc and is

358

C.C. Chow / Physica D 118 (1998) 343–370

monotone decreasing for t > tc (which applies to Type I neurons) then condition iii) in Proposition 1 is equivalent to T ? φT > tc . The decay rate of (t) is controlled to some extent by the decay rate of the synaptic current as well as the intrinsic membrane decay rate (see Eq. (2.7)). The period on the other hand is an emergent property of the network which depends on the applied current, synaptic strength, and synaptic timescales. For Type I neurons the period can be estimated from Eq. (3.9) [28]. For certain Type II neurons it may be possible that condition iii) can never be satis?ed, so stable synchronous solutions would not be possible with inhibition. Proposition 1 gives suf?cient conditions for stability. By Corollary 1, synchrony may be possible even if Proposition 1 is violated. To ?nd when locking is not possible we examine the necessary condition given by Theorem 2. From Eq. (3.6) we note that the necessary condition for stability G (0) ≡ J T ˙ (0) + 2J T
l ≥1

˙ (lT ) > 0

(5.1)

is satis?ed for Proposition 1 as expected. If ˙ (0) > ?2 l ≥1 ˙ (lT ) then synchrony is unstable. Satisfying this condition depends on ˙ (0) which is governed by the rise time of the synaptic current. For stability, the neuron must be able to ?re without being immediately and strongly inhibited by the synapse. This relates to the van Vreeswijk et al. [10] statement that the rise of the synapse must be slower than the generation of a spike for stability to be possible. For the integrate-and-?re model with nonsaturating synapses, a zero rise time guarantees instability. This can be seen by examining Eq. (3.7) where it can be shown that G (0) < 0 if α → ∞ and T > 0. This veri?es previous results of local instability of the synchronous state for zero rise time with inhibition in the integrate-and-?re model [10,12]. By Proposition 1, the anti-synchronous solution is stable as long as the synapse is fast enough so that ˙ (0.5T + lT ) < 0, for l ≥ 0. From (3.6) we ?nd that the necessary condition for stable anti-synchrony is given by G (0.5) ≡ 2J T
l ≥1

˙ (lT ? 0.5T ) > 0.

(5.2)

Condition (5.2) shows that G (0.5T ) could be negative if the onset of inhibition is very slow. For the integrate-and?re model this requires that α and β be small enough in Eq. (2.7) (see Fig. 4). This is in accordance with the results of van Vreeswijk et al. [10]. The synchronized and anti-synchronized solutions can both be stable simultaneously if the rise time is nonzero and (t) is monotone decreasing for t ≥ 0.5T . If this is the case then additional phase-locked solutions must also exist. Since G(φ) is continuous, if G (φ) is positive at 0, 0.5 and 1 then it must cross zero at least once between 0 and 0.5 and between 0.5 and 1. If it crosses only once between 0 and 0.5 and between 0.5 and 1 then these solutions must be unstable (see Fig. 4). Since G(φ) is odd about φ = 0.5, we can think of these two new solutions as one. This third solution was obtained numerically in the integrate-and-?re model with nonsaturating synapses in [10,11] where it was found to arise from a pitchfork bifurcation of the anti-synchronous solution. It approaches the synchronous solution as the synapse becomes faster compared to the period (α and β get larger). The location of the third solution with respect to the synchronous and anti-synchronous demarcates the basin of attractions for the two stable solutions. All phase-locked solutions can be unstable if the onset of the synapse is so slow that (t) is still increasing at t = T because M is no longer a positive matrix. For slow enough synapses, the necessary conditions for stability given by (5.1) and (5.2) can both be violated. 5.1.2. Excitatory coupling The case of slow excitatory coupling for two neurons has been analyzed previously by van Vreeswijk et al. [10] and Hansel et al. [11]. They both show that for Type I neurons, synchrony is generally unstable and anti-synchrony is stable for slow synapses. Our analysis agrees with their results. For two neurons, stability is not ensured since

C.C. Chow / Physica D 118 (1998) 343–370

359

˙ (lT ) < 0, violating Proposition 1. From Eq. (5.1) we ?nd that the necessary condition for synchrony in Theorem 2 is violated for excitatory coupling unless ˙ (0) is positive and large. An examination of G (0.5) for J > 0 in (5.2) shows that anti-synchrony is stable if the rise time of the synaptic current is very slow. However, stability might be possible for Type II neurons with excitation as pointed out by Hansel et al. [11]. This can be seen in our formalism because ˙ (lT ) could be positive for l ≥ 1 if the phase response curve had both positive and negative parts. We show an example in Fig. 3. For the Hodgkin–Huxley equation, Kistler et al. [27] have shown that the response of the membrane to a synaptic input can behave in a way to support synchrony. For the integrate-and-?re model with instantaneous (delta function) synaptic coupling, stable synchrony is possible. To see this consider the integrate-and-?re model (2.4) with a decay constant γ for the membrane potential restored. This implies that η(t) = ? exp(?γ t) and (t) = exp(?γ t). Thus ˙ (t) = ?γ exp(?γ t), which implies M12 in (4.13) is negative. However, these negative elements can be made arbitrarily small if γ and J are made small and stability is possible by Corollary 1. For a combination of small coupling and weak decay, the synchronous state can be stable. This is in accordance with the result of Peskin [35] who showed that the product J γ must be small for synchrony in the integrate-and-?re model. For the integrate-and-?re model, the third phase-locked solution which comes from a pitchfork bifurcation of the anti-synchronous solution is always stable [10,11]. As the synapse becomes faster and more pulse-like, the third solution approaches the synchronous solution. Thus for very brief synapses, two neurons locked in the third solution can be considered to be in a state of almost synchrony [42]. Remark 4. It is important to note that our local stability is different from the global stability of Mirollo and Strogatz [9] and related work [8,13] for pulse-coupled oscillators. In these papers, an absorption condition is assumed for synchrony. By this it is meant that if two oscillators ?re together within a certain time window, then they are considered synchronized. These papers show that for pulsatile excitation and almost all initial conditions, the oscillators will eventually get close enough to be absorbed. What we have examined here, which is in line with Refs. [10–12,35], is local stability of coupled oscillators that still in?uence each other even if completely synchronized. After the oscillators have been absorbed they may not synchronize in the local sense. They could be in a state of almost synchrony, make small oscillations around synchrony, or even be repelled away [12,42].

5.2. Heterogeneous neurons We now apply our formalism to two coupled heterogeneous neurons. In this case, phase-locked solutions need not exist. For large enough heterogeneity or weak enough synaptic strength, condition (3.5) may not be satis?ed (see Fig. 4). When there are no phase-locked solutions, several possibilities can arise. If no stable phase-locking exists and both neurons are able to ?re within a period then a possible outcome is a state of asynchrony, where each neuron is effectively decoupled and ?res independently of the other neuron. For strong inhibitory coupling, we can get a state of harmonic locking where the neurons are locked at some ratio of their periods. The locking period can be very long. If the inhibition is so strong so that the slower neuron never ?res then we call this state suppression [21]. See Fig.1 for numerical examples of these states. We will discuss harmonic locking and suppression in more detail in the following sections. Under certain conditions, stable phase-locked solutions can exist. For weakly heterogeneous neurons these solutions will be near to the homogeneous locked solutions. We concentrate on the near synchronized and near anti-synchronized solutions. From (3.5) we see that for a given level of heterogeneity, |J | must be large enough for a phase-locked solution to exist (see Fig. 4). The phase φ is determined by the heterogeneity and the coupling strength. For small heterogeneity and small phase, φ varies linearly with I1 ? I2 and inversely with J .

360

C.C. Chow / Physica D 118 (1998) 343–370

For excitatory coupling, the stability results for homogeneous neurons can be applied to the heterogeneous case. Near synchrony is generally unstable for Type I neurons. For integrate-and-?re neurons, by analyzing M(l) in (4.13), we see that stability may be possible in the case of weak coupling, weak decay of η(t), and extremely brief synapses as in the homogeneous case. The conditions for stability of the near anti-synchronous solution is possible if the heterogeneity is mild and the synapse is slow. Near synchrony can be stable for heterogeneous Type II neurons by Proposition 1. For inhibitory coupling (J < 0), stability of near anti-synchrony can be obtained directly from Proposition 1 but with φ shifted slightly from 0.5. As long as conditions (ii) and (iii) are still satis?ed then anti-synchrony is stable. Stability of the near synchronous case where φ is small but nonzero will depend on the neuron properties. If J ˙ (φT ) ≥ 0, then stability is ensured by Proposition 1. For Type I neurons this will not be generally true but ˙ (φT ) = 0 may hold as a result of delays in the onset of the synapse or from refractoriness of the neuron. (We note that J ˙ (φT ) < 0 holds for the integrate-and-?re model (2.4).) For J ˙ (φT ) < 0, stability is not guaranteed because elements of the second row of M(l) could become negative. However, by Corollary 1, negative elements do not immediately imply instability because the eigenvalues can remain within the unit circle even if M(l) takes on negative values. This leads to the following proposition on the stability of inhibitory heterogeneous neurons. Proposition 2. Stability of near synchrony for Type I inhibitory coupling is possible if the heterogeneity and strength of the applied currents are both weak enough, and the synaptic strength is strong enough. We show this by examining the elements of M(l) in (4.13) individually and then determine conditions for them to be positive which would ensure stability. As before, we assume that η(lT ˙ ) > 0, for all l . We also assume that (t) has a single maximum at t = tc and that T ? φT > tc so that ˙ (lT ± φT ) < 0, for all l . This implies that M11 (l) > 0 and M12 (l) > 0. Heterogeneous phase-locking implies that φ > 0 which affects elements M21 (l) and M22 (l). First consider M21 (l) = J ˙ (lT + φT ) J ˙ (φT )η(lT ˙ ) + . v ˙2 v ˙1 v ˙2 (5.3)

For J < 0, this element is positive if ˙ (φT ) > 0 is small enough and ˙ (lT + φT ) < 0 is negative enough and decays away as slow or slower than η(lT ˙ ). If we assume that ˙ (0) = 0 then ˙ (φT ) small implies that φT must be small. But if T is too small then ˙ (lT + φT ) can become close to zero or even positive leading to a loss of stability. Thus for M21 (l) to be positive we want the period T to be long and the phase φ to be small. The phase φ can be made smaller if the level of heterogeneity is reduced or if the synaptic strength |J | is increased. From Eq. (3.9) we see that T is made longer if the synaptic strength is increased or if the applied current is reduced. Thus in order for M21 (l) to be positive we require weak enough heterogeneity and strength of the applied currents and a strong enough synapse as stated in Proposition 2. We must also make sure that these conclusions are consistent with element M22 (l) = η(lT ˙ ) J 2 ˙ (φT ) ˙ (lT ? φT ) + . v ˙2 v ˙1 v ˙2 (5.4)

This element is positive if ˙ (φT ) > 0 is small and the decay of ˙ (lT ? φT ) < 0 is as fast or faster than η(lT ˙ ). Thus ˙ ) must have requiring M21 (l) > 0 and M22 (l) > 0 to both hold simultaneously implies that ˙ (lT ± φT ) and η(lT the same decay rate for l ≥ 1. Otherwise there will be negative elements in the second row of M(l) for some l . For saturating synapses and fast rise times, this is satis?ed (see Eq. (2.7)). However, the decay rates will in general be different for nonsaturating synapses. Even if this is the case, as long as the decay rates are fairly close or fairly fast, the negative elements can be made arbitrarily small and stability is still possible by Corollary 1. This result indicates

C.C. Chow / Physica D 118 (1998) 343–370

361

that synchronization should be easier to achieve with saturating synapses as opposed to nonsaturating synapses for heterogeneous neurons. Remark 5. Proposition 2 con?rms the numerical observations in [21] where it was found that synchrony was only observed in what we called the phasic regime. It was shown in [28] that the conditions of Proposition 2 lead to the phasic regime where βT ? 1. We can use this to estimate the boundary between synchrony and asynchrony in the J –β parameter plane. From our above analysis and previous numerical simulations [21] we ?nd that βT is approximately constant in the synchronous regime. At the boundary, for the level of heterogeneities chosen (5% in intrinsic frequencies), we found that βT ? 1 de?ned the boundary. Recall that the period of a synchronized network with saturating synapses approximately satis?es the transcendental relation (3.9) [28], which we rewrite as J =? [I (1 ? e?τ ) ? 1](1 ? τ ?1 ) , (e?1 ? e?τ ) (5.5)

where we have substituted T ? τ ≡ β ?1 . Eq. (5.5) gives the synchrony–asynchrony boundary for g as a function of β . We can obtain a simpler expression if we let that τ ? 1 + δ and expand for small δ : J ? ?[I (e ? 1) ? e + I δ ]. Thus we see that J depends linearly on δ ? τ ? 1. This was observed in numerical simulations [21]. 5.3. Suppressed states In the previous analysis for phase-locking, it was assumed that there was a phase-locked solution with period T and phase φ . However, the locking condition (3.2) may not have a solution. For inhibitory coupling, it may be that some neurons never ?re at all. We called such a state suppression in [21]. We illustrate this for two neurons. Suppose I1 > I2 and neuron 1 ?res before neuron 2. If neuron 2 never ?res (i.e. limt →∞ v2 (t) < 1) then we say neuron 2 is suppressed. As we will discuss in Section 6.3, suppression can also occur in networks of N neurons. Here we derive the conditions for suppression for two inhibitory neurons. We assume that neuron 1 has been ?ring periodically in the past at t = (n ? l)T , where l ≥ 1 and nT denotes the current time. We also assume that neuron 2 never ?res. We let J = ?g , where g > 0. Using Eqs. (2.3) and (3.2), this can be expressed as v1 (nT ) = 1 = I1 +
l ≥1 n n→∞

(5.6)

η(lT ), (lT ) < 1.

(5.7) (5.8)

lim v2 (nT ) = I2 ? g
l =0

Eqs. (5.7) and (5.8) give the conditions for suppression. Using the kernels (2.6) and (2.7) for integrate-and-?re neurons with saturating synapses and using the condition α β we can compute the sum explicitly in Eq. (5.8) and obtain the condition I2 ? g e?βT ? e?T < 1. 1 ? β (1 ? e?T ) (5.9)

We see that suppression can be achieved with large enough g or small enough T . Without self-inhibition, from Eq. (5.7) the period is approximately given by T = ln I1 I1 ? 1 (5.10)

362

C.C. Chow / Physica D 118 (1998) 343–370

If we include the effects of self-inhibition (i.e. the neuron inhibits itself with strength g when it ?res), the period T ? = I1 . We can then use Eq. (3.9) to simplify Eq. (5.9) and obtain is given by Eq. (3.9) with I I 2 ? I1 + 1 <1 1 ? e?T (5.11)

From this expression we see that if the neurons are homogeneous (I1 = I2 ) then suppression is not possible. For heterogeneous cells, if I2 < I1 then suppression will de?nitely occur for large enough T . The longer the period the more likely suppression will occur as opposed to the case without self-inhibition. We can estimate the transition boundary to suppression in g ? β space. Without self-inhibition, the period of neuron 1 is not affected by the synapse and we can assume that the boundary to suppression is given by condition (5.9) with an equality I2 ? g e?βT ? e?T ? 1, 1 ? β (1 ? e?T ) (5.12)

where the period is given by (5.10). With self-inhibition, the period now also depends on the synaptic parameters. From (5.11) we see that at a critical period suppression will occur. The boundary will then be given by the conditions specifying that period. Consider the case where T 1 but βT ? 1. This corresponds to what we have called the phasic regime [21]. This allows us to ignore factors of e?T in the period condition (3.9). The suppression boundary is then given by I1 ? g e?βT ? 1, 1?β (5.13)

for ?xed T such that condition (5.11) is minimally satis?ed. Solving for g we then obtain g ? (I1 ? 1)(1 ? β)eβT . We see that for β < 1 and T large enough, g increases as β increases as was observed numerically [21]. 5.4. Harmonic locking It is also possible that two neurons may be harmonically locked with a period ratio p1 : p2 . The ratio need not be rational. Suppression can be considered to be harmonic locking at the ratio 1 : ∞. Consider two neurons locked in the ratio of p1 : p2 so that neuron i will ?re pi times before repeating its ?ring pattern. Let the length of this periodic cycle be T , which is the same for both neurons. Fig. 1 and 2 show examples of harmonically locked states. Consider two inhibitory neurons ?ring in the past at tli = ? l T ? φlimod pi T , pi (5.15) (5.14)

where [·] indicates least integer, l = 0, 1, 2, . . ., and φlimod pi is the phase advance per ?ring which repeats after i , where m takes the values 0, 1, . . . , p ? 1. If we pi ?rings. In the next cycle, neuron i will ?re at t = T ? φm i i i substitute this into Eq. (2.3) we obtain
i ) = 1 = Ii + vi (T ? φm i

η T +
l ≥mi

l i )T T + (φlimod pi ? φm i pi l j i )T T + (φl mod pj ? φm i pj , (5.16)

?
l ≥n

g

T +

C.C. Chow / Physica D 118 (1998) 343–370
j

363

i > 0 is minimal (i.e. n gives the index of the nearest ?ring of the other where n is the index for which φn mod pj ? φm i neuron). Eq. (5.16) is a system of M ≡ p1 + p2 equations that must be solved simultaneously for the period T and i . Because of the freedom to shift globally in time we can ?x one of the phases (e.g. set φ 1 = 0), the phases φm 0 i which then leaves M equations for the remaining M ? 1 phases and the period T . For a ?xed set of parameters, there is an in?nite set of possible harmonic states with different period ratios. Stability of the harmonic state can be established by applying results of the previous section. A stability analysis can be performed as before by perturbing the time of each ?ring and constructing a mapping for the perturbations. The mapping will map all the perturbations from one cycle at t = kT to the next cycle at t = (k + 1)T , where k is the least integer of l/pi . The suf?cient condition for stability will be given by Theorem 1. We can examine the local stability of the relative ?ring times of the two neurons. As shown in Section 4.1, for inhibitory coupling, a suf?cient condition for stable phase-locking is that the phases between two consecutive spikes be either very small j i )T ) will arise in the perturbation series. The or far apart. More speci?cally terms such as g ˙ (kT + (φl mod pj ? φm i terms with k = 0 correspond to the term g ˙ (φT ) in Section 5.2. Recall that stability required g ˙ (φT ) to be small. Terms for k > 0 will be positive if the synapse decays fast enough. This then implies that stable harmonic locking will involve spiking patterns where the spikes from the two neurons are generally either coincident, close together or fairly far apart (see Fig. 1). This was observed numerically [21]. Numerically, harmonic locking was more likely to be observed if self-inhibition was included [21]. One reason may be that the locking conditions (5.16) are easier to satisfy with self-inhibition because the equations become more symmetric. Another reason may be that the transition from synchrony to suppression is slower as the parameters are changed for neuronal networks with self-inhibition compared to those without it, leading to a larger region in parameter space where harmonic locking is observed.

6. N coupled neurons We now examine the dynamics of a network of N neurons. Many of the results from two coupled neurons can be carried over to the N -neuron network. For N homogeneous neurons there are many possible phase-locked solutions. As discussed in Section 3, all symmetric combinations of the phases are solutions. Theorem 1 gives the condition for stability. We examine the existence and stability of synchronous, splay-phase, clustered, and unlocked states. Large homogeneous networks have been examined previously using different formalisms [11,12,15]. We consider the effects of weak heterogeneity. It is known that strong heterogeneity can induce various complicated states [18,19]. 6.1. Synchrony For a homogeneous network, the synchronous state is a phase-locked solution. Stability has been examined in the limit of large N by Gerstner et al. [12] (see Remark 2). They proved that as long as the in?uence of the synaptic kernel is rising at the time of the next spike then synchrony is stable, i.e. l ≥1 J ˙ (lT ) > 0, where J < 0 (J > 0) corresponds to inhibitory (excitatory) coupling. Implicit in their proof is that the synaptic kernel does not act immediately (i.e. ˙ (0) = 0). With the inclusion of heterogeneity, the phases are shifted apart for a locked solution. However, for weak heterogeneity, a stable near synchronous solution is still possible. This was observed in numerical simulations [21]. Here we illustrate that for inhibitory coupling, if the neurons are indexed so that I1 > I2 > · · · > IN then there could be a stable phase-locked solution where the phases are arranged as φ1 < φ2 < · · · < φN provided the heterogeneity is small and close to uniformly distributed. We translate time so that the phases are all positive and small.

364

C.C. Chow / Physica D 118 (1998) 343–370

First consider a homogeneous system with a period satisfying the condition ?+ 1=I
l ≥1

[η(lT ) ? g(N ? 1) (lT )].

(6.1)

? ≡ (1/N) j Ij . We now let Ii = I ? + ?i . Since the applied currents are all positive where J = ?g , g > 0, and I then ?1 > 0 and ?N < 0. If a phase-locked solution exists in the heterogeneous system, it must satisfy the condition ? + ?i + 1=I
l ≥1

η(lT ) ?
j =i,l ≥1

g (lT + (φi ? φj )T ) ?
j <i

g ((φi ? φj )T )

(6.2)

Subtracting condition (6.1) from (6.2) gives ?i =
j =i,l ≥1

g [ (lT + (φi ? φj )T ) ? (lT )] +
j <i

g ((φi ? φj )T ).

(6.3)

Consider condition (6.3) for neuron 1 (i = 1): ?1 =
j =1,l ≥1

g [ (lT + (φ1 ? φj )T ) ? (lT )].

(6.4)

Since we have assumed that φ1 < φj , then condition (6.4) can be satis?ed if (t) is monotone decreasing for t > T + (φ1 ? φN )T . Now consider condition (6.3) for neuron N : ?N =
j <N,l ≥1

g [ (lT + (φN ? φj )T ) ? (lT )] +
j <N

g ((φN ? φj )T ).

(6.5)

In this case ?N < 0 and φN > φj . Condition (6.5) can be satis?ed if (t) is monotone decreasing for t > T + (φN ? φN ?1 )T and ((φN ? φj )T ) is not too large. We can perform a similar analysis for all the other neurons. The result is that as long as (t) is monotone decreasing for t > T + (φ1 ? φN )T and it does not rise up too quickly for small t then near synchrony is possible with the phases arranged as φ1 < φ2 < · · · < φN , i.e. they will all ?re in rapid succession. We should note that the period is controlled by the applied currents, the synaptic strengths, and the synaptic time course. Stability is determined from the perturbation series for stability (4.5) which will be composed of terms with coef?cients of the form η(lT ˙ ), ?g ˙ (lT + (φi ? φj )T ), and ?g ˙ ((φi ? φj )T ), for all i , j and l . By Theorem 1, the phase-locked solution is stable if all the coef?cients are positive. The monotone decreasing condition guarantees that the ?rst two sets of coef?cients are positive. For a homogeneous system, the third set is zero. For a heterogeneous system, the third set can be negative for Type I neurons. However, by constructing the return mapping (4.8) it can be shown that stability is possible as long as they are small by using Corollary 1 and generalization of Proposition 2 to the heterogeneous N neuron network. In the case of Type I excitatory coupling and zero rise time of the synaptic current, synchrony can be globally stable but it can be locally unstable [12]. If a neuron ?res late with respect to the others then it will be pushed forward and eventually catch up to the synchronous population, but if it ?res early it will be pushed away until it catches up to the next cycle. Mirollo and Strogatz [9] considered an absorption condition so that once two neurons ?re within a given time window they are considered synchronized and remain synchronized. This absorption condition ensures global stability. 6.2. Splay-phase and clustered states For a homogeneous network, the splay-phase state also known as anti-phase, ‘rotating wave’ or the more colorful ‘ponies on a merry-go-round’ [11,15,31–34] is a possible phase-locked solution. In this state, all N neurons in the

C.C. Chow / Physica D 118 (1998) 343–370

365

network ?re with a phase separation of 1/N . This solution can be veri?ed by examining condition (3.2). With weak heterogeneity, a near splay-phase state should exist by a similar argument to the perturbed synchronous case. We can show the stability of the splay-phase state with an N neuron version of Proposition 1. We note that the coef?cients of (4.5) are composed of factors of the form η(lT ˙ ) and J ˙ (lT ± j φT ) where φ = 1/N , and j = 1, . . . , N is the neuron index. If these factors are positive then there is stability by Theorem 1. In the excitatory case, the splay-phase solution can be stable if the rise time of the synapse is very slow so that ˙ (j φT ) > 0. In this case the splay-phase is stable for similar reasons as to why the anti-synchronous solution is stable for two neurons. This has been noted previously [11,15]. With weak heterogeneity, stability of the near splay-phase state can still be stable provided the synapse is slow enough. The splay-phase state can also be stable for Type I inhibitory connections if the synapse has a fast enough rise time so that ˙ (T /N) < 0 and (t) is monotone decreasing for t > T /N . Stability may be affected if the action potentials (spikes) of the individual neurons overlap in the splay-phase state. This can occur if the network is so large that T /N is on the order of the spike width. We have veri?ed the existence of the splay-phase state numerically in our biophysical model for N = 4. The stability of the splay-phase is interesting because in the mean ?eld limit (N → ∞) it has been shown to be unstable for inhibitory coupling [26,29,43]. In the mean ?eld limit for homogeneous neurons, the splay-phase state corresponds to an ‘asynchronous’ state where the density of the ?ring times of the neurons is constant. We can see how this instability arises in our formalism. As more neurons are added, the synapse must rise faster for stability. In the limit of in?nite numbers of neurons, stability is not possible no matter how fast the synapse since (T /N) → 0. With weak heterogeneity, a stable perturbed splay-phase state is possible as long as the rise time of the synapse is shorter than the smallest phase difference. Stability should be more affected by heterogeneity as the network size increases since the phase differences become smaller. There is also the possibility of phase-locked clustered states. Clustering in neuronal networks has been investigated before [11,15,44–46]. For homogeneous networks, groups of neurons can ?re together in unison but at different phases with other groups. Clustered states can be shown to be stable if ˙ (0) = 0 and J ˙ (lT + (φa ? φb )T ) > 0, where a and b are the indices for different groups. There are many possible coexistent clustered states that depend on the initial conditions. With heterogeneity, clustered solutions of near synchronous neurons are possible as in the near synchronous solution. If the heterogeneity is nonuniform and clustered then groups may form more easily among neurons with similar intrinsic frequencies. 6.3. Asynchrony, harmonic locking and suppression Heterogeneity breaks the symmetry of the network. As we have seen above, if the heterogeneity is mild then phaselocked solutions can exist. For large enough heterogeneity there are no phase-locked solutions. When phase-locking does not exist several possibilities can arise. We have discussed asynchrony, harmonic locking and suppression for two coupled neurons. These states have their counterparts for a network of N neurons as well. If the heterogeneity is large then our analysis in Section 6.1 shows that stable synchrony can be broken. If the applied currents are very strong or the synapses are very weak then the period can become too short to support stable synchrony. In this case, the neurons will become effectively decoupled and ?re asynchronously as in the two neuron case. This was observed numerically [21]. For heterogeneous neurons, the average density of neurons at a given phase need not be constant as in the homogeneous case. If the heterogeneity is very broad then it is possible that a set of neurons can be synchronous whereas the remaining set can be asynchronous [18,19]. If the applied currents are weak and the synapses are strong then stable synchrony can again be broken. For inhibition, the possible outcomes are harmonic locking or suppression. Here the slower neurons cannot keep up with the faster ones. The dynamics will be similar to two self-inhibited neurons. The in?uence of the other neurons in the network with weakly heterogeneous applied currents acts like self-inhibition. If a neuron is too slow compared

366

C.C. Chow / Physica D 118 (1998) 343–370

to a group of other neurons, then it can drop out of the rhythm. As the synaptic strength is increased, more and more neurons will drop out until only one remains. The condition for suppression can be estimated from the combined synaptic input the neuron receives over a period. Similarly, complicated harmonic patterns can form where the slower neurons lock to sub harmonics of faster neurons. In this case a large set of coexistent patterns exist [21].

7. Discussion and conclusions We have examined the existence and stability of phased-lock solutions in a network of neurons with heterogeneous intrinsic frequency. With mild heterogeneity, we have shown that stable phase-locking is possible but it is fairly fragile. When periodic phase-locking breaks, states of asynchrony, harmonic locking or suppression may arise. In accordance with previous work [10–12], we have found that phase-locking tendencies are strongly determined by the type (I or II) of the neuron model. We also distinguish between saturating and nonsaturating synaptic types and ?nd that there can be differences in phase-locking properties. Neurons with saturating synapses tolerate the effects of heterogeneity better for stability of synchronization. Our analysis may clarify some of the apparent differences in previously derived conditions for synchronization with inhibition. We ?nd that stable synchrony requires a nonzero rise time of the synaptic current as emphasized in [10,22]. This was implicitly assumed in [12]. We also ?nd that it requires that the rise time not be too slow which was emphasized in [12] and assumed in [22]. However, Ref. [22] also requires that the synaptic decay time not be too fast which was also observed in numerical simulations in [39]. This condition may be implicit in the behavior of the kernels. For some neuron models, the kernel may become negative if the synaptic decay time becomes too fast. For example, if the threshold of the integrate-and-?re model changes with the synapse then this condition may be present [22]. Our results rely on the spike response method to characterize the dynamics of a threshold crossing variable in terms of a set of response kernels. The existence and stability of phase-locked solutions are then established by the time courses of these kernels. Critical to our analysis is that the kernels exist and can be computed. For the integrate-and-?re model the kernels can be obtained explicitly. For more complicated models, they most likely must be obtained numerically [27]. However, we can still infer general properties on how they should behave. The neuron type (I or II) governs the general characteristics of the synaptic kernels. For example, positive inputs always advance the ?ring of Type I neurons. This implies that the synaptic kernel always remains positive. As noted previously [11,23], we ?nd that Type I and Type II synapses will greatly affect the phase-locking tendencies of the neurons. Our formalism shows why this is the case from the shapes of the synaptic kernels. It should be noted that the expansion into kernels is not always guaranteed. In particular, the separation of a membrane kernel from a synaptic kernel may not always be possible if the neuron is strongly nonlinear. In such a case, an integral approximation can still to be obtained and a linearization around phase-locked solutions may be possible to examine stability. We emphasize that we have computed the local stability of the states. This is in contrast to previous strategies where the global stability is computed [8,9,13]. In these works, the class of initial conditions that will approach a synchronous solution is found. An absorption criterion is used that assumes that once the oscillators are within a given phase of each other then they remain locked and no longer in?uence each other. In our formulation, the synaptic coupling can still in?uence locked neurons and be enough to destabilize the dynamics if certain other criteria are not met. In the case of excitatory coupling, we ?nd that although oscillators may approach each other globally they may not be locally stable. It had been previously assumed that the reasons for disparate results for the stability of synchronized pulse-coupled integrate-and-?re oscillators were a result of the singular nature of the

C.C. Chow / Physica D 118 (1998) 343–370

367

delta-distribution synaptic coupling [11,12]. However, it is the assumption regarding what constitutes a locked state that separates the conclusions. Ideally, a combined analysis should be done where a global stability analysis is used to examine the class of initial conditions that approach a phase-locked solution and then a local stability analysis is performed to see if the locked state is sustainable. We only considered heterogeneity in the applied currents which is a very simple form of disorder because it is not necessary to average over different instances of the heterogeneity. For instance if we choose N random current values from a ?xed probability distribution and if N is large enough, then a single instance is equivalent to any other instance. This is because the network is insensitive to rearrangements of the current values between the neurons. We can always relabel the neurons in the network. Only the distribution of the gaps between current values are important and this is ?xed for large enough N . This property allowed us to gain insight into the N neuron network from the two neuron network. On the other hand if the heterogeneity were to be included in the coupling strengths then a single instance of the heterogeneity may not give an accurate picture of what a typical random network will do. A rearrangement of the coupling strengths could greatly affect the network dynamics. Here, an understanding of two heterogeneously coupled neurons may not give any insight into the behavior of a large network. One must average over the heterogeneity to get an idea of what a typical network will do. Then one must probe the ?uctuations around this average. For two neurons it seems that heterogeneity for synaptic strengths would have similar results to what we derived for heterogeneity in the applied current. Phase-locked solutions would be shifted away from the homogeneous solutions. Our stability results would then apply analogously. However, the behavior of large networks may be quite different from the behavior of two neurons.

Acknowledgements I would like to thank N. Kopell for many clarifying discussions and for constructive critiques of the manuscript. I thank S. Epstein, B. Ermentrout, W. Gerstner, and J. White for many helpful suggestions. I especially thank J. Ritt for performing numerical simulations, generating ?gures, and for his insightful comments on the manuscript. This work was supported by the National Science Foundation grant DMS-9631755.

Appendix A. Conductance-based neuron model For numerical simulations we used a single-compartment model neuron with inhibitory synapses obeying ?rstorder kinetics [21]. Membrane potential in each point neuron obeyed the current balance equation C dVi = Ii ? INa ? IK ? IL ? Is , dt (A.1)

4 where C = 1 ?F/cm2 , Ii is the applied current, IN a = gN a m3 ∞ h(Vi ? VN a ) and IK = gK n (Vi ? VK ) are the N gs Sj (t)(Vi ?Vs ) Hodgkin–Huxley type spike generating currents, IL = gL (Vi ?VL ) is the leak current and Is = j is the synaptic current. The ?xed parameters used were: gN a = 30 mS/cm2 , gK = 20 mS/cm2 , gL = 0.1 mS/cm2 , VNa = 45 mV, VK = ?75 mV, VL = ?60 mV, Vs = ?75 mV. The phenomena described here seem largely independent of speci?c neuronal parameters. The activation variable m was assumed fast and substituted with its asymptotic value m∞ (v) = (1+exp[?0.08(v + 26)])?1 . The inactivation variable h obeys

h∞ (v) ? h dh = , dt τh (v)

(A.2)

368

C.C. Chow / Physica D 118 (1998) 343–370

with h∞ (v) = (1 + exp[0.13(v + 38)])?1 , τh (v) = 0.6/(1 + exp[?0.12(v + 67)]). The variable n obeys dn n∞ (v) ? n = , dt τn (v) with n∞ (v) = (1 + exp[?0.045(v + 10)])?1 , τn (v) = 0.5 + 2.0/(1 + exp[0.045(v ? 50)]). The gating variable Sj for the synapse is assumed to obey ?rst-order kinetics of the form dSj = F (Vj )(1 ? Sj ) ? Sj /ts , dt (A.4) (A.3)

where F (Vj ) = 1/(1 + exp[?Vj ]), and Vj is the voltage of the pre-synaptic neuron and transmission delays are neglected. The ODEs were integrated using a fourth-order Runge–Kutta method.

Appendix B. Phase-locking of two coupled neurons with no memory As an example, we consider a model where the response kernels have no memory beyond the previous ?ring. The phase-locked solutions are then given by v1 (T + φ1 T ) = 1 = I1 + η(T ) + J (T ? φT ), v2 (T + φ2 T ) = 1 = I2 + η(T ) + J (φT ) + J (T + φT ), (B.1) (B.2)

This implies that G(φ) = J [ (φT ) + (T + φT ) ? (T ? φT )]. The perturbed dynamics equations (4.7) are given by δ(n) = M · δ(n ? 1), where ? ? η(lT ˙ ) J ˙ (lT ? φT ) ? ? δ1 (n) v ˙1 v ˙1 ?, δ(n) = , M=? (B.3) ? J ˙ (lT + φT ) J ˙ (φT )η(lT ˙ ) η(lT ˙ ) J 2 ˙ (φT ) ˙ (lT ? φT ) ? δ2 (n) + + v ˙2 v ˙1 v ˙2 v ˙2 v ˙1 v ˙2 ˙ ) + J ˙ (T ? φT ), (B.4) v ˙1 = η(T ˙ ) + J ˙ (φT ) + J ˙ (T + φT ), v ˙2 = η(T (B.5) and the prime refers to the ?rst derivative with respect to time. By de?nition, v(t) must approach the threshold from below which implies that v ˙i > 0. The eigenvalues can be solved for exactly. We can verify that the rows of M(l) in (B.3) sum to unity and so one eigenvalue is 1. The other eigenvalue is then λ = Tr M ? 1 ≡ det M . For stability, we require |λ| < 1 which implies that 0 < Tr M < 2 or ?1 < det M < 1, where η(T ˙ ) η(T ˙ ) ˙ (φT ) ˙ (T ? φT ) + + J2 , v ˙1 v ˙2 v ˙1 v ˙2 ˙ (T ? φT ) ˙ (T + φT ) η(T ˙ ) ? J2 . det M = v ˙1 v ˙2 v ˙1 v ˙2 Tr M = (B.6) (B.7)

We now apply our results of Section 4 and compare to the exact eigenvalues. If all the elements of M are positive then we immediately see that stability is achieved. Since the rows must sum to 1 and the elements are positive then each element itself is less than 1. This immediately implies that 0 < det M < 1 as required for stability. However, this condition although suf?cient for stability is certainly not necessary. An examination of Tr M or det M shows that the elements of M can become negative and still maintain ?1 < det M < 1 required for stability.

C.C. Chow / Physica D 118 (1998) 343–370

369

We next examine the necessary condition G (φ) > 0 of Theorem 2. We ?rst note that G (φ) = J T [ ˙ (φT )+ 1 (T + φT ) + ˙ (T ? φT )]. Consider the synchronous solution for synchrony (φ = 0) where G (0) = J T [ ˙ (0) + 2 ˙ (T )] > 0 is clearly necessary for stability. If in addition, ˙ (0) = 0 and η(T ˙ ) > 0 then an inspection of condition (B.6) shows that G (0) > 0 is necessary and suf?cient for stability. For the anti-synchronous solution (φ = 0.5), G (0.5) = J T [2 ˙ (0.5T ) + ˙ (1.5T )] > 0 is necessary for stability from condition (B.6) or (B.7). However, unlike the synchronous case, it cannot be made into a suf?cient condition without making additional nontrivial constraints on ˙ (0.5T ) and ˙ (1.5T ). References
[1] R. Llin? as, U. Ribary, Coherent 40-Hz oscillation characterizes dream state in humans, Proc. Natl. Acad. Sci. USA 90 (1993) 2078–2081. [2] C.M. Gray, Synchronous oscillations in neuronal systems: mechanisms and functions, J. Comput. Neurosci. 1 (1994) 11–38. [3] M.A. Whittington, R.D. Traub, J.G.R. Jefferys, Synchronized oscillations in interneuron networks driven by metabotropic glutamate receptor activation, Nature 373 (1995) 612–615. [4] J.G.R. Jefferys, R.D. Traub, M.A. Whittington, Neuronal networks for induced ‘40 Hz’ rhythms, Trends Neurosci. 19 (1996) 202–207. [5] R.D. Traub, M.A. Whittington, S. Colling, G. Buzs? aki, J.G.R. Jefferys, Analysis of gamma rhythms in the rat hippocampus in vitro and in vivo, J. Physiol. 493 (1996) 471–484. [6] A.T. Winfree, Biological rhythms and the behavior of populations of coupled oscillators, J. Theoret. Biol. 16 (1967) 15–42. [7] S. Strogatz, Norbert Wiener’s Brain Waves, Lecture Notes in Biomathematics, vol. 100, Springer, Berlin, (1993). [8] A.V.M. Herz, J.J. Hop?eld, Earthquake cycles and neural reverberations: collective oscillations in systems with pulse-coupled threshold elements, Phys. Rev. Lett. 75 (1995) 1222–1225. [9] R. Mirollo, S. Strogatz, Synchronization of pulse-coupled biological oscillators, SIAM J. Appl. Math. 50 (1990) 1645–1662. [10] C. van Vreeswijk, L. Abbott, G.B. Ermentrout, When inhibition not excitation synchronizes neural ?ring, J. Comput. Neurosci. 1 (1994) 313–321. [11] D. Hansel, G. Mato, C. Meunier, Synchrony in excitatory neural networks, Neural Comput. 7 (1995) 307–337. [12] W. Gerstner, J.L. van Hemmen, J. Cowen, What matters in neuronal locking?, Neural Comput. 8 (1996) 1653–1676. [13] S. Bottani, Pulse-coupled relaxation oscillators: from biological synchronization to self-organized criticality, Phys. Rev. Lett. 74 (1995) 4189–4192. [14] U. Ernst, K. Pawelzik, T. Geisel, Synchronization induced by temporal delays in pulse-coupled oscillators, Phys. Rev. Lett. 74 (1995) 1570–1573. [15] C. van Vreeswijk, Partial Synchronization in populations of pulse-coupled oscillators, Phys. Rev. E 547 (1996) 5522–553. [16] Y. Kuramoto, in: H. Araki (Ed.), International Symposium on Mathematical Problems in Theoretical Physics, Lecture Notes in Physics, vol. 39, Springer, New York, 1975. [17] S. Strogatz, R. Mirollo, Collective synchronisation in lattices of non-linear oscillators with randomness, J. Phys. A. 215 (1988) L699–L70. [18] D. Golomb, J. Rinzel, Dynamics of globally coupled inhibitory neurons with heterogeneity, Phys. Rev. E 484 (1993) 4810–481. [19] M. Tsodyks, I. Mitkov, H. Sompolinsky, Phys. Rev. Lett. 71 (1993) 1280–1283. [20] X.-J. Wang, G. Buzs? aki, Gamma oscillation by synaptic inhibition in an interneuronal network model, J. Neurosci. 163 (1996) 6402–641. [21] J.A. White, C.C. Chow, J. Ritt, C. Soto-Trevino, N. Kopell, J. Comput. Neurosci. 5 (1998) 5–16. [22] D. Terman, N. Kopell, A. Bose, Dynamics of two neurons coupled by mutual slow inhibition, Physica D 117 (1998) 241–275. [23] G.B. Ermentrout, Type I membranes, phase resetting curves, and synchrony, Neural Comput. 81 (1996) 979–100. [24] J. Rinzel, G.B. Ermentrout, in: C. Koch, I. Segev (Eds.), Methods in Neuronal Modelling, MIT Press, Cambridge, 1989. [25] W. Gerstner, J.L. van Hemmen, Network 34 (1992) 139–16. [26] W. Gerstner, Time structure of the activity in neural network models, Phys. Rev. E 518 (1995) 738–75. [27] W. Kistler, W. Gerstner, J.L. van Hemmen, Reduction of the Hodgkin–Huxley Equations to a single-variable threshold model, Neural Comput. 90 (1997) 1069–110. [28] C. Chow, J.A. White, J. Ritt, N. Kopell, Frequency control in synchronous networks of inhibitory neurons, J. Comput. Neurosci., in press. [29] L. Abbott, C. van Vreeswijk, Asynchronous states in networks of pulse-coupled oscillators, Phys. Rev. E 48 (1993) 1483–1490. [30] H. Tuckwell, Stochastic Processes in the Neurosciences, SIAM, Philadelphia, 1989. [31] D.G. Aronson, M. Golubitsky, J. Mallet-Paret, Ponies on a merry-go-round in large arrays of Josephson junctions, Nonlinearity 40 (1991) 903–91. [32] J.W. Swift, S.H. Strogatz, K. Wiesenfeld, Averaging of globally coupled oscillators, Physica D 55 (1992) 239–250.

370

C.C. Chow / Physica D 118 (1998) 343–370

[33] S.H. Strogatz, R.E. Mirollo, Splay states in globally coupled Josephson arrays: analytical prediction of Floquet multipliers, Phys. Rev. E 477 (1993) 220–22. [34] I.B. Schwartz, K.Y. Tsang, Antiphase switching in Josephson junction arrays, Phys. Rev. Lett. 73 (1994) 2797–2800. [35] C. Peskin, Mathematical Aspects of Heart Physiology, Courant Institute of Mathematical Sciences, New York University, New York, 1975, pp. 268–278. [36] G.B. Ermentrout, N. Kopell, Frequency plateaus in a chain of weakly coupled oscillators, SIAM J. Math. Anal. 157 (1984) 215–237. [37] Y. Kuramoto, Chemical Oscillations, Waves and Turbulence, Springer, Berlin, 1984. [38] L. Abbott, A network of oscillators, J. Phys. A 23 (1990) 3835. [39] X.-J. Wang, J. Rinzel, Alternating and synchronous rhythms in reciprocally inhibitory model neurons, Neural Comput. 4 (1992) 84–97. [40] G.B. Ermentrout, N. Kopell, Inhibition-produced patterning in chains of coupled nonlinear oscillators, SIAM J. Appl. Math. 54 (1994) 478–507. [41] R.A. Horn, C.A. Johnson, Matrix Analysis, Cambridge University Press, Cambridge, 1985. [42] P. Pinsky, SIAM J. Appl. Math. 55 (1995) 220–241. [43] A. Treves, Mean-?eld analysis of neuronal spike dynamics. Network 4 (1993) 259–284. [44] D. Golomb, D. Hansel, B. Shraiman, H. Sompolinsky, Phys. Rev. A 45 (1992) 3516. [45] D. Golomb, J. Rinzel, Clustering in globally coupled inhibitory neurons, Physica D 72 (1994) 259–282. [46] N. Kopell, G. Le Masson, Rhythmogenesis, amplitude modulation and multiplexing in a cortical architecture, Proc. Nat. Acad. Sci (USA) 91 (1994) 10 586–19 590.


赞助商链接

更多相关文章:
更多相关标签:

All rights reserved Powered by 甜梦文库 9512.net

copyright ©right 2010-2021。
甜梦文库内容来自网络,如有侵犯请联系客服。zhit325@126.com|网站地图