## Abstract

Phase response curves (PRCs) for a single neuron are often used to predict the synchrony of mutually coupled neurons. Previous theoretical work on pulse-coupled oscillators used single-pulse perturbations. We propose an alternate method in which functional PRCs (fPRCs) are generated using a train of pulses applied at a fixed delay after each spike, with the PRC measured when the phasic relationship between the stimulus and the subsequent spike in the neuron has converged. The essential information is the dependence of the recovery time from pulse onset until the next spike as a function of the delay between the previous spike and the onset of the applied pulse. Experimental fPRCs in *Aplysia* pacemaker neurons were different from single-pulse PRCs, principally due to adaptation. In the biological neuron, convergence to the fully adapted recovery interval was slower at some phases than that at others because the change in the effective intrinsic period due to adaptation changes the effective phase resetting in a way that opposes and slows the effects of adaptation. The fPRCs for two isolated adapting model neurons were used to predict the existence and stability of 1:1 phase-locked network activity when the two neurons were coupled. A stability criterion was derived by linearizing a coupled map based on the fPRC and the existence and stability criteria were successfully tested in two-simulated-neuron networks with reciprocal inhibition or excitation. The fPRC is the first PRC-based tool that can account for adaptation in analyzing networks of neural oscillators.

## INTRODUCTION

Phase response curves (PRCs) (Glass and Mackey 1988; Winfree 1980) quantify the change in cycle period of an oscillator due to an external perturbation as a function of the phase within the oscillatory cycle at which the perturbation is delivered. Previous work on neural oscillators has focused on two assumptions that allow PRCs to be used to predict the synchronization of mutually coupled oscillatory neurons. In the first case, the coupling is assumed to be weak (Ermentrout 2002; Jeong and Gutkin 2007; Mancilla et al. 2007; Netoff et al. 2005a; Preyer and Butera 2005; Rinzel and Ermentrout 1998), whereas in the second the coupling is assumed to be pulsatile (Canavier 2005; Maran and Canavier 2008; Netoff et al. 2005b; Oprisan et al. 2004; Pervouchine et al. 2006). The approaches differ in the kind of PRC that is used and in their limitations. For the weak coupling approach, the infinitesimal PRC (iPRC)—the response to an infinitesimally short duration perturbation—is required. The limitation of the weak coupling approach is that the perturbations received within a network must be sufficiently weak that the phase resetting due to the actual inputs received can be calculated as the convolution of this input with the iPRC, that is, the phase resetting should scale linearly with both input magnitude and duration. Theoretically, this method could handle inputs with effects that lasted over several cycles if they were sufficiently weak but could still be measured accurately. For the pulsatile coupling approach, the perturbation used to generate the PRCs should approximate those inputs that will actually be received in the circuit. The limitation of this approach is that the phase resetting due to each perturbation must be complete before the next is received, so for one-to-one lockings it must be complete within a single cycle. Neither approach is well suited to neurons that show slow adaptation in response to strong inputs.

Many neurons have slow currents (such as the M current) that cause spike frequency adaptation on a slower timescale than a single cycle (Benda and Herz 2003; Descalzo et al. 2005). Spike frequency adaptation (SFA) is exhibited by a majority of pyramidal neurons in neocortex and hippocampus (Avoli et al. 1994), which play important roles in advanced cognitive functions. Furthermore, these slow processes can change the synchronization properties of the neurons (Crook et al. 1998; Ermentrout et al. 2001; van Vreeswijk and Hansel 2001). Knowledge of the phase resetting characteristics of adapting neurons can therefore be crucial for the understanding of synchrony of neuron networks and populations.

In this study, we propose an improvement to phase resetting theory that allows prediction of phase locking, even in cases of strong coupling with significant adaptation. Instead of phase resetting curves (PRCs) generated with a single perturbation, we use functional PRCs (fPRCs) generated using a perturbation repetitively triggered at a fixed delay after an onset of each action potential. The intent is to mimic the natural context of a network that exhibits a 1:1 phase-locked periodic rhythm, in which case the inputs also arrive periodically and thus at a fixed delay since the last firing time in the oscillator receiving the input. Others have used a similar protocol to characterize bursting neurons (Pinsker and Kandel 1977), autorhythmic animal hearts (Kunysz et al. 1997; Levy et al. 1972), and the respiratory rhythm (Lewis et al. 1987, 1990). Foss and Milton (2000) characterized the stability of stimulation with a fixed delay greater than the firing period and tested their results in an experimental preparation. Here we focus on fixed delays shorter than the intrinsic period. We characterize adaptation in the spiking neurons in abdominal ganglion of *Aplysia californica* using the fPRC and use fPRCs to successfully predict the synchronization of inhibitory and excitatory networks composed of two adapting model neurons.

First, we address the protocols of obtaining phase response curves (PRCs) and functional phase response curves (fPRCs) and present the experimental measurement of PRCs and fPRCs from spiking neurons of *Aplysia californica*. Then we simulate the Wang and Buzsáki (1996) (WB) model neurons being driven with a fixed delay to generate the fPRC with an adaptation current (Ermentrout 1998). These model neurons exhibit type I excitatory PRCs (Ermentrout 1996) similar in shape to those of our experimental recordings. Type I PRCs consist of only advances for excitation and only delays for inhibition, whereas type II can have both advances and delays for either inhibition or excitation. Next we derive existence and stability criteria for circuits of two coupled adapting neurons using the fPRCs. Finally, we successfully test these criteria in a network consisting of two adapting WB model neurons reciprocally coupled with inhibition or excitation.

