Clustering in neural networks with heterogeneous and asymmetrical coupling strengths

Physica D 180 (2003) 210–234

Clustering in neural networks with heterogeneous and asymmetrical coupling strengths Yue-Xian Li∗ Departments of Mathematics and Zoology, University of British Columbia, Vancouver, BC, Canada V6T 1Z2 Received 20 September 2002; received in revised form 20 February 2003; accepted 24 February 2003 Communicated by Y. Kuramoto

Abstract The existence and stability of phase-clustered states have been studied previously in networks of weakly coupled oscillators with uniform coupling strengths [Physica D 63 (1993) 424]. However, several studies have shown that if the coupling is uniform and repulsive, it is hard to obtain stable phase-clustered states in networks of realistic neural oscillators when noise is present [Neural Comput. 7 (1995) 307; Phys. Rev. E 57 (1998) 2150]. This problem was avoided by introducing heterogeneity in the distribution of coupling strengths [J. Phys. Soc. Jpn. 72 (2003) 443]. It has been shown that heterogeneous coupling strengths make the occurrence of stable clustered states possible in small networks of repulsively coupled neural oscillators of all kinds [J. Comput. Neurosci. 14 (2003) 139; SIAM J. Appl. Math., submitted for publication]. The present work extends these results to large networks of N identical neurons that are globally coupled with heterogeneous and asymmetrical coupling strengths. Conditions for the existence and stability of a state of n synchronized clusters at evenly distributed phases, called the state of n splay-phase clusters, are derived. Clusters of different sizes, i.e. containing different numbers of neurons, are studied. The existence of such a state is guaranteed if the strength of the coupling originating from one neuron to other neurons is inversely proportional to the size of the cluster to which it belongs. This condition is called the rule of inverse cluster-size. At the state of n splay-phase clusters, the N -neuron network behaves like a network of n “big neurons”. Stability of this state is determined by n eigenvalues of which only one determines the stability of intra-cluster phase differences. The remaining n − 1 conditions determine the stability of inter-cluster phase differences, but only nh = (n − mod (n, 2))/2 of them have distinct real parts due to symmetry. Heterogeneous coupling makes the stability conditions depend on coupling strengths. This analysis not only reveals how clustered states occur in more general kinds of networks, but also illustrates how the stability of clustered states can be achieved in networks of repulsively coupled neural oscillators. Results on clustered states with phases that are not evenly distributed in the phase space are also presented. Potential applications of these results are discussed. © 2003 Elsevier Science B.V. All rights reserved. PACS: 07.05.Mh; 84.35.+i; 87.10.+e; 87.22.Jb Keywords: Neural networks; Synchronization; Phase clustering; Phase locking; Coupled oscillators

∗

Tel.: +1-604-822-6225; fax: +1-604-822-6074. E-mail address: [email protected] (Y.-X. Li). 0167-2789/03/$ – see front matter © 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0167-2789(03)00067-8

Y.-X. Li / Physica D 180 (2003) 210–234

211

1. Introduction A large variety of biological functions are accomplished through highly coordinated actions of cells. Thus, the origin of many biological rhythms can be attributed to synchronous rhythmic outputs from networks of coupled excitable cells or neurons. For example, heart beat results from orchestrated activities of muscle cells [7], hormonal rhythms are driven by rhythmic outputs from networks of endocrine neurons in the hypothalamus [8], respiratory rhythms has been found to originate from several hundred neurons in the preBötzinger Complex in the ventrolateral medulla of the brain stem [9], patterns of lamprey locomotion are produced by a central pattern generator (CPG) located in the spinal cord of the animal [10,11]. In many cases, synchrony between coupled neurons or cells is crucial in the occurrence of these rhythms. Therefore, understanding the mechanisms leading to robust synchrony is crucial for us to understand the origin of many biological rhythms and to find novel ways to treat diseased rhythms [12]. Of particular interest are the CPGs that generate rhythmic patterns of locomotion in multi-legged animals. Typical patterns observed in the movement of legged animals are called gaits [13–15]. Thus, trot is a symmetrical gait widely used by many four-legged animals such as horse and turtle. Each leg is in synchrony with its diagonal counter-part, but at anti-phase with its contralateral counter-part. A similar pattern is used by hexapodal runners such as cockroaches, called a tripod gait [15,16]. In tripod gait, the first and the third legs on one side of the body and the second leg on the opposite side move in-phase, while the remaining three legs also move in-phase. But the two triplets are at anti-phase, i.e. each leg moves at anti-phase with its contralateral counter-part. In the study of networks of weakly coupled neural oscillators, these states are called phase-clustered states, in which neurons in the network are divided into two or more synchronized clusters [17,18]. Conditions for the existence and stability of a variety of phase-clustered states have been obtained in small networks containing four and six neural oscillators derived from realistic neural models [5,6]. These studies revealed that allowing heterogeneity in the distribution of coupling strengths is often a necessity for the stability of these clustered states in all known neural models studied so far. The present work is to extend these results to large networks containing N neurons and to patterns containing n clusters. Influence exerted by each one of a pair of neurons on the other is determined by the nature of the coupling and the sensitivity of each neuron to the coupling. The nature of coupling is determined by three factors: (1) the coupling strength; (2) the direction of polarization; (3) the time scale [19]. When the coupling is strong, the coupling current can overwhelm the intrinsic ionic currents. The behavior of a coupled neuron can be qualitatively different from its behavior in isolation, i.e. when de-coupled from the other neurons. In many realistic networks of synaptically coupled neurons, however, the influence of each individual synapse is weak. Strong coupling can occur through temporal summation of a large number of synaptic currents. When the coupling is weak, the intrinsic properties of each coupled neuron remain unchanged. In this case, coupling plays the role of coordinating the activities of neurons in the network. In this paper, we shall focus on networks of weakly coupled neurons. In a weakly coupled network of neural oscillators, the outcome strongly depends on the direction of polarization caused by the coupling current. Gap-junctional coupling is always bidirectional, because it could be either polarizing or depolarizing depending on the voltage difference between the coupled neurons. Synaptic currents can be either excitatory or inhibitory as well as bidirectional depending on the synaptic reversal potential. Another important factor is the time scale of synaptic coupling. The classical notions that excitation leads to synchrony and inhibition leads to asynchrony were based on the assumption that the coupling is fast [20–22]. The outcome is often reversed when the coupling is slow [2,23,24]. In the present study, we shall generally assume that the coupling is slow enough so that excitation leads to asynchrony and inhibition leads to synchrony unless specified otherwise. In the study of network models of CPGs that generate rhythmic gaits observed in multi-legged animals, symmetry plays an important role [14,15]. Given the symmetry of a network determined by a given distribution of coupling strengths, symmetry analysis reveals all possible types of model-independent gait patterns that can occur. However, a couple of recent studies [5,6] of coupled phase oscillators showed that all model-independent solutions coexist.

212

Y.-X. Li / Physica D 180 (2003) 210–234

Furthermore, it was shown that these model-independent states often coexist with many other weakly modeldependent states. Thus, knowing that a phase-locked state can occur does not guarantee that this state can actually be observed both in numerical simulations of the model network and in experimental measurement of the realistic network. The state that can actually be observed must be stable even in the presence of small noise. However, if multiple coexisting states are stable under the same conditions, the occurrence of one such state must also depend on the initial phase relations between different neurons. For the network to evolve eventually into one of these states, the initial phases must be in the attraction basin of that state. These results indicate that knowing the existence of certain gait patterns in a network is important, but knowing how to make only one such pattern occur even when it coexists with all the other patterns is also important. The latter is called pattern segregation [6] and can only be achieved by analyzing the stability conditions of these patterns. The problem of pattern segregation becomes even more important as the network gets larger, since the number of coexisting patterns increases dramatically as the size of the network increases [5]. Thus, a mechanism that allows the network to selectively stabilize each one of the coexisting patterns is important for a large network to carry out its pattern-generating functions. It has been shown recently [5,6] that changing the coupling parameters in a network with heterogeneous coupling strengths can provide a sufficient mechanism for pattern segregation. However, this is achievable only when conditions for the stability of these patterns are known together with the conditions for their existence. The present work provides a theoretical framework for determining the stability of phase-clustered states in a large network of coupled neurons. In Section 2, previously published results on the existence and stability of clustered states in a network of uniformly coupled neurons are discussed. The need for extending these results to networks with heterogeneous and asymmetrical distribution of coupling strengths is presented. In Section 3, the spatial structure of a network is introduced and the phase equations describing a structured network are derived. In Section 4, the conditions for the existence of clustered states are obtained and the rule of inverse cluster-size is revealed. In Section 5, conditions that are both necessary and sufficient for the stability of phase-clustered states are derived. In Section 6, eigenvalue spectrum determining the stability of the asynchronous state is calculated in the limit when the network becomes infinitely large. In Section 7, the existence and stability conditions of clustered states with phases that are not evenly distributed in the phase space are presented. The theoretical and practical significance of the results presented here is discussed in Section 8.

2. Splay-phase clusters in networks of uniformly coupled neurons In this section, we first review results previously obtained on the existence and stability of phase-clustered states in networks of uniformly coupled oscillators. Then, we present reasons for why an extension of these results to networks of all-to-all coupled oscillators with heterogeneous and asymmetrical distribution of coupling strengths is of great importance in achieving stable phase-clustered states in networks of repulsively coupled neurons. These reasons have been demonstrated in a couple of recent studies [5,6]. We here present an overview of these results in comparison to those obtained in [1]. 2.1. Previous results on the existence and stability of clustered states A network of identical neural oscillators that are all-to-all coupled through weak and uniformly distributed coupling strengths can be described by the following system of coupled phase equations [25]: N dθi w
H (θj − θi ) = dt N j =1

(i = 1, 2, . . . , N ),

(1)

Y.-X. Li / Physica D 180 (2003) 210–234

213

