Paradoxical phase response of gamma rhythms facilitates their entrainment in heterogeneous networks

The synchronization of different $\gamma$-rhythms arising in different brain areas has been implicated in various cognitive functions. Here, we focus on the effect of the ubiquitous neuronal heterogeneity on the synchronization of PING (pyramidal-interneuronal network gamma) and ING (interneuronal network gamma) rhythms. The synchronization properties of rhythms depends on the response of their collective phase to external input. We therefore determined the macroscopic phase-response curve for finite-amplitude perturbations (fmPRC), using numerical simulation of all-to-all coupled networks of integrate-and-fire (IF) neurons exhibiting either PING or ING rhythms. We show that the intrinsic neuronal heterogeneity can qualitatively modify the fmPRC. While the phase-response curve for the individual IF-neurons is strictly positive (type I), the fmPRC can be biphasic and exhibit both signs (type II). Thus, for PING rhythms, an external excitation to the excitatory cells can, in fact, delay the collective oscillation of the network, even though the same excitation would lead to an advance when applied to uncoupled neurons. This paradoxical delay arises when the external excitation modifies the internal dynamics of the network by causing additional spikes of inhibitory neurons, whose delaying within-network inhibition outweighs the immediate advance caused by the external excitation. These results explain how intrinsic heterogeneity allows the PING rhythm to become synchronized with a periodic forcing or another PING rhythm for a wider range in the mismatch of their frequencies. We demonstrate a similar mechanism for the synchronization of ING rhythms. Our results identify a potential function of neuronal heterogeneity in the synchronization of coupled $\gamma$-rhythms, which may play a role in neural information transfer via communication through coherence.


Author Summary
The interaction of a large number of oscillating units can lead to the emergence of a collective, macroscopic oscillation in which many units oscillate in near-unison or near-synchrony. This has been exploited technologically, e.g., to combine many coherently interacting, individual lasers to form a single powerful laser. Collective oscillations are also important in biology. For instance, the circadian rhythm of animals is controlled by the near-synchronous dynamics of a large number of individually oscillating cells. In animals and humans brain rhythms reflect the coherent dynamics of a large number of neurons and are surmised to play an important role in the communication between different brain areas. To be functionally relevant, these rhythms have to respond to external inputs and have to be able to synchronize with each other. We show that the ubiquitous heterogeneity in the properties of the individual neurons in a network can contribute to that ability. It can allow the external inputs to modify the internal network dynamics such that the network can follow these inputs over a wider range of frequencies. Paradoxically, while an external perturbation may delay individual neurons, their ensuing within-network interaction can overcompensate this delay, leading to an overall advance of the rhythm.

