## Abstract

Electrical synapses are ubiquitous in the mammalian CNS. Particularly in the neocortex, electrical synapses have been shown to connect low-threshold spiking (LTS) as well as fast spiking (FS) interneurons. Experiments have highlighted the roles of electrical synapses in the dynamics of neuronal networks. Here we investigate theoretically how intrinsic cell properties affect the synchronization of neurons interacting by electrical synapses. Numerical simulations of a network of conductance-based neurons randomly connected with electrical synapses show that potassium currents promote synchrony, whereas the persistent sodium current impedes it. Furthermore, synchrony varies with the firing rate in qualitatively different ways depending on the intrinsic currents. We also study analytically a network of quadratic integrate-and-fire neurons. We relate the stability of the asynchronous state of this network to the phase-response function (PRF), which characterizes the effect of small perturbations on the firing timing of the neurons. In particular, we show that the greater the skew of the PRF toward the first half of the period, the more stable the asynchronous state. Combining our simulations with our analytical results, we establish general rules to predict the dynamic state of large networks of neurons coupled with electrical synapses. Our work provides a natural explanation for surprising experimental observations that blocking electrical synapses may increase the synchrony of neuronal activity. It also suggests different synchronization properties for LTS and FS cells. Finally, we propose to further test our predictions in experiments using dynamic clamp techniques.

**electrical synapses****conductance-based model****synchrony****neuronal network model****intrinsic currents****neocortex**

## Introduction

Electrical synapses are sites at which gap-junctions bridge the membranes of two neurons. They have long been known to exist in invertebrates (Watanabe, 1958; Furshpan and Potter, 1959), but only recently has evidence of their ubiquity been unequivocally found in the mammalian brain. Electrical synapses are present in the inferior olive (Llinas and Yarom, 1986), the hippocampus (Draghun et al., 1998; Venance et al., 2000), the cerebellum (Mann-Metzer and Yarom, 1999), the locus coereleus (Christie et al., 1989; Alvarez et al., 2002), the striatum (Kita et al., 1990), the neocortex (Galarreta and Hestrin, 1999; Gibson et al., 1999), the reticular thalamic nucleus (Landisman et al., 2002), and between motoneurons (Kiehn and Tresch, 2002).

Experiments have revealed that electrical synapses are involved in synchronizing neural activity (Draghun et al., 1998; Mann-Metzer and Yarom, 1999; Beierlein et al., 2000; Perez-Velazquez and Carlen, 2000; Tamas et al., 2000; Deans et al., 2001; Hormuzdi et al., 2001). In contrast to the results of these studies, it has been reported recently that inspiratory motoneurons may display more strongly synchronized activity in presence of carbenoxolone (CBX), a blocker of electrical synapses, than in the control situation (Bou-Flores and Berger, 2001). Therefore, in this case, electrical synapses desynchronize neural activity.

The dynamics of networks of neurons interacting via chemical synapses have been studied extensively (Golomb et al., 2001). However, only a few theoretical studies have addressed the dynamics of networks in which neurons are coupled by electrical synapses. Models for pattern generation in the lobster pyloric system have been investigated (Kepler et al., 1990; Abbott et al., 1991; Meunier, 1992). Stable antiphase locking and transitions between in-phase and antiphase locking were demonstrated in pairs of model neurons coupled by electrical synapses, and its impact on rhythmogenesis was investigated (Sherman and Rinzel, 1992; Cymbalyuk et al., 1994; Han et al., 1995). More recently, the discovery of electrical synapses in the CNS of mammals has motivated simulations of conductance-based models (Traub et al., 2001) and analytical studies in the framework of leaky integrate-and-fire (LIF) models (Chow and Kopell, 2000; Lewis and Rinzel, 2003).

Although it is clear that synchronization properties of neurons depend on their intrinsic currents, the general principles that determine this dependence are still unknown in the case of neurons interacting via electrical synapses. In this paper, we combine numerical simulations of conductance-based network models and analytical investigations of quadratic integrate-and-fire (QIF) neurons (Hansel and Mato, 2003) to discover such principles. In particular, we show that potassium currents promote synchrony, whereas sodium currents impede it. We compare our results with existing experimental data. We suggest that our findings may account for those reported by Bou-Flores and Berger (2001). We also propose experiments to further test the predictions of our work. Part of this work has been published previously in abstract form (Pfeuty et al., 2002).

## Materials and Methods

*The conductance-based model.* The membrane potential, *V*, of the neuron follows the equation:
1 where *I*_{L} = -*g*_{L}(*V* - *V*_{L}) is a leak current and Σ_{ion}*I*_{ion} is the sum of all voltage-dependent ionic currents. The currents *I*_{ext} and *I*_{noise} are a constant external current and a Gaussian white noise with a zero mean and SD, σ, respectively.

Our neuron model incorporates an inactivating sodium current, *I*_{Na} = *g*_{Na}*m*_{∞}^{3}*h*(*V* - *V*_{Na}); a delayed rectifier potassium current, *I*_{K} = *g*_{K}*n*^{4}(*V* - *V*_{K}); a slow potassium current, *I*_{Ks} = *g*_{Ks}*s*^{4}(*V* - *V*_{K}) (Erisir et al., 1999); and a noninactivating (persistent) sodium current (French et al., 1990), *I*_{NaP} = *g*_{NaP}*p*_{∞}(*V* - *V*_{Na}). The kinetics of the gating variable *h, n, s* are given by:
2 with *x* = *h, n, s* and α_{h}(*V*) = 0.21*e*^{-}^{(}^{V}^{+}^{58)/20}, β_{h}(*V*) = 3/(1 + *e*^{-}^{(}^{V}^{+}^{28)/10}), α_{n}(*V*) = 0.03(*V* + 34)/(1 - *e*^{-}^{(}^{V}^{+}^{34)/10}), β_{n}(*V*) = 0.375*e*^{-}^{(}^{V}^{+}^{44)/80}, α_{s}(*V*) = 0.07(44 + *V*)/(1 - *e*^{-}^{(}^{V}^{+}^{44)/4.6}), β_{s}(*V*) = 0.008*e*^{-}^{(}^{V}^{+}^{44)/68}. The activation functions *m*_{∞} and *p*_{∞} are given by: *m*_{∞}(*V*) = α_{m}(*V*)/(α_{m}(*V*) + β_{m}(*V*)), where α_{m}(*V*) = 0.1(*V* + 35)/(1 - *e*^{-}^{(}^{V}^{+}^{35)/10}), βm(*V*) = 4*e*^{-}^{(}^{V}^{+}^{60)/18}, and *p*_{∞}(*V*) = 1/(1 + *e*^{-}^{(}^{V}^{+}^{50)/6}).

Throughout this work, the parameters *g*_{Na} = 35 mS/cm ^{2}, *V*_{Na} = 55 mV, *V*_{K} = -90 mV, *g*_{L} = 0.1 mS/cm ^{2}, *V*_{L} = -65 mV, and *C* = 1 μF/cm ^{2} are kept constant, and we study the ways in which the network dynamics depend on the conductances *g*_{K}, *g*_{Ks}, and *g*_{NaP}.