where θi is the relative phase variable of the ith neuron and w > 0 the uniform coupling strength. H (·) is the phase-resetting function that describes the phase response of each neuron to the action of another neuron through either an excitatory or an inhibitory synapse, or a mixture of both. For more details of the phase theory see [25,26]. In summary, the H (·) function carries the information of how sensitive is each neuron to perturbations and how neurons are coupled to each other. In such a network, the existence and stability of a state of n splay-phase clusters of equal size (N/n) has been derived by Okuda [1]. Here, splay-phase clusters refer to synchronized clusters phase-locked at evenly distributed phases in the phase space which is simply a unit ring. Thus, if we use Θk to denote the phase of the kth cluster (k = 0, 1, 2, . . . , n − 1), the state of n splay-phase clusters is given in the form: Θk = f (n) t +

k n

(k = 0, 1, 2, . . . , n − 1),

(2)

where f (n) is the change in the oscillation frequency as a result of phase interactions between n different clusters. (2) is a solution of (1) provided that   n−1 w
k H f (n) = . (3) n n k=0

It is clear that f (n) is the average value of the phase-resetting function, H (·), sampled at evenly distributed points over one full cycle of oscillation. This value will determine whether the coupled neurons will oscillate with a higher (f (n) > 0) or lower (f (n) < 0) frequency than that of an isolated neuron. Linear stability analysis of this state yielded the following eigenvalues: λ

(n)

  n−1 w
 k =− H n n

(N − n multiplicity),

k=0

λ(n) p

  n−1 w
 k =− H (1 − e−i2πkp/n ) n n

(p = 0, 1, 2, . . . , n − 1).

(4)

k=0

In (4), λ(n) characterizes the stability of disturbances in the phase differences between neurons within one syn(n) chronous cluster, called intra-cluster phase disturbances [5], while λp (p = 0) characterizes the stability of inter-cluster phase disturbances. It is important to note that stability conditions determined by eigenvalues given by (4) are completely determined by model-dependent properties contained in the H (·) function and the number of clusters n. The coupling strength w (>0) plays no role in determining the sign of these eigenvalues. In a later publication [27], conditions sufficient for the stability of such clustered states were derived which had already been covered by the conditions given by (4). However, numerical simulations of realistic neuronal models coupled with uniform distribution of coupling strengths showed that the stability conditions given by (4) are rarely satisfied if the coupling is repulsive [2,3], i.e. if the coupling favors anti-phase relationship between a pair of neurons and destabilizes the in-phase relationship. As we shall explain next, repulsive coupling is often desired in order to avoid having a stable fully synchronous state in which all neurons oscillate in phase. 2.2. Clustering in networks of repulsively coupled neurons Synchrony in networks of coupled oscillators has long been a research topics interested by a broad spectrum of scientists. In particular, in the study of the CPGs that produce the periodic signals driving the movement of

214

Y.-X. Li / Physica D 180 (2003) 210–234

different legs in multi-legged animals, it has been well demonstrated that phase-clustered states of certain symmetry represent the most frequently observed patterns of animal locomotion or gaits [13–15]. The possible existence of such phase-clustered states has been studied extensively by analyzing symmetries in the distribution of coupling strengths [14,15,28]. Among the gaits observed in four-legged animals, trot, pace, and bound are the most frequently observed gaits involving synchronized movement of different legs. Several models of four coupled nonlinear oscillators have been studied extensively to show that phase-locked solutions corresponding to these gaits can be produced [14]. However, recent studies of networks of phase-coupled oscillators [5,6] showed that all such solutions coexist if the network possesses square or rectangular symmetry. They also coexist with other model-independent solutions such as pronk, characterized by periodic movement of four fully synchronized legs, and walk (when the network possesses square symmetry), which is characterized by evenly spaced phases between the four legs. In order to make the network generate only one such gait, we need to find conditions under which only one such gait solution is stable but all other coexisting solutions are unstable. To accomplish such a goal, it is necessary to study the stability conditions for these phase-clustered states. In a network of four identical coupled oscillators, trot, pace, and bound correspond to states of two anti-phase clusters, each containing two synchronized oscillators. When the coupling between each pair of neurons is synchronizing, e.g. when they are coupled through inhibitory synapses, stability conditions of these clustered states are easy to satisfy since synchrony is favored by synchronizing coupling [5]. However, such clustered states must coexist with the fully synchronized state which is also stable. Thus, perturbations in the phases of different legs can induce either a switch from trot, pace, or bound into pronk or a switch between the three clustered states. To avoid multi-stability in the coexisting gait solutions, repulsive coupling is favored when trot is the desired solution and is selected as the uniquely stable solution under certain choices of coupling parameters [5,6]. Based on the eigenvalues given by (4), the stability of a state of two anti-phase clusters (such as trot, pace, or bound) in a network of four oscillators with uniform coupling strength is guaranteed if λ(2) = − 21 w[H  (0) + H  ( 21 )] = − 21 w[g  (0) + g  ( 21 )] < 0, (2)

λ1 = −wH ( 21 ) = −wg ( 21 ) < 0,

(5) (6)

where g(φ) = [H (φ) − H (−φ)]/2 is the odd part of the phase-resetting function H (φ). Since H (φ) is 1-periodic, H  (φ) is also 1-periodic. Thus, H  (0) = g  (0) and H  (1/2) = g  (1/2). For slow excitatory coupling that is repulsive, condition (6) is satisfied since g  (1/2) > 0. However, when we consider realistic models of neurons coupled through realistic excitatory synapses, condition (5) can rarely be satisfied with typical choices of parameter values [5]. This is because g  (0) < 0 for repulsive coupling, which explains why the fully synchronous state is unstable with excitatory coupling. Thus, condition (5) is violated in excitatory networks if |g  (0)| > |g  (1/2)|. It was shown that in a number of realistic model neurons, |g  (0)| |g  (1/2)| is typically true [5]. This explains why robust clustered states have rarely been found in numerical simulations of repulsively coupled networks of realistic model neurons before the results in [5] were obtained. The only exception to our knowledge was the result in [4] in which heterogeneous distribution of coupling strengths was used. 2.3. Clustering in networks with heterogeneous and asymmetrical coupling strengths The difficulty in obtaining stable phase-clustered states in networks of repulsively coupled neurons was solved by the introduction of selective heterogeneity in the distribution of coupling strengths [5,6]. The main idea is actually rather simple. Although full synchrony in networks of repulsively coupled neurons cannot be achieved, partial synchrony in forms of clustered states is still possible provided that the stabilizing inter-cluster phase interactions are stronger than the destabilizing intra-cluster phase interactions. This is actually the physical meaning expressed

Y.-X. Li / Physica D 180 (2003) 210–234

215

by condition (5). Heterogeneity in the coupling that strengthens the inter-cluster phase interactions and/or weakens the intra-cluster phase interactions makes this condition easily achievable for all phase-coupled models of oscillators [5,6]. However, results in [5,6] were obtained only in small networks involving four and six neurons. There is an immediate interest to extend these results to large networks containing N neurons and to clustered states containing an arbitrary number of clusters, n. Furthermore, clustered states corresponding to gait patterns observed in multi-legged animals often involve clusters with different sizes, i.e. clusters containing different numbers of synchronized legs. Thus, canter is a gait in quadrupeds characterized by three “clusters”, one contains two synchronized legs and the other two contains only one leg each (see e.g. [14,15]). Also, the gait lurch in hexapods is characterized by two anti-phase clusters, one contains four legs and the other contains only two (see e.g. [15]). Thus, there is also a need to generalize results obtained in [1] and [5,6] to clusters of different sizes. As we shall demonstrate later, conditions for the existence of clustered states involving clusters of different sizes are not trivial. Asymmetrical coupling between two neurons in clusters of different sizes is one easily understandable way to achieve such a clustered state. Interestingly, the stability conditions for such asymmetric clustered states, derived without knowing the results in [1], turned out to be almost identical to those obtained in [1] when the phases of the clusters are evenly distributed. However, there is an important difference between results obtained here and those in [1]. With heterogeneous coupling strengths, the stability of such a clustered state is dependent on the coupling strengths as well as the model-dependent properties contained in H (·) and the total number of clusters, n. As shown in [5,6], the inclusion of coupling parameters in the conditions that determine the stability of clustered states has profound consequences on pattern selection in models of CPGs (see Section 8 for more discussion). Furthermore, results on clustered states with phases that are not evenly distributed previously obtained in small networks [5,6] will also be extended to large networks. This kind of clustered states was not considered in [1].

3. Structured networks with heterogeneous coupling strengths Consider a network of N identical neural oscillators that are all-to-all and weakly coupled. The phase description of this network is given by N


dθi wji H (θj − θi ) (i = 1, 2, . . . , N ), = dt

(7)

j =1

where θi is the relative phase variable of the ith neuron, and wji the strength of the synaptic coupling from the j th neuron to the ith neuron. To exclude self-coupling, we assume that wjj = 0 for all j . The only difference between (7) and (1) is that now wji is no longer assumed constant but differs for different values of i and j . However, the distribution of coupling strengths is not random but with certain symmetries that are necessary for model-independent clustered states to occur. Such a network is thus referred to as spatially structured. Symmetries in this spatial structure have strong influence on the existence and stability of clustered states in this network. Equations in (7) describe the dynamics of a nonlinear phase-coupled system on an N -dimensional torus. Since all the oscillators are identical as reflected by the independence of the phase-resetting function H (·) on the subscripts i, j , this torus is highly symmetrical. This rather strong constraint has profound consequences on the types of solutions that can occur as well as the stability of these solutions. Here is a brief summary of the approach used in this research. We first assume that neurons in the network are divided into n clusters. Neurons in each cluster are not necessarily synchronized. Such a hypothetical structure is called the desired solution structure which, at this stage, is unrelated to the spatial structure of the network. We express the phase equations based on this desired solution structure. We then assume that neurons in each cluster