Introduction
Collective oscillations or rhythms representing the coherent dynamics of a large number of coupled oscillators play a significant role in many systems. In the technological realm they range from laser arrays and Josephson junctions to micromechanical oscillators [1,2]. Among the important biological examples are the heart rhythm, the circadian rhythm generated by the suprachiasmatic nucleus [3], the segmentation clock controlling the somite formation during development [4], and brain waves [5]. One prominent brain rhythm is the widely observed γ-rhythm with frequencies in the range 30-100Hz. The coherent spiking of the neurons underlying this rhythm likely enhances the downstream impact of the neurons participating in the rhythm. The rhythmic alternation of low and high activity has been suggested to play a significant role in the communication between different brain areas [6,7]. That communication has also been proposed to be controled by the coherence of the rhythms in the participating brain areas [8][9][10][11][12][13]. For collective oscillations or rhythms to play a constructive role in a system they need to respond adequately to external perturbations and stimuli. For instance, for the circadian rhythm it is essential that it can be reliably entrained by light and phase-lock to its daily variation. Similarly, if rhythms are to play a significant role in the communication between different brain areas, their response to input from other areas represents a significant determinant of their function. Moreover, the stimulation and entrainment of γ-rhythms by periodic sensory input is being considered as a therapeutic approach for some neurodegenerative diseases [14]. Even small perturbations can affect oscillations significantly in that they can advance or delay the oscillations, i.e. they can change the phase of the oscillators. This change typically depends not only on the strength of the perturbation but, importantly, also on the timing of the perturbations and is expressed in terms of the phase response curve (PRC), which has been studied extensively for individual oscillators [15]. For infinitesimal perturbations the PRC can be determined elegantly using the adjoint method [16]. If the collective oscillation of a network of interacting oscillators is sufficiently coherent, that system can be thought of as a single effective oscillator. Consequently, the response of the macroscopic phase of the collective oscillation to external perturbations and the mutual interaction of multiple collective oscillations is of interest. The macroscopic phase-response curve (mPRC) has been obtained in various configurations, including noise-less heterogeneous phase oscillators [17,18], noisy identical phase oscillators [19,20], noisy excitable elements [21], and noisy oscillators described by the thetamodel [22], which is equivalent to the quadratic integrate-fire model for spiking neurons. Recent work has used the reduction of networks of quadratic integrate-fire neurons to two coupled differential equations for the firing rate and the mean voltage [23], which is related to the Ott-Antonsen theory [24,25], to develop a method to obtain the infinitesimal macroscopic PRC (imPRC) for excitatoryinhibitory spiking networks [26,27].
A key difference between the response of an individual oscillator to a perturbation and that of a collective oscillation is the fact that the degree of synchrony of the collective oscillation can change as a result of the perturbation, reflecting a change in the relations between the individual oscillators. Thus, the phase response of a collective oscillation to a brief perturbation consists not only of the immediate change in the phases of the individual oscillators caused by the perturbation, but includes also a change in the collective phase that can result from the subsequent convergence back to the phase relationship between the oscillators corresponding to the synchronized state, which is likely to have been changed by the perturbation [18]. Interestingly, it has been observed that the infinitesimal macroscopic phase response can be qualitatively different from the phase response of the individual elements. Thus, even if the individual oscillators have a type-I PRC, i.e. a PRC that is strictly positive or negative, the mPRC of the collective oscillation can be of type II, i.e. it can exhibit a sign change as a function of the phase [21,22,28]. Here we investigate the interplay between external perturbations and the internal interactions among neurons in inhibitory and in excitatory-inhibitory networks exhibiting γ-rhythms of the ING-and of the PING-type. We focus on networks comprised of neurons that are not identical, leading to a spread in their individual phases and a reduction in the degree of their synchrony. How does this phase dispersion affect the response of the macroscopic phase of the rhythm to perturbations? Does it modify the ability of the network to follow a periodic perturbation ? We show that the dispersion in the phase together with the within-network interactions among the neurons can be the cause of a paradoxical phase response: an external perturbation that delays each individual neuron can advance the macroscopic rhythm. We identify the following mechanism underlying this paradoxical response: external perturbations that delay individual neurons sufficiently allow the within-network inhibition generated by early-spiking neurons to suppress the spiking of less excited neurons. This results in a reduced within-network inhibition, which reduces the time to the next spike volley, speeding up the rhythm. This paradoxical phase response increases with the neuronal heterogeneity and allows the network to phase-lock to periodic external perturbations over a wider range of detuning. Thus, the desynchronization within the network enhances its synchronizabilty with other networks. The mechanism is closely related to that underlying the enhancement of synchronization of collective oscillations by uncorrelated noise [29] and the enhanced entrainment of the rhythm of a homogeneous network to periodic input if that input exhibits phase dispersion across the network [30,31]. We demonstrate and analyze these behaviors for networks of inhibitory neurons (ING-rhythm) and for networks comprised of excitatory and inhibitory neurons (PING-rhythm).

Results
We investigated the impact of neuronal heterogeneity on the response of the phase of γ-rhythms to brief external perturbations and the resulting ability of rhythms to synchronize to periodic input. As described in the Methods, we used networks comprised of minimal integrate-fire neurons that interact with each other through synaptic pulses modeled via delayed double-exponentials. To study ING-rhythms all neurons were inhibitory, while for the PING-rhythms we used excitatory-inhibitory networks. In both cases, the coupling within each population was all-to-all. Throughout, we implemented the neuronal heterogeneity by injecting a different steady bias current I bias into each neuron. Our analysis suggests that the origin of the neuronal heterogeneity plays only a minor role as long as it leads to a dispersion of their spike times [29].