## METHODS

### Experimental

##### ELECTROPHYSIOLOGY.

*Aplysia californica* were purchased from University of Miami Aplysia Resource Facility (Miami, FL) and maintained in a tank with filtered artificial sea water (ASW; Instant Ocean, Burlington, NC) at a temperature of about 19°C until used.

The animal was anesthetized by injecting about 50% of the animal's weight of isotonic MgCl_{2} [71.2 g MgCl_{2} in 1 L of 1× ASW solution containing (in mM) 460 NaCl, 10 KCl, 11 CaCl_{2}, 30 MgCl_{2}, 25 MgSO_{4}, and 10 HEPES (pH 7.6)]. The animal was opened and the abdominal ganglion was dissected out and pinned dorsal side up in a dish coated with Sylgard (Dow Corning, Midland, MI) in a saline solution of 30% isotonic MgCl_{2} and 70% 1× ASW. The abdominal ganglion was desheathed with fine scissors and fine forceps.

Desheathed abdominal ganglion was perfused with high-Mg^{2+}–low-Ca^{2+} saline (Nowotny et al. 2003), containing (in mM) 330 NaCl, 10 KCl, 90 MgCl_{2}, 20 MgSO_{4}, 2 CaCl_{2}, and 10 HEPES (pH 7.6). Intracellular recordings from the spiking neurons in the lower left quadrant (Perkel et al. 1964) of abdominal ganglion were recorded using Clampex 8.2 software with an Axoclamp 2B amplifier (Axon Instruments, Foster City, CA) in bridge mode using a microelectrode (5–12 MΩ) filled with 3 M potassium acetate. The membrane potential was amplified and digitized with a Digidata 1322A board (Axon Instruments) at a rate of 10 kHz.

##### MODEL REFERENCE CURRENT INJECTION (MRCI) DYNAMIC CLAMP.

We recorded the membrane potential (*V*) of the spiking neurons and used dynamic clamp (Sharp et al. 1993a,b) to generate the synaptic current input to the spiking neurons on abdominal ganglion of *Aplysia californica*. The excitatory or inhibitory stimulus current was generated from an alpha-shaped synaptic conductance waveform by the MRCI dynamic-clamp system (Butera et al. 2001; Raikov et al. 2004). The equations to generate the alpha shaped conductance and the injected current are given by (1) Here, *I*_{syn} represents the synaptic current injected into the neuron; α(*t*) represents an alpha-shaped waveform with a time constant τ, normalized by the term *e*/τ to a height of one, which provides a conductance of unit amplitude regardless of the value of τ; *g*_{syn} represents the maximum conductance; *i*_{trig} is set as a square pulse with amplitude 1 and duration of 1 ms; and *E*_{syn} is set at 0 mV for excitatory input and −70 mV for inhibitory input.

Use of this differential form to generate an alpha-shaped conductance can accommodate any number of inputs arriving at any time without the need for any additional bookkeeping of input history.

##### PRC MEASUREMENT IN EXPERIMENTS.

Transient phase response curves (PRCs) were obtained by measuring the change in period caused by a transient stimulus. A small amplitude and brief stimulus delivered at a particular time in the firing cycle of a spiking neuron advances or delays the next firing time. The stimulus time *ts*, indicated in Fig. 1 *A*, is the time between spike onset and perturbation onset. The phase φ is equal to the stimulus time normalized by the intrinsic period *P _{i}*, or

*ts*/

*P*. In our experiments, we used spiking neurons in the lower left quadrant (Perkel et al. 1964) of the abdominal ganglion of

_{i}*Aplysia californica*. A stimulus current (as in

*Eq. 1*) generated from an alpha-shaped synaptic conductance waveform by the MRCI dynamic-clamp system was delivered at time

*ts*after the onset of a spiking (Fig. 1

*A*). The period of the perturbed cycle and of the following one are

*P*

_{1}and

*P*

_{2}, respectively. The average of five periods of unperturbed cycles preceding the stimulus is defined as the intrinsic period

*P*. The first-order resetting was obtained as

_{i}*f*

_{1}(φ) = (

*P*

_{1}−

*P*)/

_{i}*P*, and the second-order resetting as

_{i}*f*

_{2}(φ) = (

*P*

_{2}−

*P*)/

_{i}*P*.

_{i}##### FPRC MEASUREMENT IN EXPERIMENTS.

The functional PRC was measured by repeatedly perturbing the spiking neuron of *Aplysia californica* at a fixed-delay time after each action potential until the interspike intervals (ISIs) reach a stable value (Fig. 1*B*), defined as the absolute difference in two successive ISIs being <20 ms. The repeated stimuli were turned off for a period of time for the restoration of the intrinsic period. The fixed-delay time was swept through all phases of the cycle in a random order. The blue part with the next episode of repetitive perturbation with another fixed-delay time is continuous with the green part of the current perturbation.

The temporal change of the ISIs is illustrated by the distance between the dots above the spiking in Fig. 1*B*. Averages of the periods of the last five unperturbed cycles preceding the stimuli and of the last 20 perturbed cycles are defined as the intrinsic period *P _{i}* and the network period

*P*, respectively. The fPRC is calculated as (

_{n}*P*−

_{n}*P*)/

_{i}*P*. The fixed-delay time δ is normalized by the intrinsic period

_{i}*P*to calculate the stimulus phase φ = δ/

_{i}*P*. If second-order resetting is present, then φ +

_{i}*f*

_{2}(φ) = δ/