216

Y.-X. Li / Physica D 180 (2003) 210–234

are synchronized and work backward to determine the conditions for the existence of such a state. These conditions are satisfied only when the coupling parameters are structured in a specific way, i.e. with desired symmetries that allow a such clustered state to exist. Thus, we obtain the spatial structure or symmetry group that is necessary for the existence of this clustered state. We then determine the conditions for the stability of this state by linearizing the phase equations in (7) near such a solution. We consider the general case in which the size (i.e. the number of neurons) in each cluster is different. Let us assume that there are np neurons in the pth cluster (p = 0, 1, 2, . . . , n − 1). Thus, n−1


np = N

(8)

p=0

for n = 1, n0 = N. The whole network itself becomes a single cluster network. For n > 1, i.e. for a network of multiple clusters, two subscripts are required in the definition of the phase variables. Thus, θpi denotes the phase of the ith neuron in the pth cluster, where i = 1, 2, . . . , np . Therefore, in terms of the new notation for the phase variables, the phase equation for the ith neuron in the pth cluster (p = 0, . . . , n − 1, i = 1, . . . , np ) becomes intra-cluster interactions

dθpi = dt

 np








wpj pi H (θpj − θpi ) +

j =1

inter -cluster interactions nq


n−1






wqj pi H (θqj − θpi ),

(9)

q=0,q=p j =1

in which the first term on the right-hand side specifies inputs from neurons within the same cluster (called intra-cluster interactions), and the second term describes inputs from neurons in other clusters (called inter-cluster interactions). Again, we assume that wpj pj = 0 for all p and j to exclude self-coupling. Phase-locked states are characterized by fixed phase differences between each pair of neurons in the network. Thus, they correspond to steady-state solutions of the equations that describe the phase differences in the network. Phase-clustered states are special phase-locked states in which neurons in each cluster are synchronized. We introduce φpi qj = θpi − θqj

(p, q = 0, 1, . . . , n − 1, i = 1, 2, . . . , np , j = 1, 2, . . . , nq )

(10)

as the phase difference between the ith neuron in the pth cluster and the j th neuron in the qth cluster. The equations that govern the time evolution of the phase differences are

φ˙ pi qj =

 np


intra-cluster interactions



wpk pi Hpk pi −

k=1

nq




wqk qj Hqk qj +

 ns n−1



inter -cluster interactions



wsk pi Hsk pi −

s=p k=1

k=1

ns n−1



 wsk qj Hsk qj ,

(11)

s=q k=1

where p, q = 0, 1, . . . , n − 1, i = 1, 2, . . . , np , j = 1, 2, . . . , nq . The argument of the H (·) function is abbreviated with the subscripts of the argument, e.g. H (θpi − θqj ) = H (φpi qj ) = Hpi qj . Let Φpi qj be a steady state of (11), i.e. the right-hand side vanishes for all p, q, i, j at this state:

0=

 np
k=1

intra-cluster interactions



wpk pi Hpk pi −

nq


k=1



wqk qj Hqk qj +

 ns n−1

s=p k=1

inter -cluster interactions



wsk pi Hsk pi −

ns n−1



 wsk qj Hsk qj .

(12)

s=q k=1

The spatial structure or symmetry required for the existence of this solution is determined by the set of coupling strengths {wpi qj , p, q = 0, . . . , n − 1, i = 1, . . . , np , j = 1, . . . , nq } that satisfy (12). Solutions of (12) can be

Y.-X. Li / Physica D 180 (2003) 210–234

217

divided into two general types: (i) model-independent solutions that exist for all possible 1-periodic H (·) functions, and (ii) model-dependent solutions that exist only for some 1-periodic H (·) functions with certain special properties. As shown in [5,6], this classification can be quite subtle since there could be solutions that exist for all possible 1-periodic H (·) functions, but their existence conditions have explicit dependence on H (·) and the coupling parameters. In the following sections of this paper except Section 7, we shall concentrate on studying the model-independent solutions. To determine the stability of this state, we introduce perturbations, δpi qj , to all the phase differences, so that φpi qj = Φpi qj + δpi qj . Linearizing (11) in δpi qj , we obtain the linear system:

δ˙pi qj =

 np
k=1

intra-cluster interactions



wpk pi Hp k pi δpk pi −

nq


k=1



wqk qj Hq k qj δqk qj +

 ns n−1

s=p k=1

intra-cluster interactions



wsk pi Hsk pi δsk pi −

ns n−1

s=q k=1

 wsk qj Hsk qj δsk qj , (13)

where H  (φpi qj ) is the derivative of H evaluated at φpi qj and is abbreviated by Hp i qj . 4. The rule of inverse cluster-size and the existence of n clusters For a network of N neurons, there exists a total of N (N − 1)/2 phase differences of which only N − 1 are independent. The reason is easy to understand if we number the neurons from 1 to N and pick the phase differences between nearest neighbors only: φi(i+1) = θi − θi+1 (i = 1, . . . , N − 1). Thus, there are a total of N − 1 phase differences between nearest neighbors. Any phase difference between neurons that are not nearest neighbors, e.g. φi(i+3) , can be obtained by propagating the phase differences between the nearest neighbors that form a path connecting the two, φi(i+3) = φi(i+1) + φ(i+1)(i+2) + φ(i+2)(i+3) . Thus, we only need to study the equations of N − 1-independent phase differences. This is because if you place N neurons on a unit ring, the sum of all phase differences must be one, i.e. N


φi(i+1) = 1 (where N + 1 = 1).

(14)

i=1

However, we will study the equations for N phase differences instead. This, as we shall demonstrate later, will make the analysis easier because of the symmetry in the matrix problem of the linearized system. The choice of the N phase differences is crucial for the symmetry of the related matrix problem. Our experience in studying the stability of small networks with up to six clusters [5] leads to the following way of choosing these variables that retains the required symmetry for later analysis. First, we pick the phase differences between different clusters. Since there are many neurons in each cluster, we need to pick a reference neuron in each cluster, e.g. the neuron that is numbered the first in each cluster. Now the phase difference between two clusters is defined as the phase difference between the first neurons in each cluster. Phase difference between two clusters that are nearest neighbors is φ(p+1)1 p1 = θ(p+1)1 − θp1

(p = 0, 1, . . . , n − 1),

(15)

where φ(p+1)1 p1 is the phase difference between the first neurons of the (p + 1)th and the pth clusters. These phase differences will be called the inter-cluster phase differences. There is a total of n inter-cluster phase differences between nearest neighbors. Only n − 1 are independent since one of them can always be expressed in terms of the other n − 1 phase differences due to the constraint specified by (14). The reason for including the nth linearly dependent phase difference is to retain the symmetry of the related matrix problem.

218

Y.-X. Li / Physica D 180 (2003) 210–234

Then, we pick the phase differences between neurons within each cluster. Since we have already picked the first neuron in each cluster as the reference, we continue to use it as a reference in defining the intra-cluster phase differences. Thus, we pick the phase difference between each neuron and the reference neuron in each cluster as the intra-cluster phase difference: φpj p1 = θpj − θp1

(j = 2, 3, . . . , np ).

(16)

It is easy to check that there are np − 1 intra-cluster phase differences in the pth cluster, thus the total number of intra-cluster phase differences is n−1


(np − 1) =

p=0

n−1


np −

p=0

n−1


1 = N − n.

(17)

p=0

These N − n phase variables are all independent. Differential equations for the N phase differences that have been selected are as follows. For the N − n intra-cluster phase differences (j = 2, 3, . . . , np , p = 0, 1, . . . , n − 1): φ˙ pj p1 =

np


[wpk pj Hpk pj − wpk p1 Hpk p1 ] +

ns n−1



(wsk pj Hsk pj − wsk p1 Hsk p1 ).

(18)

s=p k=1

k=1

For the n inter-cluster phase differences (p = 0, 1, . . . , n − 1): φ˙ (p+1)1 p1 =

n(p+1)




w(p+1)k (p+1)1 H(p+1)k (p+1)1 −

k=1

wpk p1 Hpk p1

k=1

ns n−1



+

np


wsk (p+1)1 Hsk (p+1)1 −

ns n−1



wsk p1 Hsk p1 .

(19)

s=p k=1

s=p+1 k=1

Note that whenever the subscript p + 1 for the clusters is larger than n − 1, we take mod (p + 1, n) instead owing to the cyclic nature of the indices. Thus, the index n is identical to the index 0. All other phase differences can be expressed as linear combinations of the N selected phase differences. Now, we study the n-cluster state in which neurons in each cluster are synchronized and maintain a finite phase difference with neurons in a different cluster. Thus, φpj p1 = 0 for all p and j , and φ(p+1)1 p1 = Φ(p+1)p (p = 0, 1, . . . , n − 1) is the fixed phase difference between the (p + 1)th and pth clusters. Substituting this steady state into (18) and (19), we obtain 0 = H (0)

0 = H (0)

np


n−1


[wpk pj − wpk p1 ] +

H (Φsp )

k=1

s=p

n(p+1)


np


k=1

w(p+1)k (p+1)1 −

k=1

ns


[wsk pj − wsk p1 ],

(20)

k=1

wpk p1 +

n−1
s=p+1

H (Φs(p+1) )

ns
k=1

wsk (p+1)1 −

n−1
s=p

H (Φsp )

ns


wsk p1 ,

k=1

(21) where wpi pi = 0 for all p and i, since self-coupling is excluded in the networks we study here. For an arbitrary choice of inter-cluster phase differences, Φpq (p = q), these conditions cannot be satisfied in a way that is independent of the function H , which is determined by specific model and coupling properties. To study phase-locked states that are model-independent, we study a case in which the phases of the n clusters are

Y.-X. Li / Physica D 180 (2003) 210–234

219

