Definition and Methodology

Pulse coupled oscillators are assumed to interact with each other
at discrete times such that the coupling can be described as an
effect multiplied by a delta function. The effect is usually
captured as a phase response curve (PRC). This formulation allows
the time evolution of the coupling to be described by a map from
one cycle to the next. The assumption of instantaneous phase
resetting requires that the time for the resetting to be complete
is short compared to the time between input pulses. The
perturbation used to tabulate the PRCs in this context is a change
in synaptic conductance evoked by a presynaptic action potential
or burst of action potentials, depending upon the type of
oscillators under study (Canavier, 2005). Although the pulse coupled methods and weak coupling methods (Ermentrout and Kopell, 1991; Hansel et al., 1995 and Netoff et al., 2005b) both utilize the phase response curve, the weak coupling method convolves a perturbation with the infinitesimal PRC whereas the pulse coupled method simply uses the perturbation itself to generate the phase response curve.

Application of the PRC to a Periodically Forced Oscillator

Perkel et al. (1964) measured the phase resetting function \(F(\phi)\ ,\) in invertebrate neurons that act as pacemakers in that they spontaneously and repetitively fire action potentials in a periodic manner. They induced precisely timed inhibitory or excitatory synaptic potentials in the pacemaker cell by evoking an action potential in an excitatory or inhibitory presynaptic neuron. \(F(\phi)\) tabulates the change in cycle length as a
function of the phase at which the perturbation was applied. This function is often normalized by the intrinsic unperturbed period (see Phase Response Curve). Note that we use the opposite sign convention than the one used in that article in order to be consistent with
the literature on this subject. Perkel et al. (1964) then used the PRC to predict stable \(m:1\) entrainment of the pacemaker neuron by a periodic train of evoked synaptic potentials. In a \(1:1\) phase-locked entrainment, the period of the forcing (\(P_F\)) must be equal to the intrinsic period of the forced neuron (\(P_i\)) plus any phase resetting at the equilibrium phase (\(\phi^*\)) of the perturbation\[P_F = P_i \{ 1 + F(\phi^*) \}\ .\] For \(m:1\) lockings in which only every \(m^{th}\) cycle contains a perturbation the appropriate expression is \( P_F = (m-1) P_i + P_i \{ 1 + F(\phi^*) \}\ .\) Therefore, all possible \(m:1\) lockings can be determined from the PRC alone. Perkel et al. (1964) also determined which lockings were stable by using a map based on the PRC. The stimulus times
\(ts[n]\) indicate the time elapsed between an action potential in the pacemaker and the arrival of a presynaptic potential due to the periodic forcing. For \(1:1\) forcing, the map is generated as follows:

A related idea is to force an oscillator with a feedback input that has a fixed delay from each time that the oscillator fires (Foss and Milton, 2000). The stability result in this case is \(-1/k< F'(\phi) < 1\ ,\) where \(k\) is the number of feedback inputs that
have been received between the time a neuron fires at time \(ts\) and the feedback input that is due to the spike at time \(ts\ .\)

Application of the PRC to Two Coupled Oscillators

The next level of complexity is two reciprocally coupled oscillators. Peskin (1975) considered synchronization in a network of two mutually excitatory pulse-coupled leaky integrate and fire neural oscillators.
\[
dV/dt = S_O - \gamma V(t)
\]
The pulsatile coupling resulting from the firing of one oscillator caused the voltage of the second oscillator to instantaneously become more depolarized by a fixed amount.
\[
V_i(t) = 1 \, \Rightarrow \, V_j(t)=V_j(t)+\epsilon \,\, \forall j \neq i.
\]

This implicitly defines a phase resetting curve because if a perturbation is received at time \(t\ ,\) the next spike will be advanced by an amount that can be calculated using the explicit solution for the differential equation that describes the leaky integrate and fire oscillator.
\[
F(\phi) = ln(1 - \gamma \epsilon S_O^{-1} e^{A \phi})/A
\]
where \(A = ln(S_0/(S_0 - \gamma))\ .\) Peskin derived a return map for the phase of the non-firing oscillator immediately after its partner fired, with the oscillators reversing roles
on each map iteration. This mapping was proven to have a unique, unstable fixed point at a phase of 0.5. This fixed point repelled trajectories toward synchronization at a phase of zero or one. Once the oscillators synchronize, the coupling term drops out because each neuron is already at threshold when its partner fires, hence synchronization results for every set of initial conditions.

Mirollo and Strogatz (1990) generalized these results to any two coupled oscillators in which the state variable (i.e. the membrane potential \(V\)) is a smooth monotonically increasing and concave down function, as it is for the leaky integrator. The pulse coupling when one oscillator fired depolarized the other by an amount \(\epsilon\) or pulled it up to the firing threshold, whichever was less. That is,
\[
V_i(t) = 1 \, \Rightarrow \, V_j(t)=min(1,V_j(t)+\epsilon) \,\, \forall j \neq i.
\]
A return map was again formulated in terms of the phase of the non-firing oscillator immediately after its partner fired. This mapping was proven to have a unique fixed point which was unstable. Once again the firing is globally repelled toward synchrony where the coupling term disappears allowing the synchronized state to exist.

The work presented so far on two identical coupled oscillators focused exclusively on excitation that leads to synchrony, that is, a phase locking with no phase difference between the oscillators. For identical oscillators in which the coupling term drops out at a phase of both \(0\) and \(1\ ,\) synchrony is always a solution. However, other solutions are possible, and for systems in which the phase resetting does not disappear at \(0\) and \(1\) synchrony may not be a solution. Therefore criteria are required for both existence and stability. The criteria given below (Dror et al., 1999 and Oprisan et al., 2004) are formulated for oscillators that are not necessarily identical. The existence criteria may be obtained from the algebraic relationships between the intervals at steady state, whereas a mapping similar to the one presented in the section on forced oscillators is required.

Figure 2: Stability Analysis for two coupled oscillators.

For \(1:1\) phase locking, the map based on the PRCs for two pulse coupled
oscillators is generated as follows:
\[
P_1 \{ 1 + F_{1}(\phi_1[n])\} + ts_1[n+1] = ts_1[n] + P_2 \{ 1 + F_{2}(\phi_2[n])\}
\]
\[
ts_1[n+1]-ts_1[n] = P_2 \{ 1 + F_{2}(\phi_2[n])\}-P_1 \{ 1 + F_{1}(\phi_1[n])\}.
\]
Here, \(F_{i}(\phi_i[n])\) corresponds to the first order phase resetting (see Phase response curve) of the \(i^{th}\) oscillator in response to the input received in the form of a pulse from another oscillator. In the neighborhood of \(\phi_i^*\ ,\) \(F_{i}(\phi_i[n])= F_{i}(\phi_i^*) + F_{i}'(\phi_i^*) \Delta \phi_i[n]\) where \(\Delta \phi_i[n] = \phi_i[n] - \phi_i^*\) and \(F_{i}'(\phi_i^*)\) is the slope of the first order PRC at the locking point \(\phi_i^*\ .\) The stimulus times can also be rewritten\[ts_i[n]=P_i \{ \phi_i^* + \Delta \phi_i[n] \}\ .\] At steady state both the oscillators are phase-locked at a common period (\(P_e\)). Therefore, \(P_i \{ 1 + F_{i}(\phi_i^*)\}=P_e\) for \(i=1,2\ .\) Using this and substituting for the stimulus times in terms of phase with further simplification produces:
\[
\Delta \phi[n+1] = \{ F'_{1}(\phi_1^*)-1 \}\{
F'_{2}(\phi_2^*)-1 \} \Delta \phi[n]
\]
If \(|{\{ F'_{1}(\phi_1^*)-1 \}\{ F'_{2}(\phi_2^*)-1 \}}| < 1\) then \(\Delta \phi[n+1]\) goes to zero and the locking at \(\phi^*=(\phi_1^*,\phi_2^*)\) is stable. Delays can easily be incorporated into this formalism by including them in the definition of the appropriate interval (Oprisan et al., 2004).

The above derivation did not apply to exact synchrony because the assumed firing order could not be guaranteed for a small perturbation from synchrony. Goel and Ermentrout (2002) obtained analogous results with a slightly different approach, and extended them to exact synchrony
in identical oscillators by making additional assumptions.
The phase transition curve \(f(\phi)=\phi - F(\phi)\) was used to define the return map
\[
G(\phi) = 1 - f(1-f(\phi))
\]
Linearizing the map and taking the derivative with respect to \(\phi\) yields
\[
G'(\phi) = [1 - F'(\phi)][1 - F'(1 - \phi + F(\phi))]
\]
They assumed that \(F(0)=F(1)=0\) so that synchrony is always a solution. Then the assumption was made that \(f(\phi)\) is monotonically increasing so that \(f'(\phi) > 0 \) and \([1 - F'(\phi)] > 0\ .\) This ensures that the firing order cannot change so that one oscillator cannot fire twice in a row. The map then has a stable fixed point at synchrony if \([1 - F'(0^+)][1 - F'(1^-)] < 1 \ .\) The use of the left and right slopes acknowledges the possibility that the first order PRC is discontinuous at \(0\) and \(1\ .\) The source of this discontinuity is often (Maran and Canavier, 2007) the existence of second order resetting (the effect on the second cycle after the perturbation see Phase Resetting Curve). In a network with all to all connections and a monotonically increasing phase transition function, the firing of one neuron cannot cause any changes in the firing order because if the old phase of one neuron was greater than the old phase of another neuron, the new phase will also be greater. The assumption of a monotonic phase transition implies that \(F'(\phi) < 1\ ,\) therefore \((1 - F'(\phi)) < 1\ ,\) and taking the absolute value of the product is no longer necessary. For modes other than synchrony, the assumption of a monotonically increasing phase transition function can be dropped because the firing order is constant under small perturbations without this assumption.