*P*, so the calculation of φ is not as straightforward. In the presence of adaptation the phase is no longer clearly defined in terms of the intrinsic period and thus the normalization is solely to facilitate comparison with the single-pulse PRCs described earlier.

_{i}### Theoretical

##### MODEL NEURON NETWORKS.

A network consisting of two single-compartment WB model neurons was used. We chose this model for two reasons: *1*) it is a relatively simple but sufficient model to produce quite similarly shaped PRCs to the *Aplysia* neurons that we record from and *2*) it is easily modified to include adaptation. It has been shown that these *Aplysia* spiking neurons exhibit type I PRCs (Preyer and Butera 2005) for a wide range of excitatory stimulus magnitudes. The WB model produced type I PRCs not just of similar shape but of similar changes in shape as the conductance stimulus amplitude is increased, as revealed previously (Preyer 2007). The methods presented herein are general and not dependent on any particular model or on any particular PRC shape.

The model neurons were reciprocally coupled by identical inhibitory or excitatory synapses. The equations for each neuron are given by where the capacitance *C* is 1 μF/cm^{2}; *V* is the cell membrane voltage in mV; *t* is time in unit of ms; *h* is the inactivation of the sodium current; *n* is the activation of the potassium current; *I _{z}* is the adaptation current; and

*m*

_{∞}is the steady-state activation of sodium current. Reversal potentials,

*E*

_{Na},

*E*

_{K}, and

*E*

_{L}, are 55, −90, and −60 mV, respectively. The maximum conductances for sodium

*g*

_{Na}, potassium

*g*

_{K}, and leak

*g*

_{L}are 35, 9, 0.1 mS/cm

^{2}, respectively. The scale factor γ is set to 5.

The adaptation current *I*_{z}, based on the formalism described by Ermentrout (1998), is given by where *a* = 0.005 ms^{−1}, *b* = −5 mV, *k* = 2.0 mV, and *g*_{z} = 1 mS/cm^{2}. Because of the steep threshold for activation at −5 mV, this type of adaptation is induced only during spikes, similar to a calcium-dependent afterhyperpolarizing current (*I*_{AHP}) but not to a muscarinic current (*I*_{M}), which allows adaptation to be produced at subthreshold voltages (Ermentrout 1998; Prescott et al. 2006).

The synaptic current *I*_{syn} is given by where *g*_{syn} is the maximum synaptic conductance; *E*_{syn} is the reversal potential, which is set to −75 mV for inhibitory perturbation and 0 mV for excitatory perturbations; *s* is the gating variable; *V*_{pre} is the voltage of the presynaptic cell; *k*_{s} = 6.25 ms^{−1} is the rate constant of the synaptic activation; τ_{syn} is the synaptic decay time constant; and *I*_{app} represents the applied current to the model cell to control the spiking frequency. In our model networks, the two cells have identical values of *I*_{app}. The differential equations were integrated using a driver program written in C that called the Fortran subroutine radau5 (Hairer and Wanner 1991), using a variable step implicit fifth-order Runge–Kutte method on a Beowulf cluster running CentOS.

##### SIMULATED PRCS AND FPRCS.

The protocols of obtaining PRCs for the adapting model neuron (Maran and Canavier 2008) were similar to the experimental method. The fPRCs were calculated at delays in increments of one hundredth of the period. Fifty stimuli at a fixed delay were applied at each delay, then the period was allowed to stabilize for 2 s with no perturbations before another stimulus at a fixed delay was applied.

## RESULTS

In this section, we present experimental fPRCs for both inhibitory and excitatory perturbations and obtain similar fPRCs using the modified WB model with adaptation. For strong perturbations, the assumption that the resetting adds linearly is violated, so the sum of the first-, second-, and higher-order PRCs does not converge to fPRCs. We derive the existence and stability criteria for a two-neuron network based on information acquired using the fPRC protocol and use these criteria to predict phase-locking behavior in networks consisting of two adapting model neurons coupled with either excitatory or inhibitory synapses. We show that in the presence of adaptation, the predictions made by fPRCs are much better than those obtained using existing methods based on standard PRCs.

### Experimental

When the fPRC protocol was applied as in Fig. 1*B*, the most dramatic changes in ISI length occurred when the train of fixed-delay stimuli was turned on or off. The ISI between the first and second red dots (after the pulse train was turned on) is much smaller than the intrinsic free-running period, as measured between the last blue dot and the first red dot; this change is primarily due to the effects of the phase resetting that produces an advance and therefore a faster frequency. On the other hand, the ISI between the last red dot and the first green dot (when the pulse train was turned off) is longer than the intrinsic period, presumably due to the increase in activation of an adaptive current caused by the increase in frequency. The corresponding time sequence of ISIs in the process of measuring this particular fixed delay (a single point in the fPRC) is illustrated in Fig. 1*C*. The ISI before the repetitive stimuli were turned on (blue dots) was almost steady at 420 ms; the ISI dropped to 220 ms (first red dot) once the pulse train was switched on at a fixed-delay time of 160 ms (indicated by the cyan star) and gradually reached a steady value around 280 ms (red dots); the ISI increased to 500 ms when the repetitive perturbations were turned off (green dots) and gradually reached the intrinsic period. The fixed-delay time was swept through values from zero to that of a single intrinsic period in a random order to minimize any systematic error in the experiment.