distributed evenly over one period. We call it the state of n splay-phase clusters. In this case, Φpq = (p − q)/n or Φ(p+1)p = 1/n. Therefore, (20) and (21) become 0 = H0

np


[wpk pj − wpk p1 ] +

n−1


k=1

0 = H0

n(p+1)


n(p+s)

Hs

s=1

w(p+1)k (p+1)1 −

k=1

np





[w(p+s)k pj − w(p+s)k p1 ],

k=1



wpk p1 +

n−1


Hs

n(p+s+1)


s=1

k=1

(22)

n(p+s)

w(p+s+1)k (p+1)1 −

k=1




w(p+s)k p1 ,

(23)

k=1

where H0 = H (0) and Hs = H (s/n). Therefore, the state of n splay-phase clusters is a model-independent solution if the following conditions are satisfied: (a)

np


wpk pj =

np


wpk p1 for all p and j. k=1

np Note that Ipj = k=1 wpk pj is the sum of all intra-cluster inputs to the j th neuron in the pth cluster. Thus, Ipj = Ip1 = Ip (for all p) requires that, within each cluster, the sum of all intra-cluster inputs to each neuron k=1

be equal. This condition can be satisfied if the coupling strength in each cluster is uniform, i.e. wpk pj = wpp for all p, k, j . In this case, Ip = (np − 1)wpp after excluding self-coupling. n(p+1)

(b)




w(p+1)k (p+1)1 =

k=1

np


wpk p1

for all p, i.e. I(p+1)1 = Ip1 .

k=1

Combined with condition (a), this implies that Ip+1 = Ip = wI (for all p), where wI denotes the sum of all intra-cluster coupling strengths. Thus, conditions (a) and (b) together require that the sum of all intra-cluster inputs each neuron receives be equal in all clusters. In the simple case in which the coupling strengths between neurons within each cluster is uniform, this condition implies wpp =

wI np − 1

for all p,

(24)

where w I = Ip (for all p) can also be considered as the total intra-cluster coupling strength. n(p+s)

(c)




n(p+s)

w(p+s)k pj =

k=1




w(p+s)k p1 for all s, p, and j.

k=1

nq wqk pj is actually the sum of all inter-cluster inputs from the qth Note that for all q = p, Oqpj = k=1 cluster to the j th neuron in the pth cluster. Thus, O(p+s)pj = O(p+s)p1 = O(p+s)p requires that the sum of all inter-cluster inputs from one cluster to each neuron in another cluster be equal. This condition can be satisfied if the inter-cluster coupling strength between neurons in two different clusters is uniform, i.e. wqk pj = wqp (for all q, p). n(p+s+1)

(d)


k=1

n(p+s)

w(p+s+1)k (p+1)1 =




w(p+s)k p1 for all s and p, i.e. O(p+s+1)(p+1)1 = O(p+s)p1 .

k=1

The simplest way to satisfy (c) and (d) is to require that for all p, O(p+1)p = w O . w O is the sum of all inter-cluster inputs from one cluster to each neuron in another cluster. This condition requires that the strength of the summed inter-cluster inputs from any cluster to any neuron in another cluster be equal. When the inter-cluster coupling strengths from neurons in one cluster to all neurons in another are identical, this condition

220

Y.-X. Li / Physica D 180 (2003) 210–234

becomes wqp =

wO nq

for all q = p,

(25)

where wO will be referred to as the total inter-cluster coupling strength. It is obvious that asymmetrical coupling is required if nq = np because wqp = wO /nq = wO /np = wpq . Therefore, the state of n splay-phase clusters is a solution if the total intra-cluster strength is identical in all clusters, and if the total inter-cluster strength from one cluster to another is identical between all clusters. There are multiple ways these conditions can be satisfied. Here is one simple, model-independent way to make it occur. We only need to make the output coupling strengths uniform within each cluster with a magnitude that is inversely proportional to the size of the cluster. Let wpq be the strength of synaptic coupling from all neurons in the pth cluster to neurons in the qth cluster, then the state of n splay-phase clusters exists if  I w   if p = q,   np wpq = (26)  wO    if p = q, np where wI and w O are, respectively, the total intra-cluster and the total inter-cluster coupling strengths. Note that when p = q, wI /np instead of wI /(np − 1) is used. The reason will become clear in (38). This result leads to the following rule in determining the existence of the n-cluster state. Theorem 4.1 (The rule of inverse cluster-size). The existence of the state of n splay-phase clusters is guaranteed in a network if neurons are organized into n clusters such that the strength of synaptic projections from each neuron is inversely proportional to the size of the cluster to which it belongs. This rule actually requires that neurons in a larger cluster project with weaker synaptic strength to other neurons. In other wards, if one neuron belongs to a larger cluster, then it has to work less harder to maintain the clustered state. This rule only provides a sufficient but not necessary condition for the existence of the state of n splay-phase clusters. The fact that if np = nq then wpq = wqp implies that the coupling strength is not symmetric between neurons in two clusters of different sizes. Thus, this rule suggests that the asymmetry in the sizes of the clusters can be compensated by a counter-balancing asymmetry in the coupling strengths. For a network that obeys the rule of inverse cluster-size, the N phase difference equations previously described by (9) now becomes inter -cluster interactions intra-cluster interactions       np nq n−1 I O


dθpi w w = H (θpj − θpi ) + H (θqj − θpi ), dt np nq j =1

q=0,q=p

(27)

j =1

where p = 0, 1, . . . , n − 1, i = 1, 2, . . . , np . For any phase-clustered state given by Θp = f (n) t + Φp0 ,

(28)

where p = 0, 1, . . . , n − 1, and Φp0 is the locked phase difference between the pth cluster and the 0th cluster. Note that Φpp = 0 for all p. Substituting (28) into (27), we can calculate the change in the frequency resulting from the

Y.-X. Li / Physica D 180 (2003) 210–234

221

phase interactions between neurons in different clusters: f (n) =

n−1


n0 − 1 I w H (0) + w O H (Φq0 ), n0

(29)

q=1

where self-coupling is excluded. If self-coupling is included or if n0 1, f (n) is can be expressed by f (n) =

n−1


wk H (Φk0 ),

(30)

k=0

where wk = wI for k = 0 and wk = wO for k = 0. Therefore, f (n) is the weighted average of the phase-resetting function H (·) at all the phase differences that appear in this phase-locked state. This result is reduced to (3) if wk = w/n and Φk0 = k/n for all k. By applying the rule of inverse cluster-size to (18) and (19), we obtain the N equations that describe the dynamics of the N selected phase difference variables to (18) and (19):   np np ns n−1 O


w I 
w φ˙ pj p1 = Hpk pj − Hpk p1  + [Hsk pj − Hsk p1 ], (31) np ns k=j

φ˙ (p+1)1 p1 =

s=p

k=2

n(p+1)




wI

H(p+1)k (p+1)1 −

n(p+1) k=2  n−1
wO + n(p+s+1)

k=1

np wI
Hpk p1 np k=2

n(p+s+1)




s=1

H(p+s+1)k (p+1)1 −

k=1

wO n(p+s)



n(p+s)




H(p+s)k p1 ,

(32)

k=1

where p = 0, 1, . . . , n − 1, and j = 2, 3, . . . , np . 5. The stability of the n-cluster state To determine the stability of the state of n splay-phase clusters, Φpi qj = (p −q)/n (for all p, q, i, j ), we linearize the phase difference equations near this state. Let δpi qj = φpi qj − Φpi qj be the perturbation to the n-cluster state. The linearized equations for the time evolution of these perturbations are   np np (p+s) n−1 O  n


wI H0
w Hs  δ˙pj p1 = δpk pj − δpk p1 ) + [δ(p+s)k pj − δ(p+s)k p1 ], (33) np n(p+s) k=j

 δ˙(p+1)1 p1 = H0

wI

n(p+1)  n−1
+ Hs s=1

k=2

s=1

n(p+1)




δ(p+1)k (p+1)1

k=2

wO n(p+s+1)

k=1

np wI
− δpk p1 np



k=2

n(p+s+1)


k=1

δ(p+s+1)k (p+1)1 −

wO n(p+s)

n(p+s)




δ(p+s)k p1 ,

(34)

k=1

where j = 2, 3, . . . , np , p = 0, 1, . . . , n − 1, and H0 = H  (0) and Hs = H  (s/n). This is a system of linear equations that can be written in the matrix form: δ˙ = J δ,

(35)

222

Y.-X. Li / Physica D 180 (2003) 210–234

where δ T = (δ02 01 , . . . , δn1 (n−1)1 ) is an N-dimensional vector, and J is an N × N matrix that is the Jacobian of (33) and (34). 5.1. Stability of intra-cluster phase differences Lemma 5.1. The first N − n rows of J that correspond to perturbations of the N − n intra-cluster phase differences is diagonal with identical diagonal entries of multiplicity N − n:        1 2 n−1 λ = −w∗ = −w O sH (0) + H  + H + ··· + H , (36) n n n where s = wI /w O . The proof is straightforward. Since δij = −δji and δik + δkj = δij by definition: n(p+s)




n(p+s)

[δ(p+s)k pj − δ(p+s)k p1 ] = −

k=1 np
k=j

δpk pj −




[δpj (p+s)k + δ(p+s)k p1 ] = −n(p+s) δpj p1 ,

(37)

k=1 np


δpk p1 = −[δpj p1 + δpj p1 ] −

k=2

np


[δpj pk + δpk p1 ] = −np δpj p1 .

(38)

k=2,k=j

It is interesting to note that in (38), we obtain np times δpj p1 when we combine two sums each includes only np − 1 terms. Substituting these equalities into (33), we obtain  n−1
I  O  δ˙pj p1 = − w H0 + w (39) Hk δpj p1 . k=1