Paradoxical Phase Response of Heterogeneous Networks: ING-Rhythm
In the absence of external perturbations the all-to-all inhibition among the neurons lead to rhythmic firing of the neurons. Due to their heterogeneity they did not spike synchronously but sequentially, as shown in Fig.1A, where the neurons are ordered by the strength of their bias current. The dependence of the phase dispersion on the coefficient of variation of the heterogeneity in the bias current (CV) is shown in Suppl. Figure S1. For sufficiently large heterogeneity some neurons never spiked: while the weak bias current they received would have been sufficient to induce a spike eventually, the strong inhibition that was generated by the neurons spiking earlier in the cycle suppressed those late spikes. Neurons with strong bias current could spike multiple times. A brief, inhibitory external input delivered to all neurons (green dashed line in Fig.1B) delayed each neuron. The degree of this individual delay depended on the timing of the input, as is reflected in the PRC of the individual neurons. If the perturbation was applied during the time between the spike volleys, the delay of each neuron had no further consequence and the overall rhythm was delayed. However, if the same inhibitory perturbation arrived during a spike volley (dashed green line in Fig.1B), it could advance the overall rhythm. As illustrated in Fig.1B, only the spiking of the late neurons was delayed by the perturbation. Importantly, with this delay some neurons did not spike before the within-network inhibition triggered by the early-spiking neurons (dashed blue line in Fig.1B) became strong enough to suppress the spiking of the late neurons altogether. With fewer neurons spiking, the  all-to-all inhibition within the network was reduced, allowing all neurons to recover earlier, which lead to a shorter time to the next spike volley. If the speed-up was larger than the immediate delay induced by the external inhibition, the overall phase of the rhythm was advanced by the delaying inhibition.
As the example in Fig.1B shows, the paradoxical phase response requires proper timing of the perturbation. We therefore determined quantitatively the macroscopic phase-response curve (PRC) of the rhythm. To do so we measured computationally the amount a brief current injection shifted the phase of the rhythm ( Fig.2A). We defined the phase as the normalized time since the first spike in the most recent volley of spikes. Reflecting the strictly positive PRC of the individual integrate-fire neurons, without heterogeneity (CV = 0) external inhibition always delayed the rhythm, independent of the timing of the pulse. In contrast, in heterogeneous networks the rhythm could be advanced if the same inhibitory perturbation was applied shortly after the first spikes in the spike volley (φ inh > 0). Increasing the neuronal heterogeneity enhanced this phase advance, since it shifted the within-network inhibition driven by the leading neurons to earlier times, while it delayed the lagging neurons. As a result, for the same external perturbation, a larger fraction of neurons that would spike in the absence of the external inhibition was sufficiently delayed to have their spikes be suppressed by the withinnetwork inhibition (cf. Fig.1B), reducing the within-network inhibition and with it the time to the next spike volley. To keep the frequency of the unperturbed network fixed in Fig.2A, we reduced the tonic input with increasing heterogeneity, which enhanced the phase advance. However, even if the tonic input was kept fixed, the phase advance increased with heterogeneity (Suppl. Fig.S2).
For weak heterogeneity the paradoxical phase response occurred only for sufficiently strong perturbations, i.e. it did not arise in the infinitesimal macroscopic PRC (imPRC). Thus, the phase response changed qualitatively as the amplitude of the perturbation was strong enough to delay the spikes of sufficiently many slow neurons until the self-inhibition of the network set in and suppressed their spikes (Fig.2B). As the CV of the neurons was increased, the dispersion was large enough that the spikes of the lagging neurons were suppressed by the self-inhibition of the network even in the absence of an external perturbation. Above that threshold value of CV the paradoxical phase response occurred even for infinitesimal perturbations (Fig.2C).
The paradoxical phase response was robust with respect to changes in the natural frequency of the network, the coupling strength, and the effective synaptic delay, as long as the rhythm persisted. The paradoxical phase advance increased with decreasing natural frequency of the network, since . With decreasing τ I 1 , a subharmonic peak emerged and eventually the rhythm disintegrated. Parameters as in D. In (A), (B) and (D), perturbations were made with a square-wave inhibitory current pulse with duration 0.1 ms to each interneuron. In (A) and (B), the amplitude of the current was 1600 pA, resulting in a 2 mV rapid hyperpolarization. In (D), the amplitude of the current was 400 pA, resulting in a 0.5 mV rapid hyperpolarization.  the inhibition had a stronger effect for lower mean input strength (Fig.3A). Changing the withinnetwork coupling strength by a factor of 2 up or down did not substantially affect the paradoxical phase response (Fig.3B) nor the strength of the rhythm (Fig.3C). Even without explicit synaptic delay (τ d = 0), the effective delay given by the double-exponential synaptic interaction was sufficient to render a paradoxical response (Fig.3D). However, when this effective delay was reduced by decreasing the rise time τ I 1 of the synaptic current, the rhythm itself developed a strong subharmonic component and eventually disintegrated (Fig.3E). In the subharmonic regime the paradoxical phase advance alternated in consecutive cycles of the rhythm.
In [13,27] the exact reduction of all-to-all coupled heterogeneous networks of quadratic integrate-fire neurons to 2 coupled ordinary differential equations for each network [23] has been used to obtain the infinitesimal macroscopic phase-response curve (imPRC) for ING and PING networks. They obtained biphasic response only if the excitatory perturbation was applied to the population of inhibitory neurons; for perturbations to the excitatory neurons they found only monophasic response (type-I). This is presumably due to the lack of a delay in the single-exponential synaptic interactions used in [13,27].

