Impact of Adaptation Currents on Synchronization of Coupled Exponential Integrate-and-Fire Neurons

The ability of spiking neurons to synchronize their activity in a network depends on the response behavior of these neurons as quantified by the phase response curve (PRC) and on coupling properties. The PRC characterizes the effects of transient inputs on spike timing and can be measured experimentally. Here we use the adaptive exponential integrate-and-fire (aEIF) neuron model to determine how subthreshold and spike-triggered slow adaptation currents shape the PRC. Based on that, we predict how synchrony and phase locked states of coupled neurons change in presence of synaptic delays and unequal coupling strengths. We find that increased subthreshold adaptation currents cause a transition of the PRC from only phase advances to phase advances and delays in response to excitatory perturbations. Increased spike-triggered adaptation currents on the other hand predominantly skew the PRC to the right. Both adaptation induced changes of the PRC are modulated by spike frequency, being more prominent at lower frequencies. Applying phase reduction theory, we show that subthreshold adaptation stabilizes synchrony for pairs of coupled excitatory neurons, while spike-triggered adaptation causes locking with a small phase difference, as long as synaptic heterogeneities are negligible. For inhibitory pairs synchrony is stable and robust against conduction delays, and adaptation can mediate bistability of in-phase and anti-phase locking. We further demonstrate that stable synchrony and bistable in/anti-phase locking of pairs carry over to synchronization and clustering of larger networks. The effects of adaptation in aEIF neurons on PRCs and network dynamics qualitatively reflect those of biophysical adaptation currents in detailed Hodgkin-Huxley-based neurons, which underscores the utility of the aEIF model for investigating the dynamical behavior of networks. Our results suggest neuronal spike frequency adaptation as a mechanism synchronizing low frequency oscillations in local excitatory networks, but indicate that inhibition rather than excitation generates coherent rhythms at higher frequencies.


Introduction
Synchronized oscillating neural activity has been shown to be involved in a variety of cognitive functions [1,2] such as multisensory integration [3,4], conscious perception [5,6], selective attention [7,8] and memory [9,10], as well as in pathological states including Parkinson's disease [11], schizophrenia [12], and epilepsy [13]. These observations have led to a great interest in understanding the mechanisms of neuronal synchronization, how synchronous oscillations are initiated, maintained, and destabilized.
The phase response curve (PRC) provides a powerful tool to study neuronal synchronization [14]. The PRC is an experimentally obtainable measure that characterizes the effects of transient inputs to a periodically spiking neuron on the timing of its subsequent spike. PRC based techniques have been applied widely to analyze rhythms of neuronal populations and have yielded valuable insights into, for example, motor pattern generation [15], the hippocampal theta rhythm [16], and memory retrieval [10]. The shape of the PRC is strongly affected by ionic currents that mediate spike frequency adaptation (SFA) [17,18], a prominent feature of neuronal dynamics shown by a decrease in instanta-neous spike rate during a sustained current injection [19][20][21]. These adaptation currents modify the PRC in distinct ways, depending on whether they operate near rest or during the spike [18]. Using biophysical neuron models, it has been shown that a low threshold outward current, such as the muscarinic voltagedependent K z -current (I m ), can produce a type II PRC, characterized by phase advances and delays in response to excitatory stimuli, in contrast to only phase advances, defining a type I PRC. A high threshold outward current on the other hand, such as the Ca 2z -dependent afterhyperpolarization K z -current (I ahp ), flattens the PRC at early phases and skews its peak towards the end of the period [18,22,23]. Both changes of the PRC indicate an increased propensity for synchronization of coupled excitatory cells [22], and can be controlled selectively through cholinergic neuromodulation. In particular, I m and I ahp are reduced by acetylcholine with different sensitivities, which modifies the PRC shape [23][24][25].
In recent years substantial efforts have been exerted to develop single neuron models of reduced complexity that can reproduce a large repertoire of observed neuronal behavior, while being computationally less demanding and, more importantly, easier to understand and analyze than detailed biophysical models. Two-dimensional variants of the leaky integrate-and-fire neuron model have been proposed which take into consideration an adaptation mechanism that is spike triggered [26] or subthreshold, capturing resonance properties [27], as well as an improved description of spike initiation by an exponential term [28]. A popular example is the adaptive exponential leaky integrate-and-fire (aEIF) model by Brette and Gerstner [29,30]. The aEIF model is similar to the twovariable model of Izhikevich [31], such that both models include a sub-threshold as well as a spike-triggered adaptation component in one adaptation current. The advantages of the aEIF model, as opposed to the Izhikevich model, are the exponential description of spike initiation instead of a quadratic nonlinearity, and more importantly, that its parameters are of physiological relevance. Despite their simplicity, these two models (aEIF and Izhikevich) can capture a broad range of neuronal dynamics [32][33][34] which renders them appropriate for application in large-scale network models [35,36]. Furthermore, the aEIF model has been successfully fit to Hodgkin-Huxley-type neurons as well as to recordings from cortical neurons [29,37,38]. Since lately, this model is also implemented in neuromorphic hardware systems [39].
Because of subthreshold and spike-triggered contributions to the adaptation current, the aEIF model exhibits a rich dynamical structure [33], and can be tuned to reproduce the behavior of all major classes of neurons, as defined electrophysiologically in vitro [34]. Here, we use the aEIF model to study the influence of adaptation on network dynamics, particularly synchronization and phase locking, taking into account conduction delays and unequal synaptic strengths. First, we show how both subthreshold and spike-triggered adaptation affect the PRC as a function of spike frequency. Then, we apply phase reduction theory, assuming weak coupling, to explain how the changes in phase response behavior determine phase locking of neuronal pairs, considering conduction delays and heterogeneous synaptic strengths. We next present numerical simulations of networks which support the findings from our analysis of phase locking in neuronal pairs, and show their robustness against heterogeneities. Finally, to validate the biophysical implication of the adaptation parameters in the aEIF model, we relate and compare the results using this model to the effects of I m and I ahp on synchronization in Hodgkin-Huxley-type conductance based neurons. Thereby, we demonstrate that the basic description of an adaptation current in the low-dimensional aEIF model suffices to capture the characteristic changes of PRCs, and consequently the effects on phase locking and network behavior, mediated by biophysical adaptation currents in a complex neuron model. The aEIF model thus represents a useful and efficient tool to examine the dynamical behavior of neuronal networks.

