Synchronization of Electrically Coupled Resonate-and-Fire Neurons

2.1. Single-variable integrate-and-fire models.

The leaky integrate-and-fire model consists of a single voltage variable with linear subthreshold dynamics between spikes [67, 13, 1]. External current input I is integrated through changes in the neuron’s membrane potential (voltage) subject to a “leak,” or linear decay over time (governed by the membrane capacitance C, conductance g_L, and reversal potential v_L). 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 $v_{eq}=v_L+\frac{I}{g_L}$ is pushed over the threshold 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 classic model have been proposed. In a single dimension, nonlinear variants can approximate the effects of voltage-dependent ion channel conductances on subthreshold dynamics, more accurately modeling the approach to threshold. Prominent examples include the quadratic [35, 12] and exponential [42] integrate-and-fire models, as well as arbitrary nonlinearities derived via model reduction or optimized to match recorded spike trains [58]. For the analysis of synchrony, a number of significant theoretical discoveries have been made using networks of single-variable integrate-and-fire neurons coupled by chemical or electrical synapses [36, 78, 63, 21, 44, 69, 89, 85, 70, 73].

Despite these successes, some firing properties of real cells (and of more realistic models) are not captured by single-variable models, particularly fluctuations in subthreshold voltage following the spike. The diversity of post-spike dynamics found in different cell types includes fast (Figure 2.1a left) and slow (Figure 2.1a center) AHP, plateau potentials (Figure 2.1a right), and damped oscillations. Although variations in the level of fast AHP can be modeled to a limited degree in single-variable models [89], the other examples shown cannot be represented due to nonmonotonic voltage trajectories. To capture these requires increasing the dimensionality of the model.

2.2. Multivariable linear integrate-and-fire models.

The idea of linear feedback between two factors governing a neuron’s approach to the firing threshold was developed independently by some of the earliest mathematical neuroscientists, Rashevsky [91], Hill [49], and Young [112]. Adding a second “adaptation/facilitation” current variable [105, 92, 77] allows a model to reflect slow recovery processes in the neuron (typically the gating dynamics of ion channels) that can create resonance, adaptation, and a full range of post-spike voltage dynamics.

A multivariable “generalized integrate-and-fire” model can also be derived from a Hodgkin–Huxley-type biophysical model by linearizing and exploiting separation of timescales to eliminate variables [92, 60, 56]. In the simplest example, a single additional variable captures the most significant slow current, often a calcium, potassium, or hyperpolarization-activated h-current, or a combination of channels with similar time constants. (For higher-dimensional linear models, also see [77, 29, 57, 103].) The subthreshold dynamics of a two-variable linear model with slow current z can be described by [92, 60, 56]

This reduced model, with a suitable reset map included, has been shown to capture realistic spiking dynamics in a range of biologically relevant scenarios [92, 11, 60, 6]. We discuss the details of reset maps in two dimensions below (section 2.4), but for now we simply consider the properties of a limit cycle trajectory with reset point x₀ = (v₀, z₀).

In a two-dimensional model, the reset voltage v₀ can have a strong effect on the shape of the voltage trajectory, reproducing a range of post-spike dynamics as shown in Figure 2.1b. The voltage trace for a model reset below the firing threshold resembles a fast AHP (Figure 2.1b, left), as in the leaky integrate-and-fire model. In contrast, when the subthreshold dynamics have an oscillatory component, v₀ can be reset above threshold to resemble a plateau potential (Figure 2.1b, right) or at the threshold to resemble a slow AHP (Figure 2.1b, center). Additional comparison with these biophysical models is shown in section A.1.

Despite the relevance of these distinct voltage dynamics to electrical synaptic interactions, only a few studies have addressed the synchronization of electrically coupled hybrid models with more than a single variable [23, 22], and none have considered the effects of post-spike fluctuations specifically. We also note that nonlinear integrate-and-fire models in multiple dimensions can be remarkably effective at reproducing diverse spiking patterns [100, 54, 55, 10, 82] and have been applied to study synchronization through chemical synapses [26, 20, 79, 80, 38, 64, 65]. However, such nonlinear models are in general not analytically tractable. For this reason, we choose to explore the full range of post-spike dynamics in a linear two-variable integrate-and-fire model that permits direct analysis of synchronization.