The time sequences of the ISIs reordered by the fixed-delay time for inhibitory and excitatory perturbations are illustrated in Fig. 2, *A* and *B*, respectively. The blue, green, and red dots are color-coded the same as in Fig. 1*C*. For excitation, the dramatic changes in the ISI, when the train of fixed-delay perturbations was turned on and off, are more obvious at smaller phases than at larger phases, whereas the converse is true for inhibitory perturbations. The duration of the transient period prior to converging to a stable ISI was longer at earlier phases for excitation and later phases for inhibition. This was true both for the transient at the onset of pulse stimulation and for the transient where the stimulation was turned off.

Before measuring the fPRC with a repetitive perturbation, we also measured a single-pulse PRC with the same experimental conditions. For a larger magnitude of perturbation [Fig. 3, *A1* (inhibition) and *B1* (excitation)] fPRCs (stars) deviated dramatically from the first-order PRCs (dots). For a smaller magnitude of perturbation [Fig. 3, *A2* (inhibition) and *B2* (excitation)] fPRCs are close to the first-order PRCs because both the second-order resetting and any adaptation are negligible. For excitatory cases, we found the results of Fig. 3*B* to be representative of all neurons measured. Specifically, in over 20 neurons measured, the first-order single-pulse PRC (*f*_{1}) exhibited advances at all phases; the second-order single-pulse PRC (*f*_{2}) exhibited delays at all phases; the fPRC was similar to the summation of *f*_{1} and *f*_{2}, (*f*_{1} + *f*_{2}), for sufficiently weak inputs; and the fPRC was distinct from *f*_{1} + *f*_{2} for strong inputs.

### Computational and theoretical

Motivated by the experimental demonstration of adaptation and its effect on the functional PRC, we then attempted to reproduce similar results in a model, to analyze the stability of the protocol, to use the data to predict the synchronization behavior of a circuit comprised of two adapting neurons, and to test these predictions.

The simulated measurements of fPRCs for inhibitory and excitatory perturbations are shown in Fig. 4, *A* and *B*, respectively, which capture a trend similar to that of the experimental observations as shown in Fig. 2, *A* and *B*. The simulations also capture the longer convergence times for excitation for shorter delays and for inhibition for longer delays. If the currents are adapting, the protocol of measuring fPRCs allows the currents to adapt completely. The net result is a periodic steady state in which the neuron receives periodic inputs and is engaged in a 1:1 phase-locked periodic rhythm. Thus the fPRC accounts for the degree of adaptation that normally occurs during a periodic rhythm in a network context.

The fPRC obtained in Fig. 4, *A* and *B* is plotted as the dot-dashed curve in Fig. 4, *C* and *D*. Note that the fPRC does not exist at late phases due to stability problems described in the following text. Figure 4*C* shows that the first-order resetting (*f*_{1}, solid line) for inhibition consists of delays, whereas the second-order (*f*_{2}, dotted curve) and the third-order (*f*_{3}, long dashes) phase responses consist of advances. The sum of the three resetting components (*f*_{1} + *f*_{2} + *f*_{3}, short dashes) approaches the functional PRC and might approach it more closely if higher-order terms were included, but would require a fine resolution and the analysis of many cycles. On the other hand, Fig. 4*D* shows that the first-order resetting for excitation consists of advances, whereas the second- and third-order resettings consist of delays. Both the first-order resetting and the fPRC are limited by causality. In other words, the amount of advance produced is limited by the time remaining until the next spike would have occurred anyway. The causality limit falls along the line *P _{i}*(1 −

*f*

_{1}). For both the simulated and experimental cases of excitatory input (Figs. 4

*D*and 3

*B1*), the fPRC differs from the first-order PRC at early phase, but converges on it at later phases when the causality limit is reached.

For weak perturbations, the assumption that the resetting in subsequent cycles adds linearly is valid, so the sum of the first-, second-, and higher-order PRCs converges to fPRCs. The fPRCs are very close to the first-order PRCs because second- (and higher-) order resetting is small (Fig. 3, *A2* and *B2*). The fPRCs are different from the first-order PRCs when the perturbation is strong in both the excitation and inhibition cases, as shown in Figs. 3, *A1* and *B1* and 4, *C* and *D*. For strong perturbations, first-order PRCs advance (delay) phase, whereas second- and third-order PRCs delay (advance) phase for excitation (inhibition). However, the second- and third-order resettings are not necessarily manifested under repetitive stimulation for strong perturbations because the neuron does not return to the original limit cycle between perturbations. In other words, the assumption of pulsatile coupling is violated. Therefore including the second- and third-order resettings can produce a worse match to the functional PRC at most phases than ignoring it (Fig. 4*D*). For less severe perturbations (experimental Fig. 3), adding second-order resetting can slightly improve the correspondence with the fPRCs.

Adaptation can be interpreted as changing the intrinsic period of the oscillator due to that activation of the adaptation current. The adaptation current used in the model summates temporally only during a spike and decays in between spikes. Therefore the greater the change in number of spikes per unit time, the greater the change in the level of activation of the adaptation current. An estimate of the amount of adaptation that has occurred can be obtained by observing the length of the ISI immediately after the fixed-delay protocol is terminated. We refer to this estimate as the modified intrinsic period (*P _{m}* in Fig. 4,

*A*and

*B*) and to the steady period achieved during the fixed-delay protocol as the network period (

*P*in Fig. 4,

_{n}*A*and

*B*). The actual phase at which an input is received cannot be determined by using the delay normalized by

*P*if

_{i}*P*differs significantly from

_{m}*P*; instead one should normalize by