aEIF neuron model
The aEIF model consists of two differential equations and a reset condition, if V §V cut then V~V r w~wzb: ð3Þ The first equation (1) is the membrane equation, where the capacitive current through the membrane with capacitance C equals the sum of ionic currents, the adaptation current w, and the input current I. The ionic currents are given by an ohmic leak current, determined by the leak conductance g L and the leak reversal potential E L , and a Na z -current which is responsible for the generation of spikes. The Na z -current is approximated by the exponential term, where D T is the threshold slope factor and V T is the threshold potential, assuming that the activation of Na zchannels is instantaneous and neglecting their inactivation [28]. The membrane time constant is t m :~C=g L . When I drives the membrane potential V beyond V T , the exponential term actuates a positive feedback and leads to a spike, which is said to occur at the time when V diverges towards infinity. In practice, integration of the model equations is stopped when V reaches a finite ''cutoff'' value V cut , and V is reset to V r (3). Equation (2) governs the dynamics of w, with the adaptation time constant t w . a quantifies a conductance that mediates subthreshold adaptation. Spike-triggered adaptation is included through the increment b (3).
The dynamics of the model relevant to our study is outlined as follows. When the input current I to the neuron at rest is slowly increased, at some critical current the resting state is destabilized which leads to repetitive spiking for large regions in parameter space [34]. This onset of spiking corresponds to a saddle-node (SN) bifurcation if at w vg L t m , and a subcritical Andronov-Hopf (AH) bifurcation if at w wg L t m at current values I SN and I AH respectively which can be calculated explicitly [33]. In the former case a stable fixed point (the neuronal resting state) and an unstable fixed point (the saddle) merge and disappear, in the latter case the stable fixed point becomes unstable before merging with the saddle. In the limiting case at w~gL t m , both bifurcations (SN and AH) meet and the system undergoes a Bogdanov-Takens (BT) bifurcation. The sets of points with dV =dt~0 and dw=dt~0 are called V -nullcline and w-nullcline, respectively. It is obvious that all fixed points in the two-dimensional state space can be identified as intersections of these two nullclines. Spiking can occur at a