Remark 5.1. Lemma 5.1 suggests that the linearized time evolution of the perturbation to each intra-cluster phase difference is independent of the time evolution of perturbations of all other phase differences. Remark 5.2. Although the eigenvalues λ given by (36) and λ(n) given by (4) look very similar, there is an important difference between the two. In λ(n) , the coupling strength w plays no role in determining the sign of λ(n) . In λ, the ratio between the two coupling strengths s = wI /w O plays an important role in determining the sign of λ. Introducing the coupling parameters into the stability condition has important consequences as has been shown in small networks of coupled phase oscillators [5,6]. Lemma 5.1 shows that this can generally be achieved in a large network of N neurons with heterogeneous and asymmetrical coupling strengths, and for an arbitrary number of clusters. This is one of the most important reasons for introducing heterogeneity in the coupling strengths. Remark 5.3. Since H  is a 1-periodic function, H  ((n − k)/n) = H  (−k/n). Thus, H  (k/n) + H  ((n − k)/n) = H  (k/n) + H  (−k/n) = 2g  (k/n) which implies that   k = Wintra + Winter n k=1      m
O   k  1 sg (0) + 2 =w g + mod (n − 1, 2)g , n 2

w∗ = w I H  (0) + w O

n−1


k=1

H

(40)

Y.-X. Li / Physica D 180 (2003) 210–234

223

n−1 O

 where m = [n − 1 − mod (n − 1, 2)]/2. Wintra = w I H  (0) and Winter = w k=1 H (k/n) are, respectively, the    net intra- and inter-cluster phase interactions. The identities H (0) = g (0) and H (1/2) = g  (1/2) were used.

Lemma 5.1 implies that w∗ > 0 is a necessary condition for the stability of the n-cluster state. It has been shown previously [5,6,25] that the sign of g  (k/n) determines whether the interaction between two neurons at a phase difference k/n is attractive (if g  (k/n) > 0) or repulsive (if g  (k/n) < 0). Thus, w ∗ > 0 means that the summarized, net influence on each neuron from neurons in all clusters, including influences from both intra- and inter-cluster phase interactions, is attractive. In networks of repulsively coupled neurons, the intra-cluster phase interaction is repulsive, i.e. Wintra < 0 since H  (0) = g  (0) < 0. However, if the net inter-cluster interaction is attractive or stabilizing, i.e. Winter > 0, w ∗ > 0 can still be achieved. As you can see in (40), the explicit dependence of w ∗ on s makes it much easier to achieve such a situation. This result was first obtained in small networks containing four or six neurons [5,6]. Here, Lemma 5.1 shows that it actually applies to much more general cases. 5.2. Stability of inter-cluster phase differences Lemma 5.2. The n × n block matrix, Jn , located at the lower-right corner of J is given below. This matrix is identical to the Jacobian at the splay-phase state, Js , of the phase difference equations of an n-neuron network with uniform and symmetric coupling strength:      −Σ H1 H2 · · · Hn−3 Hn−2 Hn−1       H  H1 · · · Hn−4 Hn−3 Hn−2  n−1 −Σ          Hn−2 Hn−1  −Σ · · · H H H n−4 n−3  n−5   O (41) Jn = Js = w  · · · ··· ··· ··· ··· ··· ··· ,     H4 H5 · · · −Σ H1 H2   H3    H  H3 H4 · · · Hn−1 −Σ H1   2       H1 H2 H3 · · · Hn−2 Hn−1 −Σ

 where Σ = n−1 s=1 Hs . The matrix Jn can be obtained by simplifying (34) into the following form:  n(p+s+1) n(p+s) n−1


1 1  ˙δ(p+1)1 p1 = ws Hs δ(p+s+1)k (p+s+1)1 − δ(p+s)k (p+s)1 n(p+s+1) n(p+s) s=0 k=2 k=2  n−1 n−1

O  w Hs δ(p+1)1 p1 + w O Hs δ(p+s+1)1 (p+s)1 , − s=1

(42)

s=1

where ws = wI for s = 0 and ws = wO for s ≥ 1. To derive this formula, we used the following identities: δ(p+s)k p1 = δ(p+s)k (p+s)1 + δ(p+s)1 p1 ,

(43)

δ(p+s+1)k (p+1)1 = δ(p+s+1)k (p+s+1)1 + δ(p+s+1)1 (p+1)1 ,

(44)

δ(p+s+1)1 (p+1)1 − δ(p+s)1 p1 = δ(p+s+1)1 (p+s)1 − δ(p+1)1 p1 .

(45)

The proof that Jn = Js is given in Appendix A.

224

Y.-X. Li / Physica D 180 (2003) 210–234

Remark 5.4. Lemmas 5.1 and 5.2 together illustrate that the N × N Jacobian matrix J in (35) has the following form:   Λ(N−n) 0(N−n)×n , J = (46) Cn×(N−n) Jn where Λ(N−n) is an (N − n) × (N − n) diagonal matrix with identical diagonal entry λ = −w ∗ , 0(N−n)×n an (N − n) × n zero matrix, Cn×(N−n) a nonzero n × (N − n) matrix, and Jn the n × n matrix given by (41). Remark 5.5. Lemma 5.2 not only gives the matrix that determines the conditions for the stability of the inter-cluster phase differences but also shows that these stability conditions are identical to those for the stability of the splay-phase state in an n-neuron network. Remark 5.6. We were ignorant of the results in [1] when Lemma 5.2 was obtained. However, we need to point out that Lemma 5.2 is identical to the result obtained by Okuda [1] about a decade ago. The only difference is that we showed his results are actually applicable to more general cases considered in this study. 5.3. Stability of the splay-phase state in an n-neuron network Lemma 5.2 reduces the problem of determining the stability of the state of n splay-phase clusters partly into that of determining the stability of the splay-phase state in an n-neuron network. Therefore, we first need to solve the stability of the splay-phase state. The analysis becomes exactly the same as was already done by Okuda in [1]. Here, we simply present the result in the form of a theorem while leaving the details of the proof in Appendix A. Theorem 5.1. The stability of a splay-phase state in a network of n identical oscillators that are weakly coupled with uniform coupling strength is determined by the following n eigenvalues: λl = −w O

n−1
k=1

Hk [1 − ei(2π/n)kl ],

l = 0, 1, . . . , n − 1,

(47)

where λ0 = 0 for l = 0 was expected since only n − 1 rows of Jn are independent. This leads to the following proposition. Proposition 5.1. The real and the imaginary parts of all the eigenvalues involve only the values of g  (φ) and h (φ), respectively, at all the possible phase differences that appear in the splay-phase state, i.e.      m
Re{λl } 2π  k  1 = − 2g 1 − cos kl − mod (n − 1, 2)g [1 − cos lπ ], (48) wO n n 2 k=1

m

Im{λl }
 = 2h wO k=1

  k 2π sin kl, n n

(49)

where l = 1, . . . , nh with nh = m + mod (n − 1, 2) and m = [n − 1 − mod (n − 1, 2)]/2. Thus, m = (n − 2)/2 if n is even and m = (n − 1)/2 if n is odd. g(φ) and h(φ) are, respectively, the odd and even parts of H (φ).  = H . The proof of this proposition is straightforward. Note that H  (φ) = g  (φ)+h (φ) is 1-periodic, thus H−k n−k Now, g  (φ) = (1/2)[H  (φ) + H  (−φ)] is even and h (φ) = (1/2)[H  (φ) − H  (−φ)] is odd. On the other hand,

Y.-X. Li / Physica D 180 (2003) 210–234

225

ei(2π/n)kl = e−i(2π/n)(n−k)l . Note also that h (1/2) = 0. Substituting these identities into (47), we obtain (48) and (49). 5.4. Stability of the state of n splay-phase clusters Now, conditions for the stability of the state of n splay-phase clusters are identical to those for the stability of the splay-phase state in a network of n “big” neurons except that now we need an extra condition to guarantee the stability of intra-cluster phase differences. The stability of the splay-phase state is determined by nh conditions given in Proposition 5.1. Thus, the stability of the state of n splay-phase clusters is determined by nh + 1 conditions. These conditions are listed in Table 1 for n = 2, . . . , 7. These results are summarized in the following theorem. Theorem 5.2. Let us consider a network of N identical, weakly coupled neurons. Suppose that the network is structured and that the rule of inverse cluster-size is satisfied. At the state of n splay-phase clusters, this network behaves like n giant neurons with evenly distributed phases. The stability conditions are w∗ = wI H  (0) + w O

n−1
s=1

H

s  n

> 0,

    m −Re{λl }
  s  2π  1 = 2g 1 − cos kl + mod (n − 1, 2)g [1 − cos lπ ] > 0, wO n n 2

(50)

l = 1, . . . , nh ,

k=1

(51) where nh = m + mod (n − 1, 2) and m = [n − 1 − mod (n − 1, 2)]/2. The conditions given in Theorem 5.2 have been verified for n = 2, 3, 4, 5, 6, 8 by calculating the eigenvalues directly from the Jacobian matrices of smaller networks involving 2–8 neurons [5]. Extensive numerical simulations have been carried out to confirm these results in smaller networks of realistic model neurons that are coupled through realistic models of synapses [5,6]. Table 1 shows that as the number of clusters increases, the number of possible ways to satisfy the stability conditions increases dramatically. For n = 2, 3, the sufficient condition is also necessary. But for n > 3, stability can occur at weaker conditions than the sufficient condition. For n = 7, there are many possible ways to satisfy the conditions for stability.

6. The stability of the asynchronous state Stability of the asynchronous state in neural networks has been an interesting topics since the firing rate models work best if the population of neurons settles quickly to an asynchronous state [28]. Stability conditions of asynchronous states have been investigated in a number of previous works [29–31] using mean field analysis. Results derived in the previous section can be readily applied to the stability of the asynchronous states in a phase model. Such a state corresponds to the splay-phase state of a network of n neurons when n → ∞. For this purpose, it is sufficient to study the eigenvalue spectrum of the linearized matrix as n → ∞. Note that the coupling strength in (47) must scale with the size of the network n, thus wO = WO /n. As n → ∞, (47) turns into the following integral:    1 n−1
k 1 H H  (θ )[1 − ei2πlθ ] dθ. (52) [1 − ei2πl(k/n) ] = −WO λ(l) = −WO lim n→∞ n n 0 k=1