Enhancing entrainment of ING-rhythms through network heterogeneity
In order to allow communication by coherence [11,32], the rhythms in different brain areas need to be sufficiently phase-locked with each other. As a simplification of two interacting γ-rhythms, we therefore investigated the ability of the rhythm in a network to be entrained by a periodic external input, particularly focusing on the possibly facilitating role of neuronal heterogeneity. Motivated by the paradoxical phase response induced by the heterogeneity, we addressed, in particular, the question whether an ING network can be sped up by inhibition to entrain it with a faster network.
The network considered here was the same as that used to analyze the fmPRC. The within-network interaction was an all-to-all inhibition with synaptic delay τ d , resulting in a rhythm with natural frequency f natural , Each neuron received heterogeneous input I bias and inhibitory periodic pulses with frequency f clock . The latter can be considered as the output of another ING-network and were, in fact, generated that way (Fig.4). We refer to this external input as the 'clock'. The detuning ∆f = f clock − f natural was a key control parameter.
For periodic input the fmPRC allows the definition of an iterated map describing the response of the network. For periodic δ-pulses that map is shown in Fig.5A. For positive detuning, i.e. when the clock is faster than the network, entrainment requires that the phase response is paradoxical in order for the rhythm to be sped up by the inhibition. If the heterogeneity and the resulting phase response are sufficiently large, the maximum of the iterated map crosses the diagonal, generating a stable and an unstable fixed point. The former is the desired entrained state.
As the detuning is increased the iterated map is shifted downward. This can decrease the slope of the iterated map at the fixed point below -1, destabilizing the fixed point in a period-doubling bifurcation. Hz. The distance between the diagonal and subdiagonal line represents the detuning between the network and periodic input. In (A), the fmPRC was determined for a δ-pulse perturbation, in (B) for a double-exponential inhibitory current (cf. (2,3)) was used as in Fig.6.
For periodic pulses comprised of double-exponential inhibitory currents (cf. (2,3)) a rich bifurcation scenario emerged (Fig.5B). Note that the full map is not continuous and not unimodal (cf. first bottom panel of Fig.5B). Nevertheless, for ∆f < 7.17 Hz the attractor remains near the unstable fixed point and displays a period-doubling cascade to chaos and multiple periodic windows. For ∆f > 7.28 Hz, however, the attractor includes points on both sides of the discontinuity (cf. third bottom panel in Fig.5B).
Having clarified the role of the fmPRC in the network's synchronizability and ability to phase-lock, we investigated the role of neuronal heterogeneity in more detail (Fig.6). To do that, we adjusted for each value of the input heterogeneity the mean input strength I (I) so as to keep the natural frequency f network constant (f network = 44 Hz). Then we determined the extent of synchronization and phaselocking of the network under the influence of periodic inhibitory input as a function of the detuning ∆f and network heterogeneity CV . As shown above, the fmPRC of a heterogeneous network could be biphasic with the amplitude of the paradoxical phase response increasing with neuronal heterogeneity. Expecting that for sufficiently large heterogeneity an ING-rhythm could be accelerated by a faster periodic inhibition, we tested phase-locking predominantly for positive detuning, corresponding to f clock > f network . We first investigated how neuronal heterogeneity affected the synchronization by comparing the dominant frequency f dom in the Fourier spectrum of the network's LFP with f clock . In Fig.6A, the color hue indicates the ratio f dom : f clock . For small heterogeneity, f dom was a rational multiple of f clock that depended on the detuning, while for sufficiently large CV the network became synchronized in the sense that f dom = f clock (yellow). The range of ∆f allowing synchronization became wider with increasing neuronal heterogeneity, implying that the neuronal heterogeneity enhanced the synchronization of the ING-rhythm. However, note that f dom = f clock did not imply a perfectly  Δf=f clock -f network (Hz) synchronized or a 1:1 phase-locked state. In fact, various different subharmonic responses arose: example 2 shows a period-4 state, while in example 3 the dynamics were actually chaotic (Fig.6A) even though f dom = f clock . Motivated by these observations, we divided the states into three types: • Type 1: f dom = f clock , not synchronized, not phase-locked ( Fig.6 example 4).
The phase diagram Fig.6A does not differentiate between types 2 and 3. It only shows that neuronal heterogeneity enhanced the synchronization of the network by shifting f dom to f clock . Therefore, we studied whether neuronal heterogeneity also enhanced the synchronization by weakening the subharmonic response and changing the synchronized state from type 2 to type 3, as well as whether the dynamics of the fmPRC shown in the bifurcation diagram Fig.5B could predict the phase relationship between the network and the clock. Using the same simulation setup as in Fig.6A, the subharmonic response is shown in Fig.6B. The color hue indicates the ratio f sub : f clock , where f sub is the frequency of the dominant peak of the LFP power spectrum that satisfies f sub < f clock . The color saturation gives the ratio of the powers at f sub and f clock (capped at 1). Thus, over most of the range of positive detuning and neuronal heterogeneity tested, the fading-away of the color with increasing heterogeneity reveals that the neuronal heterogeneity weakened the subharmonic response. Over a small range of positive detuning, increasing neuronal heterogeneity from small values induced perfect synchronization (type 3) by weakening the subharmonic response with frequency ratio f sub : f clock = 1 : 2; the system traversed a continuous period-doubling bifurcation in reverse, with type 2 (red) giving way to type 3 (white). Together with Fig.6A, this showed that neuronal heterogeneity could enhance the synchronization both by making f dom = f clock (from type 1 to type 2) and by weakening the subharmonic response (from type 2 to type 3). The range of detuning where increasing heterogeneity induced a type 3 synchronization became wider for larger synaptic delay within the network (Suppl. Figure S3). Note that the bifurcation diagram (Fig.5B) based on the fmPRC agrees well with the subharmonic response marked along the dashed line at CV = 0.1 in Fig.6B, suggesting that the fmPRC can well predict the subharmonic response and persistent phase response of the network. In addition to enhancing the frequency synchonization, neuronal heterogeneity was also able to increase the tightness of the phase-locking. Over most of the parameter regime investigated, the variance of the phase of the network relative to the periodic input (var(Φ inh )) decreased with increasing heterogeneity, as indicated by the decrease in the color saturation in Fig.6C. In fact, for detuning between 0 Hz and 2 Hz the heterogeneity reduced var(Φ inh ) to 0 (white), corresponding to the 1:1 phase-locked state. Even for the 1:2 phase-locked state (cf. the red area in Fig.6B) var(Φ inh ) was very small for a range of heterogeneity and detuning (2 Hz to 4 Hz), indicating tight phase locking. Except for type-3 synchronized states the size of the spike volleys varied between clock cycles. In fact, over wide ranges of the parameters the network did not spike in each of the clock cycles, as indicated by the color hue in Fig.6C, which gives the fraction of cycles with no network spikes (e.g. Fig.6 example 4).

