{'fraction_non_alphanumeric': 0.03219660659509203, 'fraction_numerical': 0.008950824386503067, 'mean_word_length': 4.802071746384872, 'pattern_counts': {'":': 1, '<': 9, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Synchronization of electrically coupled resonate-and-fire neurons

Electrical coupling between neurons is broadly present across brain areas and is typically assumed to synchronize network activity. However, intrinsic properties of the coupled cells can complicate this simple picture. Many cell types with strong electrical coupling have been shown to exhibit resonant properties, and the subthreshold fluctuations arising from resonance are transmitted through electrical synapses in addition to action potentials. Using the theory of weakly coupled oscillators, we explore the effect of both subthreshold and spike-mediated coupling on synchrony in small networks of electrically coupled resonate-and-fire neurons, a hybrid neuron model with linear subthreshold dynamics and discrete post-spike reset. We calculate the phase response curve using an extension of the adjoint method that accounts for the discontinuity in the dynamics. We find that both spikes and resonant subthreshold fluctuations can jointly promote synchronization. The subthreshold contribution is strongest when the voltage exhibits a significant post-spike elevation in voltage, or plateau. Additionally, we show that the geometry of trajectories approaching the spiking threshold causes a"reset-induced shear"effect that can oppose synchrony in the presence of network asymmetry, despite having no effect on the phase-locking of symmetrically coupled pairs.


Introduction
Synchronization of activity between neurons has been hypothesized to contribute to a variety of brain functions [74], including motor control [18], memory [30], and coordination between brain regions [11]. This synchrony can be supported by either electrical or chemical synapses, or some combination of the two. Because electrical synapses (gap junctions) diffusively couple the voltages of connected cells, their effect is typically thought to be synchronizing, an idea with support from both theoretical and experimental studies [5,14,55]. However, their effect is potentially more complex, in part because the coupling combines effects at the dramatically different timescales of 1 spiking activity and subthreshold fluctuations of membrane potential [14,50,13,67]. Observations have shown that many cell types with strong electrical coupling exhibit resonant properties, which can create distinctive voltage fluctuations between spikes [37]. The same intrinsic properties that determine a frequency-selective spiking response to input also cause subthreshold dynamics such as transient oscillations, hyperpolarization followed by rebound, or a depolarized plateau following the spike; any of these effects can potentially contribute to synchronization. Our goal is to explore the interaction between electrical coupling and resonant intrinsic dynamics of spiking neurons, to understand both the dynamical mechanisms involved and their relevance to the function of neural systems.
We study the effect of electrical coupling on the synchronization of resonant spiking neurons, by applying the theory of weakly coupled oscillators to reduce the complexity of the synchronization problem and gain analytical insight [4]. This technique relies on a perturbative approximation to derive a reduced phase model for limit cycle oscillators [72]. Synchronization of the phase model is determined by the interaction function, which captures the effect of coupling as a function of the phase of each oscillator along its periodic limit cycle. Determining how the interaction function depends on a property of the oscillator, such as subthreshold resonance, spike size, or post-spike behavior, shows how that property contributes to synchronization. Common challenges in phase reduction analysis are that it may not be possible to independently vary the dynamical properties of interest, or to analytically compute the interaction function.
For the phase reduction to be analytically tractable, we use a minimal hybrid neuron model for the dynamics of resonant spiking. Hybrid neuron models, such as the integrate-and-fire model and its generalizations, idealize spiking as a threshold crossing with discrete post-spike reset, combined with continuous subthreshold dynamics between spikes. We focus on the resonate-and-fire model [38], which has linear damped oscillations as its subthreshold dynamics. The simplification of spiking allows the model to remain analytically tractable, while the discrete reset map helps to create complex resonant or integrator-like spiking dynamics. With the separation of discrete and continuous dynamics, we can independently vary the subthreshold and spiking properties of the model and determine their effects on the synchrony of small model networks.
On the other hand, the discrete reset map complicates the application of weakly coupled oscillator theory for the analysis of synchrony. The reset creates a discontinuity of the limit cycle, breaking a standard assumption in one step of the phase model reduction: calculation of the phase response curve. The phase response curve (PRC) measures the phase shift resulting from a perturbation to the oscillator at any point along the limit cycle. The discontinuous spiking dynamics lead to discontinuity of the PRC. In certain cases, discontinuous PRCs can be calculated directly [60,17,27], but the discontinuity in general necessitates an extension of standard methods for calculating the PRC. Shirasaka et al. [73] recently proved that for general hybrid models the PRC can be calculated by a variation to the standard "adjoint method," as previously suggested by Ladenbauer et al. [46]. We present an alternative, intuitive derivation of this result, elucidating the connection to the geometry of the threshold and reset in the context of hybrid neuron models, and apply this understanding to our resonate-and-fire analysis.
The paper is organized as follows. In section §2, we describe the general properties and history of hybrid models in neuroscience, and define our generalization of the resonate-and-fire model. In section §3, we review the theory of weakly coupled oscillators and present our approach to calculating the PRC for hybrid models. The remainder of the paper contains our analysis of synchronization in the resonate-and-fire model. In section §4, we apply our adjoint method approach to obtain an analytical expression for the PRC and interaction function of the electrically coupled resonate-and-fire neuron. To explore the dependence of the interaction function on model parameters, we focus separately on the even-and odd-symmetric components of the interaction function, which lead to distinct effects on synchrony. In section §5, we show that the spike itself always promotes robust synchrony through its contributions to the odd component, while the subthreshold fluctuations can additionally strongly promote synchrony in a "plateau potential" regime with strong resonant dynamics. We also show, in section §6, that the threshold and reset can phase shift the interaction function through "reset-induced shear" arising from the geometry of trajectories crossing the threshold, leading to an even component that can have complex effects on synchronization. In example three-cell networks, the presence of this even component leads to often-ignored effects on network synchronization when heterogeneity of frequencies or coupling breaks the symmetry of the interactions.

.1 Hybrid models
We first introduce some fundamental examples of hybrid models in neuroscience, and provide notation that will be used in the following sections. Hybrid models have a central role in the history of mathematical neuroscience. Well before the detailed processes generating action potentials were understood, Lapicque [49] postulated that inputs to a neuron accumulate in a continuous process of integration, eventually triggering a spike. This idea led to the leaky integrate-and-fire model (2.1) and a number of variations that are still widely used [25,38,39,7,1]. The separation of the dramatically different timescales of spiking and subthreshold dynamics into distinct mechanisms gives these models both computational efficiency and analytic tractability. Generalized integrateand-fire models are remarkably effective at reproducing diverse spiking behaviors [39,40,57] and can even be fit directly to spike trains [41].
Hybrid models have surprisingly rich and complex dynamics, inspiring active study from both neuroscience and dynamical systems perspectives [77,71,12]. For the analysis of synchrony, a number of studies have considered networks of coupled integrate-and-fire neurons and singlevariable variants [26,58,45,15,32,61,50,67,64,51,54]. However, despite the importance of synchronization in neural dynamics, only a few studies have addressed the synchronization of more complex hybrid models with more than a single variable, which is necessary to exhibit resonance [19,59,60,16,27,46,47].
The leaky integrate-and-fire model consists of a single voltage-like variable with linear subthreshold dynamics between spikes [49,10,1]. External current input I (t) is integrated through changes in the neuron's membrane potential (voltage) subject to a "leak," or linear decay over time. When voltage crosses a threshold v T from below, a spike occurs and the voltage is reset to v R .
For sufficiently large constant current input, the equilibrium voltage is pushed over the threshold, I (t) = v eq > v T , and the model exhibits a regular spiking state, a limit cycle with periodic firing. Note that the existence of a limit cycle in one dimension is enabled by the discontinuous reset.
A number of variations on this model have been proposed. Some models add nonlinear subthreshold dynamics in the form of a quadratic [25,9], exponential [31], or other function to more accurately model the approach to threshold. Adding a second "adaptation" variable [78,69] allows the model to reflect various slow recovery processes in the neuron (typically the gating dynamics of ion channels) that can create resonance and adaptation. With the presence of a second variable, not only can the subthreshold dynamics be more complex, but the threshold voltage v T becomes a threshold manifold, and the reset voltage v R can generalize to a reset map from this threshold manifold to a corresponding reset manifold. While models with a "hard reset," like (2.1), map all trajectories to a single point, assuming the adaptation process saturates during the spike and eliminates any history in the dynamics, models with a "soft reset," like the adaptive exponential integrate-and-fire model [7], instead map to a line of constant voltage. Generally, the soft reset increments the adaptation variable, modeling a rapid change in the state of ion channels which, unlike the hard reset, allows changes in the adaptation variable to persist over multiple cycles.
Although a hybrid neuron model can in general consist of any number of dimensions with complex threshold manifolds and reset maps [57], a two-dimensional model is sufficient to illustrate the essential properties. For our explanations, we restrict models to two dimensions x = (v, w) with threshold manifold T , reset manifold R, and corresponding reset map R : T → R.

Resonate-and-fire model
The basic elements of the resonate-and-fire model were developed independently by some of the earliest mathematical neuroscientists, Rashevsky [68] and Hill [36]. Izhikevich [38] introduced a modern version of the model and coined the "resonate-and-fire" name. The model's response dynamics have been shown to reproduce essential features of the dynamics of resonant cells, including a frequency-selective firing response to periodic input and a rebound of voltage following hyperpolarizing input [38]. These properties have been widely observed in many classes of neurons, and can be caused by a variety of voltage-dependent ion channels [37]. Note that since rebound from hyperpolarization requires crossing of subthreshold voltage trajectories, our model must have more then one dimension. The same holds for non-monotonic voltage trajectories with a post-spike plateau potential. The linear two-dimensional resonate-and-fire model is thus a minimal model for studying these phenomena. The resonate-and-fire model has previously been studied with periodic forcing [43] and in the context of pulse-coupled pairs of neurons [59] and globally pulse-coupled populations [60]. The subthreshold dynamics of the resonate-and-fire model are linear in two dimensions, consisting of linear decay of both the voltage variable and the adaptation variable towards an equilibrium voltage v eq , combined with coupling between the variables to create linear damped oscillations. 3) The resonate-and-fire hard reset rule specifies that if a threshold v T in the voltage coordinate is .1: Spiking limit cycle of the resonate-and-fire model with hard reset, shown in the (v, w) phase plane. Two nearby trajectories (blue) converge to the limit cycle (red) immediately following the instantaneous reset (dotted lines) at the reset point x R . Isochrons (gray) connect points of the same asymptotic phase (including along the spiking threshold at v T = 0). Parameters: crossed from below, the state is instantaneously reset to a single point We set the threshold voltage to an arbitrary value v T = 0, rescale time such that ω = 1, and unless otherwise specified choose a fixed value for the decay parameter, λ = 0.1. Our results are qualitatively similar as long as λ is sufficiently small. In the opposite extreme, taking the limit λ → ∞ with the product λω fixed recovers the leaky integrate-and-fire model. 1 With these parameters fixed, the reset parameters and equilibrium voltage determine the existence and properties of spiking in the model. Varying the equilibrium voltage v eq relative to the threshold reflects a combined effect of altering the biophysical resting potential and the tonic (constant) current input to the cell. If a trajectory starting from the reset point x R crosses the spiking threshold and is again reset, it forms a spiking limit cycle corresponding to regular spiking, shown in figure 2.1. We explore more detailed existence conditions in section 4.1, but a limit cycle will always exist for sufficiently large constant current input, moving the equilibrium voltage v eq over the spiking threshold just as in the integrate-and-fire model. The reset parameters, especially the reset voltage v R , also have a strong effect on the shape of the voltage trajectory for this limit cycle, as shown in figure 2.2. The voltage trace for an oscillator reset near the peak of its oscillation (large positive reset v R , figure 2.2 right) resembles a plateau potential, a sustained post-spike elevation in voltage. The trace for an oscillator reset near its trough (large negative reset v R , figure 2.2 left) resembles a fast after-hyperpolarization (AHP), a post-spike dip in voltage. We will refer to these strong negative and positive reset regimes simply as "plateau" and "AHP" below. An intermediate case with v R = 0 is shown in figure 2.2 center. As several studies have pointed out, these distinct after-potential regimes can significantly affect synchronization mediated by electrical synapses [13,67,55], which we will investigate further using the resonate-and-fire model. The regular spiking solution for the resonate-and-fire model is defined to start from the reset point at time t = 0 and is valid for times up to the period T , when the trajectory crosses threshold.

Extended resonate-and-fire model: soft reset, spikes, and coupling
From the perspective of a general hybrid model, the spike-and-reset dynamics restrict the definition of both the threshold and reset: with the hard reset map R (w) = x R from (2.4), the threshold manifold is a line of constant voltage T = {(v, w) : v = v T }, and the reset manifold is a single point R = {x R = (v R , w R )}. In order to expand the possible dynamics of resonate-and-fire models while retaining simplicity, we also consider the "soft reset" variation commonly used in other hybrid models, which increments the adaptation variable by a constant value. The reset manifold for the soft reset is thus a line of constant voltage More generally, the threshold and reset manifolds could have arbitrary orientations and locationswe briefly discuss considerations for the general case in section 4.2.
To build a network model of resonate-and-fire oscillators, we define the electrical coupling as direct exchange of current, coupling the voltage but not the adaptation variables [51]. Specifically, we assume the coupling to be linear in the voltage difference and let k i j be the coupling strength between neuron j and neuron i. For the cells in the network, we introduce heterogeneity of intrinsic frequencies ω i about a mean frequencyω = 1, while keeping other parameters fixed for the population. The resulting network model is as follows, where I c x i , The spiking dynamics of the resonate-and-fire model (2.3-2.6) describe the subthreshold dynamics, threshold crossing, and post-spike reset but not the spike itself. In order to model the effects of electrical coupling, which transfers current based on voltage fluctuations both during and between spikes, we must supplement the model by a description of the transient voltage spike. We choose to define the spike as a δ-function at the instant of crossing threshold [51], The spike magnitude parameter M determines the integral of the spike over time, or current exchanged. When M is large relative to the subthreshold fluctuations, the effect of electrical coupling is essentially pulse coupling, as studied by Miura and Okada [59]. (However, their results are not directly comparable, as they consider pulse coupling to the adaptation (w) variable rather than the voltage (v) variable.) We focus on a regime where the spike effect is comparable to or smaller than the subthreshold fluctuations, and use the theory of weakly coupled oscillators to analyze the spike and subthreshold contributions to synchrony independently.

Phase mapping for weakly coupled oscillators
The theory of weakly coupled oscillators is a method of model reduction with the goal of creating simpler models for the dynamics of interacting limit-cycle oscillators. The state of each oscillator is mapped to a single variable: the phase (or timing) of oscillations. This process is referred to as "phase reduction" and the result as a "phase model." The oscillators' intrinsic dynamics and the coupling between oscillators together determine the interaction function, which captures the modulation of instantaneous frequency caused by coupling. This scalar function defines the coupling between the phase oscillators, and thus determines the dynamics of the coupled system (together with any heterogeneity of intrinsic frequencies). Dynamical properties such as the stability of phase-locked states can be directly assessed from features of the interaction function, as we will show in section §6. Although strictly derived in the limit of weak coupling and heterogeneity, predictions from the phase model often remain valid even at moderate levels of coupling. Here we first provide a brief derivation of the phase model for a continuous coupled oscillator system. We then explain the specific challenges presented by hybrid systems and our approach to overcoming these challenges. The dynamics of a general system of coupled oscillators are described by The state space of each individual oscillator is x i ∈ R n , where n is the dimension of the oscillator model. The intrinsic dynamics of each oscillator are defined by f + g i , where f gives an average of the intrinsic dynamics across all oscillators, and g i captures the heterogeneity of the oscillators. The average dynamics must have an asymptotically stable T -periodic limit cyclex (t) defined by The pairwise coupling is defined by the coupling function I c , which is weighted by connection strengths k i j and summed over all pairs to give the total coupling. The phase reduction requires the assumption of weak coupling and weak heterogeneity, meaning that the heterogeneity g i and the total coupling term must both be small compared to the average intrinsic dynamics f . To reduce the model for this coupled system, we define a phase mapping Θ : R n → R, from the state x to a scalar phase variable θ that uniquely identifies points on the limit cycle. The dynamics of the coupled system are translated through this mapping to define the phase model, The dynamics of this reduced system depend only on the interaction function H and the heterogeneity of frequencies Ω i . Below, we show how these are derived from the coupling I c and heterogeneity g i of the full model. We first define the phase map for points on the limit cycle, giving a periodic "time" coordinate, Phase is unique modulo T , with θ = 0 determined by our choice of initial conditionx (0) for the reference limit cycle. (Note that phase is sometimes rescaled to a period of 1 or 2π, but we choose to keep the natural units of time.) The phase map can then be extended beyond the limit cycle to give the "asymptotic phase" of all points in the basin of attraction. Trajectories that eventually converge are assigned the same phase, i.e., Θ (x (t)) = Θ (x (θ)) = θ iff lim Although typically not calculated explicitly, the full phase mapping can be visualized by its level surfaces, referred to as isochrons, linking points that approach the limiting trajectory with the same timing (figure 2.1). 8

Phase response curve and the adjoint method
The full phase mapping is difficult to find analytically in all but the simplest contexts. Fortunately, the weak coupling assumption allows the phase reduction to proceed with a linear approximation of the mapping about the limit cyclex (t), which is much easier to compute. For a trajectory close to the limit cycle, phase can be approximated as linearly dependent on the deviation away from the limit cycle, ∆x (t).
The infinitesimal phase response curve Z (iPRC or PRC, also called the "phase sensitivity function") is defined as the proportional shift in phase caused by infinitesimal perturbations to the limit cycle, Note that the PRC is naturally defined as a vector-valued function, giving the effect of perturbations to each state variable. In the context of neural synchrony, however, the voltage component is usually most important, because perturbations tend to be currents and thus directly affect only the current balance equation for the dynamics of voltage.
A direct method to approximate the PRC either experimentally or computationally is simply to measure phase shifts caused by many small but finite perturbations. We instead follow the "adjoint method," which derives and solves a differential equation for the PRC tied to the limitcycle dynamics. Below we provide a brief exposition that captures the essence of this method and its proof.
From the definition of asymptotic phase (3.4), the phase difference between the limit cyclex (t) and a nearby trajectory x (t) must be constant in time. That is, the separation ∆x (t) = x (t) −x (t) must satisfy To first order, the deviation from the limit cycle ∆x (t) evolves according to the Jacobian matrix of derivatives of the system dynamics, D f , evaluated on the limit cycle given by (3.2).
Substituting this expression into (3.7), we obtain Because this holds for any ∆x, Z must satisfy the following T -periodic linear differential equation known as the adjoint equation, defined by the adjoint of the linearized limit-cycle dynamics, The PRC is the unique periodic solution of (3.8), given a normalization constraint that follows from our definition of phase on the limit cycle (3.3).
Note that if this constraint is satisfied at any single time, (3.8) ensures it will remain satisfied for all time.

The phase model and interaction function
Using the PRC result derived via the adjoint method, we can complete the derivation of the phase reduction of the coupled oscillator system (3.1). The effect of any weak time-dependent perturbation on the phase of an oscillator, including the effect of coupling, is governed by the PRC. Specifically, if the dynamics of the limit cycle are continuously (weakly) perturbed according to Since the perturbation is weak, its effects occur on a slow timescale, O (ǫ), which can be separated from the faster dynamics of unperturbed phase, dθ dt = 1. If the perturbation is also periodic with the intrinsic period T , this separation allows us to eliminate the explicit time-dependence p (t) by averaging the slow effect of coupling over a full period of the fast phase dynamics.
Because the perturbations that define the coupled population model (3.1) are close to periodic, their effects can be approximated by the method of averaging. 2 The heterogeneity and coupling perturbations to an oscillator are functions of its own trajectory x i and those of the other oscillators x j . Each trajectory is approximated by the T -periodic average limit cycle, x i ≈x (θ i ) (where θ i = Θ (x i )); therefore we can approximate these perturbations as periodic functions of phase, The final result of the phase reduction is the phase model, i.e.
where The frequency term Ω i arises from the intrinsic heterogeneity in cellular properties, and the interaction function H from the coupling. We note that the form of H as a function of phase differences, φ j−i = θ j − θ i , arises from the application of averaging.

PRC for hybrid models
For hybrid models, the PRC as well as the trajectory may be discontinuous at the threshold crossing. Without the periodicity constraint acting as a boundary condition, the adjoint equation (3.8) with normalization (3.9) no longer has a unique solution, and a naive application of the adjoint method fails to find the PRC. Intuitively, the resolution is to find the appropriate "reset" map or boundary condition for the PRC. Shirasaka et al. [73] present a boundary condition linked to the saltation matrix, a correction to the linearized dynamics of a hybrid system to account for discontinuity across a boundary [6]. They prove that the solution to the adjoint problem with their boundary condition gives the PRC for the asymptotic phase of hybrid systems. Special cases of this result have also been presented by Ladenbauer et al. [46] for a specific hybrid neuron model and by Park et al. [65] for nonsmooth systems with discontinuous boundaries but no reset map. We present a brief heuristic derivation for an equivalent adjoint boundary condition that follows directly from the existence of an appropriate differentiable phase mapping, and show that the condition has an intuitive form tied to the geometry of the threshold. Consider a limit cycle trajectoryx (t), reaching the threshold manifold atx (T ) = x T = (v T , w T ), as shown in figure 3.1. The phase on this trajectory cannot change across the instantaneous reset. That is, since Θ (x (t)) must be continuous in t, where the right limit of the reset discontinuity is given by the reset point R (x T ). We introduce a unit tangent vector u along the threshold manifold T at x T (figure 3.1). Assuming the PRC is well-defined before and after reset, we can apply a directional derivative along this tangent vector (D u ) to both sides of (3.13).
This equation gives a boundary condition for the PRC, replacing the standard assumption of periodicity. Together with the normalization condition (3.9), it determines the unique solution to the adjoint problem for the discontinuous limit cycle of the hybrid model. With the reset map only defined on the threshold manifold, the derivative is only defined in the tangent direction, and there is no corresponding constraint on the perpendicular component of the PRC. In N dimensions the threshold manifold is (N − 1)-dimensional, so rather than a single vector u, we enforce (3.14) for each of N − 1 vectors u i spanning the tangent space. 3 This boundary condition expresses an intuitive fact about perturbations to the limit cycle: the phase difference between the limit cycle and a trajectory perturbed along the threshold (given by Z u (T − )) must be the same as the difference after both trajectories are reset. The difference after reset is expressed by the PRC at the reset point, Z (0 + ), with the perturbation transformed by the reset map to a distinct perturbation along the reset manifold, approximated by D u R (figure 3.1).
4 Phase reduction of the resonate-and-fire model 4

.1 Existence and stability of spiking limit cycles
Before applying the theory of phase reduction to any model, we must ensure the system exhibits a stable limit cycle. We begin our analysis of the resonate-and-fire model by finding the existence and stability conditions for spiking limit cycles. These conditions define boundaries of the stable limit cycle regime in parameter space, bifurcations of the model.

Hard reset
We first discuss the existence conditions of the spiking limit cycle with hard reset, (2.4). The spiking limit cycle exists whenever a trajectory starting from the reset point crosses the threshold. The limit cycle is lost in a "grazing bifurcation" when the trajectory becomes tangent to and then fails to cross the firing threshold. Beyond this bifurcation, trajectories show decaying subthreshold oscillations, approaching rest at the equilibrium voltage. The hard reset map ensures that the spiking limit cycle is always stable, as the reset erases all effects of small perturbations to the cycle by projecting to the single reset point.
We can visualize the spiking regime boundaries in a two-dimensional parameter space that captures the most important dimensions of variability of the model dynamics. We fix the frequency ω = 1 without loss of generality, and choose a small decay parameter λ = 0.1 to give slowly decaying subthreshold oscillations. The dynamics then depend on the reset parameters v R and w R -2 Rest Spiking In both rest and depolarization block regimes, the system has a stable fixed point (quiescent state); this state is below threshold at rest and above threshold for depolarization block.
as well as the equilibrium v eq , but because the model is invariant to uniform rescaling of the (v, w) phase space, we can rescale to |v R | = 1 without loss of generality. Therefore, in figure 4.1 we explore two distinct two-dimensional (v eq , w R ) parameter spaces for positive and negative reset voltage, fixed at v R = ±1 (see figure 2.2). In the both negative and positive reset regimes, a single grazing bifurcation occurs for low equilibrium voltage v eq (negative tonic input current), below which the model is quiescent (at rest). In the positive reset regime, a second grazing bifurcation occurs for high v eq , when the voltage fails to dip below threshold, corresponding to depolarization block. The bifurcations shown in figure 4.1 curve away from the origin because any increase in the magnitude of w R increases the radius of the orbit, moving the system away from the bifurcation. We can express the tangency condition for the grazing bifurcation in terms of the orientation of the trajectory when crossing threshold, θ H = tan −1 v(T ) w(T ) . At the bifurcation, this orientation matches that of the threshold,

Soft reset
Although the subthreshold dynamics always lead to decay of perturbations, the soft reset map can in some cases amplify perturbations, making the limit cycle unstable. In addition to the grazing bifurcation boundaries, the stable spiking regime can also be lost in a saddle-node bifurcation of limit cycles, where the stable and unstable cycles collide and annihilate. In addition, the dynamics with soft reset allow for limit cycles with multiple spikes, and we will show that the spiking limit cycle can also lose stability in a period-doubling bifurcation.
To determine existence and stability conditions for the soft reset limit cycles, we reduce the dynamics to a Poincaré return map on the reset manifold (also referred to as an adaptation map [77]). This map takes a value of the adaptation variable at the k th reset, w k , to the value at the following reset, w k+1 = P (w k ). A limit cycle corresponds to a fixed point of the map,w = P (w), and the cycle is asymptotically stable if the fixed point satisfies dP dw (w) < 1. The derivative of the return map characterizes the degree of attraction to the limit cycle and is therefore also tied to the validity of the asymptotic approximation of the phase reduction. We begin by deriving the return map for a general reset map from the threshold (v T , w) to R (w) = (v R , R w (w)). We then evaluate the fixed points and their stability for the soft reset map.
We define the flow where the trajectory x (t) satisfies the subthreshold dynamics from equation (2.3) starting from an initial condition x (0) = (v R , w k ) T on the reset manifold. The flow evaluates to We then define the spike time map τ (w k ), giving the time it takes such a trajectory to reach the threshold, A trajectory starting at w k crosses threshold at the point F (w k , τ (w k )), and is reset to R (F w (w k , τ (w k ))).
The w-component of this reset point, R w , gives the desired return map, A fixed pointw = w 0 of the return map corresponds to a limit cycle with period T = τ (w 0 ). The stability of a limit cycle with soft reset is assessed by evaluating the derivative of the return map for the soft reset rule R w (w) = w + ∆w (2.6).
where we recall the trajectory's orientation at threshold θ H = θ 0 + T + π 2 + tan −1 λ. For the soft reset rule then, the limit cycle can be stable or unstable depending on the decay λ and the geometry of the threshold and reset manifolds. Loss of stability occurs when which implies that bifurcations occur in the full resonate and fire model when The negative slope instability corresponds to a subcritical (proof not shown) period doubling bifurcation, where an unstable period-two limit cycle collides with a stable period-one limit cycle. The positive slope instability corresponds to a saddle-node bifurcation of limit cycles, with a stable and magenta the positive slope (saddle-node) bifurcation. Solid red/blue indicates unstable cycles past the bifurcation, slopes dP dw > 1. Black lines show the grazing bifurcation of the unstable limit cycle. Parameters λ = 0.1 and v R = ±1 fixed while varying v eq and ∆w, plotted using w 0 to facilitate comparison with hard reset parameter space.
cycle coalescing either with a finite unstable cycle, or at infinity (return map slope approaching unity as w 0 → −∞ for finite parameter values).
In figure 4.2, we plot the stability and grazing bifurcations for soft reset together in v eq , w 0 parameter space. Since w 0 corresponds to the parameter w R for hard reset, the grazing bifurcations (solid lines in figure 4.2) are identical to the hard reset grazing bifurcations in v eq , w R coordinates from figure 4.1. The loss of stability bifurcations and the grazing bifurcations form the two boundaries of a narrow unstable limit cycle regime. These bifurcations are both related to threshold crossing and necessarily lie close together. Perturbations are amplified by the reset map, causing instability, because of a mismatch between the angles of incidence with the threshold and reset manifolds, which increases as the trajectory approaches tangency to the threshold (the grazing bifurcation condition). We note that, near the saddle-node bifurcations, the return map can have two fixed points, representing stable and unstable cycles. The representation in figure 4.2 unfolds the bifurcation so that a single point in the parameter space v eq , ∆w corresponds to two points in v eq , w 0 space, stable and unstable cycles on opposite sides of the dashed bifurcation line.

Resonate-and-fire PRC
The first step in proceeding with the phase reduction of the resonate-and-fire model is to evaluate its phase response curve (PRC), expressing the effect of perturbations to the limit cycle as a proportional phase shift (section 3.2). In general, the PRC can be evaluated either by direct calculation or by the adjoint method. Although directly calculating the effect of perturbations typically must be carried out computationally, the hard reset of the original resonate-and-fire model allows the PRC to be calculated directly from an analytical expression for phase [60]. Since all points on the threshold are mapped onto the single point x R , the threshold itself is an isochron and can serve as a reference point to define the phase map for all trajectories. However, we take a different approach deriving the PRC from the adjoint equation for a general reset rule, following the theory described in section 3.4. This approach allows us to see the hard reset as a special case and to capture the geometric intuition in the relationship between the dynamics, the adjoint equation, and the reset rule.
We proceed by first evaluating the adjoint equation (3.8) for the subthreshold dynamics. In general, the adjoint equation evaluates the dynamics linearized about the limit cycle. Since the resonate-and-fire dynamics (2.3) are linear, the adjoint equation is simply defined by the negative transpose of the time-independent linear operator, The PRC solution exhibits exponentially growing oscillations, The PRC is defined as written for times 0 < t < T , and extends periodically to all t modulo T , with a discontinuity at t = 0 (due to the discontinuous reset map skipping over dynamics during the spike). The amplitude A = The general form for the boundary condition from (3.14) simplifies given our assumption that the threshold and reset manifolds are in the w-direction (constant v), with reset map R w (w) For the hard reset (2.4), mapping to a constant reset point, the derivative of the reset map is zero. Thus (4.8) reduces to the terminal condition Z w (T − ) = 0, corresponding to a phase α = 0. This result is equivalent to the geometric constraint that the PRC must be perpendicular to the threshold, or oriented in the v-direction at time T . That is, perturbations along the threshold have no effect after the reset, and the threshold is an isochron, as shown in figure 2.1. For the soft reset rule (2.6), an increment of w, the boundary condition (4.8) mandates continuity of the w-component of the PRC, Z w (T − ) = Z w (0 + ) . Intuitively, this tells us that a perturbation to w immediately before the spike has the same effect as a perturbation after the spike; i.e., perturbations tangent to the threshold are unchanged by the soft reset map. This continuity boundary condition leads to a phase shift of Z v , α = arctan sin T cos T − e λT .

16
For soft reset, ∆w was set to match the hard reset limit cycle, so that w 0 = w R .
Note that this implies that the threshold manifold is not an isochron. This phase shift is typically small for positive v R but can grow more significant in parts of the negative reset regime. (See example in figure 4.3, and full calculation in section A.2.) We note that in the more general case of a nonlinear reset map, the result depends on dR w dw , for which the hard and soft reset are the special cases dR w dw = 0 and 1 respectively. A nonlinear reset map with derivative close to either extreme would lead to small corrections to the corresponding α value. Similarly, small variations in the orientation of the threshold or reset manifold lead to minor adjustments to (4.8) and to the resulting phase shift α.
We also note that the amplitude A has a singularity when This singularity is associated with the high sensitivity to perturbations near bifurcations of the limit cycle. For the hard reset α = 0, θ H = ± π 2 is the grazing bifurcation condition (4.1). For the soft reset, the singularity occurs when tan θ H = − cot α = e λT −cos T sin T , which is condition (4.5) for the saddle-node bifurcation of the limit cycle (positive slope instability of the return map). Interestingly, the negative slope loss of stability (period-doubling bifurcation) is not reflected in the PRC. It seems that in the positive slope instability, phase shifts continue to accumulate progressively on each cycle, while in the negative slope case, positive and negative phase shifts alternate as they grow, leading the PRC to reflect an averaged finite phase shift that fails to capture the loss of stability.

Resonate-and-fire phase model
We now proceed to construct the interaction function and phase model for electrically-coupled resonate-and-fire neurons. In the weak coupling regime, this reduced model with a single phase variable for the state of each oscillator captures the full synchronization properties of the resonateand-fire network (2.7). The interaction function H and heterogeneity Ω of the phase model are calculated according to equation (3.10), from the electrical coupling and heterogeneity of frequencies in the full model along with the PRC from (4.7). The heterogeneity of frequencies in the phase model follows directly from the resonate-and-fire model's frequency heterogeneity. Evaluating the integral (3.11) and applying the normalization condition (3.9) shows that Ω i = ω i (the subthreshold angular frequency). We note that weak heterogeneity in any other parameter of the resonate-and-fire model (e.g., v R ) would create an additional additive contribution to the phase model heterogeneity Ω i proportional to that parameter heterogeneity.
The interaction function H defines the nonlinear coupling between cells in the phase model. It can be expressed as a convolution integral of the coupling current and the PRC, according to (3.12). In the case of electrical coupling, the coupling current I c = v (t + φ) − v (t) depends on the voltage component of the limit cycle. The resulting interaction function, depicted in figure 4.4, is Since the interaction function convolves the PRC and the voltage limit cycle, with a linear dependence on both, we can separate the subthreshold and spiking components, H = H sub + H spike , corresponding to the subthreshold and spiking components of the limit cycle from equation (2.8), The spike δ-function component of the limit cycle determines the spike interaction function H spike , representing the effect of this voltage transient through the electrical coupling. This is essentially a pulse-coupling interaction, as used in simple models of excitatory chemical synapses [58,45,59], and the effect is entirely determined by the PRC and the amplitude of the spike. With zero-width or δ-function pulses, the spike interaction function can be discontinuous at the origin, as shown in figure 4.4. We return to the effects of this interaction in section 5.3.
The subthreshold interaction function H sub captures the effect of subthreshold fluctuations of the limit cycle, v sub (t) = r 0 e −λt cos (t + θ 0 ). Analysis of this component will allow us to determine how the resonant properties of the model contribute to synchronization. Each parameter of the model affects synchronization both through its effect on the limit cycle and on the PRC, making the combined effect (encoded by H sub ) potentially complex. The only general constraints on the subthreshold interaction function are that it must be continuous and pass through the origin; H sub (0) = 0 because the gap junction coupling is diffusive, proportional to the difference of voltages.
To simplify the analysis, we split the subthreshold interaction function into three terms with distinct parameter dependence (calculation in section A.3). (4.11) where θ T = θ 0 + T is the angular coordinate of the trajectory at threshold. We note that the Sfunction closely resembles sine, while the two C-functions resemble a vertically shifted cosine, as shown in the example in figure 4.5. This approximate odd and even symmetry about the origin means they contribute in distinct ways to synchronization of simple networks. In the following sections, we consider these effects by studying networks of two and three cells.

Synchronization of two electrically coupled resonate-and-fire neurons
Our primary goal with the phase reduction of the resonate-and-fire model is to provide insight into the synchronization of networks of electrically coupled resonant neurons. Even after phase reduction, the analysis of large systems with realistic network architecture is hindered by the nonlinear phase coupling and number of degrees of freedom, and in general must be carried out numerically. Therefore, we focus on minimal networks of two or three cells and save the analysis of large-scale networks for future work. In this idealized context, we can explain how the cellular properties of the resonate-and-fire oscillators determine network synchrony through the distinct components of the interaction function (in terms of slopes, amplitudes and discontinuities). Although large networks cannot be completely understood by their two-and three-cell subnetworks, in many cases the intuition built on these minimal networks holds [2].

General considerations
We begin by examining the synchronization of a symmetrically coupled pair of oscillators, the simplest and most commonly analyzed network [26,44,50,55]. The assumption of symmetry implies that the equation governing the phase difference between the oscillators (5.2) isolates a component of the interaction function with odd symmetry, simplifying the analysis. Experiments have shown that electrical coupling is often close to symmetric, especially when between the same cell type and formed by the common symmetric channel type connexin-36 [5]. Nonetheless, asymmetry can and does arise, either from rectification (favoring current flow in one direction) in the gap junction channels that make up the electrical synapse, or from asymmetries in size or gap junction location. We will address the effects of asymmetry briefly in section §6. The phase model (3.10) for a pair of symmetrically coupled oscillators (k 12 = k 21 = K) is given Expressed in terms of the phase difference φ i− j = θ i − θ j and frequency difference ω 1−2 , (5.1) reduces toφ where H odd (φ) = 1 /2 (H (φ) − H (−φ)) = H (φ) − H even (φ). (Some analyses refer to the G-function, G (φ) = −2H odd (φ) [72]). Fixed points of (5.2) correspond to phase-locked states of the coupled pair, the existence and stability of which are determined by properties of H odd , as depicted in figure 5.1a and described below. For identical oscillators (ω 1−2 = 0), the odd symmetry of equation (5.2) implies a pair of fixed points at φ 1−2 = 0 and φ 1−2 = T /2, for synchronous and antiphase states respectively. If H odd has only a single local maximum, which is typical for the resonate-and-fire model, only one of these two fixed points is stable and no additional fixed points exist. The synchronous state is stable and the antiphase state unstable when the slope H ′ odd (0) is positive, and the reverse holds for negative slope. We note that in larger networks a similar dependence on the slope can be shown: a strong positive slope leads to global synchrony, while a negative slope leads to global incoherence [2]. (For simplicity, we shorten references to the slope evaluated at the origin to "slope".) For hybrid models, the slope of the full interaction function H ′ (0) may be undefined, with different right and left limits, but the odd symmetry forces H ′ odd (0) to always be well-defined. As the frequency heterogeneity of the pair increases, the pair of fixed points shift progressively in the relative phase of the oscillators. We refer to these states as near-synchronous and near-antiphase. For small frequency heterogeneity ω 1−2 > 0, the phase difference in the nearsynchronous state is approximately inversely proportional to the slope, φ 1−2 ≈ ω 1−2 2KH ′ odd (0) . Note that this phase difference is in units of time, while to compare phase-locked states across oscillators with different periods we should evaluate phase in radians. To account for this we rescale both the phase and the slope H ′ odd , For larger heterogeneity, the phase difference continues to increase until ω 1−2 is greater than the amplitude of H odd (red level in figure 5.1a). At this point the fixed points of (5.2) are lost in a saddle-node bifurcation, and the phase-locked state transitions to "phase walk-through," with the phase of the cells slipping past each other (red curve in figure 5.1b right).
We can thus quantify the robustness of synchrony to heterogeneity in two ways: the slopê H ′ odd (0) gives the strength of synchrony for small heterogeneity, while the amplitude of H odd gives the critical heterogeneity at which the near-synchronous state is lost. Below, we consider the robustness and stability of near-synchronous phase-locking in more detail for the subthreshold and spiking components of the resonate-and-fire interaction function.

Subthreshold contribution
For the analysis of synchrony in the resonate-and-fire model, we begin by evaluating the contribution to the odd component of the subthreshold interaction function H sub , and return to the , so its slope and amplitude, our two measures of robustness, are roughly proportional. We first identify the condition for stable synchrony (positive slope), and show that the subthreshold contribution is virtually always synchronizing, with negative slope only possible along the spiking regime boundary . We then explore the robustness of synchrony as measured by the magnitude of the positive slope, identifying a dependence on the shape of the voltage trajectory through the reset voltage v R . We use the decomposition (4.11) of H sub into component functions with approximate symmetry (see figure 4.5) to simplify our calculation of the slope. We consider the limit of small decay (λ ≪ 1) by expanding to first order in the decay parameter (see calculation in section A.4).
The O (1) contribution toĤ ′ odd (0) from S is the primary subthreshold factor determining of the strength and stability of synchrony. The sign of this slope is determined by A S , the coefficient of the S component. Specifically, the condition for stable synchrony is approximated to order λ by Note that this condition differs only to order λ from the PRC singularity condition (4.9), which determines the location of both the grazing bifurcation for hard reset and the saddle-node bifurcation for soft reset. Near these boundaries, a negative slope can result from either A S < 0 or from A S ≈ 0 and A C1 < 0 or A C2 < 0. We show the slopeĤ ′ odd (0) for the full parameter space in figure 5.2, with the negative slope region highlighted by zooming in near the spiking boundary (figure 5.2b). Additionally, because the slopeĤ ′ odd scales with the diverging PRC amplitude A, both negative and positive slopes near these boundaries can grow extremely large. We note, however, that this result should be interpreted with caution, as the assumption of weak coupling also breaks down approaching these boundaries.
The other significant trend in the slope, determining the robustness of synchrony, is the difference between the positive and negative reset regime. The slope is roughly unit magnitude in the positive reset regime, sharply contrasting with the negative reset regime where the slope is uniformly small (figure 5.2). By explicitly plotting the slope against the reset point v R in figure 5.3, we can see clearly the presence of two regimes with a distinct transition in between. As shown in figure 2.2 these two regimes have characteristic voltage waveforms. A large positive v R is characterized by a plateau potential in the voltage trajectory, while large negative v R is characterized by after-hyperpolarization (AHP). We conclude that for the resonate-and-fire model in the plateau potential regime, the subthreshold dynamics significantly contribute to synchrony. In the AHP regime, with little to no subthreshold contribution, the resonate-and-fire model can only be We can further investigate this trend through our slope approximation (5.3), focusing on the S component as the dominant slope contribution. As v R is increased, the coefficient A S varies little, whileŜ ′ odd ≈ T − sin T increases sharply with an increase in the period T . Note that withω = 1 fixed, T represents two distinct factors: the extent of the cycle in radians and its duration. However, varying the duration alone (via ω) has no effect on the phase-reduced dynamics (see section A.3), so the effects on H odd are entirely due to the extent of the limit cycle covering a larger fraction of a cycle of continuous oscillation.

Spike interaction function
The spiking component of the interaction function also has a significant odd component, and thus contributes to the synchronization of the coupled pair. We show here that this contribution always promotes synchrony, reinforcing previous results for excitatory pulse coupling of resonant neurons [34].
Because we model the spike as a δ-function (with magnitude M), the interaction function convolution equation (3.12) has a simple form for the spiking component, i.e., a time-reversed copy of the PRC.
The term δ (t) leads to an additional constant term Z v (0), the value of the PRC at the spiking discontinuity. Although the left and right limits of this point are defined by (4.7), Z v (0) is unconstrained. In experiments and biophysically detailed models, perturbations during a spike typically have little to no effect, due to the many ion channels open during an action potential lowering the input resistance of the cell [33]. Thus, we assume Z v (0) = 0, although we note that this choice determines only a constant offset and does not affect the discontinuity of H spike . Since the spike interaction function is discontinuous at zero, its contribution to the stability and robustness of the synchronous state depends primarily on the size and direction of the discontinuity. 4 The jump discontinuity can stabilize the fully synchronous state even with nonzero heterogeneity, if |ω 1−2 | < H odd (0 + ) = ∆H spike [20,51]. (For more in-depth analysis of the limit approaching discontinuity in the interaction function, see Shirasaka et al. [73].) We evaluate this discontinuity directly from the PRC: The condition ∆H spike > 0 is always satisfied for hard reset (α = 0), and in the soft reset case, ∆H spike = 0 evaluates to exactly the stability boundary (4.5), with ∆H spike > 0 for all stable cycles. The spike interaction function thus always promotes the stability of the synchronous state, potentially synchronizing heterogeneous oscillators even when the subthreshold contribution is not synchronizing. As we showed above, the subthreshold contribution to the slope H ′ odd (0) is near zero or negative for a significant portion of the parameter space, including most of the AHP regime, thus requiring this spike contribution in order to synchronize.

Synchronization of a three-cell network: effect of the even component
Due to the symmetry of the coupled pair, the odd component of the interaction function alone determines the evolution of the phase difference, and the even-symmetric component has no effect on synchronization. However, the even component has the potential to strongly affect synchronization in both larger resonate-and-fire model networks and in actual biological networks. Although many studies of synchronization in model neurons ignore this possibility by focusing on symmetrically coupled pairs, the effect of coupling with an even component has been explored in large, sparse neuronal networks [32] and in large regular networks of phase oscillators (Kuramoto model and generalizations) [2]. Here, we show that similar effects can be seen in minimal networks, specifically asymmetrically coupled pairs and three-cell networks, the study of which can help us understand the more complex properties of larger networks. We first provide intuition on the effects of the even component by revisiting the dynamics of a coupled pair, in terms of the phase difference φ 1−2 and the mean phaseθ.
In the symmetrically coupled pair (∆k = k 12 −k 21 =0), the even component has no effect on the phase difference, but it does shift the frequency dθ dt of the phase-locked state. In the case of asymmetric coupling strength k 12 k 21 , the term ∆kH even (φ) can be interpreted as a differential shift in the instantaneous frequencies of the two oscillators. This effective frequency shift can either promote or oppose synchrony depending on whether the sign of the product adds to or cancels with the intrinsic frequency heterogeneity ω 1−2 .
In a three-cell network, an even component term in the phase difference dynamics can also arise from coupling to a third oscillator with a different intrinsic frequency, even if the coupling is fully symmetric. The phase difference equation for a symmetric three-cell network is While the first two H terms partially cancel to isolate the odd component, the latter two terms depend essentially on the even component. We will first demonstrate the effects of these additional even component terms on phase-locking, along with the effective frequency shifts from coupling asymmetry, in the context of three-cell networks. We will then return to the resonate-and-fire model, investigating the size of the even component of the subthreshold interaction function, its potential effects on neuronal networks, and its origin in the model dynamics.

Phase-locking of three cells
A network of three cells is a minimal case that allows network asymmetry (of coupling or frequency) to trigger an effect of the even component even when the pairwise coupling is symmetric. In this context, the amplitude of the odd component (relative to the heterogeneity) is still the primary factor determining the existence of synchronous phase-locking, but the addition of an even component can modify the outcome dramatically, especially when the even component grows large relative to the odd component. The phase model for the three-cell network is assumed to take the general form of (3.10). We additionally assume pairwise symmetry of the coupling, k i j = k ji , and expand the phase difference equations for the network, To simplify calculations in this analysis, we restrict the interaction function to its first Fourier components, 5Ĥ φ = A odd sinφ + A even 1 − cosφ , (6.3) whereφ andĤ indicate phase in radians (note that we drop the hat notation below). We also impose the constraint that cosine be accompanied by a constant offset, H even ∝ 1 − cos φ, from the diffusive coupling condition H (0) = 0. We parametrize the amplitude of the even component by fixing the odd component and varying the amplitude ratio, Examples of the two-dimensional phase plane from equations (6.2) are shown in figure 6.1a, for the phase reduction of the resonate-and-fire model on a symmetric three-cell network (approximated in the form (6.3)). The intersections of the nullclines are fixed points of the system, corresponding to phase-locked states. For β = 0, the odd coupling strength A odd is set at the critical value for phase locking of the network. Oscillators 1 and 2 are closely locked, while 1 and 3 are locked at φ 1−3 ≈ π 2 . For small changes in the even component, β > 0 shifts the nullclines together to promote phase-locking, while β < 0 shifts them apart. Corresponding simulations of the full resonate-and-fire model network are shown in figure 6.1b. For β < 0, oscillators 1 and 2 remain entrained, but oscillator 3 slips past in relative phase.
In a general three-cell network with arbitrary connection weights, this shift of the nullclines combines with the synchronizing effect of the odd component and with the effective heterogeneity from coupling asymmetry to determine the presence or absence of phase locking. In the limit of small β, we can show analytically that the shift of the nullclines described above is the generic first-order effect on symmetric networks, and separates additively from the frequency and coupling heterogeneity terms when both are present. To do this, we first simplify the form of the nullclines by rewriting the coupling as a phase lag accompanied by an offset, Taking one of the nullcline equations from (6.2), we then combine all the coupling terms into a single sine function, capturing parametric dependence of the nullcline in the offset, amplitude, and phase of this effective interaction. where f α ≈ arctan −k 23 sin (φ 1−2 + β) 2k 31 + k 23 cos (φ 1−2 + β) .  The odd component affects only the term ω 1−3 A odd , decreasing the effect of the intrinsic frequency heterogeneity. This can be shown to increase the extent of the nullclines, promoting a fixed point as with the coupled pair. The frequency shift term ω e f f gives the effective heterogeneity from coupling asymmetry. Similar to the asymmetrically coupled pair (6.1), this supports or opposes phase locking depending on its sign relative to the intrinsic heterogeneity ω 1−3 .
The remaining terms represent the additional effects of the even component through their joint dependence on β and φ 1−2 , which for small β reduces to a dependence on the sum (φ 1−2 + β). In a network with symmetric coupling (ω e f f = 0), this dependence on (φ 1−2 + β) causes the φ 1−3 nullcline to shift along the φ 1−2 axis with a change in β (and likewise for the φ 1−2 nullcline with respect to φ 1−3 ). In the absence of frequency heterogeneity, the nullclines are straight lines φ 1−3 = 0 and φ 1−2 = 0 (see (6.2)), and thus are unaffected by these shifts. For larger frequency heterogeneity, the nullclines form closed curves and a phase-locked fixed point can be lost in a saddle-node bifurcation if the shifts move the nullclines apart. This can occur near the bifurcation for a small change in β of the proper sign (relative to the frequency heterogeneity) as shown in figure 6.1, or for larger changes in β regardless of sign. As β grows larger, the shift of the nullclines increases proportionally (accompanied by distortion from higher order terms not included in 6.5). An example where this larger shift eliminates the synchronous fixed point regardless of the sign of β is shown in figure 6.2.
To generalize from these minimal network examples to dynamics on larger networks, we investigate which effects of the even component occur generically for random networks versus contingent on the specific coupling and frequencies. The global synchrony of a large network will in some sense average the synchrony of random subnetworks, so we explicitly average over random three-cell networks. In figure 6.3 we plot the Kuramoto order parameter R 2 = 1 or coupling heterogeneity) are averaged out; we see no effect for small β of either the nullcline shifts or coupling asymmetry effects. For large β, we see a significant decrease in synchrony with both types of heterogeneity, as the even component effects begin to dominate over the intrinsic frequency heterogeneity (as in figure 6.2). This result resembles analytical results for chains of oscillators with spatially constrained coupling, where an increase in the critical coupling occurs for β ≥ π 4 [63,79]. It should be noted that in large networks the effect of the even component on synchronization can depend dramatically on the probability distribution from which the frequencies are drawn [62], whereas this dependence for random three-cell networks is minimal.

Even component of subthreshold interaction function
In the previous section we showed that the even component of the interaction function can have significant effects on phase locking. We now proceed to assess the magnitude of the even component for the resonate-and-fire oscillator. In figure 6.4, we plot the amplitude ratio β of the even component for the subthreshold interaction function. We note a dramatic difference between the positive and negative reset regimes (as seen with the odd component slope in figure 5.2). β is consistently large and negative in the negative reset regime (except for a small positive region for soft reset where the odd component is negative), and varies significantly across the positive reset regime. Although both the odd and even amplitudes decrease with v R (towards v R = −1), the odd component remains smaller, explaining the large amplitude ration in the negative reset regime. However, since the overall amplitude of the subthreshold interaction function is small, any effects of the even component here are likely to be dominated by the spike interaction, as discussed in section 5.3. In the positive reset regime, both the overall and relative amplitudes can be significant. β increases from negative to positive with increases in v eq , and is largest in magnitude at the spiking regime boundaries. This trend is driven by the component function C1 from (4.11), for which the coefficient A C1 ≈ A 2 v eq (for hard reset). The component function C2 also contributes to the even component, but to a lesser degree.
Where the even component is small in magnitude (figure 6.4, pale color region) the effects depend on the sign of β and cancel when averaged over random networks. This is unlikely to significantly affect large biological networks. For local adaptation to overcome this averaging and support or oppose synchrony, the even component of a connection would need to adapt based on frequencies and coupling strengths of both coupled cells, a possibility that seems biologically unrealistic. However, the even component is sufficiently large in a significant portion of the parameter space to potentially oppose synchrony in random networks (as in figure 6.3). In the neural context, for systems in which synchrony supports biological function, cells may need to adapt their dynamical properties to keep the even component small. This can occur most directly through shifting the equilibrium v eq , either by slower changes in conductances or synaptic weights, or by faster shifts in the tonic input to the cell. This mechanism could potentially enable rapid adaptive shifts (up or down) in the level of synchrony. Finally, it is possible that a large even component could instead promote specific functional states, such as a chimera state [48], rather than simply opposing synchrony.

Origin of the even component
In some of the mathematical literature on coupled oscillators, strictly odd-symmetric coupling is assumed as an easier to analyze default case [e.g., 21,24]. Our analysis shows that, for the resonate-and-fire model, strictly odd coupling is certainly not generic. To understand the factors determining the variation in the even component, it helps to consider the general case as arising from phase shifts relative to the purely odd special case. Consider an approximation in which the limit cycle, PRC, and resulting interaction function are all sinusoidal, differing only in phase shifts. (This is strictly true only in the small decay and large period limit.) The phase of the interaction function follows from the relative phase of v and Z v , determined by the boundary condition (3.14) for the adjoint equation. For the hard reset, if the limit cycle crosses the threshold at θ T ≈ 3π /2 (v eq ≈ 0), v and Z v are out of phase by ≈ π 2 , resulting in an odd interaction function H ≈ sin. Any phase shift in the PRC away from this produces a corresponding shift of the interaction function, introducing a nonzero even component. This can occur if v eq shifts away from 0, or from the soft reset boundary condition phase shift α.
A more geometrical perspective on the relative phase of the limit cycle and PRC is through the concept of dynamical shear, variation of angular velocity with radial displacement from a limit cycle. In a model without shear, perturbations perpendicular to the cycle do not cause phase shifts, so the limit cycle is normal to its isochrons and parallel to the (vector) PRC. This is close to the condition that components of the limit cycle and PRC be π 2 out of phase, and tends to lead to an odd interaction function for diffusive coupling. For instance, in the Stuart-Landau oscillator, a minimal model for shear about a limit cycle, the shear term in the dynamics directly scales a cosine term in the interaction function [3]. In the resonate-and-fire model, although there is no shear in the linear subthreshold dynamics, the effect of the threshold shifts the orientation of the isochrons relative to the limit cycle exactly like shear in the Stuart-Landau oscillator. This effective "reset-induced shear" is caused by perturbed trajectories on one side of the limit cycle crossing threshold earlier than those on the other side. In the soft reset case, this effect also depends on the geometry of the reset manifold. We can assess the validity of this explanation by evaluating a condition for no shear in the resonate-and-fire model, In figure 6.4, we show condition (6.6) for zero reset-induced shear as a dotted line. We see this intuition breaks down in the negative reset regime, where the odd component is extremely small. Otherwise, the condition closely approximates β = 0, with small variations that result from effects of the discontinuities in v and Z v not accounted for in this analysis.

Discussion
We applied the theory of weakly coupled oscillators to study the synchronization of resonateand-fire neurons coupled by electrical synapses. The use of a minimal hybrid model to capture the resonant dynamics allowed much of this analysis to proceed analytically. We calculated the phase reduction of the resonate-and-fire model neuron using the adjoint method for the PRC of hybrid models, following Shirasaka et al. [73]. We also presented a simplified derivation of their technique. We found that the resonant properties give rise to a potentially strong contribution of the subthreshold fluctuations to synchronization, in addition to the synchronizing effect of the spike. We also showed that, despite having no effect on coupled pairs, effects arising from the reset (i.e. reset-induced shear) have the potential to impair synchronization in certain network configurations.

Synchronization of Resonate-and-fire oscillators
Our analysis focused on the resonate-and-fire oscillator as an idealized model to study synchronization of electrically coupled resonant neurons. The phase reduction technique allowed us to separate the interaction function into components from spiking and subthreshold voltage fluctuations, dissecting their distinct contributions to synchronization. We showed that, in the resonate-and-fire model, the effect of a fast voltage spike transmitted through electrical coupling is always synchronizing, reinforcing previous work on pulse coupling of resonant neurons [34,46,23,76]. The contribution from subthreshold fluctuations generally promotes synchrony as well, but can actively oppose synchrony in a small region of parameter space. Where the spike and subthreshold synchronizing effects oppose one another, our model predicts that the net effect will depend linearly on the relative magnitude of the subthreshold fluctuations and the spike. Although our model does not explore the factors determining spike size, this can be inferred from more detailed biophysical models and applied to a more quantitative analysis of the resonate-and-fire spike effects.
In coupled pairs and other networks with extensive symmetry, synchronization is solely determined by the odd component of the interaction function. For the resonate-and-fire model, the subthreshold contribution to this odd component generally has a positive slope, promoting stable synchrony. This subthreshold contribution is small when the reset voltage v R is strongly negative (corresponding to after-hyperpolarization). A stronger contribution to the odd component supporting synchronization occurs when v R is well above threshold (corresponding to a plateau potential). The only significant departures from this rule are strong subthreshold effects near the boundaries of the spiking regime, including the small region with desynchronizing effects. However, the assumption of weak coupling breaks down near these bifurcations, so this conclusion should be verified by different methods of analysis.
Finally, we showed that in networks with less symmetry, significant reset-induced effects on synchronization can appear. The even component of the interaction function is often ignored, both because analyses focus on symmetrically coupled pairs (e.g., [51]) and because, as observed by Sakaguchi [70], it has complex "ambivalent effects on mutual entrainment." We analyze three-cell networks to show how the even component has the potential to oppose synchrony, especially when large enough to dominate over intrinsic frequency heterogeneity. In the resonate-and-fire model, the even component varies strongly with the equilibrium voltage, potentially interfering with the subthreshold synchronizing effect in parts of the positive reset regime. In general, any phase shift of the interaction function will introduce an even component; our derivation of the adjoint method for hybrid model PRCs clarifies a mechanism for such phase shifts linked to the post-spike reset. The boundary condition for the hybrid PRC determines the phase shift, dependent on the reset map (hard or soft reset) and on the geometry of the trajectory, threshold, and reset manifold. We characterize this effect as "reset-induced shear": a phase shift results when trajectories on one side the limit cycle cross threshold and are reset ahead of the limit cycle trajectory.

Comparison of resonator and integrator neurons
Taken as a whole, our results show that subthreshold resonance of model neurons can have a significant synchronizing effect in electrically coupled networks, contrasting with typical observations of integrator neurons. Previous work on single-variable integrate-and-fire models has found that the subthreshold effect of electrical coupling tends to oppose synchrony [50,67]. Because the reset voltage must be below threshold, the effect of the reset is desynchronizing [51] and tends to dominate the small synchronizing effects of other subthreshold fluctuations. Thus, in simple integrator models synchronization must rely on transmission of the spike only, functionally similar to pulse coupling from fast excitatory chemical synapses. In contrast, electrically coupled resonator neurons may combine pulse and continuous coupling to synchronize. This may help explain experimental observations of a loose correlation across brain regions between resonant properties of neurons and the prevalence of electrical synapses [37,66].
We found that the subthreshold contribution to synchrony is strongest in the plateau potential regime. Since the PRC is not dramatically different between the plateau and AHP regimes, our analysis suggests that the subthreshold synchronizing effect of resonance is mediated primarily by the temporal extent of the voltage fluctuations. In the plateau regime the subthreshold voltage waveform extends close to a full sinusoidal cycle, providing greater opportunity for exchange of current. This synchronizing effect likely extends beyond our resonate-and-fire analysis to the electrical coupling of other resonator neurons. Experimental results show plateau potentials in resonant neurons with widespread electrical coupling in the inferior olive [53,56], suggesting a potential synchronizing effect of the plateau. Synchronization of subthreshold oscillations in the absence of spiking [52] may also rely on a similar mechanism. Our predictions concerning resonance and subthreshold effects are directly testable experimentally, using pharmacological manipulation of resonant properties or dynamic clamp techniques to perturb and test single neurons and circuits, supplemented by analysis of detailed biophysical models.
We note, however, that integrator versus resonator is not a strict classification and does not always correspond directly with synchronization properties, despite the general trends observed. Although type I excitability (associated with a SNIC bifurcation), type I PRCs (Z v strictly positive), and the integration of input are often taken as loosely equivalent properties, Ermentrout et al. [28] clarified that systems near a SNIC bifurcation can have type II PRCs, with strong negative lobes. Additionally, Dodla and Wilson [20], analyzing synchrony based only on the shapes of the PRC and voltage fluctuations, emphasize that the type of PRC alone is insufficient to determine the synchronization of electrically coupled oscillators. Our work reinforces these results, showing that a resonator model can in certain regimes have integrator-like properties, both in the PRC and in the interaction function and synchrony.

Synchronization of hybrid model neurons
Our use of an idealized hybrid model neuron for this study necessitated modification of the PRC analysis techniques typically applied to biophysically detailed continuous models (section 3.4). The change in perspective from continuous to discontinuous dynamics may help provide new insights into basic questions of synchronization, such as whether resonance or other properties of neurons support the synchronizing effects of the spike. The discontinuous hybrid model PRC leads to a spike interaction function that is discontinuous at the origin, creating an especially strong syn-chronizing effect when the PRC jump is negative (positive for the interaction function). On the other hand, estimates of the (infinitesimal) PRCs from real neurons or biophysical models are continuous and approximately zero at the instant of spiking. If we smooth a hybrid neuron PRC to match these observations, the negative jump translates to a smooth peak skewed "rightward," to the latter portion of the PRC closely preceding the spike [51]. Realistic PRCs generally show this rightward skew, which gives a synchronizing positive slope to the interaction function, matching the strictly positive resonate-and-fire discontinuity. The skew has been shown to vary with adaptation in a range of models and experiments [23], including in hybrid models [67,46]. Future work could further link these adaptation skew effects and the resonance effects that we study, as well as bridging the gap between hybrid and continuous models by explicitly considering specific spike shapes along with the hybrid model dynamics.
For the subthreshold effects of electrical coupling, the interaction function depends on both the PRC and the limit cycle. This allows for significant variation in both odd and even components, even for the simple resonate-and-fire model. In other hybrid models, the odd component subthreshold effects have been shown to vary widely between models and with variation of parameters [67,17], consistent with the diversity of electrical coupling effects on synchrony (from combined spike and subthreshold effects) in many biophysical models [42,13,55]. Our work reinforces this observation for the odd component and also emphasizes similar variability in the even component, which likely generalizes to other hybrid models. Specifically, reset-induced shear is a newly identified factor that introduces a variable even component to the interaction function. Our analysis techniques allowed us to link this to effects of the reset map and the geometry of threshold and reset, and can be applied more generally to disentangle complex odd and even component effects in other hybrid models.
Despite the many advantages of hybrid models, little is formally known about the bounds on their validity. Ideally one should have a rigorous understanding of the hybrid model as a suitable asymptotic limit of related continuous models, tying the reset map to the detailed dynamics of spiking [41]. Some work has touched on this for simple cases: comparing input response [75,8], spiking transitions [22], or network spiking dynamics [35] between hybrid and biophysical models. Still, any direct comparison between idealized hybrid model dynamics and more complex biophysical models is challenging, whether at the single-cell or population level. Phase reduction provides a possible locus for such comparison, since detailed models can be phase-reduced computationally, translating them into the same "language" as our study of the resonate-and-fire model. The geometric insight into the PRC from our derivation of the adjoint method for hybrid models facilitates the interpretation of such comparisons. Our conclusions regarding the synchronization of resonant neurons can thus be verified and extended by comparisons with the computational phase reduction of detailed biophysical models and with the empirical phase response analysis of real neurons. A A.1 Connection to Shirasaka et al. [73] Here we will demonstrate that the boundary condition (3.14) for the hybrid model PRC across the reset discontinuity, which we derived in section 3.4, matches the result derived by Shirasaka et al. [73] following techniques from nonsmooth dynamical systems theory. The primary difference between the two results is that we present N − 1 boundary conditions, for a reset map R defined on the (N − 1)-dimensional threshold manifold, while Shirasaka et al. defined a reset map Φ in N dimensions (on an open neighborhood of the threshold) and presented N distinct boundary conditions for the adjoint problem. We show here that the N − 1 conditions corresponding to our result match exactly, and that the remaining condition simply enforces the normalization (3.9) (regardless of the definition of Φ off the threshold manifold).
Their result is formulated in terms of the saltation matrix C, where f v is the v-component of the dynamics,v is a v-direction unit vector, and T − and T + are the left and right limits of the boundary crossing. The Jacobian DΦ corresponds to our directional derivatives D u R along threshold. The row space of matrix M is thev direction only (perpendicular to the threshold), so for any component Z i along the threshold, M does not contribute to the boundary condition and (A.1) reduces to our result (3.14).
The remaining v-component boundary condition from (A.1) can be shown to simply enforce the normalization condition (3.9). We first evaluate this final component, We then expand (DΦ f (T )) · Z (T + ) as a sum over components Σ N i=1 f i (T − ) D i Φ · Z (T + ). If the Nth term is the v-component f v (T − ) D v Φ · Z (T + ), the remaining N − 1 components along the threshold reduce to Σ N−1 i=1 f i (T − ) Z i (T − ) according to (A.3).
Thus, we see that the final boundary condition does not depend on the specific definition of the reset map Φ off the threshold manifold, and simply enforces that the normalization condition f · Z = 1 holds across the reset.

A.2 PRC phase shift for soft reset
The general form of the resonate-and-fire PRC is Z v (t) = A r 0 e λt cos (t − T + α) , Z w (t) = A r 0 e λt sin (t − T + α) .
We show the phase shift evaluated over the full resonate-and-fire parameter space in figure A.1.

A.4 Slope of interaction function components
Here we evaluate the slope of each component of the resonate-and-fire interaction function, its contribution to the slope of the odd component of the interaction function, and expand the final result to first order in the decay parameter λ.