*The network architecture.* The network consists of *N* neurons randomly connected by bidirectional electrical couplings. The average number of connections per neuron is denoted by *M*. We assume that all existing synapses have the same conductance, *g*, and we define a connectivity matrix by *J*_{ij} = *g* if neuron *i* and *j*, (*i, j* = 1..., *N*), are connected and *J*_{ij} = 0 otherwise. The connection between neuron *i* and *j* adds a contribution, -*g*(*V*_{i} - *V*_{j}), to the total synaptic current received by neuron *i*, where *V*_{i} and *V*_{j} are the membrane potential of neurons *i* and *j*, respectively. Therefore, to include the effect of the electrical synapses in the dynamics of the network, an additional current *I*_{i}^{ES} = -Σ_{j}*J*_{ij}(*V*_{i} - *V*_{j}) is added to the right side of the differential equation satisfied by the membrane potential of neuron *i* (see Eqs. 1, 9). In our simulations, the network size is *N* = 1600, the average connectivity is *M* = 10, and the conductance is the same for all of the synapses: *g* = 0.005 mS/cm ^{2}.

*Definition of synchrony and measure of the synchrony level in numerical simulations.* By synchrony of the activity of a pair of neurons, we mean the tendency of these neurons to fire spikes at the same time. In this sense, a pair of neurons will be said to fire in synchrony if the cross-correlation of their spike trains displays a central peak that is above chance. The degree of synchrony can be characterized by the amplitude of this peak. According to this definition, if two neurons fire action potentials periodically and in antiphase, they can be considered to fire asynchronously.

In a large network of neurons, the activity is asynchronous if at any time the number of action potentials fired in the network is the same up to some random fluctuations. This implies that in a large network, the asynchronous state is unique. When this state is unstable, the network activity is necessarily synchronous. In this case, neurons tend to fire preferentially in some windows of time.

Note that a stable asynchronous state can coexist with a stable synchronous state if the network dynamics display multistability. In this work we focus on the conditions under which asynchronous states in large networks and antiphase locking states in pairs of identical neurons become unstable because of electrical synaptic interactions.