Author Summary
Synchronization of neuronal spiking in the brain is related to cognitive functions, such as perception, attention, and memory. It is therefore important to determine which properties of neurons influence their collective behavior in a network and to understand how. A prominent feature of many cortical neurons is spike frequency adaptation, which is caused by slow transmembrane currents. We investigated how these adaptation currents affect the synchronization tendency of coupled model neurons. Using the efficient adaptive exponential integrate-and-fire (aEIF) model and a biophysically detailed neuron model for validation, we found that increased adaptation currents promote synchronization of coupled excitatory neurons at lower spike frequencies, as long as the conduction delays between the neurons are negligible. Inhibitory neurons on the other hand synchronize in presence of conduction delays, with or without adaptation currents. Our results emphasize the utility of the aEIF model for computational studies of neuronal network dynamics. We conclude that adaptation currents provide a mechanism to generate low frequency oscillations in local populations of excitatory neurons, while faster rhythms seem to be caused by inhibition rather than excitation.
constant input current lower than I SN or I AH depending on whether the sequence of reset points lies exterior to the basin of attraction of the stable fixed point. This means, the system just below the bifurcation current can be bistable; periodic spiking and constant membrane potential are possible at the same input current. Thus, periodic spiking trajectories do not necessarily emerge from a SN or AH bifurcation. We determined the lowest input current that produces repetitive spiking (the rheobase current, I rh ) numerically by delivering long-lasting rectangular current pulses to the model neurons at rest. Note that in general I rh depends on V r , such that in case of bistability, I rh can be reduced by decreasing V r [33].
We selected realistic values for the model parameters (C~0:1 nF, g L~0 :01 mS, E L~{ 70 mV, D T~2 mV, V T{ 50 mV, t w~1 00 ms, V r~{ 60 mV) and varied the adaptation parameters within reasonable ranges (a[½0,0:1 mS, b[½0,0:2 nA). All model parametrizations in this study lead to periodic spiking for sufficiently large I, possibly including transient adaptation. Parameter regions which lead to bursting and irregular spiking [34] are not considered in this study. V cut was set to {30 mV, since from this value, even without an input current, V would rise to a typical peak value of the action potential (v50 mV) within less than 1 ms while w essentially does not change due to its large time constant. Only in Fig. 1A-C we used V cut~2 0 mV to demonstrate the steep increase of V past V T .

Traub neuron model
In order to compare the effects of adaptation in the aEIF model with those of I m and I ahp in a biophysically detailed model and with previously published results [18,22,40] we used a variant of the conductance based neuron model described by Traub et al. [41]. The current-balance equation of this model is given by where the ionic currents consist of a leak current , and the slow K z -currents I m~gm v(V {E K ), and I ahpg The gating variables m, h and n satisfy first-order kinetics with a m~0 :32(V z54)=(1{exp({(V z54)=4)) and b m~0 :28 where , and the intracellular Ca 2z concentration ½Ca 2z is described by Units are mV for the membrane potential and ms for time. Note that the state space of the Traub model eqs. (4)-(9) is sixdimensional.
The dynamics of interest is described below. Starting from a resting state, as I is increased, the model goes to repetitive spiking. Depending on the level of I m , this (rest-spiking) transition occurs through a SN bifurcation for low values of I m or a subcritical AH bifurcation for high values of I m , at input currents I SN and I AH , respectively. The SN bifurcation gives rise to a branch of stable periodic solutions (limit cycles) with arbitrarily low frequency. Larger values of I m cause the stable fixed point to lose its stability by an AH bifurcation (at I AH vI SN ). In this case, a branch of unstable periodic orbits emerges, which collides with a branch of stable limit cycles with finite frequency in a fold limit cycle bifurcation at current I FLC vI AH . The branch of stable periodic spiking trajectories extends for currents larger than I AH and I SN . This means that in the AH bifurcation regime, the model exhibits hysteresis. That is, for an input current between I FLC and I AH a stable equilibrium point and a stable limit cycle coexist. On the contrary, I ahp does not affect the bifurcation of the equilibria, since it is essentially nonexistent at rest.

Network simulations
We considered networks of N coupled neurons with identical properties using both models (aEIF and Traub), driven to repetitive spiking with period T, where the vector x i consists of the state variables of neuron i for the Traub model), f governs the dynamics of the uncoupled neuron (according to either neuron model) and the coupling function h ij contains the synaptic current I syn (received by postsynaptic neuron i from presynaptic neuron j) in the first component and all other components are zero. I syn was modeled using a bi-exponential description of the synaptic conductance, where g ij denotes the peak conductance, s the fraction of open ion channels, d ij the conduction delay which includes axonal as well as dendritic contributions, and E syn the synaptic reversal potential. c is a normalization factor which was chosen such that the peak of s equals one. The spike times t j of neuron j (at the soma) correspond to the times at which the membrane potential reaches V cut (in the aEIF model) or the peak of the action potential (in the Traub model). t r and t d are the rise and decay time constants, respectively. For excitatory synapses the parameters were chosen to model an AMPA-mediated current (E syn~0 mV, t r~0 :1 ms, t d~1 ms), the parameters for inhbitory synapses we set to describe a GABA A -mediated current (E syn~{ 80 mV, t r0 :5 ms, t d~5 ms). We simulated the aEIF and Traub neuron networks, respectively, taking F :~T {1~4 0 Hz, homogeneous all-to-all connectivity without self-feedback (g ii~0 ), and neglecting conduction delays (d ij~0 ). We further introduced heterogeneities of several degrees w.r.t. synaptic strengths and conduction delays to the computationally less demanding aEIF network. Specifically, g ij (i=j) and d ij were sampled from a uniform distribution over various value ranges. The neurons were weakly coupled, in the sense that the total synaptic input received by a neuron from all other neurons in the network (assuming they spike synchronously) resulted in a maximal change of ISI (T) of less than 5%, which was determined by simulations. As initial conditions we used points of the spiking trajectory at times that were uniformly sampled from the interval ½0,T, i.e. the initial states were asynchronous. Simulation time was 20 s for each configuration of the aEIF networks and 10 s for the Traub neuron networks. All network simulations were done with BRIAN 1.3 [42] applying the secondorder Runge-Kutta integration method with a time step of 1 ms for coupled pairs and 10 ms for larger networks.
We measured the degree of spike synchronization in the simulated networks using averaged pairwise cross-correlations between the neurons [43], where s k i~1 if neuron i spikes in time interval k, otherwise s k i~0 , for k~1, . . . ,T k =t. h:i indicates the average over all neuronal pairs (i,j) in the network. Calculation period T k was 1 s and time bin t was 2:5 ms. k assumes a value of 0 for asynchronous spiking and approaches 1 for perfect synchronization.
In order to quantify the degree of phase locking of neurons in the network we applied the mean phase coherence measure s [44,45] where Q k ij is the phase difference between neurons i and j at the time of the k th spike t k is the largest spike time of neuron j that precedes t k i , t kz1 j is the smallest spike time of neuron j that succeeds t k i . K is the number of spikes of neuron i in the calculation period T K .
and h:i denotes the average over all pairs (i,j). s~0 means no neuronal pair phase locks, s~1 indicates complete phase locking. s was calculated using for T K the last 10 s (aEIF networks) or 5 s (Traub networks) of each simulation.

PRC calculation
The PRC can be obtained (experimentally or in simulations) by delivering small perturbations to the membrane potential of a neuron oscillating with period T at different phases q and calculating the change of the period. The PRC is then expressed as a function of phase as PRC(q)~T{T pert (q), where T pert (q) is the period of the neuron perturbed at q. Positive (negative) values of PRC(q) represent phase advances (delays). An alternative technique of determining the PRC is to solve the linearized adjoint equation [22,[46][47][48][49] subject to the normalization condition q(0) T f( x x(0))~1 (see Text S1 ). x, f are as described above (cf. eq. (10)) and D x f is the Jacobian matrix of f. x x denotes the asymptotically stable Tperiodic spiking trajectory as a solution of the system of differential equations and a reset condition in case of the aEIF model. Eq. (16) together with the reset condition describe the dynamics of an uncoupled neuron.
x x is an attractor of this dynamical system and nearby trajectories will converge to it. To obtain x x, we integrated the neuron model equations for a given set of parameters and adjusted the input current I, such that the period was T. Analysis was restricted to the regular spiking regime (cf. [34] for the aEIF model). Parameter regions where bursting and chaotic spiking occurs were avoided.
For Traub model trajectories, the peak of the action potential is identified with phase q~0, for aEIF trajectories q~0 corresponds to the point of reset. The first component q V of the normalized Tperiodic solution q of eq. (15) represents the PRC, also called infinitesimal PRC, which characterizes the response of the oscillator to a vanishingly small perturbation (cf. Text S1). For continuous limit cycles x x, as produced by the Traub model, q can be obtained by solving eq. (15) backward in time over several periods with arbitrary initial conditions. Since x x is asymptotically stable, the T-periodic solution of the adjoint system, eq. (15), is unstable. Thus, backward integration damps out the transients and we arrive at the periodic solution of eq. (15) [48][49][50]. In case of the aEIF model with an asymptotically stable T-periodic solution x x, that involves a discontinuity in both variables V V (t), w w(t) at integer multiples of T, we treated the adjoint equations as a boundary value problem [18]. Specifically, we solved the adjoint system subject to the conditions where q V ,q w denote the two components of q, and q w (T { ) :~lim t8T q w (t) is the left-sided limit. Eq. (19) is the normalization condition. Eq. (20) is the continuity condition, which ensures T-periodicity of the solution (see Text S1, derivation based on [51][52][53]). From the fact, that the end points of T-periodic aEIF trajectories differ, i.e.
, which in turn leads to q(0)=q(T { ). Perturbations of the same strength, which are applied to V just before and after the spike, have therefore a different effect on the phase, leading to a discontinuity in the PRC.
The PRCs presented in this study were calculated using the adjoint method. For validation purposes, we also simulated a number of PRCs by directly applying small perturbations to the membrane potential V V of the oscillating neuron at different phases and measuring the change in phase after many cycles -to ensure, that the perturbed trajectory had returned to the attractor x x. The results are in good agreement with the results of the adjoint method.

Phase reduction
In the limit of weak synaptic interaction, which guarantees that a perturbed spiking trajectory remains close to the attracting (unperturbed) trajectory x x, we can reduce the network model (10) to a lower dimensional network model where neuron i is described by its phase q i [48][49][50]54,55] as follows.
where q V i is the PRC of neuron i and V V i the first component (membrane potential) of the spiking trajectory x x i (see previous section and Text S1). H d ij is the T-periodic averaged interaction function calculated using I syn with conduction delay d ij (11). Note that d ij simply causes a shift in the interaction function: . H d ij only depends on the difference of the phases (in the argument) which is a useful property when analyzing the stability of phase locked states of coupled neuronal pairs. In this case (without self-feedback as already assumed) the phase difference Q :~q 2 {q 1 evolves according to the scalar differential equation whose stable fixed points are given by the zero crossingsQ Q of H D for which lim e:0 dH D (Q Q{e)=dQv0 and lim e:0 dH D (Q Qze)=dQ v0. If H D is differentiable atQ Q, these left and right sided limits are equal and represent the slope. Note however that H D is continuous, but not necessarily differentiable due to the discontinuity of the PRC of an aEIF neuron. Therefore, the limits might not be equal in this case. The case where H D is discontinuous atQ Q, which can be caused by d-pulse coupling, i.e. I syn is replaced by a d-function, is addressed in the Results section. We calculated these stable fixed points, which correspond to stable phase locked states, for pairs of identical cells coupled with equal or heterogeneous synaptic strengths and symmetric conduction delays, d :~d 12~d21 , using PRCs derived from the aEIF and Traub neuron models, driven to 40 Hz periodic spiking. Periodic spiking trajectories of both models and PRCs of Traub neurons were computed using variable order multistep integration methods, for PRCs of aEIF neurons a fifth-order collocation method was used to solve eqs. (17)- (20). These integration methods are implemented in MATLAB (2010a, The MathWorks). Bifurcation currents of the Traub model were calculated using MATCONT [56,57].

PRC characteristics of aEIF neurons
We first examine the effects of the adaptation components a and b, respectively, on spiking behavior of aEIF neurons at rest in response to (suprathreshold) current pulses ( Fig. 1A-C). Without adaptation (a~b~0) the model produces tonic spiking (Fig. 1A). Increasing a or b leads to SFA as shown by a gradual increase of the inter spike intervals (ISI) until a steady-state spike frequency F is reached. Adaptation current w builds up and saturates slowly when only conductance a is considered (Fig. 1B) in comparison to spike-triggered increments b (Fig. 1C). Fig. 1D,E depicts the relationship between F and the injected current I for various fixed values of a and b. Increased subthreshold adaptation causes the minimum spike frequency to jump from zero to a positive value, producing a discontinuous F -I curve (Fig. 1D). A continuous (discontinuous) F -I curve indicates class I (II) membrane excitability which is typical for a SN (AH) bifurcation at the onset of spiking respectively. An increase of a causes this bifurcation to switch from SN to AH, thereby changing the membrane excitability from class I to II, shown by the F -I curves. An increase of b on the other hand does not produce a discontinuity in the F -I curve, i.e. the membrane excitability remains class I (Fig. 1E). Furthermore, increasing a shifts the F -I curve to larger current values without affecting its slope, while an increase of b decreases the slope of the F -I curve in a divisive manner. When b is large, the neuron is desensitized in the sense that spike frequency is much less affected by changes in the driving input.
In Fig. 2A,B we show how a and b differentially affect the shape of the PRC of an aEIF neuron driven to periodic spiking. The PRCs calculated using the adjoint method (solid curves) match well with those obtained from simulations (circles). While nonadapting neurons have monophasic (type I) PRCs, which indicate only advancing effects of excitatory perturbations, increased levels of a produce biphasic (type II) PRCs with larger magnitudes, which predict a delaying effect of excitatory perturbations received early in the oscillation cycle. An increase of b on the other hand flattens the PRC at early phases, shifts its peak towards the end of the period and reduces its magnitude. The type of the PRC however remains unchanged (type I). Indeed, if a~0 the PRC must be type I, since in this case the component q V of the solution of the adjoint system, eqs. (17)-(20), can be written as where c(s) is given by the right-hand side of eq. (17). Thus, q V cannot switch sign.
To provide an intuitive explanation for the effects of adaptation on the PRC, we show the vector fields, V -and w-nullclines, and periodic spiking trajectories of four aEIF neurons (Fig. 2C-F). One neuron does not have an adaptation current (a~b~0), two neurons possess only one adaptation mechanism (a~0:1 mS, b~0 nA and a~0 mS, b~0:2 nA, respectively) and for one both adaptation parameters are increased (a~0:1 mS, b~0:2 nA). An excitatory perturbation to the non-adapting neuron at any point of its trajectory, i.e. at any phase, shifts this point closer to V cut along the trajectory, which means the phase is shifted closer to T, hence the advancing effect (Fig. 2C). The phase advance is strongest if the perturbing input is received at the position along the trajectory around which the vector field has the smallest magnitude, i.e. where the trajectory is ''slowest''. In case of subthreshold adaptation (Fig. 2D), the adapted periodic spiking trajectory starts at a certain level of w which decreases during the early part of the oscillation cycle and increases again during the late part, after the trajectory has passed the w-nullcline. A small transient excitatory input at an early phase pushes the respective point of the trajectory to the right (along the V -axis) causing the perturbed trajectory to pass through a region above the unperturbed trajectory, somewhat closer to the fixed point around which the vector field is almost null. Consequently, the neuron is slowed down and the subsequent spike delayed. An excitatory perturbation received at a later phase (to the right of the dashed arrow) causes phase advances, since the perturbed trajectory either remains nearly unchanged, however with a shorter path to the end of the cycle, compared to the unperturbed trajectory, or it passes below the unperturbed one where the magnitude of the vector field (pointing to the right) is larger. Note that for the parametrization in Fig. 2D, both, the resting state as well as the spiking trajectory are stable. In this case, a strong depolarizing input at an early phase can push the corresponding trajectory point into the domain of attraction of the fixed point, encircled by the dashed line in the figure, which would cause the resulting trajectory to spiral towards the fixed point and the neuron would stop spiking. On the other hand, increasing I would shrink the domain of attraction of the fixed point and at I~I AH , it would be destabilized by a subcritical AH bifurcation. When a~0 and bw0, we obtain a type I PRC (Fig. 2E), as explained above. The advancing effect of an excitatory perturbation is strongest late in the oscillation cycle, indicated by the red arrow, where the perturbation pushes a trajectory point from a ''slow'' towards a ''fast'' region closer to the end of the cycle, as shown by the vector field. When a as well as b are increased, the PRC exhibits both adaptation mediated features (type II and skewness), see Fig. 2F. A push to the right along the corresponding trajectory experienced early in the cycle brings the perturbed trajectory closer to the fixed point and causes a delayed next spike. Such an effect persists even if the fixed point has disappeared due to a larger input current. In this case, the region where the fixed point used to be prior to the bifurcation, known as ''ghost'' of the fixed point, the vector field is still very small. This means that type II PRCs can exist for larger input currents IwI SN . Note that differences of the vector fields and the shift of the nullclines relative to each other in Fig. 2C,D as well as Fig. 2E,F are due to different input current values (as an increase of I moves the V -nullcline upwards). The maximal phase advances, indicated by solid arrows in Fig. 2A,B, are close to the threshold potential V T (where the Vnullcline has its minimum) in all four cases.
We next investigate how the changes in PRCs caused by either adaptation component are affected by the spike frequency. Bifurcation currents, rheobase currents and corresponding frequencies, in dependence of a and b, as well as regions in parameter space where PRCs are type I and II, are displayed in Fig. 3A-D. Fig. 3E,F shows how individual PRCs are modulated by spike frequency (input current). Both PRC characteristics, caused by a and b, respectively, are more pronounced at low frequencies. Increasing I changes a type II PRC to type I and shifts its peak towards an earlier phase. The input current which separates type I and type II PRC regions (in parameter space) increases with both, a and b (Fig. 3A,B). That is, an increase of b can also turn a type I into a type II PRC, by bringing the spiking trajectory closer to the fixed point or its ''ghost''. This is however only possible if the system is in the AH bifurcation regime (awC=t w ) or close to it. Spike-triggered adaptation thereby considerably influences the range of input currents for which the PRCs are type II. The spike frequency according to the input current, at which a type II PRC turns into type I increases substantially with increasing a, but only slighly with an increase of b (Fig. 3C,D). The latter can be recognized by the similarity of the respective (green) curves in the subfigures C and D. Type II PRCs thus only exist in the lower frequency band whose width increases with increasing subthreshold adaptation.

Phase locking of coupled aEIF pairs
In this section, we examine how the changes in phase response properties due to adaptation affects phase locking of coupled pairs of periodically spiking aEIF neurons. Specifically, we first analyze how the shape of the PRC determines the fixed points of eq. (23) and their stability, and then show how the modifications of the PRC mediated by the adaptation components a and b change those fixed points. Finally, we investigate the effects of conduction delays and heterogeneous coupling strengths on phase locking in dependence of adaptation.
Relation between phase locking and the PRC. In case of identical cell pairs and symmetric synaptic strengths, g 12~g21 , the interaction functions in eq. (23) are identical, H d   Note that the region in I-a space where the PRCs are type II is very shallow in A compared to B, the corresponding regions in F -a space shown in C and D however are rather similar. This is due to the steep (flat) F -I relationship for b~0 nA (b~0:2 nA) respectively (see Fig. 1D which in turn is equivalent to PRC(T { )wPRC(0) for excitatory coupling and PRC(T { )vPRC(0) for inhibitory coupling. Hence, it is the discontinuity of the PRC which determines the stability ofQ Q in this case. A synaptic current with finite rise and decay times causes an additional rightwards shift and a smoothing of the interaction function. The stability of the fixed pointQ Q is then determined by the slope of the PRC and its discontinuity on the interval (d{Q Q, d{Q Qze), where ew0 is on the order of the synaptic timescale. If the PRC slope is negative on this interval and its discontinuity (if occurring in the interval) is also negative, i.e. PRC(T { )wPRC(0), thenQ Q is stable for excitatory coupling and unstable for inhibitory coupling. In Fig. 4A we show the effect of the synaptic timescale, i.e. t r and t d , on the interaction function for a given PRC. Fig. 4B,C illustrates how the stability of the synchronous state of a neuronal pair is given by the slope of the PRC, for three different delays. The slope of the PRC is positive at q~d z 1 , q~d 2 and negative at q~d 3 and remains positive (negative) until I syn has decayed to a small value. Therefore, synchrony is unstable for delays d 1 , d 2 and stable for d 3 , indicated by the slope of H d at Q~0, which is negative for the first two and positive for the third delay.
Effects of adaptation on phase locking of coupled aEIF pairs. First, consider pairs of identical aEIF neurons with the PRCs shown in Fig. 2A,B, symmetrically coupled through instantaneous synapses (t r :0 and t d :0) and without conduction delays (d~0). When the coupling is excitatory, the in-phase locked state (synchrony) is unstable in case of type I PRCs, since they have a positive ''jump'' at q~0, i.e. PRC(T { )vPRC(0). Synchrony is stable for pairs with type II PRCs however, as PRC(T { )wPRC(0). The anti-phase locked state on the other hand is unstable because of the positive PRC slopes at q~T=2. In case of inhibitory coupling, synchrony is stable for type I pairs and the anti-phase locked state is stable for all pairs. This means, bistability of in-phase and anti-phase locking occurs for inhibitory neurons with type I PRCs.
Next, we consider pairs that are coupled through synaptic currents I syn with finite rise and decay times, as described in the Methods section. In Fig. 5 we show how the stable (and unstable) phase locked states of pairs of neurons with symmetric excitatory (A, B) and inhibitory (C, D) synaptic interactions and without conduction delays change, when the PRCs are modified by the adaptation components a and b. For excitatory pairs, stable fixed points shift towards synchrony, when a or b is increased. The phase differences become vanishinly small, when the PRCs switch from type I to type II due to subthreshold adaptation. Perfect synchrony is stabilized, where the PRC slopes at q~e for small ew0 become negative, due to even larger values of a (not shown) or lower spike frequency (see Fig. 3C-F). Neurons that have type I PRCs with a pronounced skew, as caused by spike-triggered adaptation, lock almost but not completely in-phase, if the adaptation is sufficiently strong. Inhibitory pairs on the other hand show stable synchrony independent of PRC type and skewness. Larger values of a or b lead to additional stabilization of the anti-phase locked state. That is, strong adaptation in inhibitory pairs mediates bistability of in-phase and anti-phase locking. All phase locking predictions from the phase reduction approach are in good agreement with the results of numerically simulated coupled aEIF pairs.
Phase locking of aEIF pairs coupled with delays. We next investigate how phase locked states of excitatory and inhibitory pairs are affected by synaptic currents that involve conduction delays, considering the PRC of a neuron without adaptation, and two PRCs that represent adaptation induced by either a or b. Neurons symmetrically coupled through excitatory synapses with a conduction delay do not synchronize irrespective of whether adaptation is present or not (Fig. 6A-C). Instead, stable states shift towards anti-phase locking with increasing mutual delays. Inhibitory pairs on the other hand synchronize for all conduction delays (Fig. 6D-F), but the anti-phase locked states of coupled inhibitory neurons with type II PRCs or skewed type I PRCs are destabilized by the delays. The bistable region is larger in case of spike-triggered adaptation compared to subthreshold adaptation (Fig. 6E,F). Again, all stable phase locked states obtained using phase reduction are verified by numerical simulations. Fig. 7 illustrates the phenomenon that synchronous spiking of excitatory pairs is destabilized by the delay, while obtained for synaptic conductances with three different sets of synaptic time constants: t r~0 :01 ms, t d~0 :1 ms (blue), t r~0 :25 ms, t d~2 :5 ms (green); t r~0 :75 ms, t d~7 :5 ms (magenta), and d~0. The synaptic current I syn associated with each pair of time constants (center) illustrates the three synaptic timescales relative to the period T~25 ms. Note that I syn shown here is received by the neuron at the beginning of its ISI. B: PRC (solid black) of an aEIF neuron spiking at 40 Hz and excitatory synaptic currents I syn with t r~0 :1 ms, t d~1 ms (dashed blue) received at three different phases. Assuming the input comes from a second, synchronous neuron, these phases represent three different conduction delays d 1~0 ms, d 2~1 0 ms, and d 3~2 0 ms. Note that synaptic input received at an earlier phase causes a larger peak of I syn , due to the smaller value V of the membrane potential which leads to a larger difference E syn {V to the synapse's reversal potential E syn . C: Interaction functions H di (Q) for pairs of neurons with the PRC shown in B, coupled by excitatory synapses with t r~0 :1 ms, t d~1 ms, and delays d 1 ,d 2  synchrony remains stable for inhibitory pairs. Consider two neurons oscillating with a small phase difference Q~q 1 {q 2 w0 (neuron 1 slightly ahead of neuron 2). Then, a synaptic input received by neuron 2 at a delay QvdvT=2 after neuron 1 has spiked, arrives at an earlier phase (q 2~d {Q) compared to the phase at which neuron 1 receives its input (q 1~d zQ). Consequently, if the synapses are excitatory and the PRCs type I, the leader neuron 1 advances its next spike by a larger amount than the follower neuron 2 (Fig. 7A). In case of excitatory neurons and type II PRCs, depending on Q and d, the phase of neuron 1 is advanced by a larger amount or delayed by a smaller amount than the phase of neuron 2, the latter of which is shown by the changed spike times in Fig. 7B. It is also possible that the phase of the leader neuron is advanced while that of the follower neuron is delayed. Hence, for either PRC type, Q increases due to delayed excitatory coupling, that is, synchrony is destabilized. For inhibitory synapses and type I PRCs, the leader neuron 1 delays its subsequent spike by a larger amount than the follower neuron 2 (Fig. 7C). In case of type II PRCs, neuron 1 experiences a weaker phase advance or stronger phase delay than neuron 1, or else the phase of neuron 1 is delayed while that of neuron 1 is advanced, depending on Q and d (Fig. 7D). Thus, delayed inhibitory coupling causes Q to decrease towards zero for either PRC type, that is, synchrony is stabilized.
Phase locking of aEIF pairs coupled with delays and unequal synaptic strengths. In the following we analyze phase locking of neuronal pairs with unequal synaptic peak conductances g 12 =g 21 . Due to the linearity of the integral in eq. (21) we can substitute H d ij~: g ijH H d ij in eq. (23), which yields By setting eq. (25) to zero, we obtain the condition eq. (26) for the existence of phase locked states, Phase locked states therefore only exist if the ratio of conductances g 12 =g 21 is not larger than the maximum of the periodic function . For a type II PRC on the other hand, this minimum is zero (unless the negative lobe of the PRC is small and the synapse slow), from which follows that max Q R(Q)??. The effects of heterogeneous synaptic strengths on phase locking of neuronal pairs without adaptation, as well as either adaptation parameter increased, are shown in Fig. 8. For excitatory pairs coupled without a conduction delay it is illustrated, how the right hand side of eq. (25) changes when the coupling strengths are varied (A-C). In addition, stable phase locked states of excitatory and inhibitory pairs coupled through synapses with various mutual conduction delays (d~0, 3, or 6 ms) are displayed as a function of g 12 =g 21 (D-I). When the ratio of conductances g 12 =g 21 is increased, the zero crossings of dQ=dt given by eq. (25), i.e. phase locked states, disappear for neurons with type I PRCs (through a SN bifurcation). Q then continuously increases (or decreases) (mod T) as shown by the dashed curves (without roots) in Fig. 8A,C and indicated by the arrows in Fig. 8D,F,G,I. This means, the spike frequency of one neuron becomes faster than that of the other neuron. Neurons with type II PRCs on the other hand have stable phase locked states even for diverging coupling strengths. Bistability of two phase locked states can occur for a ratio g 12 =g 21 close to one (equal coupling strengths), depending on the PRC and the delay. Synchronization of excitatory-inhibitory pairs is not considered in this paper. It should be noted however, that if both neurons have type I PRCs, phase locking is not possible, irrespective of the ratio of coupling strengths. In this case, one interaction function is strictly positive and the other strictly negative and thus, the condition (26) for fixed points of eq. (25) cannot be fulfilled.

Synchronization and clustering in aEIF networks
In order to examine how the behavior of pairs of coupled phase neurons relates to networks of spiking neurons, we performed numerical simulations of networks of oscillating aEIF neurons without adaptation and with either a subthreshold or a spiketriggered adaptation current, respectively, and analyzed the network activity. The neurons were all either excitatory or inhibitory and weakly coupled. Fig. 9 shows the degree of synchronization k (A, C) and the degree of phase locking s (B) for these networks considering equal as well as heterogeneous conduction delays and synaptic conductances. An increase of either adaptation parameter (a or b) leads to increased k in networks of excitatory neurons with short delays. It can be recognized however, that k increases to larger values and this high degree of synchrony seems to be more robust against heterogeneous synaptic strengths, when the neurons are equipped with a subthreshold adaptation current (Fig. 9A,C). These effects correspond well to those of the adaptation components a and b on synchronization of pairs, presented in the previous section. Parameter regimes (w.r.t. a,b,d ij and g ij ) that cause stable in-phase or near in-phase locking of pairs, such as subthreshold adaptation in case of short delays or spike-triggered adaptation for short delays and coupling strength ratios close to one (Fig. 6A-C and Fig. 8D-F), lead to synchronization, indicated by large k values, in the respective networks. Networks of non-adapting excitatory neurons remain asynchronous as shown by the low k values. For equal synaptic strengths, these networks settle into splay states where the neurons are pairwise phase locked, with uniformly distributed phases (Fig. 9B,D). When the delays are large enough and the synaptic strengths equal, splay states also occur in networks of neurons with large b, indicated by low k and high s values in Fig. 9A,B. As far as inhibitory networks are concerned, non-adapting neurons synchronize, without delays or with random delays of up to 10 ms. Furthermore, synchrony in these networks is largely robust against heterogeneities in the coupling strengths (Fig. 9A). Networks of inhibitory neurons with subthreshold adaptation only show synchronization and pairwise locking for larger delays (i.e. d ij random in ½0,5 ms or larger). Spike-triggered adaptation promotes clustering of the network into two clusters, where the neurons within a cluster are in synchrony, as long as the delays are small. These cluster states seem to be most robust against heterogeneous synaptic strengths when the delays are small but not zero. For larger delays, inhibitory neurons of all three types (with or without adaptation) synchronize, in a robust way against unequal synaptic strengths. The behaviors of inhibitory networks are consistent with the phase locked states found in pairs of inhibitory neurons (Fig. 6D-F). Particularly, stable synchronization of pairs with larger conduction delays and the bistability of inphase and anti-phase locking of pairs with spike-triggered adaptation for smaller delays, nicely carry over to networks. In the former case, synchrony of pairs relates to network synchrony, in the latter case, bistability of in-phase and anti-phase locking of individual pairs can explain the observed two cluster states. Note that bistability of in-phase and anti-phase locking is also shown for inhibitory pairs with subthreshold adaptation and d~0 ms. In this case however, the slope of H D (Q) at Q~T=2 is almost zero (not shown), which might explain why the corresponding networks do not develop two-cluster states. The behavior of all simulated networks does not critically depend on the number of neurons in the network, as we obtain qualitatively similar results for network sizes changed to N~50 and N~200 (not shown). The numerical simulations demonstrate that stable phase locked states of neural pairs can be used to predict the behavior of larger networks. To understand the biophysical relevance of the subthreshold and spike-triggered adaptation parameters, a and b, in the aEIF model, we compare them with the adaptation currents I m and I ahp in a variant of the Hodgkin-Huxley type Traub model neuron. Specifically, in this section we investigate the effects of the low-and high-threshold currents I m and I ahp , respectively, on spiking behavior, F -I curves and PRCs of single neurons, and on synchronization of pairs and networks, using the Traub model, and compare the results with those of the previous two sections. It should be stressed, that the aEIF model was not fit to the Traub model in this study. Therefore, the comparison of how adaptation currents affect SFA, PRCs and synchronization in both models, are rather qualitative than quantitative.
PRC characteristics of Traub neurons. Without adaptation, g m~gahp~0 (hence I m~Iahp~0 ), the model  (25), as a function of Q for pairs of excitatory aEIF neurons coupled with different ratios of synaptic conductances g 12 =g 21 (d~0). Zero crossings with a negative slope indicate stable phase locking and are marked by black dots. Adaptation parameters of the neurons and PRCs are shown in the top row. D-I: Stable phase locked states of excitatory (D-F) and inhibitory (G-I) pairs as a function of the synaptic conductance ratio, for three different conduction delays d~0, 3 and 6 ms (black, brown, green). Unstable states are not shown for improved clarity. Dashed lines denote equal synaptic strengths, grey arrows indicate a continuous increase or decrease of Q (mod T) for ratios g 12 =g 21 at which phase locked states do not exist (see main text). doi:10.1371/journal.pcbi.1002478.g008 exhibits tonic spiking in response to a rectangular current pulse (Fig. 10A). When either adaptation current is present, that is the conductance g m or g ahp is increased to 0:1 mS, the membrane voltage trace reveals SFA. Note that I ahp causes stronger differences in subsequent ISIs after stimulus onset, when comparing the V -traces of neurons with either adaptation conductance set to 0:1 mS. The F-I curves in Fig. 10B indicate that the presence of I m predominantly has a subtractive effect on the neuron's F -I curve and gives rise to class II excitability. The presence of I ahp on the other hand flattens the F -I curve, in other words its effect is divisive. Furthermore, an increase of I m changes a type I PRC to type II, whereas increased I ahp reduces its amplitude at early phases and skews its peak to the right (Fig. 10C).
Evidently, the effects of I m and I ahp on SFA, F -I curves and PRCs of Traub neurons are consistent with the effects of the adaptation parameters a and b in aEIF neurons (Figs. 1, 2).
We further show how the PRC characteristics caused by the adaptation currents depend on the injected current I, hence the spike frequency F , and the bifurcation type of the rest-spiking transition (Fig. 10D-I). An increase of I reduces the effects of I m and I ahp on the PRC. That means, at higher frequencies F , larger levels of I m and I ahp are required to obtain type II and skewed PRCs, respectively. This frequency dependence of adaptation current-mediated changes of the PRC is similar in both neuron models (Figs. 3, 10D-I). Note, that in the Traub model a rather low value of g m (25 nS) is sufficient to guarantee a type II PRC for spike frequencies of up to 100 Hz (Fig. 10F,G), compared to the aEIF model, where a much larger value of a (w0:1 mS) would be necessary (Fig. 3C,D). As far as the bifurcation structures of both models are concerned, an increase of the low-threshold adaptation parameters g m and a has a comparable effect in the Traub and the aEIF models, respectively, changing the transition from rest to spiking from a SN via a BT to an AH bifurcation. The exact conductance values at which this change, i.e. the BT bifurcation, occurs, differ (g m~0 :02 mS for the Traub model and a~0:001 mS for the aEIF model).
Synchronization of coupled Traub neurons. We show the effects of the adaptation currents I m and I ahp on phase locked states of pairs of Traub neurons symmetrically coupled without conduction delays in Fig. 11A-D. Excitatory pairs of neurons without adaptation phase lock with a small phase difference. Low levels of I m are sufficient to stabilize in-phase locking, by turning the PRC from type I to II (Fig. 11A), while an increase of I ahp reduces the locked phase difference to almost but not exactly zero, that is, near in-phase locking, by skewing the PRC (Fig. 11B). Inhibitory synaptic coupling produces bistability of in-phase (synchrony) and anti-phase locking (anti-synchrony) for pairs of neurons without adaptation or either adaptation current increased (Fig. 11C,D). Note that the domain of attraction of the antisynchronous state grows with increasing I m or I ahp , while that of the synchronous state shrinks. In contrast to the aEIF model, this bistability also occurs for neurons without an adaptation current (compare Figs. 5C,D, 11C,D). The effects of I m or I ahp on synchronization of networks of Traub neurons coupled without conduction delays and equal synaptic strengths, are shown in Fig. 11E,F. In correspondence with the effects on pairs, I m and I ahp promote synchronization of excitatory networks, shown by the course of network synchronization measure k over time (Fig. 11E). The mean values of phase locking measure s are 0.26 for nonadapting neurons and 0.98 for networks where either adaptation current is increased. An increased adaptation current I m leads to larger k values, compared to an increase of I ahp , which is similar to the aEIF networks where increased a causes larger k values than an increase of b (compare Figs. 9C, 11E). In contrast to networks of excitatory aEIF neurons without adaptation, which develop splay states, k values of nonadapting excitatory Traub neuron networks increase to about 0.5, while low s values indicate poor phase locking, hence splay states do not occur (Fig. 11F). Networks of inhibitory neurons organize into clusters, indicated by k values that converge to 0.5 (Fig. 11E) and large s values (0.96 without adaptation, 0.94 for either I m or I ahp increased). Particularly, clustering into two clusters was revealed by the raster plots, see Fig. 11F. These twocluster states of networks can be explained by the bistability of synchrony and anti-synchrony of individual pairs. Clustering emerges for all three types of Traub neurons, with and without adaptation, as opposed to networks of inhibitory aEIF neurons, Figure 11. Influence of adaptation on synchronization properties of Traub model neurons. A-D: Stable (solid black) and unstable (dashed grey) phase locked states of coupled pairs of Traub neurons with identical PRCs, as a function of conductances g m and g ahp , respectively. Corresponding changes in PRCs are displayed in the top row. The neurons are coupled through excitatory or inhibitory synapses as indicated by the diagrams on the left, with equal synaptic strengths, g 12~g21 and d~0. E: Network synchronization k over time, of N~50 coupled excitatory (solid) and inhibitory (dashed) Traub neurons without, g m~gahp~0 mS (black) or with adaptation, g m~0 :1 mS, g ahp~0 mS (blue) and g m~0 mS, g ahp~0 :2 mS (red), driven to 40 Hz spiking. The neurons are all-to-all coupled with equal synaptic conductances, g ij~0 :06 nS (black and blue), g ij~0 :18 nS (red), but without self-feedback, g ii~0 , and conduction delays, where cluster states only occur in case of spike-triggered adaptation (Fig. 9). Considering the collective behavior of coupled excitatory neurons, the synchronizing effects of I m and I ahp in the Traub model are comparable to those of the adaptation components a and b in the aEIF model.

Discussion
In this work we studied the role of adaptation in the aEIF model as an endogenous neuronal mechanism that controls network dynamics. We described the effects of subthreshold and spiketriggered adaptation currents on the PRC in dependence of spike frequency. To provide insight into the synchronization tendencies of coupled neurons, we applied a common phase reduction technique and used the PRC to describe neuronal interaction [48,55]. For pairs of coupled oscillating neurons we analyzed synchrony and phase locking under consideration of conduction delays and heterogeneous synaptic strengths. We then performed numerical simulations of aEIF networks to examine whether the predicted behavior of coupled pairs relates to the activity of larger networks. Finally, to express the biophysical relevance of the elementary subthreshold and spike-triggered adaptation mechanisms in the aEIF model, we compared their effects with those of the adaptation currents I m and I ahp in the high-dimensional Traub neuron model, on single neuron as well as network behavior.
Conductance a, which mostly determines the amount of adaptation current in absence of spikes, that is, subthreshold, qualitatively changes the rest-spiking transition of an aEIF neuron, from a SN to an AH via a BT bifurcation as a increases. Thereby the neuron's excitability, as defined by the F -I curve, and its PRC, are turned from class I to class II, and type I to type II, respectively. A similar effect of a slow outward current that acts in the subthreshold regime on the PRC has recently been shown for a two-dimensional quadratic non-leaky integrate-and-fire (QIF) model derived from a normal form of a dynamical model that undergoes a BT bifurcation [18,48]. The relation between the PRC and the bifurcation types has further been emphasized by Brown et al. [47] who analytically determined PRCs for bifurcation normal forms and found type I and II PRC characteristics for the SN and AH bifurcations, respectively. A spike-triggered increment b of adaptation current does not affect the bifurcation structure of the aEIF model and leaves the excitability class unchanged. When a is small such that the model is in the SN bifurcation regime, an increase of b cannot change the PRC type. In the AH bifurcation regime, b substantially affects the range of input current for which the PRC is type II but causes only a small change in the corresponding frequency range. Furthermore, spike-triggered adaptation strongly influences the skew of the PRC, shifting its peak towards the end of the ISI for larger values of b. Such a right-skewed PRC implies that the neuron is most sensitive to synaptic inputs that are received just before it spikes. Similar effects of spike-triggered negative feedback with slow decay on the skew of the PRC have been reported for an extended QIF model [18,22,48,58].
PRCs determine synchronization properties of coupled oscillating neurons. When the synapses are fast compared to the oscillation period, the stability of the in-phase and anti-phase locked states (which always exist for pairs of identical neurons) can be ''read off'' the PRC for any mutual conduction delay, as we have demonstrated. A similar stability criterion that depends on the slopes of the PRCs at the phases at which the inputs are received has recently been derived for pairs of pulse-coupled oscillators [59]. Under the assumption of pulsatile coupling, the effect of a synaptic input is required to dissipate before the next input is received. In principle, the synaptic current can be strong, but it must be brief such that the perturbed trajectory returns to the limit cycle before the next perturbation occurs [14].
We have shown that, as long as synaptic delays are negligible and synaptic strengths equal, excitatory pairs synchronize if their PRCs are type II, as caused by a, and lock almost in-phase if their PRCs are type I with a strong skew, as mediated by b. Inhibitory pairs synchronize in presence of conduction delays and show bistability of in-phase and anti-phase locking for small delays, particularly in case of skewed PRCs. Conduction delays and synaptic time constants can affect the stability of synchrony in a similar way, by producing a lateral shift of the interaction function H d (Q), as shown in Fig. 4. Note however, that the synaptic timescale has an additional effect on the shape of H d (Q), smoothing it for slower synaptic rise and decay times. We have further demonstrated that heterogeneity in synaptic strengths desynchronizes excitatory and inhibitory pairs and leads to phase locking with a small phase difference in case of type II PRCs and small delays. While neurons with type II PRCs have stable phase locked states even for large differences in synaptic strengths, pairs of coupled neurons with type I PRCs are only guaranteed to phase lock when the synaptic strengths are equal. Similar effects of heterogeneous synaptic conductances have recently been observed in a computational study of weakly coupled Wang-Buszaki and Hodgkin-Huxley neurons (with class I and II excitability, respectively) [60].
The activity of larger aEIF networks, simulated numerically, is consistent with the predictions of the behavior of pairs. In fact, knowledge on phase locking of coupled pairs helps to explain the observed network states. Both adaptation mediated PRC characteristics, i.e. a negative lobe or a pronounced right skew, favor synchronization in networks of excitatory neurons, in agreement with previous findings [17,22,61]. This phenomenon only occurs when the conduction delays are negligible. It has been shown previously that synchrony in networks of excitatory oscillators becomes unstable when considering coupling with delays [62,63]. We have demonstrated that increased conduction delays promote asynchrony in excitatory networks, with or without adaptation currents. Inhibitory neurons on the other hand are able to synchronize spiking in larger networks for a range of conduction delays. This provides support to the hypothesis that inhibitory networks play an essential role in generating coherent brain rhythms, as has been proposed earlier [43,64], [2] for review. Inhibition rather than excitation has been found to generate neuronal synchrony particularly in case of slow synaptic rise and decay [40,61,65], and in the presence of conduction delays as has recently been shown experimentally [66]. In regimes that lead to bistability of in-phase and anti-phase locking according to our analysis of pairs, the simulated networks break up into two clusters of synchronized neurons. Recently it has been shown that a stable two cluster state of pulse coupled neural oscillators can exist even when synchrony of individual pairs is unstable [67]. Such cluster states have been invoked to explain population rhythms measured in vitro, where the involved neurons spike at about half of the population frequency [68].
Spike frequency has been shown to affect the skewness of PRCs, using type I integrate-and-fire neurons with adaptation [58], and to modulate the negative lobe in type II PRCs of conductance based model neurons [45]. Using the aEIF model we have demonstrated that the spike frequency strongly attenuates the effect of either adaptation mechanism on the PRC. At high frequency, unphysiologically large adaptation parameter values are necessary to produce a negative lobe or a significant right-skew in the PRC. This means, for a given degree of adaptation in excitatory neurons, synchronization is possible at frequencies up to a certain value. The stronger the adaptation, the larger this upper frequency limit. It has been previously suggested that the degree of adaptation can determine a preferred frequency range for synchronization of excitatory neurons, based on the observation (in vitro and in silico) that the neurons tend to spike in phase with injected currents oscillating at certain frequencies [69]. This preferred oscillation frequency increases with increasing degree of SFA. According to our results, at low frequencies synchronization of local circuits through excitatory synapses is possible, provided that the neurons are adapting and delays are short. At higher frequencies, adaptation much less affects the synchronization tendency of excitatory neurons and inhibition may play the dominant role in generating coherent rhythms [43,64].
The adaptation currents I m and I ahp have previously been found to influence the phase response characteristics of the biophysical Traub neuron model, turning a type I PRC to type II (through I m ) and modulating its skew (through I ahp ) [18,22]. We have shown that these changes of the PRC are reflected in the aEIF model by its two adaptation parameters and that in both models (aEIF and Traub) these changes are modulated by the spike frequency. As a consequence, the adaptation induced effects on synchronization of pairs and networks of oscillating neurons are qualitatively similar in both models. Quantitative differences with respect to these effects may well be reduced by fitting the aEIF model parameters to Traub neuron features.
Our analysis of phase locked states is based on the assumption that synaptic interactions are weak. Experimental work lending support to this assumption has been reviewed in [14,50]. Particularly for stellate cells of the entorhinal cortex, synaptic coupling has been found to be weak [70]. Another assumption in this study is that the neurons spike with the same frequency. Considering a pair of neurons spiking at different frequencies, equation (23) needs to be augmented by a scalar v, which accounts for the constant frequency mismatch between the two neurons [71]: dQ=dt~vzH d 21 ({Q){H d 12 (Q). In this case, the condition for the existence of phase locked states is D(Q) :~H d 21 ({Q){H d 12 (Q)~v. Due to the assumption of weak synaptic strengths however, max Q DD(Q)D must be small, which means that the above condition can only be met if v is small. In other words, in the limit of weak coupling phase locking is only possible if the spike frequencies are identical or differ only slightly. The phase reduction technique considered here, and PRCs in general, are of limited applicability for studying network dynamics in a regime where individual neurons spike at different frequencies, or even irregularly. How adaptation currents affect network synchronization and rhythm in such a regime nevertheless remains an interesting question to be addressed in the future.

Supporting Information
Text S1 Supplementary Methods. A) Calculation of the PRC using the adjoint method. B) Phase reduction. (PDF)