2.3. Resonate-and-fire model.

The resonate-and-fire model is a two-dimensional linear integrate-and-fire model introduced by Izhikevich [53] with the goal of capturing essential features of the dynamics of resonant cells. Several studies have examined dynamics of the resonate-and-fire model in biologically relevant scenarios, including spiking triggered by periodic and noisy input [11, 92, 106, 60] and synchronization of pulse-coupled networks [79, 80, 81], and a variation on the model has been used to fit realistic spike train data [17]. Given the rich range of possible dynamics, we choose to analyze synchronization in the resonate-and-fire model rather than a more general two-dimensional linear integrate-and-fire form.

The subthreshold dynamics of the resonate-and-fire model are a slightly restricted canonical form of (2.2) in the damped oscillation regime (complex eigenvalues). Additionally, the model restricts the parameter space by assuming that a single time constant λ governs the decay of both variables, $\lambda =\frac{g_L}{C}=\frac{1}{\tau }$.¹ To simplify our notation, the slow current variable z is replaced by an effective adaptation/facilitation variable $w=\frac{1}{\Omega C}z$ with $\Omega =\sqrt{\frac{k}{\tau C}}$, voltage is measured relative to v_z, and the tonic current input is subsumed into the equilibrium voltage $v_{eq}=v_L-v_z+\frac{I}{g_L}$;

Without loss of generality, we set the threshold voltage to v_T = 0 and rescale time such that W = 1 (though allowing variations about this mean frequency in the network setting). Unless otherwise specified, we choose a fixed value for the decay parameter, λ = 0.1; our results are qualitatively similar as long as λ is sufficiently small (relative to W = 1). In the opposite extreme, taking the limit λ → ∞ with W fixed recovers the leaky integrate-and-fire model.² With these parameters fixed, the equilibrium voltage and details of reset determine the existence and properties of spiking in the model. If a trajectory starting from the reset point x₀ crosses the spiking threshold and is again reset, it forms a spiking limit cycle corresponding to regular spiking, shown in Figure 2.1b and Figure 2.2. We explore more detailed existence conditions in section 4.1, but a limit cycle will always exist for v_eq over the spiking threshold just as in the integrate-and-fire model (large constant current input).

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.³

where (r₀, θ₀) are the polar coordinates of the reset point of the limit cycle relative to equilibrium, such that (v₀, w₀) = (v_eq + r₀ cosθ₀, r₀ sinθ₀).

To build a network model of electrically coupled resonate-and-fire oscillators, we introduce diffusive coupling of the voltage variables, i.e., make the coupling linear in the voltage difference between the oscillators. Experiments have shown that electrical coupling current is typically linear and symmetric over a wide range of voltages, particularly for the common channel type connexin-36 [101, 7, 19]. Thus, we include pairwise coupling current terms I_ij (x_i, x_j) = k_ij (v_j − v_i) in the voltage dynamics, with coupling strengths k_ij = k_ji. However, note that exceptions to this form can arise as a result of gap junction channel rectification or asymmetry in cell properties. We briefly discuss the effects of asymmetry in section 6. We also allow parameters of the intrinsic dynamics to vary across cells, introducing heterogeneity that will oppose synchronization. The resulting network model is

We take the single-oscillator dynamics (2.3) with W = 1 as the mean dynamics of the (uncoupled) population, which is the focus of our analysis until we consider network structure in section 5. To complete our definition of the electrically coupled resonate-and-fire model, we need to fill in two details: the reset map in the two-dimensional phase space and the voltage spike itself.

2.4. Details of the reset map.

In multiple dimensions, the reset map need not map all trajectories crossing threshold to a single point as in (2.1). We distinguish a “hard reset,” which directly generalizes (2.1) by mapping to a single point x_R [30, 9], from a “soft reset,” which allows perturbations to persist over multiple cycles [77, 92]. We consider both hard and soft resets for the resonate-and-fire model. The resonate-and-fire hard reset rule specifies that if the threshold v_T = 0 is crossed from below, the state is instantaneously reset to the point x_R = (v_R, w_R).⁴