Application of the PRC to a Unidirectional Ring of Coupled Oscillators

Dror et al. (1999) considered the case of a unidirectional ring. Each oscillator receives exactly one input per cycle. Stable phase-locked modes must satisfy the following criterion:
\[
\sum_{i=1}^N{P_i\phi_i^*} = j*P_e
\]
where \(P_i \{ 1 + F_{i}(\phi_i^*)\}=P_e\) and \(ts_i[n] = P_i \{ \phi_i^* + \Delta \phi_i[n] \}\) as before and \(j \in [0,N-1]\ .\) Each value of \(j\) produces a different firing order. \(N+j-1\) cycles elapse before
the feedback due to a spike in a particular oscillator propagates back to that oscillator,
similar to the integer \(k\) in Foss and Milton (2000). Let \(\Delta \phi_i[n+1]\) denote the change in phase of the \(i^{th}\) oscillator in the \((n+1)^{st}\) cycle. \(\Delta \phi_1[n+1]\) in a ring of \(N\) coupled oscillators using PRCs is fully described by the \(N-1\) state variables \(\Delta \phi_q[n]\) where \(q \in [1,N-1]\ .\) Here, unlike the case of the two coupled oscillators, we have
\[\tag{1}
\Delta \phi_1[n+1] =A \,\Delta \phi_1[n]
\]

where \(A\) is the Jacobian matrix of the discrete-time system determined by the PRCs (Dror et al., 1999). The matrix \(A\) has been worked out in (Dror et al., 1999) where time was used instead of the phase variable. The matrix \(A\) is unchanged for identical oscillators and is given in (Oprisan and Canavier, 2001) for heterogeneous rings of three. The matrix \(A\) is distinct for each assumed firing order.

In the case of \(N\)-ring coupled oscillators the eigenvalues of the matrix \(A\) determine the stability of the system given the values of \(F_{i}(\phi_i^*)\) for \(i=1,\ldots,N\ .\) If \(\lambda_{max}\) denotes the maximum eigenvalue and if \(|\lambda_{max}| < 1\) stability is guaranteed. These existence and stability criteria predict the activity of model neural networks accurately as long as the duration of the coupling is short compared to the period (Canavier et al., 1997 and 1999; Luo et al., 2004).

Goel and Ermentrout (2002) used PRCs to study similar geometries such as a bidirectional ring or chain of \(N\) coupled oscillators. The firing order cannot be assumed to be resistant to all small perturbations in such a system so the issue of stability for this general case remains as an open question.

Application of the PRC to N all to all Coupled Oscillators