_{i}*P*. Furthermore, the single-pulse PRC measured in the absence of adaptation may be quite different from that measured in its presence, if that were possible. Therefore single-pulse PRCs are not strictly applicable in the presence of significant adaption, as we will show.

_{m}Note that the inhibitory PRC (Fig. 4*C*, solid line) is larger at late phases, whereas the excitatory PRC (Fig. 4*D*, solid line) is larger at early phases. Slow convergence to a stable phase relationship results because phase resetting changes the period, which causes adaptation that changes the effective intrinsic period (*P _{m}*) in the opposite direction. For inhibition, resetting increases the period, which decreases the activation of the adaptation current, which decreases the effective intrinsic period. Decreasing the effective intrinsic period increases the effective phase, which generally further increases the observed network period, slowing and opposing the effect of adaptation until the modified intrinsic period stops changing (adaptation is complete). A complementary process occurs for excitatory inputs. The interplay between phase resetting and adaptation is complex. Each pulse has an effect on the timing of the next spike due to resetting the phase of the ongoing oscillation, but it also affects the timing of the next spike by altering the amount of time available for the adaptation current to decay—and these two effects usually oppose each other. It is not easy to predict whether the opposing effects will cancel each other out or whether one will dominate.

### Existence and stability criteria for two-neuron network

The theoretical method for predicting the existence and stability of phase-locked modes using single-pulse PRCs for networks of two bursting neurons has been previously developed (Oprisan et al. 2004; Sieling et al. 2009). We describe it here to show how the new methods based on the fPRC are an improvement. The assumed firing pattern in a two-cell network with a 1:1 locking is given in Fig. 5 *A1*. Neuron 1 (2) sends input to neuron 2 (1) when it fires. The recovery interval (*tr*) is defined as the time elapsed between when a neuron receives an input and when it next fires. The stimulus interval (*ts*) in a network is the time between when a neuron fires and when it receives an input. By definition of 1:1 phase locking, *tr*_{1}[*n* − 1] = *ts*_{2}[*n*] and *tr*_{2}[*n*] *= ts*_{1}[*n*].

The assumptions of the single-pulse PRC method can be summarized as follows: *1*) each neuron in the network can be represented as a limit cycle oscillator, *2*) the trajectory of each neuron returns to its limit cycle between inputs (this is the assumption of pulsatile coupling), and *3*) the input received by each neuron has the same effect in the closed-loop circuit as in the open-loop circuit used to generate the phase response curves. The key idea is that under the assumption of pulsatile coupling, the first-order PRC (*f*_{1}) and the second-order PRC (*f*_{2}) can be used to calculate the values of the intervals as follows: *ts* = *P _{i}*[φ +

*f*

_{2}(φ)],

*tr*=

*P*[1 − φ +

_{i}*f*

_{1}(φ)], such that

*tr*is a function of

*ts*,

*tr*=

*g*(

*ts*). When phase-locking occurs,

*tr*

_{1}=

*ts*

_{2}and

*tr*

_{2}=

*ts*

_{1}, so the intersection of the

*tr*

_{1}=

*g*(

*ts*

_{1}) and

*ts*

_{2}=

*g*

^{−}

^{1}(

*tr*

_{2}) curves yields the fixed points. The stability is calculated using the slopes of the first- and second-order phase resetting curves as in Oprisan et al. (2004).

In contrast, the fPRC requires only the third assumption just cited, that the fPRC is generated with an input similar to the one it will receive in the network. The key idea of the functional PRC is to characterize the response of each oscillator to inputs arriving periodically with a fixed delay to characterize the response of the steady phase-locked mode when the delay times within the network are perturbed away from the steady value. The functional PRC protocol (Fig. 5*A2*) allows the tabulation of the steady-state recovery interval (*tr*) in a given neuron as a function of the delay (δ) between when the neuron fires and when it next receives an input, such that *tr = g*(δ). Figure 5*A3* illustrates the assumption that in a closed-loop circuit the recovery intervals for each neuron can be calculated as *g*(δ) by assuming that at steady state the stimulus interval is the fixed delay δ. Thus at steady state in a periodic locking δ_{2} = *g*_{1}(δ_{1}) and δ_{1} = *g*_{2}[*g*_{1}(δ_{1})]. The modes that satisfy the existence criteria can be located graphically (Fig. 5*B*) at the intersection of *tr*_{2}[∞] *= g*_{1}(*ts*_{1}[∞]) and *ts*_{2}[∞] = *g _{2}^{−1}*(

*tr*

_{2}[∞]), where the minus 1 indicates the inverse function and

*tr*

_{1}[∞]

*= ts*

_{2}[∞]

*=*δ

_{2}and

*tr*

_{2}[∞]

*= ts*

_{1}[∞]

*=*δ

_{1}. Although the delays are not fixed in the closed-loop circuit, the stability of the locking can be approximated by assuming a perturbation of two neurons that are actually coupled by fixed delays. Linearizing the map δ

_{1}[

*n*]

*= g*

_{2}{

*g*

_{1}(δ

_{1}[

*n*− 1])} provides a stability criterion that −1 <

*g**(

_{2}*ts**)

_{2}*g**(

_{1}*ts**) < 1, where the asterisk indicates the value at the fixed point and the prime indicates the slope at the fixed point. Note that all intervals are defined in terms of time rather than phase. The modes that satisfy the existence criteria are located graphically by finding the intersection of

_{1}*tr*

_{1}=

*g*

_{1}(

*ts*

_{1}) and

*ts*

_{2}=

*g*(

_{2}^{−1}*tr*

_{2}), where