Most examples of soft reset models, like the adaptive exponential integrate-and-fire model [10], instead reset voltage to a fixed value while incrementing the adaptation variable. We define the resonate-and-fire soft reset rule as follows, with a fixed increment of the adaptation variable Δw:

The contrasting modeling choices reflected by these two reset rules can be understood with reference to ion channel gating variables in the Hodgkin–Huxley formalism, which approach a sigmoidal function of voltage at steady state. Gating variables with fast time constants during the spike approach their saturated limits as in the hard reset, whereas slow variables are incremented by an approximately fixed amount that depends on their respective time constants.

2.5. Details of voltage spikes.

The definition of an integrate-and-fire model (2.3)–(2.7) typically describes the subthreshold dynamics, threshold crossing, and post-spike reset, but not the suprathreshold portion of the spike itself. Thus, in order to capture the full effects of electrical coupling, with current transfer driven by fluctuations both during and between spikes, we must supplement the model by a description of the effects of the spike. We focus on narrow spikes and with minor assumptions are able to model the spike effects in a parsimonious form that cleanly separates the subthreshold and spike effects in our analysis (section 4.3). Our narrow spike limit is sufficient for most examples of fast sodium action potentials but would need to be extended to account for broader calcium spikes or extremely high frequency spiking.

We assume a cell is insensitive to input during and immediately following the spike, and thus the spike has a characteristic waveform v_spike (t) (0 < t < τ_spike), regardless of the state of any coupled cells. This assumption reflects observations in experiments and biophysically detailed models of a refractory period where perturbations typically have little to no effect, due to the many ion channels open during an action potential significantly lowering the input resistance of the cell [46].⁵

If we also assume the voltage during the spike is significantly greater than the voltage of any (nonspiking) coupled cells, the integrated current transfer from electrical coupling during the spike is approximated by the constant

If the time interval of the spike is small relative to the period (τ_spike ≪ T), the effect of this current on coupled cells is well approximated as an instantaneous transfer of electrical charge, or increment of voltage by Δv = q_spike [69, 89, 70]. We can equivalently represent this effect as a δ-function term in the coupling current for each time a coupled cell spikes (defining S_j as the set of spike times of cell j):

2.6. Notation for general hybrid models.

In order to illustrate the broad applicability of our analysis techniques, we introduce our methods using general notation for an arbitrary hybrid model. The threshold voltage v_T becomes a threshold manifold $\mathcal{T}$, with an arbitrary nonlinear reset map $R:\mathcal{T}\mapsto \mathcal{R}$ mapping to a corresponding reset manifold $\mathcal{R}$. Although a hybrid neuron model can in general consist of any number of dimensions with complex threshold manifolds and reset maps [77], a two-dimensional model is sufficient to illustrate the essential mathematical properties, so we restrict our analysis to models in two dimensions x = (v, w).

Note that for the resonate-and-fire model, the threshold manifold and soft reset manifold are both lines of constant voltage $\mathcal{T}=\left\{(v,w):v=v_T\right\}$ and $\mathcal{R}=\left\{(v,w):v=v_R\right\}$, while the reset manifold for hard reset is degenerate, a single point $\mathcal{R}=\left\{x_R=\left(v_R,w_R\right)\right\}$. We address the possibility of arbitrary orientations and locations of the threshold and reset manifolds in the resonate-and-fire model in section 4.2.

3.1. 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 5. Although a formal derivation assumes weak coupling and weak 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.

We note that the phase reduction approach can also be augmented by considering a second coordinate to capture motion transverse to the limit cycle, following a construction closely resembling that which we present here for phase [5, 107, 110, 98, 109], and which has recently been extended to nonsmooth dynamics [108]. In the case of the simple models we study here such a reduction would do little to simplify the dynamics, but in other cases a second coordinate can allow reduction under fewer conditions and may generate more realistic predictions.

The dynamics of a general system of coupled oscillators are described by