226

Y.-X. Li / Physica D 180 (2003) 210–234

Table 1 Conditions for the stability of the state of n splay-phase clustersa n

nh

Conditions for the stability of inter-cluster phase differences

2

1

g

3

1

4

2

5

2

6

3

7

3

  1 >0 2   1 g >0 3       1 1 1 > 0, g  + g >0 g 4 4 2     1 2 + g > 0, 5 5   √     1 2 3+ 5 g + ag >0 a= = 2.618 5 5 2 ag

        1 1 1 1 2g  + g > 0, g  + g > 0, 6 2 6 3       1 1 1 g + 3g  + 2g  >0 6 3 2

      1 2 3 + bg + cg > 0, 7 7 7       1 2 3 bg + cg + g > 0, 7 7 7       1 2 3 cg + g + bg >0 7 7 7   1 − cos (4π/7) 1 − cos (6π/7) b= = 3.247, c = = 5.049 1 − cos (2π/7) 1 − cos (2π/7) g

Stabilizing phase difference(s) (additional conditions if required are listed in brackets) 1 (i) 2 1 (i) 3        1 1  1 (i) g > g  4 4 2  1 1 (ii) & 4 2        2  1  1 g > a g  (i) 5 5 5         2 2 1  (ii) g > a g  5 5 5  2 1 (iii) & 5 5             1  1 1   1 g  g (i) > 3 g  + 2  6 6 3  2             1 1 1 1   1  1 (ii) & 2g  > g  , 3g > g 6 3 6 2  3 6                1 1 1 1 1 1  , g  (iii) & g > g  > g  6 2 6 3  2 3  1 1 1 (iv) & & 6 3 2            1 2  3  1 g > b g  + c g  (i)  7 7 7 7             2 1  3  2 g > c g  + b g   7 7 7 7              2  3 3 1  g  (iii) g > b g  + c  7 7 7  7           2 1 2 3  1 g + bg > c g  (iv) & 7 7 7 7 7           1 3 2  3 1 bg + g > c g  (v) & 7 7 7 7 7           2 2 3 1  3 (vi) & g + bg > c g  7 7 7 7 7  1 2 3 (vii) & & 7 7 7 (ii)

a The conditions listed in this table are identical to those for the stability of the splay-phase state in a network of n neurons. The additional

 condition that is not listed in the table is w ∗ = w I H0 + w O n−1 s=1 Hs > 0, which guarantees the stability of intra-cluster phase differences. In the last column, a stabilizing phase difference is defined as a phase difference, φ, such that g  (φ) > 0. This implies that interaction between neurons at a phase difference φ is stabilizing. Distinct ways to satisfy the stability conditions are numbered. Note that g  (φ) is an even, 1-periodic function, thus g  (φ) = g  (−φ) = g  (1 − φ). For example, if g  (1/5) > 0, then g  (4/5) = g  (−1/5) = g  (1/5) > 0.

Y.-X. Li / Physica D 180 (2003) 210–234

227

Since H (θ ) is a 1-periodic smooth function, we can express it as a Fourier series ∞

a0
[al cos 2πlθ + bl sin 2π lθ ], + 2

H (θ) =

(53)

l=1

 al = 2

H (θ ) cos 2πlθ dθ,

l = 0, 1, . . . ,

(54)

H (θ ) sin 2πlθ dθ,

l = 1, 2, . . . .

(55)

0

 bl = 2

1

1

0

Its derivative is H  (θ ) = 2π

∞


m[−al sin 2πlθ + bl cos 2π lθ ] =

l=1

∞


[αl sin 2π lθ + βl cos 2π lθ ],

(56)

l=1

where αl = −2π lal and βl = 2π lbl . H  (θ ) is also a 1-periodic function but it does not contain a constant term. Substitute into (52), we obtain    1 βl αl π WO l H  (θ )[ cos 2πlθ + i sin 2π lθ ] dθ = WO (57) +i = [bl − ial ]. λ(l) = WO 2 2 2 0 Therefore, the stability of the asynchronous state is determined by the signs of the Fourier coefficients of the odd part of the function H (θ ), i.e. bl . If H is even (i.e. H  is odd), then the asynchronous state is neutrally stable (unstable) since the real parts of all eigenvalues are zero. If any one of the Fourier coefficients bl is positive, then the asynchronous state is unstable. Requiring that bl < 0 for all l is equivalent to requiring that each Fourier mode of the odd part of the H function (i.e. g) stabilizes anti-phase interactions. Even in case this is satisfied, the stability is usually fragile anyway as n → ∞. This is because if bl = O(l −γ ) and γ > 1, then Re{λ(l)} =

πWO lbl = O(l −(γ −1) ) → 0, 2

as l → ∞.

(58)

Therefore, there will be infinitely many eigenvalues with infinitely small real parts. When this occurs, the splay-phase state becomes irregularly asynchronous in the presence of noise. For the integrate-and-fire model, it was shown that bl = O(1/ l 3 ) [32]. For most neuronal models with synaptic coupling, the H function is continuous and smoother which usually results in a more rapid convergence of bl to zero as l → ∞.

7. Clustered states with phases that are not evenly distributed Results in previous sections are focused on clustered states with evenly distributed phases, called splay-phase states. Splay-phase relation between different clusters imposes a cyclic symmetry in the clustered solution. With this symmetry, we showed that the stability of these states is determined only by the slopes of the odd part of the phase-resetting function, i.e. g  (ψ) = [H  (ψ) + H  (−ψ)]/2, sampled evenly in the phase space at n different phases: ψk = k/n (k = 0, 1, . . . , n−1). However, for clustered states with phases that are not evenly distributed, the stability depends on slopes of both the odd and the even parts of the function H (·) [5,6]. To make the analysis simple to explain, we shall focus only on the case for n = 2 when we extend results on clustered states with nonuniformly distributed phases previously obtained in small networks [5,6] to large networks with N neurons and n clusters. Thus, we shall only study the state of two clusters phase-locked at an arbitrary phase difference ψ other than 0 and 1/2.

228

Y.-X. Li / Physica D 180 (2003) 210–234

It is known for a long time that every zero of the function g(ψ) gives rise to a phase-locked solution in two symmetrically coupled neurons [25]. This result can be extended to states of two clusters at an arbitrary phase difference ψ in networks of four and six symmetrically coupled neurons [5,6]. It was also shown that, in networks of generic limit cycle oscillators with a g(ψ) that has only two zeros (one at ψ = 0 and the other at ψ = 1/2), such a state can still occur with appropriately structured heterogeneity in coupling strengths. Such a state is called a weakly model-dependent state. This is because the model-dependent function H (ψ) appears in the condition that determines its existence, but this dependence does not prevent a network from satisfying this condition for any generic limit cycle oscillator. It only makes the exact combinations of coupling strengths in this condition differ for different models with different forms of H (·) function. Here, we show how to extend this result to large networks of N neurons with heterogeneous and asymmetrical coupling strengths. For n = 2, the sizes of the two clusters, n0 and n1 , are related by n1 = N − n0 .

(59)

Since now we are dealing with two clusters only, the phase difference equations (18) and (19) are simplified to the following system of equations: φ˙ pj p1 =

np


n(p+1)

[wpk pj Hpk pj − wpk p1 Hpk p1 ] +

k=1

φ˙ 11 01 =

n1





[w(p+1)k pj H(p+1)k pj − w(p+1)k p1 H(p+1)k p1 ],

(60)

k=1

w1k 11 H1k 11 −

k=1

n0


w0k 01 H0k 01 +

k=1

n0


w0k 11 H0k 11 −

k=1

n1


w1k 01 H1k 01 ,

(61)

k=1

where p = 0, 1 (p+1 = mod (p+1, 2)), and j = 2, 3, . . . , np . Note that there is only one independent inter-cluster phase difference since n = 2. Thus, φ11 01 = φ10 = −φ01 . The existence of such a state of two clusters phase-locked at φ10 = ψ or φ01 = −ψ = 1 − ψ is determined by 0 = H (0)

np


n(p+1)

[wpk pj − wpk p1 ] + H(p+1)p

k=1

0 = H (0)

n 1
k=1

w1k 11 −

n0
k=1






[w(p+1)k pj − w(p+1)k p1 ],

(62)

k=1

w0k 01 + H (−ψ)

n0


w0k 11 − H (ψ)

k=1

n1


w1k 01 ,

(63)

k=1

where φ00 = φ11 = 0 was used, and H(p+1)p = H (φ(p+1)p ). If we assume that wpq specifies the coupling strength from each neuron in cluster p to any neuron in cluster q (p, q = 0, 1), equations (62) and (63) are satisfied by the following choice of these parameters:   wI   if p = q,   np wpq = (64) O  w  pq   if p = q.  np The rule of inverse cluster-size still holds but the inter-cluster coupling strengths between the two clusters can be asymmetrical and satisfy O O H (ψ) = w01 H (−ψ) w10

or g(ψ) = (s − 1)H (−ψ),

O /w O is the ratio between the two asymmetrical coupling strengths. where s = w01 10

(65)

Y.-X. Li / Physica D 180 (2003) 210–234

229

Stability of such a clustered state can be obtained by solving the linearized system of (60) and (61) at this two-cluster state using the coupling parameters given by (64). The linearized system can be obtained in a similar way as we derived Eqs. (33) and (34) in Section 5, it reads O δ˙pj p1 = −[w I H  (0) + w(p+1)p H  (φ(p+1)p )]δpj p1 ,

(66)

O  O  δ˙11 01 = −[w10 H (ψ) + w01 H (−ψ)]δ11 01 + off diagonal terms,