*tr*

_{1}=

*ts*

_{2}and

*tr*

_{2}=

*ts*

_{1}

_{.}Then the stability criterion is applied to determine whether the mode is stable, that is, whether the inevitable perturbations away from the presumed fixed-delay times in the network will grow or shrink. Figure 5

*B*shows the intersection of the

*tr*

_{1}=

*g*

_{1}(

*ts*

_{1}) and

*ts*

_{2}=

*g*(

_{2}^{−1}*tr*

_{2}) curves for the WB model with adaptation and mutual excitation. The reciprocal fixed-delay paradigm described here does not completely characterize the circuit because it neglects the feedback interaction between oscillators, although the results in Fig. 5,

*C*and

*D*described in the following text confirm that this is a reasonable approximation. The idea is that if each neuron can lock to an input given at a fixed delay, when they are reciprocally coupled with the same effective (but not fixed) delays, the resultant locking should be stable.

The intersection method given earlier does not find exact synchrony because it presumes an alternating firing order. A separate check is run to determine whether the recovery intervals at a phase of zero are measurable, indicating that the fixed-delay protocol produces a stable locking and that they are equal. Period two modes that satisfy the mapping, *tr*_{2}[*n*] *= g*_{2}〈*g*_{1}{*g*_{2}[*g*_{1}(*tr*_{2}[*n* − 2])]}〉, could theoretically be found using the zeros of the firing time map Δ = *g*_{2}〈*g*_{1}{*g*_{2}[*g*_{1}(*tr _{2}*[

*n*− 2])]}〉 −

*tr*

_{2}[

*n*] (but see next section).

For perturbations from exact synchrony, the criterion −1 < *g* _{2}*(

*ts**)

_{2}*g**(

_{1}*ts**) < 1 should be applied at

_{1}*ts**= 0 and

_{1}*ts**=

_{2}*g*(0). In other words, there are two possible perturbations from synchrony, one in which a given neuron slightly leads it partner and one in which it lags. When it leads, it will receive an input at a phase slightly greater than zero (to the right of zero). The network period is given by

*g*(0); therefore when the neuron lags its partner it will receive an input at a time just less than

*g*(0) [to the left of

*g*(0)]. In some cases a stable locking could not be obtained for the fPRC protocol at a phase of

*g*(0). In that case, there is no way to measure this slope from the left; thus in practice the criterion cannot be applied. However, the absence of a stable locking implies that synchrony in the two-neuron network would be vulnerable to perturbations and probably unstable. The inability to rigorously apply the criterion is of concern in only a narrow region on the border of stability in which the attraction of the stable locking at zero delay compensates for the unstable locking at a delay of one network period.

The choice of a type I PRC for the simulated example was purely phenomenological to match the experimental neuron. Neither the conventional PRC nor the functional PRC methods depend on the shape of the PRC for their implementation, although the shape of the PRC or fPRC determines which modes can exist as well as their stability.

### Testing the predictions with networks of adapting model neurons

To determine any advantage of using the fixed-delay protocol, we compared the predictions made with multiple-pulse fPRCs and single-pulse PRCs. For a network consisting of two adapting WB model neurons reciprocally coupled by inhibitory and excitatory synapses, the predictions based on the fPRC protocol (green circles) gives a better agreement with the observations (red crosses) than the prediction by standard PRCs (blue circles) for both inhibitory (Fig. 5*C*) and excitatory (Fig. 5*D*) coupling. As indicated in Fig. 5*C*, in an inhibitory network, antiphase synchrony was perfectly predicted by the fPRC method (green circles) for all the parameters tested. In an excitatory network (Fig. 5*D*), the nearly synchronous mode, the nearly antiphase mode, and the perfect antiphase mode are correctly predicted by the fPRC method. The standard PRC method correctly predicted for only very limited ranges of small conductance for both inhibition and excitation networks.

With mutual inhibition (Fig. 5*C*) the antiphase mode is always observed and predicted for the range of values of synaptic conductance. In addition, in-phase synchrony (not shown) is observed in the range of conductance values 0.01 to 0.06 mS/cm^{2}, but is predicted from only 0.01 to 0.04 mS/cm^{2} with multiple-pulse fPRCs. This is because at conductance values of 0.05 to 0.06 mS/cm^{2} no stable locking fixed-delay protocol existed at a delay of one network period corresponding to the network period at a delay of zero (see above-cited methods). Therefore synchrony was predicted to be unstable because there was no way to calculate the required slope for a perturbation from synchrony. At conductance values from 0.07 to 0.18 mS/cm^{2}, near-leapfrog modes in which the order of firing changes every cycle (Maran and Canavier 2008; Oh and Matveev 2008) were observed to be bistable (Skinner et al. 2005) with the antiphase mode. These leapfrog modes are not shown because they could not be predicted as the fixed points of the map described in the preceding session. The stable lockings under the fixed-delay protocol at the appropriate delays required to measure the *g* functions used in the maps did not exist. On the other hand, the method based on single-pulse PRCs performed much more poorly, predicting stable synchrony at conductance values from 0.01 to 0.16 mS/cm^{2}.