The state space of each individual oscillator is $x_i\in {ℝ}^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. We assume the average dynamics to have an asymptotically stable T-periodic limit cycle $\overline{x}(t)$ defined by

We have also restricted the pairwise coupling I_ij (x_i, x_j) to a common functional form defined by the coupling function I_c, which is weighted by connection strengths k_ij 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 $\Theta :{ℝ}^n\mapsto S^1$, 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 condition $\overline{x}(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.,

3.2. 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 cycle $\overline{x}(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 (iPRC or PRC, also called the “phase sensitivity function”) Z 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 limit-cycle 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 cycle $\overline{x}(t)$ and a nearby trajectory x(t) must be constant in time. That is, the separation $\Delta x(t)=x(t)-\overline{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, Df, 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.

3.3. 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 weakly perturbed according to $\frac{dx}{dt}=f(x)+ϵp(t)(\text{for }ϵ\ll 1)$, the perturbed phase θ= Θ (x (t)) satisfies

Since the perturbation is weak, its effects occur on a slow timescale, O(ϵ), which can be separated from the faster dynamics of unperturbed phase, $\frac{d\theta }{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.⁶ 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\approx \overline{x}\left({\theta }_i\right)\text{ (where }{\theta }_i=\Theta \left(x_i\right))$; therefore, we can approximate these perturbations (to order ϵ) as periodic functions of phase, $g_i\left(x_i\right)\approx {\tilde{g}}_i\left({\theta }_i\right)\equiv g_i\left(\overline{x}\left({\theta }_i\right)\right)$ and likewise for ${\tilde{I}}_c$.

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.

3.4. 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, Kurebayashi, and Nakao [98] 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 [8, 68, 3]. They prove that the solution to the adjoint problem with their boundary condition gives the PRC for the asymptotic phase of hybrid systems. Related results have been independently presented by Park et al. [87] for nonsmooth systems with discontinuous boundaries but no reset map and by Ladenbauer et al. [64] for a specific hybrid neuron model. We present a brief heuristic derivation for an equivalent adjoint boundary condition and show that the condition has an intuitive form tied to the geometry of the threshold. We assume the existence of a phase mapping in a neighborhood about the limit cycle which is differentiable everywhere except in the transverse direction across the threshold and reset manifolds. This is equivalent to assuming the existence of a well-defined PRC. For additional discussion of this assumption, see [98].

Consider a limit cycle trajectory $\overline{x}(t)$, reaching the threshold manifold at $\overline{x}(T)=x_T=\left(v_T,w_T\right)$, as shown in Figure 3.1. The phase on this trajectory cannot change across the instantaneous reset. That is, since $\Theta (\overline{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 $\mathcal{T}$ at x_T (Figure 3.1). Given the assumed differentiability of the phase mapping, 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. As long as the limit cycle is not tangent to the boundary (the transverse flow condition of [87]), (3.14) is independent of the normalization condition (3.9), and together they determine the unique solution to the adjoint problem for the discontinuous limit cycle of the hybrid model. Note that the derivative of the reset map is only defined in directions tangent to the threshold, so 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.⁷

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_uR (Figure 3.1).

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, i.e., bifurcations of the resonate-and-fire model.

4.2. 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 resonate-and-fire model with hard reset allows the PRC to be calculated directly from an analytical expression for phase [80]. Since all points on the threshold are mapped onto the single point x_R, the threshold can serve as a reference point to define the phase map for all trajectories. However, we choose instead to derive the PRC from the adjoint equation with 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=\frac{1}{\sqrt{1+{\lambda }^2}\text{cos}\left({\theta }_H-\alpha \right)}$ is determined by the normalization condition,

Based on the reset map, the boundary condition for the adjoint equation (3.14) links the left and right limits of the discontinuity, determining the phase shift α. Examples of the v-component PRC for both positive and negative v_R are shown in Figure 4.3.

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.6), mapping to a constant reset point, the derivative of the reset map is zero. Thus (4.9) 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.

For the soft reset rule (2.7), an increment of w, the boundary condition (4.9) 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,

This phase shift is typically small for positive v_R but can grow more significant in parts of the negative reset regime. (See an example in Figure 4.3 and the full calculation in section A.3.)

We note that in the more general case of a nonlinear reset map, the result depends on $\frac{d{R}_w}{dw}$, for which the hard and soft resets are the special cases $\frac{dRw}{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.9) 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, ${\theta }_H=\pm \frac{\pi }{2}$ is the grazing bifurcation condition (4.1). For the soft reset, the singularity occurs when $\text{tan}{\theta }_H=-\text{cot}\alpha =\frac{e^{\lambda T}-\text{cos}T}{\text{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, whereas 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.

4.3. Phase model of electrically coupled resonate-and-fire neurons.

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 resonate-and-fire network (2.5). We explore how the phase model changes as we vary the equilibrium voltage and reset parameters, representing different cell types and classes of firing dynamics.

The interaction function H and heterogeneity ω of the phase model are calculated according to (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 represents a combination of all parameter heterogeneity in the resonate-and-fire network model, as described by the integral (3.11). We assume a fixed level of frequency heterogeneity (relative to the mean ω = 1) and thus focus on how parameter changes affect the interaction function H.

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 function I_c = v (t + ϕ)−v (t)+Mδ (t + ϕ) (cf. (2.8)) contains terms for the subthreshold voltage difference and the spike of the coupled cell. The resulting interaction function, depicted in Figure 4.4, is

Since the interaction function convolves the PRC and the coupling function, with a linear dependence on both, we can additively separate the corresponding subthreshold and spiking components, H = H_sub + H_spike:

The spike (δ-function) component of the coupling function 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 [78, 63, 79]. In the narrow spike limit that we study, the effect is entirely determined by the PRC and the magnitude of the spike, and 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. Analysis of this component allows 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.

The interaction function integral for the subthreshold component results in a complex expression, in part because the integral must be split due to the discontinuity of the orbit (full calculation in section A.4). To simplify the analysis, we split the subthreshold interaction function into three terms with distinct parameter dependence:

where θ_T = θ₀+T is the angular coordinate of the trajectory at threshold. We note that the S-function 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.

5.1. General considerations.

We begin by examining the synchronization of a symmetrically coupled pair of oscillators, the simplest and most commonly analyzed network [36, 62, 69, 74]. 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 [7]. 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₁₂ = k₂₁ = K) is given by

Expressed in terms of the phase difference ϕ_i−j = θ_i − θ_j and frequency difference Δω_1–2, we can rewrite (5.1) in terms of the odd component of the H function, H_odd(ϕ) = ½(H(ϕ) − H(−ϕ)) = H(ϕ) −H_even(ϕ), as follows:

(Some analyses refer to the G-function, G(ϕ) = −2H_odd (ϕ) [96]). 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 (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 electrically coupled 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}^{\prime }(0)$ is positive, and the reverse holds for negative slope. (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}^{\prime }(0)$ to always be well defined.

We note that in larger networks a similar dependence on the slope exists. Typically, a strong positive slope leads to global synchrony, while a negative slope leads to global incoherence [2]. Nonetheless, in some cases global features of the interaction function H (ϕ) can have significant effects of synchronization; the even component effects we discuss in section 6 are one example.

As the frequency heterogeneity of the pair increases, the pair of fixed points shifts 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 near-synchronous state is approximately inversely proportional to the slope, ${\varphi }_{1-2}\approx \frac{\Delta {\omega }_{1-2}}{2K{H}_{odd}^{\prime }(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}^{\prime }$,

(Note that here we use prime to denote the slope of both functions, so H′ is a time derivative and ${\widehat{H}}^{\prime }$ is a phase derivative.)

For larger heterogeneity, the phase difference continues to increase until Δω_1–2 is greater than the amplitude of 2KH_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 phases 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: for fixed coupling K, the slope ${\widehat{H}}_{odd}^{\prime }(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. Due to the noise and heterogeneity present in biological systems, we consider robustness to be an important criterion for synchronization in addition to the stability of the synchronous state for identical oscillators. Below, we consider the robustness and stability of near-synchronous phase-locking in more detail for the subthreshold and spiking components of the interaction function for electrically coupled resonate-and-fire neurons.

5.2. Subthreshold contribution.

For the analysis of synchrony in the electrically coupled 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 spiking component in section 5.3. For the subthreshold interaction function the odd component is typically dominated by the first Fourier component ${\widehat{H}}_{odd}(\widehat{\varphi })\propto \text{sin}(\widehat{\varphi })$, so its slope and amplitude, our two measures of robustness, are roughly proportional (see section A.6). 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 split the calculation of the slope of H_sub into components, following the decomposition (4.13) of component functions with approximate symmetry (see Figure 4.5). We consider the limit of small decay (λ ≪ 1) by expanding to first order in the decay parameter (see calculation in section A.5).

The O(1) in λ contribution to ${\widehat{H}}_{odd}^{\prime }(0)$ from S is the primary subthreshold factor determining the strength and stability of synchrony (as shown in Figure 4.5). 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 the boundary of this stable region differs only to order λ from the PRC singularity condition (4.11),

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 either from A_S < 0 or from A_S ≈ 0 and A_C1 < 0 or A_C2 < 0.

We show the slope ${\widehat{H}}_{odd}^{\prime }(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 ${\widehat{H}}_{odd}^{\prime }$ 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 regimes. 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.1a, 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 through electrical coupling. In the AHP regime, on the other hand, synchrony is stable for identical intrinsic frequencies but minimally robust to heterogeneity. Because there is little to no subthreshold contribution to synchrony in this regime, the model can only synchronize in the presence of moderate frequency heterogeneity through the effect of the transmitted spike (see section 5.3).

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 ${\widehat{S}}_{odd}^{\prime }\approx T-\text{sin}T$ increases sharply with an increase in the period T. That is, increasing the reset voltage v_R increases the angular extent of the limit cycle as it covers more of the (v, w) plane. Note that with W = 1 fixed, T represents two distinct factors: the extent of the cycle in radians and its duration. However, varying the duration alone by simply rescaling time has no effect on the strength or robustness of synchrony. The effects represented by H_odd are entirely due to the extent of the limit cycle covering a larger fraction of a cycle of continuous oscillation.

5.3. 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 [47].

In the narrow spike limit implied by the δ-function spike, the interaction function convolution equation (3.12) has a simple form for the spiking component, namely a time-reversed copy of the PRC.⁸

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.⁹ The jump discontinuity can stabilize the fully synchronous state even with nonzero heterogeneity if |Δω_1–2| < H_odd (0⁺) = ΔH_spike [28, 70, 98]. (For more in-depth analysis of the limit approaching discontinuity in the interaction function, see [98, 52].) We evaluate this discontinuity directly from the PRC:

The condition ΔH_spike > 0 is always satisfied for stable spiking limit cycles in both hard and soft reset conditions. Specifically, e^λT cos α − cos (T − α) > 0 holds for both hard and soft reset (α = 0 and the phase shift (4.10), respectively), while the scaling factor A changes sign with cos (θ_H − α) at the stability boundary (4.5). We evaluate the value of ΔH_spike over the full parameter space in section A.7. 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}^{\prime }(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.

We note that without the assumption that coupling effects are blocked during the spike (discussed in section 2.5), the spike interaction function would have an additional constant term from Z_v (0), the sensitivity to the coupling current that leaves the cell during its spike. The resulting constant offset to the interaction function would only affect synchronization in asymmetrically coupled networks, shifting the effective frequency of each cell by an amount proportional to the total strength of its connections.

6.1. 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 amplitude of the even component is large relative to that of the odd component.

The phase model for the three-cell network takes the general form of (3.10). We assume pairwise symmetry of the coupling, k_ij = 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,¹⁰

where $\widehat{\varphi }$ and $\widehat{H}$ indicate phase in radians (note that we drop the hat notation below). We also impose the constraint that the 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 (6.2) are shown in Figure 6.1a for the phase reduction of the electrically coupled 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 ${\varphi }_{1-3}\approx \frac{\pi }{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.

with the parameters approximated to order β as

The amplitude of the odd component appears only in the term $\frac{\Delta {\omega }_{1-3}}{A_{odd}}$, decreasing the effect of the intrinsic frequency heterogeneity (i.e., larger coupling is equivalent to smaller heterogeneity). When the coupling asymmetry term is small, decreasing the frequency heterogeneity term can be shown to increase the extent of the nullclines, promoting a fixed point as with the coupled pair. The frequency shift term ω_eff 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 (ω_eff = 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.

We next investigate which effects of the even component occur generically for random networks versus contingent on the specific coupling and frequencies. This may provide insight into larger network dynamics, to the extent that global synchrony of a large network can act to average the synchrony of random subnetworks. In Figure 6.3 we plot the Kuramoto order parameter $R^2=\frac{1}{N_t}{\sum }_t{|\frac{1}{N_j}{\sum }_j{e}^{i{\varphi }_j(t)}|}^2$, a measure of the strength of synchrony, averaged over time and over instantiations of random frequencies and coupling heterogeneity. Effects that depend on the sign of β (relative to the frequency or coupling heterogeneity) are averaged out; we see no effect for small β of either the nullcline shifts or the coupling asymmetry effects. For sufficiently 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 $\beta \ge \frac{\pi }{4}$ [84, 111]. 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 [83], whereas this dependence for random three-cell networks is minimal.

6.2. 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 electrically coupled 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 ratio 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.13), for which the coefficient $A_{C1}\approx \frac{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, green region) the effects depend on the sign of β and cancel when averaged over random networks. Small even components are thus unlikely to significantly affect large biological networks. For the tuning of local connection properties to overcome this averaging and support or oppose synchrony, the even component of a connection would need to be tuned based on the 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 [66], rather than simply opposing synchrony.

6.3. Origin of the even component.

Due to a predominant focus on odd component effects and symmetric coupling contexts, the literature on neural synchrony has little discussion of the factors determining variation in the even component. Thus, with the goal of illuminating this broader question, we present a brief analysis for the electrically coupled resonate-and-fire model, explaining the origin of phase shifts of the interaction function relative to the purely odd case. We 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 $\approx \frac{\pi }{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 parallel to the (vector) PRC. This is close to the condition that components of the limit cycle and PRC be $\frac{\pi }{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 [4]. 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 vector PRC 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. Radial perturbations have no effect on phase when the PRC is perpendicular to the radial direction, or parallel to the PRC,

In Figure 6.4, we show condition (6.6) for zero reset-induced shear as a dotted line. We see that 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.

7.1. Synchronization of electrically coupled 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 with a wide range of post-spike dynamics. 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 focus primarily on the subthreshold contribution, which varies strongly with the parameters of the reset. In coupled pairs and other networks with extensive symmetry, synchronization is solely determined by the odd component of the interaction function. For the electrically coupled 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 AHP. A stronger synchronizing contribution 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.

In our analysis of the resonate-and-fire model spike contribution, we showed that the effect of a fast voltage spike transmitted through electrical coupling is always synchronizing, reinforcing previous work on pulse coupling of resonant neurons [47, 64, 33, 102]. 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, a more quantitative analysis of the resonate-and-fire spike effects could integrate spike size measurements obtained from experiments or from detailed biophysical models.

Finally, we showed that in networks with sufficient asymmetry, significant reset-induced effects on synchronization can appear. The even component of the interaction function is often ignored in the neural context, both because analyses focus on symmetrically coupled pairs (e.g., [70]) and because, as observed by Sakaguchi, Shinomoto, and Kuramoto [94], it has complex “ambivalent effects on mutual entrainment.” (However, see [61, 86].) We analyze three-cell networks of generic phase oscillators to show how the even component has the potential to oppose synchrony, especially when large enough to dominate over intrinsic frequency heterogeneity. In the electrically coupled resonate-and-fire model, the even component varies strongly with the equilibrium voltage or applied current, 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 of the limit cycle cross threshold and are reset ahead of the limit cycle trajectory.

7.2. 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, by enabling significant post-spike voltage fluctuations not typically found in integrator neuron dynamics. Previous work on single-variable integrate-and-fire models has found that the subthreshold effect of electrical coupling tends to oppose synchrony [69, 89]. Because the reset voltage must be below threshold, the effect of the reset is desynchronizing [70] 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, as with the resonate-and-fire model in the more integrator-like AHP regime. The observation that electrically coupled resonator neurons may additionally rely on continuous subthreshold coupling to synchronize may help explain the loose correlation across brain regions between resonant properties of neurons and the prevalence of electrical synapses found in experiments [51, 88].

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 [72, 75], suggesting a potential synchronizing effect of the plateau. Synchronization of subthreshold oscillations in the absence of spiking [71] 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.

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 saddle-node on invariant circle (SNIC) bifurcation), type I PRCs (Z_v strictly positive), and the integration of input are often taken as loosely equivalent properties, Ermentrout, Glass, and Oldeman [39] clarified that systems near a SNIC bifurcation can have type II PRCs, with strong negative lobes. Additionally, Dodla and Wilson [28], 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.

7.3. Synchronization of hybrid model neurons.

The shift in perspective to discontinuous hybrid model dynamics enables insights not possible with biophysical models, cleanly separating distinct components of the electrical coupling effects. However, relating these results back to phenomena in more realistic models is critical for understanding the broader implications, especially where the two representations seem to diverge. The discontinuous hybrid model PRC leads to a spike interaction function that is discontinuous at the origin, whereas estimates of the (infinitesimal) PRCs from real neurons or biophysical models are continuous and approximately zero at the instant of spiking. However, if we smooth a hybrid neuron PRC (or measure the finite-perturbation direct PRC; see Figure A.1a), the negative jump translates to a continuous peak skewed “rightward,” to the latter portion of the PRC closely preceding the spike [70]. Realistic PRCs generally show this rightward skew, which gives a synchronizing positive slope to the interaction function, matching the synchronizing effect of the positive resonate-and-fire discontinuity. The skew of the PRC has been shown to vary with adaptation in a range of models and experiments [33], including in hybrid models [89, 64], related to the effects of resonance and reset on the discontinuity in our study. Future work can bridge the gap between hybrid and continuous models by further exploring shared features such as skew in the PRC, as well as considering the effects of wider spikes along with the hybrid model dynamics. To accomplish this will require a careful treatment of the small spike width and instantaneous reset limits. Our assumption of no current exchange during the spike, for example, could be replaced by the more realistic case of low susceptibility following from the model dynamics.

For the subthreshold effects of electrical coupling, the interaction function depends on both the PRC and the limit cycle, allowing for significant variation in both odd and even components to be captured in the simple resonate-and-fire model. Our work reinforces observations that subthreshold electrical coupling effects can vary widely based on intrinsic properties in both hybrid [89, 23] and biophysical models [59, 18, 74], and also emphasizes additional variability in the even component through the reset-induced shear effect. Although the hard versus soft reset distinction does not have a direct analogue in biophysical models, the fact that our resonate-and-fire results are unchanged for these two limiting cases suggests that they hold for other reset maps in this range–-even nonlinear reset maps. Our analysis techniques could be directly extended to explore these general reset properties in more depth, potentially linking them back to features of biophysical models in scenarios where they do have direct effects.

Given the many advantages of hybrid models for both analytical understanding and computational efficiency, knowing their limits is critical, but attempts at derivations from more complex biophysical models are challenging and typically give results only for narrow cases [58, 52]. Several studies instead address direct validation of dynamic properties, comparing input response [100, 11], spiking transitions [31], or network spiking dynamics [48] between hybrid and biophysical models. Phase reduction can also provide a 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. Our separation of spike, subthreshold, and reset effects can also provide insight into biophysically derived phase models, in terms of sets of ion channels that tend to be active in specific segments of the limit cycle, e.g., sodium and potassium for spiking dynamics, slow potassium currents and the hyperpolarization-activated HCN channel for subthreshold dynamics, and calcium-activated currents (KCa or CAN) for post-spike reset dynamics. 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.