Paradoxical phase response and entrainment of PING rhythms
Many γ-rhythms involve not only inhibitory neurons, but arise from the mutual interaction of excitatory (E) and inhibitory (I) neurons (PING rhythm) [33]. The key elements to obtain a paradoxical phase response and the ensuing enhanced synchronization are self-inhibition within the network, neuronal heterogeneity and effective synaptic delay. Since in PING rhythms the connections from E-cells to I-cells and back to the E-cells form an effective self-inhibiting loop, we asked whether PING-rhythms can exhibit behavior similar to the behavior we identified for ING-rhythms. Considering a PING-rhythm generated by an E-I network comprised of integrate-fire neurons, we first studied its fmPRC. To avoid that all I-cells receive identical input and therefore spike as a single unit, the I-cells received, in addition to the excitation from the E-cells, heterogeneous, tonic, Gaussiandistributed subthreshold input with mean I (I) = 36 pA and CV (I) = 0.167. The phase response of the network was probed by applying an identical external excitatory perturbation to all E-cells and recording the resulting phase shift of the LFP (cf. eqs. (7,8)) of the E-population, averaged across 500 realizations of the subthreshold input to the I-cells (Fig.7A). More specifically, the perturbations consisted of a square-wave excitatory current pulse with amplitude 76 pA and duration 0.1 ms to each E-cell, resulting in a 2 mV rapid depolarization. Without neuronal heterogeneity the external excitation always advanced the phase of the rhythm resulting in an fmPRC that was strictly positive. In the heterogeneous case, however, the PING rhythm exhibited a paradoxical phase response, whereby the collective rhythm was delayed while the individual neurons were advanced by the excitation. The delay was caused by the increase of self-inhibition within the network that was generated by the additional spikes in the E-population, which in turn drove additional spikes in the I-population. In contrast to the fmPRC of the ING-rhythm, this paradoxical phase response was not monotonic in the heterogeneity. While weak heterogeneity resulted in strong delay, the delay decreased with increasing intermediate CV-values and only increased again for larger CV (Fig.7B left top). This non-monotonicity arose because we kept the frequency of the network constant as we increased its heterogeneity. This required a decrease in the tonic input to the E-cells with increasing heterogeneity. For the stronger tonic input used for weak heterogeneity the same external perturbation elicited more additional spikes than it did for strong heterogeneity where the tonic input was weaker (cf. titles of subpanels of Fig.7B). The total number of spikes occurring in each cycle of the unperturbed network also decreased with increasing heterogeneity. Consequently, the relative change in the number of spikes induced by the perturbation was non-monotonic in the heterogeneity. As a result, the relative change in the inhibitory synaptic conductance resulting from the perturbation and with it the phase delay was also non-monotonic. As for the ING rhythm, we investigated the role of neuronal heterogeneity in the synchonizability and the ability of phase-locking of coupled PING rhythms. In analogy to the ING-case, we considered the case of the E-population of a PING network receiving periodic excitation generated by a clock PING network (Fig.4B). As before, we adjusted the tonic input strength to the E-population to keep the natural frequency of the network constant as we changed its heterogeneity (f network = 41Hz).
To probe the impact of the paradoxical phase response on the synchronization we focused on negative detuning for which the periodic external excitation needed to slow down the network in order to achieve phase-locking. Indeed, with increasing heterogeneity the network could become synchronized with the slower clock over a larger ranger of the detuning as indicated by the fading saturation of the color in Fig.7C. Here the color hue indicates the ratio f super : f clock , where f super was determined as the frequency with the most power among the frequencies higher than f clock in the Fourier spectrum of the E-population's LFP. The color saturation indicates the ratio of the power at the frequencies f super and f clock . Thus, a color hue closer to green (f super : f clock = 1 : 1) or with a lower saturation implies better synchronization. By observing how the width of the range of detuning allowing synchronization varied with neuronal heterogeneity, we concluded that generally, the synchronizability of PING rhythm was enhanced by the neuronal heterogeneity by weakening subharmonic response. Note that for CV ∈ [0, 0.1] the synchonizability of the PING rhythm decreased slightly with neuronal heterogeneity. This was consistent with the nonmonotonicity exhibited by the fmPRC seen in Fig.7B. The neuronal heterogeneity played a similar role in the tightness of the phase-locking as in the synchronizability (Fig.7D).