Figure 5*D* compares the predictions made by fPRC (green circles) and standard PRC (blue circles) for excitation networks. The predictions generated using standard PRC protocol are qualitatively incorrect beyond about 0.08 mS/cm^{2}. At smaller values of the synaptic conductance, both methods predict a nearly synchronous mode in which the firing intervals between the two neurons alternate between short and long. Above this value of conductance, both the observed modes and the modes predicted using the fPRC protocol gradually change from near synchrony to near antiphase and finally to perfect antiphase at *g*_{syn} = 0.16 mS/cm^{2} in which all firing intervals are equal (see also Fig. 5*B*). The intrinsic period of the uncoupled oscillators is 60 ms, so the phase resetting greatly shortens the period. For example, the network period is 21.75 ms at *g*_{syn} = 0.16 mS/cm^{2}. As the conductance is increased beyond 0.3 mS/cm^{2}, the network period is reduced to less than a sixth of its uncoupled value and the multiple-pulse protocol no longer accurately predicts the antiphase mode. This is because the gaps in the curves similar to those plotted in Fig. 5*B* prevent an intersection. Therefore the individual fixed-delay lockings that constitute the observed antiphase mode are unstable and are stabilized by the coupling in a manner that we cannot at this time explain or predict.

## DISCUSSION

The effect of adaptation on synchronization properties has inspired many theoretical and computational studies (Crook et al. 1998; Fuhrmann et al. 2002; Latham et al. 2000a,b; van Vreeswijk and Hansel 2001). Herein, we have proposed a novel method for studying the effects of adaptation on synchronization properties that does not require the assumption of pulsatile coupling or of weak coupling and can be applied directly to experimental data. Here we discuss issues related to our method, such as the history of the use of a fixed-delay protocol, the relationship between adaptation and PRC types, and the limitations of our method. We conclude with a discussion of significance and possible extensions of the fPRC method to larger networks.

### History of fixed-delay phase resetting experiments

The phase resetting effects of a single stimulus have been studied extensively (Glass and Mackey 1988; Winfree 1980), but the effects of repeated stimulation at a fixed delay following a marker event such as a spike have been studied much less. An early study (Pinsker and Kandel 1977) used the term “contingent” to describe pulses applied at a fixed delay after a marker event. Stimulation of a presynaptic cell was contingent on the beginning of a burst in the regularly bursting of left upper quadrant bursters in the abdominal ganglion of *Aplysia*. For some early and late phases, they found an adaptation effect for multiple-pulse PRCs. Early fixed-delay studies using repetitive burst perturbations were performed in the hearts of anesthetized dogs (Levy et al. 1972; Michaels et al. 1984). The cardiac pacemaker was driven to a synchronized frequency with the stimuli. The study of fixed-delay stimulation on spontaneous beating of chick heart cell aggregates (Kunysz et al. 1997), an experimental model of reentrant tachycardia, showed that the complex patterns of bursting behaviors arose from the interaction of phase resetting (thus the excitability) and overdrive suppression, which altered with the amplitude of the stimuli. The authors confirmed the finding with the theoretical studies on a nonlinear model and an ionic model. A study considering respiratory rhythm applications assumed a radial isochronic clock model in which a fixed delay implied that the input was received at a fixed phase, but at variable amplitude such that complicated dynamics ensued (Lewis et al. 1987). The trajectory was not constrained to return to the original limit cycle during the stimulation protocol.

More recently, single-pulse phase resetting curves (Foss and Milton 2000; Foss et al. 1997) were used to determine the existence and stability criteria for phase locking to a repetitive input at a fixed delay that could be much longer than the oscillator period; inputs at fixed delays shorter than the oscillator period were assumed to occur at a fixed phase and to be stable because the phase was not assumed to be able to change during the fixed-delay interval (second-order resetting was ignored, as was adaptation). In no case were the results of the fixed-delay protocol used to predict 1:1 phase-locking network activity.

### Relationship of adaptation and PRC types

Previous studies, such as that of Mancilla et al. (2007) and Tsubo et al. (2007), focused on changing the frequency with an applied current, allowing adaptation to stabilize, then measuring a single-pulse PRC. This method is more suited to weak coupling in which the network frequency is assumed not to be very different from the intrinsic frequency of the component neurons. There are various biophysical mechanisms that underlie spike frequency adaptation. The most commonly studied are the calcium-dependent potassium current (*I*_{AHP}) and the slow voltage-dependent potassium current (*I*_{M}) and, less frequently, sodium-activated potassium currents. As mentioned previously, there are two basic PRC types, type I and type II; they have different shapes and are associated with different synchronization properties (Ermentrout 1996; Rinzel and Ermentrout 1998). Increasing *I*_{M} can switch the PRC of hippocampal pyramidal cells (Prescott et al. 2008) from type I to type II. This result was predicted (Prescott and Sejnowski 2008; Prescott et al. 2006) by dynamic analysis of a Morris–Lecar model with both adaptation currents of *I*_{AHP} and *I*_{M.} Conversely, a decrease in the amount of *I*_{M} is predicted to change the shape of the PRC from type II to type I. This switch was demonstrated in cortical pyramidal neurons (Stiefel et al. 2008, 2009) by using a cholinergic neuromodulator to down-regulate *I*_{M}. Modeling studies (Ermentrout et al. 2001; Jeong and Gutkin 2007) suggest that in type I neurons *I*_{AHP} flattens the PRC at the beginning of the firing cycle. Our methods do not require us to determine how the shape of the single-pulse PRC changes as a function of adaptation or of frequency, but rather depend on the simple and direct measurement of the network period on the time at which an input is received.