One can characterize the degree of synchrony in a population of *N* neurons by measuring the temporal fluctuations of macroscopic observables, as for example, the membrane potential averaged over the population (Hansel and Sompolinsky, 1992; Golomb and Rinzel, 1994). The quantity
3 is evaluated over time, and the variance
of its temporal fluctuations is computed, where
denotes time-averaging. After normalization of σ_{V} to the average over the population of the single-cell membrane potentials,
, one defines χ(*N*):
4 which varies between 0 and 1. The central limit theorem implies that in the limit *n* → ∞, χ(*N*) behaves as:
5 where *a* > 0 is a constant and *O*(1/*N*) means a term of order 1/*N*. In particular, χ(*N*) = 1 if the activity of the network is fully synchronized (i.e., *V*_{i}(*t*) = *V*(*t*) for all *i*), and χ(*N*) = *O*(1/√*N*) if the state of the network activity is asynchronous. In the asynchronous state, χ(∞) = 0. More generally, the larger χ(∞), the more synchronized the population. Note that this measure of synchrony is sensitive not only to the correlations in the spike timing but also to the correlations in the time course of the membrane potentials in the subthreshold range.

*Measure of the irregularity of spike trains.* The irregularity of the spike trains is characterized by the coefficient of variability (CV) (i.e., the ratio between the SD of the interspike intervals and their mean values averaged over the entire network). For periodic spike trains, CV = 0, whereas CV = 1 for spike trains randomly distributed with Poisson statistics.

*Quadratic integrate-and-fire model.* It can be shown that near the onset of periodic firing, the detailed subthreshold dynamics of a large class of neurons can be reduced to a simple one-dimensional model in which the dynamical variable, *v*_{red}, evolves in time (Ermentrout, 1996; Hansel and Mato, 2003):
6 The constants, *C*_{red}, *A, V*^{*}, *I*_{c}^{red}, the reduced external current, *I*_{ext}^{red}, and the reduced noise, *I*_{noise red}, can be computed as functions of the parameters of the full model as shown by Ermentrout (1996) and Hansel and Mato (2003). For *I*_{red} > *I*_{c}^{red}, the solution of Equation 6 diverges in finite time. This corresponds to the firing of an action potential. If one supplements Equation 6 with the condition that the variable *v*_{red} is reset to *V*_{r} < *V*^{*} immediately after *v*_{red} reaches some threshold, *V*_{T}, from below, this yields a reduced model that accurately describes the dynamics of the neuron in the limit *I*_{ext}^{red} → *I*_{c}^{red}. In particular, the current–frequency relationship (*I*–*F* curve) of this model and the neuron behave similarly in this limit, for any value of the parameters *V*_{r} and *V*_{T}. When *I*_{ext}^{red}-*I*_{c}^{red} is not small, the reduced model is no longer an exact description of the neuronal dynamics. However, one can fit the parameters so that the reduced model provides a good approximation of the *I*–*F* curve of the neuron. The difference, *V*_{T} - *V*_{r}, will be called the reset depth.

It is convenient to rewrite the subthreshold dynamics of the reduced model in terms of the dimensionless variables *ṽ*_{red} and *Ĩ*^{red} defined by:
7
8 where τ_{0} has the dimension of a time. This yields:
9 The variable *ṽ*_{red} is reset to *ṽ*_{r} = *A*(*V*_{r} - *V*^{*})τ_{0} /*C*_{red} whenever it reaches the threshold value: *ṽ*_{T} = *A*(*V*_{T} - *V*^{*})τ_{0} /*C*_{red}. Note that *ṽ*_{r} < 0. A similar model has been used to study networks of excitatory and inhibitory neurons (Latham et al., 2000; Hansel and Mato, 2003).

The subthreshold dynamics (Eq. 9) need to be supplemented with a model for the suprathreshold part of the membrane potential time course. Assuming that the width of the action potentials is much smaller than the interspikes, we represent their time course by a δ-function at each time a spike is fired. Therefore, the (dimensionless) reduced membrane potential of the neuron can be written: 10 where θ measures the integral over time of the suprathreshold part of an action potential.

The model defined by Equations 9 and 10 will be called the QIF model. To simplify notations of the reduced model, we will drop the index (“red”) and the tildes.

*Phase reduction in the weak coupling limit.* Let us consider a neuron firing action potentials periodically with an interspike *T*. A small and instantaneous perturbation applied at time Δ*t* after a spike induces a small change in the timing of the subsequent spikes. This change, which depends on Δ*t*, or equivalently on ϕ = 2πΔ*t*/*T*, can be characterized by a function, *Z*(ϕ), which measures the delay or the advance induced in the firing times after the perturbation. A positive value of *Z*(ϕ) indicates that the perturbation advances the subsequent spikes. A negative value of *Z*(ϕ) indicates a delay. The effect of noninstantaneous weak perturbations (such as spikelets) can be estimated by convolving the phase-response function (PRF), *Z*, with the perturbation. It can be shown that the dynamic behavior of a network of weakly interacting neurons can be completely described in terms of the response function in the framework of the phase reduction approach (Kuramoto, 1984; Ermentrout and Kopell, 1986; Hansel et al., 1993, 1995; Golomb and Hansel, 2000; Neltner et al., 2000; Lewis and Rinzel, 2003). In this approach, a phase variable is associated with each neuron in the network. This variable, ϕ_{i}(*i* = 1,..., *N*), measures the time elapsed since the last action potential fired by neuron *i*. For a network of identical neurons coupled with electrical synapses, one can show (see Appendix) that the phase variables follow a set of *n* first-order coupled differential equations:
11 Γ(ϕ_{i} - ϕ_{j}) the phase coupling between neurons *i* and *j* given by:
12*J*_{ij} is the connectivity matrix, and η_{i}(*t*) is a white noise with zero mean and variance
with σ, the SD of the noise of the nonreduced dynamics.

*Numerical integration.* In the simulations of the conductance-based model, the differential equations were integrated using the second-order Runge-Kutta scheme with fixed-time step: Δ_{t} = 0.01 msec. Averaged quantities (firing rate, CV, χ) were computed over a time period of 1 sec after discarding a transient of 500 msec.

## Results

### Single neuron and coupling properties of the conductance-based model

#### Firing properties of single neurons

In the absence of persistent sodium and slow potassium currents (the “control model”), our conductance-based model neuron fires tonically for large enough external currents, I_{ext} > *I*_{c} = 0.16 μA/cm^{2}. The frequency of the discharge in response to a step of current is plotted as a function of the step amplitude (*I*–*F* curve) in Figure 1*A.* For *I*_{ext} larger than but close to *I*_{c}, the minimum value of the external current required for the neuron to fire, the firing frequency can be arbitrarily small, because action potentials appear at *I*_{c} through a saddle-node bifurcation [for a definition of a saddle-node bifurcation, see Strogatz (1994) and Rinzel and Ermentrout (1998)]. Decreasing *g*_{K} does not significantly modify the resting potential and the rheobase of the neuron, but it increases the gain of the *I*–*F* curve (Fig. 1*B*). The slow potassium (*I*_{Ks}) and persistent sodium (*I*_{NaP}) currents modify the onset of firing, because these currents are activated at rest. As should be expected, the persistent sodium current increases the excitability of the neuron and reduces its rheobase, *I*_{c} (Fig. 1*C*). For *g*_{NaP} > 0.08 mS/cm^{2}, the neuron is spontaneously active. As shown in Figure 1*D*, the main effects of adding *I*_{Ks} on the *I*–*F* curve are a translation to the right, a linearization, and, if *g*_{Ks} is large enough, a suppression of the ability of the neuron to fire at low rates (because for large enough *g*_{Ks}, a subcritical Hopf bifurcation occurs at the onset of firing [for a definition of a subcritical Hopf bifurcation, see Strogatz (1994) and Rinzel and Ermentrout (1998)].

The response of the neuron to a step of external current, whose value is adjusted to give a firing frequency of 50 Hz, is plotted in Figure 2 for different combinations of the intrinsic currents. The action potentials are shown with a higher temporal resolution in Figure 3*C.* In the control case (Fig. 2*A*; and Fig. 3*C*, solid line), action potentials are narrow and are followed by a strong afterhyperpolarization (AHP) that brings the membrane potential 15 mV below its resting value. When *g*_{K} is decreased (Fig. 2*B*), the resting potential of the neuron does not change, but the AHP is reduced and the membrane potential always remains above rest. This conductance also controls the width of the spike and the depth of the AHP [Fig. 3*C*, compare the solid line (*g*_{K} = 9 mS/cm^{2}) with the dotted line (*g*_{K} = 2.5 mS/cm^{2})].

The persistent sodium and the slow potassium currents significantly modify the resting potential and the external current required to adjust the firing rate. In contrast, they only slightly affect the depth of the repolarization, measured from threshold, after the action potential (Fig. 2*C,D*). The size and width of the action potentials remain almost unchanged when the conductances of these currents are varied [Fig. 3*C*, compare the dotted line with the dashed-dotted line (effect of *g*_{Ks}) and the solid line with the dashed line (effect of *g*_{NaP})].

#### Properties of the spikelets

When a neuron fires an action potential, it induces a depolarization, called a spikelet, in the neurons that are connected to it. The spikelets generated by a neuron firing one isolated spike in response to a very brief but strong pulse of current (Fig. 3*A*) or firing tonically at 50 Hz (Fig. 3*C*), are shown in Figure 3, *B* and *D*, respectively, for different parameters of the ionic currents.

In all four cases displayed, the width of the spikelets is much larger than the width of the presynaptic spikes that generate them. This is because the spikelets are a filtered version of the presynaptic membrane potential profile.

The size of the spikelet is practically unaffected by the presence of the persistent sodium (Fig. 3*B,D*, compare the solid line with the dashed line), because this current does not greatly affect the shape of the action potential. The only noticeable effect of the persistent sodium on the spikelet is to slow down its decay time course. This is because it reduces the input conductance of the postsynaptic neuron (by ∼50% at -76 mV).

The modulation of the spikelets by the delayed rectifier potassium current is primarily caused by changes in the presynaptic neuron voltage time course. When *g*_{K} is reduced, the action potential is broader, and this induces an increase in the amplitude and the width of the spikelets. Although the slow potassium current does not greatly affect the width of the action potential, it contributes to bringing the membrane of the presynaptic neuron transiently below its holding potential. This explains why in Figure 3*B* the spikelet is narrower and displays a faster and deeper repolarization in the presence of *I*_{Ks}.

The spikelets generated by tonically firing neurons at low firing rates are similar to those in Figure 3*B* (data not shown). For high-enough firing rates, the slow potassium saturates, affecting the dynamics primarily as would a constant hyperpolarizing current. Therefore, its effect on the voltage traces and on the spikelets becomes less pronounced (Fig. 3, compare *B* and *D*).

### Synchrony in the conductance-based network model

The goal of this section is to show how synchrony properties of our conductance-based network model depend on the intrinsic currents of the neurons. Further understanding of the effects described in this section will be provided in the next section by investigating analytically a simplified and more abstract model.

#### Synchrony is modified when intrinsic conductances are changed A pair of conductance-based neurons

The effect of intrinsic currents on neuronal network dynamics can first be demonstrated by studying the dynamics of a pair of neurons. Here, as an example, we show that a persistent sodium current tends to promote antiphase locking, whereas a potassium current promotes in-phase locking. In Figure 4*A*, the traces of two identical neurons (*g*_{K} = 9 mS/cm^{2}; *g*_{Ks} = *g*_{NaP} = 0) are plotted over a time interval of 100 sec. Both neurons receive a noisy input that makes them fire at ∼50 Hz. These traces suggest that the neurons fire action potentials in synchrony. This is confirmed by the cross-correlation of the traces computed over a longer time interval (100 sec), which displays a strong peak centered around *t* = 0 msec, as shown in Figure 4*B.* In contrast, when the persistent sodium conductance is large enough (*g*_{NaP} = 0.4 mS/cm^{2}), neurons tend to fire in antiphase (Fig. 4*C*), and the peak of the cross-correlogram is shifted to 10 msec (Fig. 4*D*). If a slow potassium current is added (*g*_{Ks} = 0.15 mS/cm^{2}) while keeping the persistent sodium conductance the same, antiphase locking becomes unstable and neurons tend to fire in-phase (Fig. 4*E, F*). Therefore, a slow potassium current is a promoter of synchronous activity.

#### Large networks of conductance-based neurons

In Figure 5 the synchrony measure, χ, is plotted against the conductances *g*_{K}, *g*_{Ks}, and *g*_{NaP}. Examples of raster plots of network activity are also displayed in this figure, for particular values of the intrinsic conductances (Fig. 5, insets). In each of the panels, the external input (average deviation and SD) is kept constant. Therefore, when the intrinsic conductances vary, the average firing rate of the neurons and the CV of their interspike intervals change. In particular, the firing rate decreases (respectively increases) when the potassium (respectively the persistent sodium) conductances increase (top right figures in each of the panels in Figure 5).

In Figure 5*A*, the effect of *g*_{K} on the synchrony level is shown (*g*_{Ks} = *g*_{NaP} = 0). When *g*_{K} is small enough (*g*_{K} < 4.5 mS/cm^{2}), χ is very close to zero. A detailed study shows that in this region, χ is on the order of 1/√*N*, where *N* is the number of neurons in the network, and vanishes in the limit of a very large network (data not shown): non-zero values of χ are attributable to finite size effects, and the network is asynchronous. For *g*_{K} of >4.5 mS/cm^{2}, χ increases rapidly up to χ ≈ 0.35. For example, in the raster plot shown in the right inset of Figure 5*A* (control situation, *g*_{K} = 9 mS/cm^{2}), χ = 0.34. This corresponds to a state in which the neurons fire together within time windows of approximately one-third of the period.

Figure 5*B* shows the effect of the slow potassium current for *g*_{NaP} = 0 and *g*_{K} = 2.5 mS/cm^{2}. For *g*_{Ks} = 0, the neurons fire asynchronously. At *g*_{Ks} ≈ 0.06 mS/cm^{2}, χ starts to increase. For large *g*_{Ks}, χ saturates at a value, χ ≈ 0.55, that is substantially larger than in the control situation. This value of χ corresponds to a tight synchrony of the action potentials fired by the neurons (Fig. 5*B* inset on the right).

In Figure 5*C*, χ is plotted against *g*_{NaP} for *g*_{K} and *g*_{Ks} at their control values (*g*_{K} = 9 mS/cm^{2}, *g*_{Ks} = 0). Clearly, χ is a decreasing function of *g*_{NaP}. When *g*_{NaP} is large enough (*g*_{NaP} > 0.1 mS/cm^{2}), the neurons fire asynchronously.

In the results presented in Figure 5, the firing rate of the neurons is not controlled. Therefore, the origin of the variation of χ with the intrinsic current conductances is not clear. It may be attributable to the fact that, in general, synchrony depends on the firing rate for given intrinsic currents. It might also be a result of the fact that dynamic properties of neurons with different intrinsic currents may be different even if they have the same firing rate.

#### Dependence of synchrony on intrinsic currents at fixed firing rates

To clarify further how potassium and sodium currents affect the synchrony of the network, we performed another set of numerical simulations in which the external input and the noise level were tuned to keep constant both the average firing rate of the neurons and the CV of their spike trains when the conductances of the intrinsic currents were changed. The results are shown in Figure 6 (three panels on the left). In all of the simulations, the average firing rate of the network was 50 Hz and the CV ≈ 0.12 (Fig. 6, panels on the right).

The level of synchrony increases with the potassium conductances. For small-enough *g*_{K} or *g*_{Ks}, the activity of the network is weakly correlated. A detailed study of the dependence of χ with the network size reveals that for *g*_{K} (respectively *g*_{Ks}) smaller than *g*^{*}_{K} ≈ 4 mS/cm^{2} (respectively *g*_{Ks} smaller than *g*^{*}_{Ks} ≈ 0.075 mS/cm^{2}), χ is on the order of 1/√*N* (data not shown). Therefore in this range of conductances, the network is in the asynchronous state. At *g*^{*}_{K} and *g*^{*}_{Ks}, a transition to a synchronous state occurs, and beyond these values the network settles into a synchronous state. This occurs although the level of noise has been increased.

In contrast to the potassium conductances, which promote synchrony, the persistent sodium current impedes it, as shown in Figure 6*C.* For large-enough *g*_{NaP}, synchrony is even destroyed. A sharp transition to the asynchronous state occurs for *g*_{NaP} ≈ 0.15 mS/cm^{2}.

Relying on the way intrinsic currents affect spikelet modulation (Fig. 3), one might naively expect that potassium currents that reduce the size of the spikelets should also reduce synchrony, and that persistent sodium current, which increases this size, should promote it. As a matter of fact, we found a trend that was exactly the opposite.

#### Firing rate affects synchrony differently when different intrinsic currents are involved

We also investigated in what ways the network dynamic state depends on the firing rate of the neurons. The firing rate was changed by varying the average external input, and χ was plotted as a function of the frequency of the neurons for several combinations of potassium and sodium conductances.

Figure 7 plots the synchrony measure, χ, versus the average firing rate of the neurons for four sets of conductance parameters and two levels of noise. The qualitative behavior of χ as a function of the firing rate depends on the intrinsic currents involved.

As expected, χ is always larger for the smaller level of noise (σ = 0.3 mV/msec^{1/2}) (Fig. 7, circles) than for the larger one (σ = 0.6 mV/msec^{1/2}) (Fig. 7, squares). In particular, for σ = 0.6 mV/msec^{1/2}, the noise is strong enough to prevent synchrony in the entire range of firing frequencies for *g*_{K} = 3.5 mS/cm^{2} (Fig. 7*B*) and for *g*_{NaP} = 0.15 mS/cm^{2} (Fig. 7*C*). This level of noise is sufficient to destroy synchrony as well in the control situation, but only at <20 Hz (Fig. 7*A*).

In the control case, χ increases monotonously with the firing rate (except for a small decrease followed by a slight increase between 50 and 100 Hz for σ = 0.3 mV/msec^{1/2}). In contrast, for smaller *g*_{K} (e.g., *g*_{K} = 3.5 mS/cm^{2}), we find a nonmonotonous variation of synchrony: χ increases at low firing rates but decreases again for firing rates of >20 Hz. The network state becomes asynchronous for firing rates of >60 Hz (Fig. 7*B*, squares). If the slow potassium current is added (*g*_{Ks} = 0.15 mS/cm^{2}), χ again varies monotonously with the firing rate (Fig. 7*D*). Finally, with enough persistent sodium current and not overly strong noise (*g*_{K} = 9 mS/cm^{2}, with *g*_{NaP} = 0.15 mS/cm^{2}, *g*_{Ks} = 0) (Fig. 7*C*, squares), χ varies nonmonotonously with the firing rate, as in the reduced *g*_{K} case (Fig. 7*B*).

### Synchrony in networks of quadratic integrate-and-fire neurons depends on the phase response function

The PRF, which characterizes how neurons respond to small perturbations, is a key concept to understand the relationship between intrinsic properties of neurons and their collective dynamics (Kuramoto, 1984; Hansel et al., 1993; van Vreeswijk et al., 1994; Kopell and Ermentrout, 2002). This function depends on the excitability properties of the neurons and therefore is determined by the intrinsic currents involved in their dynamics. One expects that the differences in the dynamic behavior of our conductance-based network model for different sets of parameters reflect the changes in the excitability properties of the neurons. However, relying on numerical simulations alone, it would very difficult to establish general principles to relate the excitability properties and the neuronal PRF to the synchronization properties. Below, we consider a network of QIF neurons (Eqs. 9, 10) fully connected by electrical synapses. As we show below, the qualitative properties of the excitability of the QIF neurons crucially depend on the parameters (the threshold and the reset voltages and the external current). Thanks to its relative simplicity, this model can be investigated using analytical techniques. It reveals how, in the framework of this model, one can relate the stability of antiphase locking of a pair of neurons and the stability of the asynchronous state of a large network to the shape of the PRF of the neurons. Subsequently, we show that similar rules can be applied to conductance-based models, and that they provide a unified framework explaining the diversity of behavior found in our numerical simulations of the conductance-based model.

#### The membrane potential and the PRF of QIF neurons

In the Appendix, it is shown that between two successive spikes the subthreshold membrane potential of a QIF neuron reads:
13 where *t* measures the time elapsed since the first spike. The firing period, *T*, is determined by the condition: *v*(*T*) = *V*_{T}. Therefore:
14 In the weak coupling limit, the dynamics of the network are completely determined by the phase-response function, Z(ϕ) (see Materials and Methods), which depends on three parameters, namely, the firing frequency of the neurons (or equivalently the external current *I*_{ext}), the threshold *V*_{T}, and the reset potential *V*_{r}. The response function of the QIF model can be derived analytically as shown in the Appendix. One finds:
15

Changing the ratio -*V*_{r}/*V*_{T} for a fixed reset depth, *V*_{T} - *V*_{r} has a strong influence on the shape of the trajectory of the phase-response function of the neuron, as shown in Figure 8. In Figure 8*A*, the membrane potential of a QIF neuron firing at 50 Hz is plotted for three values of the ratio *V*_{r}/*V*_{T} and a fixed reset depth: *V*_{T} - *V*_{r} = 3. For -*V*_{r}/*V*_{T} = 0.1, the concavity of the subthreshold trajectory is always directed upward. In contrast, for -*V*_{r}/*V*_{T} = 10, it is always downward. For -*V*_{r}/*V*_{T} = 1, the concavity changes from upward to downward.

It is easy to see that *Z*(ϕ) > 0 for all ϕ, and that it is with one maximum. The location of the maximum of *Z* depends on *V*_{r}, *V*_{T}, and *I*_{ext}. This is depicted in Figure 8, where the voltage trace and the response function are plotted for different values of these parameters. For -*V*_{r}/*V*_{T} = 1, Z(ϕ) is symmetric and its maximum is always at ϕ = π. For -*V*_{r}/*V*_{T} < 1, the maximum of *Z* is located in the first half of the firing period, whereas for -*V*_{r}/*V*_{T} > 1, it is located in the second half. The value of Z at the maximum also depends on *V*_{r}, *V*_{T}. It has smaller values when -*V*_{r}/*V*_{T} = 1. The larger the firing rate, the stronger the dependence of *Z* on -*V*_{r}/*V*_{T} (Fig. 8, compare *B*, for 50 Hz, and *C*, for 10 Hz).

#### A pair of identical QIF neurons without noise

When, in the absence of noise, two identical neurons interact in a symmetric way, they reach, at large time, a phase-locked state in which the two neurons fire action potentials with a fixed phase shift, Δ. In the Appendix, it is shown that Δ is a solution to the equation:
16 where the function Γ_{-}(ϕ) ≡ 1/2(Γ(ϕ) - Γ(-ϕ)), and Γ(ϕ) is the effective phase coupling (see Materials and Methods). Because of the symmetry of the system and the 2π periodicity of the function Γ_{-}, there are always at least two solutions to this equation: Δ = 0 and Δ = π. The first solution corresponds to the state in which the two neurons fire in-phase (in-phase locking), the second one to the state in which they fire in antiphase (antiphase locking). Other solutions with 0 < Δ < π may also exist. However, only solutions that are stable can be reached at large time. Therefore, out of all of the solutions to Equation 16, only those that satisfy the condition (see Appendix):
17 correspond to phase-locked states that the two neurons can eventually reach starting from appropriate initial conditions. In the following, we focus on the stability of the antiphase state of a pair of QIF neurons.

#### Stability conditions of antiphase locking

The voltage trajectory consists of three parts: (1) a subthreshold part, during which the membrane potential increases with time from a value *V*_{r} to a threshold value *V*_{T}; (2) a spike, modeled by a δ-function of amplitude θ, and (3) a resetting, where potential is instantaneously brought back to V_{r}. Then, as explained in the Appendix, antiphase locking of a pair of QIF neurons is stable if:
18 with:
19
20
21

The first term, *S*_{sp}, is the contribution of the presynaptic spikes to the destabilization of antiphase locking. The second term, *S*_{r}, is the contribution of the instantaneous reset of the membrane potential. The last term, *S*_{sub}, corresponds to the effect of the coupling between the two neurons when the presynaptic neuron is subthreshold. These three terms are plotted in Figure 9*A* as a function of the ratio -*V*_{r}/*V*_{T}, for θ = 1, *V*_{T} - *V*_{r} = 3, and a frequency *v* = 50 Hz.

The terms S_{r} and S_{sub} (Fig. 9, dotted line and dashed-dotted line, respectively) are invariant under the transformation -*V*_{r}/*V*_{T} → *V*_{T}/-*V*_{r}. That is why the corresponding curves in Figure 8*B* are symmetric around -*V*_{r}/*V*_{T} = 1 (note the logarithmic scale of the *x*-axis). S_{r} is always positive and therefore destabilizes antiphase locking. In contrast, *S*_{sub} is always negative and leads to stabilized antiphase locking. The sum of the two terms is also plotted in Figure 9*A* (double-dotted-dashed line) to show that their overall contribution is destabilizing except for very small or very large values of -*V*_{r}/*V*_{T}.

Because the maximum of *Z* moves from the first part of the period to the second part of the period when -*V*_{r}/*V*_{T} crosses 1 (Fig. 8), the sign of the derivative of *Z* at ϕ = π also changes from negative to positive, and the term *S*_{sp} (dashed line) tends to destabilize antiphase locking for -*V*_{r}/*V*_{T} > 1. The sum of all three terms of the right side of Equation 18 is plotted in Figure 9*A* (solid line). Antiphase locking is stable only if the phase at which *Z* reaches its maximum is small enough.

The phase diagram for the stability of antiphase locking as a function of -*V*_{r}/*V*_{T} and the firing frequencies is displayed in Figure 9*C.* The stability of antiphase locking requires that -*V*_{r}/*V*_{T} be small enough. Moreover, the range of values of -*V*_{r}/*V*_{T} where it is stable increases with the frequency of the neurons. This is because the contribution of the term *S*_{sp} in the right side of Equation 18 grows with the firing rate and tends to dominate the contribution of the two other terms. In contrast, at low frequencies, the response function weakly depends on -*V*_{r}/*V*_{T} (Fig. 8), and *S*_{sp}, which is proportional to the frequency, is negligible. Moreover, *S*_{r} + *S*_{sub} is positive except for very small or very large values of -*V*_{r}/*V*_{T}. Therefore, for small frequencies, the domain of stability of antiphase locking is very narrow.

This analysis shows that for a pair of neurons, the stability of antiphase locking depends critically on the shape of the phase-response function and in particular on the sign of the derivative of the response function at the mid-period. The more skewed toward the right (the second part of the firing period) the response function of the neurons, the more unstable antiphase locking becomes.

#### Large networks of weakly coupled QIF neurons

The stability analysis of the asynchronous state of a large network of neurons coupled all-to-all receiving a noisy input is also simplified if one assumes weak coupling. Using the phase-reduction approach, it can be shown (see Appendix) that the asynchronous state is stable if:
22 For all integers, n > 0, with *Z*_{n} and *v*_{n} the *n*th-Fourier components of the function *Z*(ϕ) and *v*(ϕ), respectively (Eqs. 13, 15) defined by
and a similar equation for *v*_{n}. If for some *n*, μ_{n} is negative (respectively positive), it is said that the mode of order n is stable (respectively unstable). For the asynchronous state to be stable, the modes at any order must be stable. It is unstable when at least one mode is unstable.

The first two terms in Equation 22 correspond to the effect of the interaction. The first term represents the effect of the spikes and can be written:
23 The second term, μ ′ _{n} = *ng Im*(*nZ*_{n}*v*_{-}_{n}), corresponds to the combined effect of the reset of the membrane potential and its subsequent subthreshold evolution.

The sign of their sum depends on the parameters, *V*_{r}, *V*_{T} and on the frequency of the neurons, *v.* The last term in Equation 22 corresponds to the effect of the noise. It is always positive. Therefore, as should be expected, noise increases the stability of the asynchronous state (because of the negative sign in front of this term in Eq. 22). The stability of the asynchronous state depends on the competition between the first two terms and the last term. In particular, the sign of μ_{n} depends on the ratio σ^{2}/*g*. Note that because of the factor *n*^{2} in the last term of Equation 22, the stability of the modes increases rapidly with their order.

We first consider the case of a noiseless network (i.e., with σ = 0). In Figure 9*B*, we plotted μ′ _{1} (dashed line) and μ ′ _{1} (dashed-dotted line) against -*V*_{r}/*V*_{T}. The qualitative behavior of these quantities is similar to those of *S*_{sp} and *S*_{sub} + *S*_{r}, respectively, for a pair of neurons. In particular μ ′ _{1} is symmetric (in a semilogarithmic scale) around -*V*_{r}/*V*_{T} = 1. Moreover, μ′ _{1} increases monotonically from negative values to positive values when -*V*_{r}/*V*_{T} increases, and it changes sign at -*V*_{r}/*V*_{T} = 1. This can be easily explained. Indeed, for -*V*_{r}/*V*_{T} > 1 (respectively, -*V*_{r}/*V*_{T} < 1), the function *Z*(ϕ) is skewed toward ϕ < π (respectively ϕ > π) where sin(ϕ) > 0 [respectively sin(ϕ) < 0]; therefore, for *n* = 1, the integral in Equation 23 is negative (respectively positive) and μ′ _{1} < 0 (respectively μ′ _{1} > 0). In particular, because for -*V*_{r}/*V*_{T} = 1, *Z*(ϕ)is symmetric around ϕ = π, the integral in Equation 23 vanishes for *n* = 1 [because sin(ϕ + π) = -sin(ϕ)]. The sign of μ_{1} (solid line) is primarily determined by the sign of μ′ _{1}. It is negative for small -*V*_{r}/*V*_{T} and becomes positive for -*V*_{r}/*V*_{T} ≈ 1. Therefore the mode n = 1 is stable only if -*V*_{r}/*V*_{T} is small enough. Similar analysis can be performed for the modes n > 1. It reveals that if -*V*_{r}/*V*_{T} is not too large, the mode *n* = 1 determines the stability of the asynchronous state.

The phase diagram for the stability of the asynchronous state as a function of -*V*_{r}/*V*_{T} and the firing frequency is plotted in Figure 9*D.* It is very similar to the phase diagram for the stability of the antiphase state for a pair of neurons (Fig. 9*C*). The only qualitative difference between these two phase diagrams is that for a pair of neurons, antiphase locking is stable at low firing rates and very large -*V*_{r}/*V*_{T}, whereas for a large network, the asynchronous state is unstable in that limit. Therefore, we can conclude that neurons coupled with electrical synapses will be more easily synchronized if their phase-response functions are skewed to the right than to the left.

The phase diagram for the stability of the asynchronous state is plotted in Figure 10 for different levels of noise. For large-enough firing rates, the instability lines for σ^{2}/*g* = 0 and σ^{2}/*g* = 0.01 are very close to each other. The distance between the two lines increases when the frequency decreases. At some value, *v* = *v*^{*}, the two lines separate completely. At *v*^{*}, the line for σ ^{2}/*g* = 0.01 changes direction and continues toward the right of the phase diagram (large -*V*_{r}/*V*_{T}). In contrast, the line for σ = 0 continues toward the left of the diagram. In particular, for σ^{2}/*g* = 0.01, the asynchronous state is stable for any value of -*V*_{r}/*V*_{T} for small-enough firing rates. As a matter of fact, one can show analytically that the asynchronous state is stable for any value of -*V*_{r}/*V*_{T} if the firing rate is smaller than a critical value, *v*^{*}(σ ^{2}/*g*), which vanishes with σ^{2}/g. The size of the region in which the asynchronous state is stable increases with the noise level, as expected. In particular, *v*^{*}, and the frequency range in which the asynchronous state is stable for all -*V*_{r}/*V*_{T}, increases with σ^{2}/*g*.

How synchrony depends on the ratio -*V*_{r}/*V*_{T} at fixed frequencies can be deduced from the phase diagram in Figure 10. For low and moderate noise levels (e.g., σ^{2}/*g* = 0.01 and σ^{2}/*g* = 0.05), we find three generic behaviors as a function of -*V*_{r}/*V*_{T}: (1) At low firing rates, the asynchronous state is stable for all -*V*_{r}/*V*_{T}. (2) At high firing rates, the asynchronous state is stable for small -*V*_{r}/*V*_{T} and becomes unstable when -*V*_{r}/*V*_{T} increases. (3) In some intermediate range of firing rates, the stability of the asynchronous state varies nonmonotonously with -*V*_{r}/*V*_{T}: when -*V*_{r}/*V*_{T} increases from 0, the asynchronous state is first stable, then unstable, and stable again for large -*V*_{r}/*V*_{T}. The size of the intermediate region in which the asynchronous state is unstable decreases when the noise level increases.

Figure 10 also allows us to predict how the stability of the asynchronous state depends on the frequency for fixed -*V*_{r}/*V*_{T}. For small-enough -*V*_{r}/*V*_{T}, the asynchronous state is stable at all firing rates. The range of -*V*_{r}/*V*_{T} where this happens is larger when the noise is stronger. For large-enough -*V*_{r}/*V*_{T}, the asynchronous state is stable at low firing rates and becomes unstable when the firing rate is large enough. If the noise is not too strong, the stability of the asynchronous state varies nonmonotonously in some intermediate range of -*V*_{r}/*V*_{T} < 1. In this domain, the asynchronous state is stable at low rates, loses stability when the rate increases, but becomes stable again at high firing rates.

### Interpretation of the simulation results of the conductance-based network

Our analysis of the synchrony properties of QIF networks shows that the stability of the asynchronous states depends crucially on the shape of the phase-response function, *Z*(ϕ), and in particular on the location of its maximum. This suggests that the desynchronization effect of the persistent sodium and the synchronization effect of the potassium currents revealed by our numerical simulations may be related to the way these currents shape the phase-response functions of the neurons.

The function, *Z*(ϕ), for our conductance-based neuron model is plotted in Figure 11 for different values of *g*_{K}, *g*_{Ks}, and *g*_{NaP}. In Figure 11*A–C*, Z(ϕ) > 0 for all ϕ, except in a very tiny region after the spike peak. In Figure 11*D*, the region of negative *Z*(ϕ) is larger. However, in all four cases, negative values of *Z*(ϕ) are much smaller than positive values. This means that most of the time, a small and brief depolarizing perturbation advances subsequent firing of action potentials.

In all of the cases presented in Figure 11, the function *Z*(ϕ) has only one maximum, whose location depends on the intrinsic currents. Comparison of Figure 11*A* with Figure 11*B* shows that the maximum of *Z* is shifted to the right when *g*_{K} increases, and that the global shape of the phase-response function is more skewed to the right (respectively to the left) for large (respectively small) values of *g*_{K}. The slow potassium current has an effect similar to that seen for the delayed rectifier current. It distorts the phase-response function toward the second half of the firing period. Therefore, increasing the potassium conductances has a similar effect on the phase-response function of the conductance-based neuron as increasing -*V*_{r}/*V*_{T} does for the QIF neuron. In contrast, the persistent sodium current shifts the peak of the response to the left and skews the shape of the phase-response function in that direction, as seen by comparing Figure 11*A* with Figure 11*C.* Hence, increasing *g*_{NaP} in the conductance-based model and decreasing -*V*_{r}/*V*_{T} in the QIF model have the same qualitative effect.

We are now able to understand the results of our numerical simulations. In Figure 6, *A* and *B*, the network state is asynchronous for small potassium conductances (both delay-rectifier and slow-potassium) and becomes synchronous when these conductances are increased. A similar behavior is found in the QIF network for *v* > *v*_{*}(σ): the asynchronous state is stable when -*V*^{r}/*V*^{T} is small, and the asynchronous state is destabilized when -*V*_{r}/*V*_{T} increases enough. This stems from the fact that the qualitative changes in the phase-response function of the conductance-based neuron when *g*_{K} increases and the QIF neuron when -*V*_{r}/*V*_{T} increases are qualitatively similar. Conversely, increasing the conductance of the persistent sodium current has the same effect on the phase-response function of the neuron as a decrease of -*V*_{r}/*V*_{T}. This explains the results of Figure 6*C.*

The phase diagram of Figure 10 also explains the different behaviors we found in our simulations for fixed-current conductances when the frequency varies (Fig. 7) in presence of noise. For instance, in the control case, the response function is skewed to the right. This is equivalent to -*V*_{r}/*V*_{T} > 1 for the QIF model. This explains the monotonous transition from the asynchronous state at low firing rates to synchrony when the firing rate increases, as shown in Figure 7, *A* and *D.* The nonmonotonous behavior in Figure 7, *B* and *C*, can also be explained. In Figure 7*C*, *g*_{K} is smaller than in the control model. The corresponding effective value of -*V*_{r}/*V*_{T} is in the intermediate range, where the asynchronous state changes stability nonmonotonously with the frequency. A similar argument can account for the results in Figure 7*B*, but here the change in the phase-response function is attributable to the persistent sodium current.

## Discussion

### How intrinsic currents affect synchrony in networks of neurons connected by electrical synapses

In this work, we have derived explicit rules for the stability of the asynchronous state in networks of neurons interacting via electrical synapses. We have shown analytically how the shape of the PRF, which depends on the model parameters, determines the stability of the asynchronous state in QIF networks. We have extended these results to provide a unified framework explaining the diversity of behavior found in our numerical simulations when potassium and sodium conductances are changed.

One can understand intuitively how potassium and sodium currents shape the PRF. The current, *I*_{K}, increases the refractoriness of the neuron; therefore, it reduces its responsiveness to small depolarization after a spike. A similar effect, but stronger and more lasting, occurs with *I*_{Ks}. This intuitively explains why we found that larger conductances of these currents skew the PRF toward the second half of the firing period. Similar effects have been found by Ermentrout et al. (2001) for M current and AHP potassium current. The current *I*_{NaP} is an inward current. It is already activated near rest, and depolarizing perturbations amplify this activation. Hence *I*_{NaP} increases the responsiveness of the neuron after a spike and shifts the maximum of the response function toward the first half of the period, as we have found here.

We predict that potassium currents, like *I*_{K} and *I*_{Ks}, tend to promote synchrony of neurons coupled via electrical synapses and that, in contrast, *I*_{NaP} tends to oppose it. Our other predictions concern the dependency of the synchrony level with the firing rate. We expect that these conclusions do not depend on the particular models we have chosen for the intrinsic currents.

### Assumption of weak coupling

The stability of the asynchronous state of a large network of all-to-all connected QIF neurons can be calculated analytically in the absence of noise for any coupling strength using the population density method (Abbott and van Vreeswijk, 1993; Hansel and Mato, 2003). An alternative approach is the phase-reduction method. Although it is only mathematically justified for weak interactions, it leads to nontrivial results that remain valid in a reasonable range of interaction strengths. Moreover, it is straightforward to study the effect of noise in this framework. These were the motivations to use the phase-reduction method in the framework of which all of the analytical results of this paper were derived.

At finite coupling strength, an additional instability of the asynchronous state appears at low firing rates and small -*V*_{r}/*V*_{T}. It corresponds to the impossibility of controlling low firing rates when the coupling is too strong and the AHP is too small. This is because for a small AHP, the recurrent electrical synapse interactions tend to depolarize the neurons on the average, and this increases their firing rates. This instability is similar to the rate instability in networks of excitatory neurons (Hansel and Mato, 2003).

A detailed analysis of the QIF network at finite coupling strength shows that in a broad range of coupling strengths, the phase diagrams for two coupled neurons and for large networks are similar to those derived here, except for the additional rate instability line. A full report of this analysis will be published elsewhere.

Simulation results were performed here for a coupling conductance, *g* = 0.005 mS/cm^{2}. This corresponds to coupling coefficients (Amitai et al., 2002) CC ≈5% and to a spikelet size of ∼0.5 mV (Fig. 3). We have verified that similar results have been found for conductances four times larger (data not shown). Therefore the conclusions of the present work are relevant for electrical synapse conductances in the physiological range (Galarreta and Hestrin, 2001, 2002; Amitai et al., 2002).

### Related works

Stable antiphase locking has been demonstrated in previous studies for a pair of neurons coupled via electrical synapses. Sherman and Rinzel (1992) used numerical simulations to show that this nonintuitive phenomenon occurs in a simple conductance-based model. However, they did not relate it to the intrinsic properties of their model. Antiphase locking was found by Han et al. (1995) in simulations of coupled Morris–Lecar oscillators. For this model, they computed the effective phase interaction, Γ, and showed that their simulation results could be predicted from it. Our work further clarifies the conditions of stable antiphase locking in that it gives general qualitative criteria that the PRF must satisfy, which allows us to predict how it depends on intrinsic currents.

Phase locking of neurons connected by electrical synapses has been also investigated analytically (Chow and Kopell, 2000; Lewis and Rinzel, 2003). In both studies, neurons were modeled in the subthreshold range as leaky integrators, but the spikes were described in different ways. Chow and Kopell (2000) used a phenomenological model to describe the time course of spikes, whereas Lewis and Rinzel (2003), like us, modeled them with a δ-function. Both studies found stable antiphase locking at low firing rates and destabilization for large-enough firing rates. This can be understood using Equations 18–21 as follows. For a passive integrator, *Z*(ϕ), *Z*(ϕ) = 1/[*I*_{ext}*T*] exp(ϕ*T*/2πτ_{m}), where T is the firing rate, I_{ext} is the external current, and τ_{m} is the neuronal membrane time constant (Kuramoto, 1991; Hansel et al., 1995). The value of *Z* diverges exponentially with the period at all ϕ, as do *S*_{sub} and *S*_{r}. Detailed analysis reveals that the most divergent of these terms is *S*_{sub} < 0. Hence, for low-enough firing rates, Γ′(π) < 0. The term *S*_{sp} is positive (*Z*′(π) > 0) and decreases with T, whereas *S*_{sub} and *S*_{r} increase. Hence, for a large-enough frequency, Γ′(π) changes sign. Therefore for passive integrators, antiphase locking is stable at low rates and loses stability as the firing rate increases. In contrast, we have found that the QIF model behaves differently than the LIF except in the limit of very large -*V*_{r}/*V*_{T}. Therefore our results show that the predictions of Chow and Kopell (2000) and Lewis and Rinzel (2003) may be relevant only for neurons with strong potassium currents.

### Relationship to experiments

Electrical synapses connect fast spiking (FS) as well as low-threshold spiking (LTS) interneurons (Kawaguchi and Kubota, 1997; Gibson et al., 1999). Recently, the synchrony level of pairs of LTS and FS cells (characterized by the amplitude of the peak of their spikes cross-correlation) has been studied as a function of the firing rate *v* (Mancilla et al., 2002). It has been found that for FS cells, synchrony increases with *v* but decreases for LTS cells. Our results may explain this difference if one assumes that *I*_{NaP} is stronger in LTS than in FS cells and/or that potassium conductances are stronger in FS than in LTS cells. This assumption is compatible with the fact that FS and LTS cells have different firing patterns: in response to current injection, the action potentials of FS but not LTS cells are followed by a pronounced AHP, which may reflect stronger potassium conductances in FS than in LTS cells. A test of this explanation would be to show that the PRFs of these LTS and FS cells have different shapes.

The dynamic clamp technique (Sharp et al., 1993) makes it possible to artificially manipulate the intrinsic currents of the neurons and to couple two neurons via artificial electrical synapses. It provides a nice way of further testing our predictions. Therefore, it enables systematic verification of the effects of sodium and potassium currents on the synchronization properties of a pair of neurons (Fig. 4) and how they correlate with their PRF.

Several studies in the mouse and the rat have reported the presence of electrical synapses in hypoglossal and other inspiratory brainstem and phrenic motoneurons (Mazza et al., 1992; Rekling and Feldman, 1998). Recently, it has been found that inspiratory motoneuron synchrony is modulated by electrical synapses: blocking these synapses increases their level of synchrony (Bou-Flores and Berger, 2001). An explanation of the surprising findings of Bou-Flores and Berger (2001) is suggested by our work. Indeed, since strong I_{NaP} (Powers and Binder, 2003) is present in these neurons, electrical synapses can reduce the level of synchrony of their activity compared with the case in which electrical synapses are blocked. Synchrony in the latter situation would be attributable, for example, to GABAergic synaptic interactions between the neurons or to a spatially correlated time-dependent external input. As a matter of fact, we have verified that this effect occurs in numerical simulations of our network model with GABAergic synapses added. A detailed study of this phenomenon will be reported elsewhere.

## Appendix

### The *I–F* curve, the membrane potential, and the phase-response function of the QIF neuron

The subthreshold membrane potential, *v*(*t*), of a QIF neuron satisfies the dynamic equation:
24 Integrating this differential equation, one finds that its general solution is for *I*_{ext} > *I*_{c}:
25 where α is a constant of integration. The condition that at *t* = 0 the membrane potential of the neuron is at its reset value, *V*_{r}, determines α. One finds:
26 The function *v*(*t*) increases monotonously to infinity. Therefore, after some time *T, v*(*t*) reaches the threshold, *V*_{T}. At that time, the neuron fires an action potential, and *v*(*t*) is instantaneously reset to *V*_{r}. Subsequently, *v*(*t*) starts again to increase until it again reaches the threshold after another time *T*. Therefore, the neuron is firing periodically. The condition: *v*(*T*) = *V*_{T} determines the firing period. One finds:
27 The *I–F* curve of the QIF neuron is given by *v* = 1/*T*(*I*_{ext}). When the subthreshold membrane potential, *v*(*t*) reaches *V*_{T}, the neuron fires an action potential represented by a δ-function:
28 where θ measures the integral over time of the suprathreshold part of an action potential.

The response function *Z* is the phase resetting curve of the neuron (Kuramoto, 1984; Hansel et al., 1993; Rinzel and Ermentrout, 1998) in the limit of vanishing small perturbations of the membrane potential. For the case of one-dimensional models, it can be written as
(Kuramoto, 1991; Hansel et al., 1995). Using Equation 24, one finds for the QIF model:
29

### Phase interaction for electrical synapses

The phase coupling function, Γ(ϕ), between two identical neurons, *i* and *j*, interacting synaptically depends on the phase difference, ϕ_{i} - ϕ_{j}, of the two neurons and is given by (Kuramoto, 1984; Hansel et al., 1993, 1995; Kopell and Ermentrout, 2002):
30 where *Z* is the phase-response function of the neurons and *I*_{syn}(ψ_{i}, ψ_{j}) is the synaptic current coming from the presynaptic neuron j to the postsynaptic neuron i. For electrical synapses, this current is:
31

### Phase locking of a pair of weakly interacting neurons

When, in the absence of noise, two identical neurons interact weakly in a symmetric way, they reach a phase-locked state at large time. This is easily shown in the phase-reduction framework (Hansel et al., 1993, 1995; Van Vreeswick et al., 1994; Kopell and Ermentrout, 2002). The phases, ϕ_{1} and ϕ_{2}, of the two neurons satisfy the two coupled equations:
32
33 Subtracting these equations, one finds:
34 where Δ is the phase shift between the two neurons, Δ= ϕ_{1} - ϕ_{2}, and Γ_{-}(ϕ) = ^{1}_{2}(Γ(ϕ) = -Γ (-ϕ)). At large time, the phase shift, Δ, reaches a fixed point (i.e., it is such that *d*Δ/*d*t = 0) and the two neurons become locked with a phase shift that satisfies the condition:
35 Because of the symmetry of the system and the 2π periodicity of the function Γ_{-}, there are always at least two solutions to this equation. One is the in-phase solution, Δ = 0, and the other is the antiphase solution: Δ = π. Other solutions with 0 < Δ < π may also exist. However, only those that are stable can be reached at large time. Linearizing Equation 34 around these solutions, one finds that the stability condition is:
36 where combining Equations 30 and 31:
37 Differentiating this function with respect to ϕ, one finds:
38 The membrane potential *V* can be expanded in its components (spike, reset, and subthreshold component) as:
39 which allows us to rewrite Γ′(Δ):
40 Integrating by parts the first term and using the properties of the δ-function, we finally obtain:
41

*n* neurons: stability of the asynchronous state for the QIF network with noise

The stability analysis of the asynchronous state can be investigated for a network of phase oscillators in the presence of white noise (zero mean and variance *D*). For a general phase-coupling interaction, Γ(ϕ), one can show (Kuramoto, 1984) that the asynchronous state is stable if for all *n* > 0:
42 where we have defined:
43 with Γ_{n}, the *n*th-Fourier component of the 2π periodic function Γ(ϕ):
44

For a network of QIF neurons interacting with electrical synapses (Eqs. 24, 28) the *n*th-Fourier component of phase interaction has the form:
45 where *Z*_{n} and *v*_{n} are the *n*th-Fourier components of the functions *Z*(ϕ) and *v*(ϕ) defined in Equations 29 and 25. The variance of the noise, D, depends on the average of the phase-response function,
and is given by *D* = σ^{2}*Z*_{0}^{2} (Kuramoto, 1984), where σ^{2} is the variance of the white noise in the full QIF model. This yields the stability condition for the asynchronous state:
46 Note that *Z*_{0} depends on the parameters *V*_{r}, *V*_{T} and on the firing frequency of the neurons. Therefore, the desynchronizing effect of the noise depends on these parameters. In particular, because *Z*_{0} decreases when the firing rate increases, the noise is more efficient at stabilizing the asynchronous state at low firing rates than at high firing rates.

## Footnotes

This work was supported by North Atlantic Treaty Organization Physical and Engineering Science and Technology Collaborative Linkage Grant 977683, Les Programmes Internationaux de Coopération Scientifique–Centre National de la Recherche Scientifique (number 837), L'Action Concertée Incitative “Neurosciences integratives et computationnelles” (Ministère de la Recherche, France), and project A99E01 from the Comisión Asesora Científica de la Cooperación Argentino-Francesa. Research by D.G. was supported by Israel Science Foundation Grant 657/01 and by the National Science Foundation, under Agreement 0112050. We thank Y. Loewenstein, C. van Vreeswijk, B. Connors, and Y. Yarom for stimulating discussions and B. Connors, Y. Loewenstein, and C. Meunier for careful and critical reading of this manuscript.

Correspondence should be addressed to Dr. David Hansel, Laboratoire de Neurophysique et Physiologie du Système Moteur, 45 Rue des Saints-Pères, 75270 Paris Cedex 06, France. E-mail: david.hansel{at}biomedicale.univparis5.fr.

Copyright © 2003 Society for Neuroscience 0270-6474/03/236280-15$15.00/0