Discussion
In this paper we have analyzed the response of collective oscillations of inhibitory and of excitatoryinhibitory networks of integrate-fire neurons to external perturbations. For ING-and PING-rhythms we have shown that the combination of neuronal heterogeneity and effective synaptic delay can qualitatively change the phase response compared to the phase response of the individual neurons generating the rhythm. Thus, perturbations that delay the I-cells can paradoxically advance the ING-rhythm and perturbations that advance the E-cells can delay the PING-rhythm. As a result, the macroscopic phase-response curve for finite-amplitude perturbations (fmPRC) of the rhythm can change sign as the phase of the perturbation is changed (type-II), even though the PRC of all individual cells is strictly positive (type-I). This change of the fmPRC enhances the ability of the γ-rhythm to synchronize with other rhythms.
The key element of the mechanism driving the paradoxical phase response and the enhanced synchronization is the cooperation of the external perturbation and the effectively delayed within-network inhibition. In the ING-network a suitably timed external perturbation delays the lagging -but not the early -neurons sufficiently to allow the within-network inhibition triggered by the early neurons to keep the lagging neurons from spiking. This reduces the overall within-network inhibition and with it the duration of the cycle. Thus, the perturbation modifies the internal dynamics of the rhythm, which leads to changes in the phase of the rhythm that can dominate the immediate phase change the perturbation induces. The situation is somewhat similar to that investigated in [18]. There it had been pointed out that an external perturbation of a collective oscillation can lead to changes in its phase in two stages: i) an immediate change of the phases of all oscillators as a direct result of the perturbation and ii) a subsequent slower change in the collective phase resulting from the convergence of the disturbed phases back to the synchronized state. That analysis was based on a network of phase oscillators and could therefore not include a key element of our results, which is the perturbation-induced change in the number of neurons that actually spike and the resulting change in the within-network inhibition that results in a change of the period of the rhythm. As discussed in [31,34], for ING-rhythms such a change in the number of spiking neurons underlies also the enhanced phase-locking found in [30]. Going beyond ING-rhythms, we showed that PING-rhythms can also exibit a paradoxical phase response via a mechanism that is analogous to that of ING-rhythms. For that analysis we have focused on excitatory-inhibitory networks with only connections between but not within the excitatory and inhibitory populations. For excitatory inputs to the excitatory cells to generate a paradoxical phase response it is necessary that the additional spikes of the excitatory neurons that are caused by the external perturbation induce additional spikes of the inhibitory neurons. This behavior arises if the inibitory population is also allowed to be heterogeneous. Moreover, the within-network inhibition has to be strong enough to be able to suppress the spiking of lagging excitatory neurons. This is, e.g., found in mice piriform cortex, where principal neurons driven by sensory input from the olfactory bulb arriving early during a sniff recruit inhibitory interneurons via long-range recurrent connections, resulting in the global, transient suppression of subsequent cortical activity [35]. A characteristic feature of the paradoxical phase response of the PING rhythm is the extended cycle time following enhanced activation of the excitatory cells. A strong such correlation between the cycle time and the previous LFP amplitude has been observed for the γ-rhythm in hippocampus [36]. To assess whether this rhythm exhibits paradoxical phase response would require comparing the macroscopic phase response [37] with that of indvidual participating neurons.
In order for the global perturbation to affect the various neurons differently, they have to be at different phases in their cycle. Our analysis suggests that the specific cause for this heterogeneity in the spike times does not play an important role. Indeed, as shown in [29], even fluctuating heterogeneities that are generated by noise rather than static heterogeneities reflecting intrinsic properties of neurons can enhance the synchronization of multiple γ-rhythms in interconnected networks of identical neurons. Note that the noise driving this synchronization is uncorrelated across neurons. The analysis of that system revealed the same mechanism at work as the one identified here. In various previous analytical and computational analyses it has been found that the dynamics of the macroscopic phase of a collective oscillation can qualitatively differ from that of the microscopic phase. Thus, for interacting groups of noisy identical phase oscillators the macroscopic phases of the groups can tend to lign up with each other, even if all pair-wise interactions between individual oscillators prefer the antiphase state, and vice versa [20]. An analogous result has been obtained for heterogeneous populations of noiseless oscillators [17]. Qualitative changes have also been found in the macroscopic phase response of rhythms in noisy homogeneous networks when the noise level was changed [21,22,28]. Using a Fokker-Planck approach for globally coupled excitable neurons, a type-I mPRC was obtained for weak noise, where the rhythm emerges through a SNIC bifurcation, while a type-II mPRC arose for strong noise that led to a Hopf bifurcation [21]. A similar approach was used to obtain the mPRC via the adjoint method for an extension of the theta-model that implements conductance-based synaptic interactions. Again, although individual theta-neurons have a type-I PRC, a type-II mPRC was obtained for the rhythm, which arose in a Hopf bifurcation [22]. This was also the case in an extension to networks of excitable and inhibitory neurons [28]. Thus, results reminiscent of those presented here have been obtained previously. However, the mechanism underlying them was not addressed in detail and remained poorly understood. We expect that our analysis will provide insight into those systems. The key element of the mechanism discussed here is that due to the dispersion of the spike times, which either results from neuronal heterogeneity or noise, the external perturbation enables the within-network inhibition to suppress the spiking of a larger number of neurons than without it. In our system this was facilitated by the delay with which spikes triggered the within-network inhibition, which allowed some neurons to escape its impact in the absence of the external perturbation, but not in its presence. Our analysis showed, however, that the explicit delay is not necessary; the effective delay resulting from a double-exponential synaptic interaction is sufficient. In fact, when reducing that effective delay the paradoxical phase response did not disappear until the delay was so short that the rhythm itself developed a strong subharmonic component and disintegrated.
In this paper we have focused on a specific, very simple neuronal model, the leaky integrate-fire model with conductance-based pulsatile coupling. In previous work on the enhanced synchronization among γ-rhythms via noise-induced spiking heterogeneity it was demonstrated that this result does not depend sensitively on the neuron type. Comparable results were obtained also with Morris-Lecar neurons for parameters in which the periodic spiking arises from a SNIC-bifurcation, resulting in a type-I PRC as is the case for integrate-fire neurons, but also for parameters for which the spiking is due to a Hopf bifurcation, resulting in a type-II PRC [29]. For networks comprised of heterogeneous neurons with type-II PRC the fmPRC of the collective oscillation is likely to be more complex, since the heterogeneity allows the same input to induce phase shifts with opposite signs for different neurons. However, we expect that the interplay between the within-network inhibition and the external perturbation can again substantially and qualitatively modify the fmPRC by changing the number of neurons participating in the rhythm. In [29] the results were also found to be robust with respect to significant changes in the network connectivity (random instead of all-to-all) as well as the reversal potential of the inhibitory synapses, as long as the rhythm itself persisted robustly (cf. [38]). In fact, the coupling did not even have to be synaptic; collective oscillations of relaxation-type chemical oscillators that were coupled diffusively were also shown to exhibit noise-induced synchronization. These results suggest that the paradoxical phase response found here arises in a much wider class of macroscopic collective oscillations. The strong paradoxical phase response that we demonstrated for heterogeneous networks allows their rhythm to synchronize with a periodic external input over a range of detuning that increases substantially with the neuronal heterogeneity. This is reminiscent of computational results for anterior cingulate cortex that investigated networks of excitatory neurons coupled via a common population of inhibitory neurons. There heterogeneity was also found to enhance the synchrony of rhythms, as measured in terms of coincident spikes within 10ms bins [39]. The heterogeneity-enhanced synchrony we have identified suggests that the coherence of γ-rhythms emerging in different interacting networks could also be enhanced by neuronal heterogeneity. It has been proposed that the coherence of different γ-rhythms, which has been observed to be modified by attention [8], plays an important role in the communication between the corresponding networks [11,32]. Computational studies have shown that the direction of information transfer between networks depends on the relative phase of their rhythms [12,13], which can be changed by switching between different base states [40,41]. Whether the enhanced synchrony resulting from neuronal heterogeneity enhances this information transfer is still an open question.
Disrupted γ-rhythms have been observed in multiple brain regions in neurological diseases, especially Alzheimer's disease. Optogenetic and sensory periodic stimulation at γ-frequencies has been found to entrain the γ-rhythm in hippocampus and visual cortex, respectively, and has resulted in a significant reduction in total amyloid level [42]. Similar neuro-protective effects of entrainment by external γstimulation have also been found for other sensory modalities [14,43]. This suggests that γ-stimulation by sensory input might be a feasible therapeutic approach. Our results suggest a potential role of neuronal heterogeneity in this context.
From a functional perspective, it has been shown that the noise-induced synchronization mentioned above can facilitate certain learning processes [44]. Specifically, a read-out neuron was considered that received input from neurons in two networks via synapses that exhibited spike-timing dependent plasticity. The two networks were interacting with each other and each of them exhibited a γ-rhythm, albeit at different frequencies. For low noise the two rhythms were not synchronized and the read-out neuron received inputs from the two networks at uncorrelated times. These inputs drove the plasticity inconsistently, leading only to a very slow overall evolution of the synaptic weights, if any. However, for stronger noise the two networks were synchronized, providing a more consistent spike timing that lead to substantial changes in the synaptic weights. As a result, the read-out neuron was eventually only driven by the network that had the larger natural frequency in the absence of the coupling between the networks. It is expected that synchrony by neuronal heterogeneity will have a similar impact.