In our experiments, we found that the *Aplysia* neurons that we studied exhibit type I standard PRC under all the values of synaptic conductance that were tested (see Fig. 2 of Preyer and Butera 2005). In the simulations, we were able to match the adaptation observed in *Aplysia* neurons by using an *I*_{AHP}-like current that temporally summates only during spikes and decays in between. The theoretical method based on fPRCs successfully predicted the synchronization of excitatory and inhibitory networks consisting of two type I neurons with an *I*_{AHP}-like adaptation current. Our method, however, is general and completely independent of the underlying PRC types.

### Limitations

A derivation for stability of a fixed-delay locking in the presence of adaptation was not obtained and we can obtain the recovery interval only as a function of the delay for stable lockings. It seems likely that if both neurons in the network can stably entrain in isolation then the network consisting of those intervals will be stable, but it is also conceivable that the feedback within the network could stabilize some delay values that are not stable in isolation. By analogy, Dror et al. (1999) showed a stable mutual locking of neurons that could not be periodically driven at the same phase they receive input in the mutual locking.

### Summary

In our experimental and numerical studies, functional PRCs were generated by delivering a repetitive-pulse perturbation at a fixed-delay time triggered by spike onset at delays ranging from zero to the intrinsic period. Our experimental results on spiking neurons of *Aplysia* have shown that the fPRCs are distinct from the standard PRC for relatively strong perturbations that induce spike frequency adaptation. The simulated fPRCs for an adapting WB model neuron were in good agreement with experimentally generated fPRCs for strong inhibitory and excitatory perturbations. The predictions of 1:1 phase-locking network activity generated using fPRCs were in good agreement with the observations of phase locking in networks consisting of two adapting neurons coupled with excitatory or inhibitory synapses.

Previously the application of PRC methods has been limited to limit cycle oscillators. However, some components of central pattern generators, such as postinhibitory rebound neurons, do not oscillate when uncoupled from the network. Others may oscillate chaotically (Elson et al. 1999). If these neurons entrain to a periodic input, it is still possible to describe their recovery times as a function of the stimulus time and analyze the network using the methods developed in this study. Therefore the significance is that this method can correctly account for the effect of adaptation on synchronization and has the potential to be applied even to neurons that are not endogenous oscillators but participate in a periodic rhythm.

A limitation of the fPRC method is that a given neuron is assumed to receive only a single input per cycle; thus it cannot be easily generalized to account for a triphasic rhythm in three neurons, for example. However, as currently formulated, the fPRC method can easily be extended to the analysis of global synchrony in an N neuron network using a method analogous to that shown for the conventional PRC (Achuthan and Canavier 2009).

The fPRC method can also be extended to any 1:1 phase locking between two clusters of neurons, again by analogy with published PRC methods (Achuthan and Canavier 2009). Locking between clusters of similar or disparate populations may be widespread in the nervous system. It has been suggested that a synchronous population of cortical neurons locks with a population of thalamic neurons in absence seizures (Perez Velazquez et al. 2007). Similarly, populations of inspiratory neurons in the pre-Botzinger complex and expiratory neurons in the parafacial respiratory group (Janczewski and Feldman 2006; Joseph and Butera 2005; Mellen et al. 2003) have been proposed to phase lock during respiratory rhythmogenesis. Two clusters of entorhinal stellate cells (Pervouchine et al. 2007), with each cluster firing at theta frequency in antiphase with the other cluster, have been invoked to explain a beta peak seen by Cunningham et al. (2004) in slices of entorhinal cortex.

## GRANTS

This study was supported by National Institute of Neurological Disorders and Stroke Grant NS-54281 to C. Canavier (Principal Investigator) under the Collaborative Research in Computational Neuroscience Program and National Institutes of Health Division of Research Resources Grant RR-020115, subcontracted to R. Butera.

## Acknowledgments

We thank the reviewers for suggestions, which greatly increased the clarity of the manuscript.

- Copyright © 2009 the American Physiological Society

## REFERENCES

- Achuthan 2009.↵
- Avoli 1994.↵
- Benda 2003.↵
- Butera 2001.↵
- Canavier 2005.↵
- Crook 1998.↵
- Cunningham 2004.↵
- Descalzo 2005.↵
- Dror 1999.↵
- Elson 1999.↵
- Ermentrout 1996.↵
- Ermentrout 1998.↵
- Ermentrout 2002.↵
- Ermentrout 2001.↵
- Foss 2000.↵
- Foss 1997.↵
- Fuhrmann 2002.↵
- Glass 1988.↵
- Hairer 2006.↵
- Janczewski 2006.↵
- Jeong 2007.↵
- Joseph 2005.↵
- Kunysz 1997.↵
- Latham 2000a.↵
- Latham 2000b.↵
- Levy 1972.↵
- Lewis 1987.↵
- Lewis 1990.↵
- Mancilla 2007.↵
- Maran 2008.↵
- Mellen 2003.↵
- Netoff 2005a.↵
- Netoff 2005b.↵
- Nowotny 2003.↵
- Oh 2008.↵
- Oprisan 2004.↵
- Perez Velazquez 2007.↵
- Perkel 1964.↵
- Pervouchine 2006.↵
- Pervouchine 2007.↵
- Pinsker 1977.↵
- Prescott 2006.↵
- Prescott 2008.↵
- Prescott 2008.↵
- Preyer 2007.↵
- Preyer 2005.↵
- Raikov 2004.↵
- Rinzel 1998.↵
- Sharp 1993a.↵
- Sharp 1993b.↵
- Sieling 2009.↵
- Skinner 2005.↵
- Stiefel 2008.↵
- Stiefel 2009.↵
- Tsubo 2007.↵
- van Vreeswijk 2001.↵
- Wang 1996.↵
- Winfree 1980.↵