(67)

where p = 0, 1 (p + 1 = mod (p + 1, 2)), and j = 2, 3, . . . , np . Thus, we obtain three distinct eigenvalues O  H (ψ)] (n0 − 1 multiplicity), λ1 = −[w I H  (0) + w10

(68)

O  H (−ψ)] λ2 = −[w I H  (0) + w01

(69)

(n1 − 1 multiplicity),

O  O  λ3 = −[w10 H (ψ) + w01 H (−ψ)].

(70)

7.1. Symmetrical inter-cluster coupling Let us consider the special case when the inter-cluster coupling is symmetrical, i.e. if O O = w01 w O = w10

(i.e. s = 1),

(71)

the existence condition (65) becomes g(ψ) = 0.

(72)

This is the same condition for two coupled oscillators to oscillate at a locked phase difference ψ [5,6,25]. In other words, the behavior of two mutually coupled clusters of synchronized neurons is indeed similar to that of just two coupled oscillators. The stability of such a state is determined by the three eigenvalues λ1 = −[w I H  (0) + w O H  (ψ)] (n0 − 1 multiplicity),

(73)

λ2 = −[w I H  (0) + w O H  (−ψ)] (n1 − 1 multiplicity),

(74)

λ3 = −2w O g  (ψ),

(75)

which are identical to results obtained in [5]. Thus, for each zero of g(ψ), there exists a phase-locked state. In general, g(0) = g(1/2) = 0 for all coupled oscillators since g(ψ) is an odd function of period 1. Thus, the fully synchronous solution (ψ = 0) and the state of two anti-phase clusters (ψ = 1/2) are always solutions of such a symmetrical network. But these are simply special cases of the splay-phase clusters that we studied in previous sections. For example, for ψ = 0, the stability conditions are reduced to g  (0) < 0, which simply requires that in-phase interaction between any pair of neurons be attractive or synchronizing. This the well-known result that we already know since the pioneer work by Kuramoto [25] on the phase theory of coupled oscillators. For ψ = 1/2, λ1 = λ2 = −[w I H  (0) + wO H  (1/2)] = −[w I g  (0) + wO g  (1/2)]. This is because H  (0) = g  (0) and H  (1/2) = H  (−1/2) = g  (1/2). Also, λ3 = −2w O g  (1/2) = −2w O H  (1/2). This result is the same as obtained in [5,6] for small networks. Furthermore, if wI = w O = w, i.e. the coupling is uniformly distributed, these eigenvalues become identical to those in (5) and (6) obtained by Okuda.

230

Y.-X. Li / Physica D 180 (2003) 210–234

7.2. Asymmetrical inter-cluster coupling For any value of ψ that is not 0 or 1/2, g(ψ) = 0 is a strongly model-dependent property that only occurs under rather restricted conditions (see Fig. 9 in [5]). However, the condition given in (65) implies that a state of two clusters phase-locked at a value of ψ other than 0 and 1/2 is possible for a generic phase-resetting function H (·). According to (65), this can occur if H (ψ) and H (−ψ) have the same sign and the ratio between the coupling strengths is selected appropriately. Using the theory developed in this paper, the eigenvalues that determine the stability of such a state can also be obtained as is shown in (68)–(70). This result again suggests that the asymmetry in the distribution of the phase differences between different clusters can also be compensated by a counter-balancing asymmetry in the distribution of coupling strengths. Based on this counter-balancing mechanism, the theory developed here can, in theory, be applied to more general clustered states involving phases that are not evenly distributed. Actually, the existence and stability of other kinds of phase-locked states with nonuniform phase distribution previously obtained in [5,6] can also be obtained in an N -neuron network. However, a complete classification of all possible weakly model-dependent solutions requires the analysis of all possible symmetry groups of the network. Such a task is not only beyond the desired scope of this paper but also beyond the capability of the author at the present stage.

8. Summary and discussion This work extends several results obtained previously in networks of all-to-all coupled oscillators with uniform coupling strengths [1] and in small networks with heterogeneous coupling strengths [5,6] to an arbitrarily large network of N oscillators with heterogeneous and asymmetrical coupling strengths. This extension was motivated by a number of results recently obtained in the study of coupled oscillators. First, several numerical and analytical studies have shown that, when noise is present, stable clustered states can hardly be obtained in networks of repulsively coupled neural oscillators when the coupling strength is uniform [2,3,33]. It has been shown in both numerical and analytical studies that the stability of such clustered states can be easily achieved by introducing heterogeneity in the coupling strengths [4–6]. The most important feature of a network of coupled oscillators with heterogeneous coupling strengths is that the stability of phase-locked solutions is determined by the coupling strengths as well as other model-dependent properties [5,6]. For weakly model-dependent, phase-locked solutions, even the existence is dependent on the coupling parameters. This feature does not exist in networks of uniformly coupled oscillators in which the coupling strength plays no role in determining the existence and stability of phase-locked solutions, provided that its magnitude is small. This difference between networks with uniform and nonuniform coupling strengths has profound consequences in the study of phase-locked solutions in coupled oscillators. It solved a practical difficulty encountered by some people in the study of synchrony in coupled oscillators (see e.g. [2,3,33]). This difficulty is related to the coexistence and multi-stability in the phase-locked solutions of coupled oscillators. It is often found in numerical simulations that the desired pattern often coexists with some other solution(s) that should be avoided. But when the system is perturbed, transition from a state that is desired to another state that is not desired often occurs since the two solutions coexist and are stable under the same conditions. This is partly due to the fact that the existence of synchronized clustered states with certain spatial–temporal symmetries are related to the spatial symmetries in the coupling strengths of the network [14,15,34,35]. In a network of coupled phase oscillators, high degree of symmetry in the coupling strengths results in the coexistence of most physiologically important solutions such as the gait solutions [5,6]. Thus, finding a mechanism that avoids multi-stability between these coexisting solutions, a process called pattern selection or segregation [6], is important for the design of a functioning pattern-generating network. Knowing

Y.-X. Li / Physica D 180 (2003) 210–234

231

exactly how a change in coupling parameters induces a change in the stability of these solution is very important for achieving pattern selection in a highly symmetric network. Results on the existence and stability of phase-clustered states that incorporate the relative magnitudes of coupling strengths are of great importance in realizing pattern selection. Results on pattern selection can also be used in constructing a network of current-based, realistic neural oscillators that possesses collective computational properties like those of Hopfield’s networks of two-state neurons [36,37]. Second, phase-clustered states in a small network of four or six oscillators are often used in the study of stereo-typical gait patterns observed in four- or six-legged animals [13–16]. In some frequently observed gaits in quadrupeds and hexapods, phase locking between different numbers of synchronized legs can occur. Thus, in the gait called lurch in hexapods, the two front legs and the two hind legs are in synchrony but these four synchronized legs are at anti-phase to the other two synchronized legs in the middle. Therefore, the lurch gait involves two anti-phase clusters of different sizes, one includes four legs and the other includes two. As the number of legs increases, the number of possible clustered states that involve clusters of different sizes also increases. A theory of coupled phase oscillators that predicts conditions for the existence and stability of clustered states involving clusters of different sizes is of great interest in the study of network models of animal gaits. Analysis presented in this work showed that the asymmetry in the size of the clusters can be compensated by a counter-balancing asymmetry in the coupling strengths. This mechanism leads to the rule of inverse cluster-size which suggests that a state of n synchronous clusters of different sizes can exist if the strengths of all synaptic outputs from each neuron is inversely proportional to the size of the cluster to which it belongs. The proportionality constant can be different for the coupling within each cluster (w I ) and between different clusters (w O ). This condition is sufficient for the existence of n splay-phase clusters. The analysis also showed that in the state of n synchronized clusters the network is almost identical to a network of n coupled giant “oscillators”. The main difference is that the total synaptic inputs each neuron receives must be stabilizing, i.e. w∗ > 0. The introduction of heterogeneity in the coupling strengths makes this condition much easier to satisfy in networks of repulsively coupled neural oscillators. Once this condition, which governs the stability of intra-cluster phase disturbances, is satisfied, the remaining conditions for the stability of the n-cluster state are identical to those governing the stability of the splay-phase state in a network of n neurons. Results on clustered states involving nonuniform inter-cluster phase differences previously obtained in small networks [5,6] were also extended to an arbitrarily large network of N neurons. It was shown that a mechanism similar to that underlying the rule of inverse cluster-size can also apply to clustered states with noneven distribution of inter-cluster phase differences. In this case, an asymmetry in the distribution of phase differences can be compensated by another counter-balancing asymmetry in the inter-cluster coupling strengths. This result suggests that the theory developed in this work is likely applicable to even more general cases of clustered states that are weakly model-dependent and involve clusters of different sizes that are locked at phases that are not evenly spaced in the phase space. This generalized theory not only makes results previously obtained in small and uniformly coupled networks applicable to more general kinds of networks with heterogeneous and asymmetrical coupling strengths, but also makes the more general types of clustered states achievable in network of realistic neural oscillators.

Acknowledgements The author is grateful to Drs. Robert Miura and Yuqing Wang for numerous discussions. Robert Miura is also acknowledged for reading an earlier version of this manuscript critically. This work was financially supported by NSERC (Natural Sciences and Engineering Research Council of Canada) grants to YXL.

232

Y.-X. Li / Physica D 180 (2003) 210–234

Appendix A. Stability of the splay-phase state of an n-neuron network Consider the phase equations of a network of n identical neurons with weak, uniform, symmetric coupling strength wO : n−1


dθi = wO H (θj − θi ) (i = 1, . . . , n). dt

(A.1)

j =i

Equations for the phase differences, φij = θi − θj , are   n−1 n−1 n−1


dφij H (φli ) − H (φlj ) = w O H (φ(l+i)i ) − H (φ(l+j )j ) , = wO  dt l=i

l=j

(A.2)

l=1

where l + i = mod (l + i, n) if l + i > n − 1. Let Φij = (i − j )/n be the splay-phase state. Linearizing (A.2) at this state, we obtain the following equations for the phase perturbations, δij = φij − Φij : n−1


dδij = wO Hl [δ(l+i)i − δ(l+j )j ], dt

(A.3)

l=1

where Hl = H  (l/n), and by definition, δii = 0. Not all phase differences are independent. We will only study the n phase differences between nearest neighbors, φ(i+1)i for 0 ≤ i ≤ n − 1. Thus, we only need to study equations for δ(i+1)i : n−1

n−1



dδ(i+1)i = wO Hl [δ(l+i+1)(i+1) − δ(l+i)i ] = wO Hl [δ(l+i+1)(l+i) − δ(i+1)(l+i) − δ(l+i)i ] dt l=1 l=1  n−1  n−1 n−1


O  O   =w Hl [δ(l+i+1)(l+i) − δ(i+1)i ] = w − Hl δ(i+1)i + Hl δ(l+i+1)(l+i) . l=1

l=1

(A.4)

l=1

This gives a Jacobian matrix Js that is identical to the matrix Jn given by (41). To calculate the eigenvalues of the matrix Jn , we note that the sum of the entries in each row and each column is zero, thus 0 is always an eigenvalue of Jn with the corresponding eigenvector δ = (1, . . . , 1). This is because only n − 1 out of the n equations are independent. This eigenvalue corresponds to the tangent direction of the phase-space trajectory of this state and is related to translational invariance in time. Also, all of the eigenvalues   are located in a 

 ring of radius ρ(Jn ) centered at − , 0 in the complex plane with a spectral radius ρ(Jn ) =   if H  (s/n) > 0 or H  (s/n) < 0 for all s = 1, . . . , n − 1 (i.e. if they are either all positive or all negative). If H  (s/n) > 0 for all s = 0, stability is guaranteed. This was the sufficient condition that was found in [27]. If H  (s/n) < 0 for all s = 0, instability is guaranteed. Another observation is that the matrix Jn is neither symmetric nor anti-symmetric. However, it is normal, i.e. Jn JnT = JnT Jn . Therefore, Jn can be diagonalized by a unitary matrix that also diagonalizes both the symmetric and the anti-symmetric parts of Jn . Note that any real matrix can be expressed as the sum of a symmetric and an anti-symmetric matrices. Since all eigenvalues of a symmetric matrix are real and all eigenvalues of an anti-symmetric matrix are pure imaginary, the eigenvalues of the symmetric part of Jn should give us the real parts of all the eigenvalues of Jn . Thus, stability can be determined by calculating the eigenvalues of the symmetric part of Jn . However, the best thing about the matrix Jn is that it is a circulant matrix [38]. Thus, the eigenvalues of Jn are given

Y.-X. Li / Physica D 180 (2003) 210–234

233

by (see Theorem 5.1): λl = −w O

n−1
k=1

Hk [1 − ei(2π/n)kl ],

l = 0, 1, . . . , n − 1.

(A.5)

These eigenvalues are exactly the same as those obtained in [1] as shown in (4). Here is a brief summary on how to prove this result. Since Jn is circulant, we can express Jn as a polynomial of the basic circulant matrix C: 


Jn  k = − I + H C = − Hk [I − C k ], k wO n−1

n−1

k=1

k=1

where I is the identity matrix and C is basic circulant matrix given by   0 1 0 ··· 0 0     0 1 ··· 0 0   0    0 0 0 ··· 0 0    C= . ··· ··· ··· ··· ··· ···      0 0 0 ··· 0 1    1 0 0 ··· 0 0

(A.6)

(A.7)

The eigenvalues of C are the n complex numbers evenly distributed on the unit ring starting at 1, i.e. µl = ei(2π/n)l , l = 0, 1, . . . , n − 1. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18]