Methods
Neuron model. Both E-cells and I-cells were modeled as leaky integrate-and-fire neurons, each characterized by a membrane potential V i (t) satisfying where V rest is the resting potential and τ E,I the membrane time constants of the E-and I-cells, respectively. I  (t) is a time-dependent external input that represents perturbations applied to determine the fmPRC or, in the study of synchronization, the periodic input generated by the clock network. I (bias) i denotes a tonic, neuron-specific excitatory bias current that implements the heterogeneity of the neuron properties. The corresponding conductances are denoted g syn , , g ext , and g bias . Upon the i th neuron reaching the spiking threshold V peak , the voltage V i was reset to the fixed value V reset . Parameters for the neuron were kept fixed throughout all simulations (see Table 1). The local field potential (LFP) of the network was approximated as the mean voltage across all neurons j = 1, ...N in the respective population. principal cells with all-toall interneuron-principal and principal-interneuron connections (i.e., without principal-principal and interneuron-interneuron connections). In PING, only principal cells received external input I ext (t).
To gain insight into the interaction between two ING rhythms, we considered the simplified situation in which all neurons in the network received strictly periodic input I (ext) , which was generated by another ING network ('clock'). Similarly, for PING rhythms, the E-cells of the PING network received strictly periodic excitatory input I (ext) from another PING network through all-to-all connection between their E populations.
Synaptic currents. We used delayed double-exponential conductance-based currents to model the excitatory and the inhibitory synaptic inputs from neuron j to neuron i, The periodic input ('clock') that was used to test the synchronizability of the ING-rhythm was gen-

ING network Parameter
Value τ I , membrane time constant 20 ms u rest , resting potential -55 mV V peak , spiking threshold -50 mV V reset , reset voltage -60 mV τ d , synaptic delay 3 ms N