Mirollo and Strogatz (1990) formulated a map for \(N\) all-to-all coupled identical oscillators (i.e. each oscillator coupled to all the others) with assumptions analogous to those for their two neuron circuit. As the system evolved in time oscillators formed groups that
fired at the same time. A larger group tended to absorb other oscillators until a single group remained. The proof that their map converged to synchrony was performed in two stages. The first part showed that for almost all initial conditions (i.e. state of all the oscillators), an absorption occurred in finite time. The second part of the proof showed that the measure of the sets of initial conditions that lived forever without experiencing the final absorption to synchrony was zero.

Goel and Ermentrout (2002) also studied \(N\) all-to-all pulse-coupled identical oscillators in which every oscillator was identically connected to every other oscillator. They again assumed that \(F(0)=0\) and \(F(1)=1\ ,\) which guarantees the existence of a synchronous solution, and that \(f'(\phi) > 0\ ,\) which guaranteed that the oscillators always fired in the same order. The oscillators were indexed in the order of their firing and an \((N-1)\)-dimensional map was formulated. Stability conditions were derived by linearizing the map about its fixed point. The stability conditions were derived based on whether the PRC was continuously differentiable or not. If \(F'(0^+) = F'(1^-)\ ,\) then stability was guaranteed if \([1 - F'(\phi)] < 1\ .\) Otherwise the following \(N - 1\) conditions must be satisfied \[ [1 - F'(0^+)]^a [1 - F'(1^-)]^{N-a} < 1\ ,\] where \( 1 \leq a < N\ .\) In the latter case, stability depends upon the size of the network and synchrony can be destabilized for large \(N\ .\)

Validity of the Assumptions

In general, pulsatile coupling methods can be extended to inputs that are not strictly pulsatile under the assumption that the effect of a perturbation has dissipated by the time the next one is received. The likelihood of a violation of this assumption is much greater in networks comprised of oscillators that receive multiple inputs per cycle compared to those that receive only a single input per cycle.

The condition that the firing order is always preserved as a fixed point is
approached is not always met in networks of model neurons (Maran and Canavier, 2007). Stability
in the case of changing firing order is an open problem because the return maps only work with a particular assumed firing order. Also, sometimes the effect of a perturbation does not dissipate before the neuron fires (or bursts) next. In this case the existence and stability criteria must include higher order terms. The contribution of second order resetting has been incorporated for the two neuron case (Oprisan et al., 2004) and for the unidirectional ring
(Oprisan and Canavier, 2001).

Examples of Applications to Biological Neurons

The existence and stability criteria for two nonidentical pulse-coupled oscillators were tested in a hybrid network of one bursting model neuron and one physiological bursting neuron (Oprisan et al., 2004). These networks were constructed using the Dynamic Clamp using reciprocal inhibition. The predictions were both qualitatively and quantitatively accurate within the limits set by experimental noise. For these bursting neurons, the contribution of second order resetting to the existence criterion was essential in some cases.

Netoff et al. (2005a) utilized similar methods in hybrid networks comprised of either two biological neurons or one model neuron and one spiking neuron. Stellate cells and oriens-lacunosum-moleculare interneurons (and their model counterparts) that fired single spikes
were used to construct reciprocally excitatory, reciprocally inhibitory, and heterogeneous (excitatory-inhibitory) two neuron networks. This study incorporated noise into their predictions of firing interval histograms. Their predictions were accurate without incorporating second order resetting. The method in this study and in that of Pervouchine et al. (2006) is called by the authors as spike time response curve (STRC) method, but it is identical to the pulsatile coupling techniques described in this article except that time is used rather than phase. This is possible as long as there is no second order resetting.
Pervouchine et al. (2006) extended the work of Netoff et al (2005) to examine the contribution of a specific conductance (the hyperpolarization activated cation conductance) to the PRC and hence to synchronization in similar hybrid networks. They also analyse three neuron
bidirectional chains in which two oscillators are connected by a third cell. In one case, the system was reduced to two oscillators by treating the third cell as a follower cell, and in another, symmetry was used to reduce the system. For one type of cell (the oriens-lacunosum-moleculare interneurons), the presence of second order resetting was significant at synchrony. The pulsatile coupling methods accurately predicted synchronization and phase locking in the two and three cell networks, and provided insight into larger networks.

Oprisan SA and Canavier CC. Stability analysis of rings of pulse-coupled oscillators: the effect of phase resetting in the second cycle after the pulse is important at synchrony and for long pulses. Differential Equations and Dynamical Systems, 9:242-259, 2001.