K. Okuda, Variety and generality of clustering in globally coupled oscillators, Physica D 63 (1993) 424–436. D. Hansel, G. Mato, C. Meunier, Synchrony in excitatory neural networks, Neural Comput. 7 (1995) 307–337. U. Ernst, K. Pawelzik, T. Geisel, Delay-induced multistable synchronization of biological oscillators, Phys. Rev. E 57 (1998) 2150–2162. Y.Q. Wang, Z.D. Wang, Y.X. Li, X. Pei, Synchronous phase clustering in a network of neurons with spatially decaying excitatory coupling, J. Phys. Soc. Jpn. 72 (2003) 443–447. Y.X. Li, Y.Q. Wang, R. Miura, Clustering in small networks of excitatory neurons with heterogeneous coupling strengths, J. Comput. Neurosci. 14 (2003) 139–159. Y.X. Li, A minimal network model for quadrupedal locomotion based on symmetry and stability, SIAM J. Appl. Math., submitted for publication. J. Keener, J. Sneyd, Mathematical Physiology, Springer, New York, 1998. G. Leng, Pulsatility in Neuroendocrine Systems, CRC Press, Boca Raton, FL, 1988. J.C. Smith, H.H. Ellenberger, K. Ballanyi, D.W. Richter, J.L. Feldman, Pre-Botzinger complex: a brainstem region that may generate respiratory rhythm in mammals, Science 254 (1991) 726–729. S. Grillner, Neurological bases of rhythmic motor acts in vertebrates, Science 228 (1985) 143–149. S. Grillner, P. Wallen, Central pattern generators for locomotion, with special reference to vertebrates, Ann. Rev. Neurosci. 8 (1985) 233–261. L. Glass, M.C. Mackey, From Clocks to Chaos, The Rhythm of Life, Princeton University Press, New Jersey, 1988. P.P. Gambaryan, How Animals Run: Anatomical Adaptations, Wiley, New York, 1974. J.J. Collins, I.N. Stewart, Coupled nonlinear oscillators and the symmetries of animal gaits, J. Nonlin. Sci. 3 (1993) 349–392. M. Golubitsky, I. Steward, P.-L. Buono, J.J. Collins, A modular network for legged locomotion, Physica D 115 (1998) 56–72. R.J. Full, R. Blickhan, L.H. Ting, Leg design in the hexapedal runners, J. Exp. Biol. 158 (1991) 369–390. D. Golomb, D. Hansel, B. Shraiman, H. Sompolinsky, Clustering in globally coupled phase oscillators, Phys. Rev. A 45 (1992) 3516–3530. D. Golomb, D. Hansel, G. Mato, Mechanisms of synchrony of neural activity in large networks, in: F. Moss, S. Gielen (Eds.), Handbook of Biological Physics, vol. 4, Neuro-informatics and Neural Modelling, Elsevier, Amsterdam, 2001, pp. 887–968.

234

Y.-X. Li / Physica D 180 (2003) 210–234

[19] J. Rinzel, Intercellular communication, in: C.P. Fall, E.S. Marland, J.M. Wagner, J.J. Tyson (Eds.), Computational Cell Biology, Springer, New York, 2002, pp. 140–168. [20] A. Winfree, Biological Rhythms and the behavior of populations of coupled oscillators, J. Theor. Biol. 16 (1967) 15–42. [21] Y. Kuramoto, in: H. Araki (Ed.), International Symposium on Mathematical Problems in Theoretical Physics, Springer, Berlin, 1975, pp. 68–77. [22] R.E. Mirollo, S.H. Strogatz, Sychronization of pulse-coupled biological oscillators, SIAM J. Math. Anal. 50 (1990) 1645–1662. [23] X.J. Wang, J. Rinzel, Alternating and synchronous rhythms in reciprocally inhibitory model neurons, Neural Comput. 4 (1992) 84–97. [24] C. Van Vreeswijk, L.F. Abbott, B. Ermentrout, When inhibition not excitation synchronizes neural firing, J. Comput. Neurosci. 1 (1994) 313–321. [25] Y. Kuramoto, Chemical Oscillations, Waves and Turbulence, Springer, New York, 1984. [26] B. Ermentrout, N. Kopell, Frequency plateaus in a chain of weakly coupled oscillators, SIAM J. Math. Anal. 15 (1984) 215–237. [27] C.C. Chow, Phase-locking in weakly heterogeneous neural networks, Physica D 118 (1998) 343–370. [28] P. Dayan, L. Abbott, Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems, MIT Press, Cambridge, MA, 2001. [29] L. Abbott, C. Van Vreeswijk, Asynchronous states in networks of pulse-coupled oscillators, Phys. Rev. E 48 (1993) 1483–1490. [30] A. Treves, Mean-field analysis of neuronal spike dynamics, Network 4 (1993) 259–284. [31] W. Gerstner, Time structure of the activity in neural network models, Phys. Rev. E 51 (1995) 738–758. [32] D. Golomb, D. Hansel, The number of synaptic inputs and the synchrony of large, sparse neuronal network, Neural Comput. 12 (2000) 1095–1139. [33] D. Hansel, G. Mato, C. Meunier, Clustering and slow switching in globally coupled phase oscillators, Phys. Rev. E 48 (1993) 3470–3477. [34] M. Golubitsky, I. Steward, P.-L. Buono, J.J. Collins, Symmetry in locomotor central pattern generators and animal gaits, Nature 401 (1999) 693–695. [35] M. Golubitsky, I.N. Steward, D.G. Schaeffer, Singularities and Groups in Bifurcation Theory, vol. II, Springer, New York, 1988. [36] J.J. Hopfield, Neurons with graded response have collective computational properties like those of two state neurons, Proc. Natl. Acad. Sci. USA 81 (1984) 3088–3092. [37] T. Aonishi, K. Kurata, M. Okada, Statistical mechanics of an oscillator associative memory with scattered natural frequencies, Phys. Rev. Lett. 82 (1999) 2800–2803. [38] R.A. Horn, C.A. Johnson, Matrix Analysis, Cambridge University Press, Cambridge, 1985.