# Phase-Response Curves Give the Responses of Neurons to Transient Inputs

## Abstract

Neuronal firing is determined largely by incoming barrages of excitatory postsynaptic potentials (EPSPs), each of which produce a transient increase in firing probability. To measure the effects of weak transient inputs on firing probability of cortical neurons, we compute phase-response curves (PRCs). PRCs, whose shape can be related to the dynamics of spike generation, document the changes in timing of spikes caused by an EPSP in a repetitively firing neuron as a function of when it arrives in the interspike interval (ISI). The PRC can be exactly related to the poststimulus time histogram (PSTH) so that knowledge of one uniquely determines the other. Typically, PRCs have zero values at the start and end of the ISI, where EPSPs have minimal effects and a peak in the middle. Where the peak occurs depends in part on the firing properties of neurons. The PRC can have regions of positivity and negativity corresponding respectively to speeding up and slowing down the time of the next spike. A simple canonical model for spike generation is introduced that shows how both the background firing rate and the degree of postspike afterhyperpolarization contribute to the shape of the PRC and thus to the PSTH. PRCs in strongly adapting neurons are highly skewed to the right (indicating a higher change in probability when the EPSPs appear late in the ISI) and can have negative regions (indicating a decrease in firing probability) early in the ISI. The PRC becomes more skewed to the right as the firing rate decreases. Thus at low firing rates, the spikes are triggered preferentially by inputs that occur only during a small time interval late in the ISI. This implies that the neuron is more of a coincidence detector at low firing frequencies and more of an integrator at high frequencies. The steady-state theory is shown to also hold for slowly varying inputs.

## INTRODUCTION

Information about a stimulus can be encoded by changes in the neuronal firing rate and/or changes in the timing of individual action potentials (see Mazurek and Shadlen 2002; Shadlen and Newsome 1998). A growing body of experimental evidence suggests that both rate and spike-time response modes coexist (Bialek and Rieke 1992; de-Ruyter-van-Steveninck 1997; Markam 1997; Panzeri et al. 2001; Prut et al. 1998; Reich et al. 1998; Reinagel and Reid 2000; Riehle et al. 1997; Thorpe 1996). In fact neurons can generate precisely timed spikes (e.g., Calvin and Stevens 1968; Fellous et al. 2000; Mainen and Sejnowski 1995; Nowak et al. 1997) that can depend on intrinsic biophysical properties of neurons as well as on the recent input history (Azouz and Gray 1999, 2000). Although the timescales of the rate and spike-time modes may differ, they are not independent of each other.

In general, individual synaptic inputs are relatively weak (Markram et al. 1997; Reyes and Sakmann 1999) so that many are needed to evoke repetitive firing (Oviedo and Reyes 2002). Nevertheless, single inputs produce a small but measurable effect on the timing of subsequent spikes. Reyes and Fetz (1993a,b) experimentally examined the effects of small transient inputs on the timing of spikes in a cortical neuron that was driven to fire repetitively by a constant depolarizing current. The excitatory postsynaptic potential (EPSP)–like inputs produced a clear shift in the time of the next spike that depended on when the input occurred relative to the previous spike. That is, the observed spike-time shifts strongly depended on the timing of the weak transient inputs. The curve of spike-timing shifts (called the phase-response curve [PRC]) provides a useful tool for the understanding of how small inputs affect spike times. An important implication is that the overall (mean) firing rate of the neuron and the timing of the action potentials in response to incoming synaptic inputs depends on both the average rate and the temporal correlation of the inputs.

We consider a periodically firing neuron as a nonlinear oscillator and then use methods from the theory of dynamical systems to describe the effects of small transient inputs on the timing of the spikes. Rinzel and Ermentrout (1998) showed that the onset of repetitive firing in neurons generally occurs through one of two mechanisms: *1*) a saddle-node bifurcation on a circle (see Fig. 1; called type I excitability) or *2*) a Hopf bifurcation (called type II excitability). Tateno et al. (2004) recently suggested that cortical pyramidal cells fall into the first class and inhibitory interneurons fall into the second class. We will focus exclusively on models of class I excitability because we are interested primarily in the responses of cortical pyramidal cells. Furthermore, the cortical pyramidal PRCs measured by Reyes and Fetz (1993) are nonnegative; neurons that begin repetitive firing at a Hopf bifurcation (type II) always have a substantial negative region in their PRCs (Brown et al. 2004a). In a previous paper (Ermentrout et al. 2001; Fig. 8), we fit the data from Reyes and Fetz (1993) to a conductance-based model with a calcium-dependent potassium current. Here we use a nonlinear two-dimensional integrate-and-fire model (θ-neuron) to fit the data. We then use this simplified model to illustrate the role that the firing rate and the degree of adaptation have on the shape of the PRC. Conductance-based models of cortical neurons are complex and have many compartments and voltage-gated conductances. However, because the onset of repetitive firing stereotypically falls into one of two classes, the present results can be regarded as being very general, as illustrated by the simple model of spike generation.

We first describe the simple spiking model that extends the θ-neuron to incorporate a slow negative feedback term. The θ-neuron is formally equivalent to every model with type I excitability near the onset of repetitive firing (Ermentrout and Kopell 1986; Hoppensteadt and Izhikevich 1997). It is also equivalent to the quadratic integrate-and-fire (QIF) model, which has been used by numerous authors as a model for cortical cells (Latham et al. 2000; Pfeuty et al. 2003). Izhikevich (2004) recently extended the QIF model to include a negative feedback term and showed that it reproduces the voltage traces of many types of neurons. A change of variables converts the QIF model with adaptation into the θ-neuron with adaptation; this will be our model of choice in the remainder of the paper.

We then show the equivalence between the PRC of a neuron to the poststimulus time histogram (PSTH). Because the PSTH shows a strong relationship to the background firing rate of neurons (Brown et al. 2004b; Herrmann and Gerstner 2001), this motivates our study of the PRC and how spike-frequency adaptation (SFA) and background firing rate affect it. We use the extended θ-neuron to predict how the PRCs are affected by neurons’ intrinsic firing properties. In addition, to directly examine the relationship between the effects of EPSPs on the cells’ overall firing rate and the timing of individual spikes, we document the changes in PRCs with the cells’ baseline firing rates. We show that the PRC formalism can work under noisy conditions. We further show that, although PRC theory is strictly valid only under the assumption of stationarity of the input (i.e., the neuron is firing at a constant frequency), the methods nevertheless work reasonably well with time-varying inputs. Strictly speaking the mathematical reduction is valid at low firing rates; to validate our use of the theta model we compare results with those obtained in a conductance-based model of pyramidal neurons.

## METHODS

Many types of neurons fire repetitively when injected with a constant current of sufficient amplitude. Hodgkin (1948) divided excitable membranes into three classes based on their responses to brief near-threshold stimuli and constant suprathreshold stimuli. Many types of cortical neurons (e.g., regular spiking units and fast-spiking units) fall into Hodgkin’s type I excitability: *1*) they fire trains of action potentials in response to constant current injections; *2*) they fire at arbitrarily low frequencies; and *3*) the relationship between firing rate and the magnitude of the injected current (*F*/*I* curve) can be fit with a square root (for the instantaneous firing rate) or a linear (for steady-state firing rate) function. Mathematical analysis of detailed biophysical models of these neurons shows that as current is increased, the stable resting state of the neuron collides with an unstable state (called a *saddle-node*) leading to the appearance of an oscillatory solution. Rinzel and Ermentrout (1998) were the first to note that this behavior coincided with the properties of type I excitability. Ermentrout and Kopell (1986) showed analytically that any model system that undergoes a saddle-node transition from rest to repetitive firing is mathematically equivalent to a simple one-dimensional model for action potential generation, the so-called *θ-neuron*. (See also Ermentrout 1996; Gutkin and Ermentrout 1998; Hoppensteadt and Izhikevich 1997 for derivations and applications to spiking neurons.) The quadratic integrate-and-fire (QIF) model can be written as (1) where *I*(t) is applied, synaptic, and other currents; *V*_{rest} is the resting potential of the neuron; *V*_{thr} is the threshold for the spike to occur; τ_{s} is a time constant; and R_{m} is the membrane resistance (Latham et al. 2000). Near the onset of repetitive firing, all type I neurons satisfy dynamics identical to *Eq. 1*, indicating that the model is very general. Because the nonlinearity is quadratic, once *V* exceeds *V*_{thr}, the voltage can expand in finite time (for *I* ≥ 0). Most modelers truncate the voltage to a finite value and then reset *V* to a different finite value, as is done in the linear integrate-and-fire model. Hansel and Mato (2003) fit the *F*/*I* curve for the QIF to that of a standard cortical neuron model. Mathematically, the truncation is not a valid approximation for the dynamics; the spike should be at +∞ and the reset at −∞. Ermentrout and Kopell (1986) showed that a simple change of coordinates makes it possible to map the full dynamics of *Eq. 1* onto the circle. Let Then, the θ-neuron satisfies (2) where γ = R_{m}/(*V*_{thr} − *V*_{rest}). θ lies on a circle so that θ = 0 and θ = 2π are the same point in the action potential cycle. Figure 1*A* summarizes the physical interpretation of θ for the case when the applied current *I*(*t*) is zero. The resting state of the neuron occurs at θ = −2 arctan (1/2) ≈ −0.92 and the spike occurs when θ = π. If constant *I* of sufficient strength *I* > (*V*_{thr} − *V*_{rest})/(4R_{m}) is injected into the θ-neuron, θ travels counterclockwise around the circle and fires at a rate Ermentrout (1998) used this formula to fit the instantaneous *F*/*I* curve for a cortical neuron reported in Stafstrom et al. (1984).

Izhikevich (2000) extended the QIF model to include an additional term analogous to spike-frequency adaptation and showed that a version of this simple model could fit a wide range of neural firing patterns (Izhikevich 2004). Singular perturbation methods can be used to justify the additional term (Izhikevich 2000), resulting in a model of the form where *V*_{a} is the strength of the adaptation (in millivolts), τ is the decay rate, f(*V*) is set to be quadratic, and δ(*V* − *V*_{spike}) means that each time the neuron spikes, the adaptation is incremented. We include this in the θ-neuron *Eq. 2*, which becomes (3) where g_{z} = *V*_{a}/(*V*_{thr} − *V*_{rest}) and we use the fact that “firing” occurs when θ = π. To facilitate analysis, we take two steps: *1*) we “smooth” the delta function with a continuous pulselike function of θ and *2*) limit the adaptation at high frequencies by multiplying by (1 − *z*). These steps do not alter the dynamics of the model, and yet make the negative feedback correspond more closely to the behavior of a biophysical version of adaptation. In particular the adaptation term qualitatively describes the influence of slow voltage-dependent K-currents (e.g., muscarine-sensitive potassium current or “M-current”). We therefore transform the extended QIF model to a more tractable and general model for adaptation (4) This simple model is sufficient to produce SFA. Figure 2 shows the activity of the model under a DC injection.

To visualize the spikes we can plot an auxiliary quantity ν = 1 − cos θ. Note that this graphical “trick” takes care of the fact that in the original QIF formulation of the model the voltage *V* increases to infinity. On the other hand θ does not expand when the phase traverses the spike. Moreover in the θ-neuron, 1 − cos θ governs the intrinsic time evolution of the phase (as opposed to the phase-response function 1 + cos θ) and the activation of the slow feedback variable z. We also plot the time course of the negative feedback variable z and the instantaneous firing frequency (1/ISI). Figure 2 shows simulations for two different values of g_{z} (*A*, low; *B*, high). We see that as the negative feedback activates, the firing slows down so that the model shows firing rate adaptation. The absolute decrease in the firing rate is determined by the strength of the feedback (g_{z}) and the time course is set by its time constant. Note that even after the model has adapted to a steady-state firing rate, time-dependent fluctuations in the adaptation z are apparent. Characteristics of these fluctuations are determined by the above parameters, by the form of the activation function, and by the threshold for activation. The amplitude of the fluctuations of the slow feedback relative to the mean level of the feedback depends on the firing rate of the neuron: given a particular value of g_{z}, at higher firing rates the average amplitude of the z-current is high, yet the size of the transients is relatively low. On the other hand, at a lower firing rate the overall amplitude of the current is low, and relative spike-dependent deviations are strong. As we shall argue later this is the key to understanding the activity-dependent changes in the PRC shape and thus in neural response to transient inputs.

We further note that the extended θ-neuron is similar to a recent QIF model presented by Izhikevich (2004), where the voltage dependency of the negative feedback was (heuristically) assumed to be linear; here the dependency is sigmoidal, a better representation of the channel kinetics that also allows us to model both voltage-dependent and spike-dependent adaptation. In this study we have set parameters to ensure the adaptation feedback to be partially activated with “depolarization ” in the “subspike ” phase region. Voltage-dependent adaptation is modeled with a low-threshold outward current (as opposed to Ca^{2+}-dependent afterhyperpolarization [AHP] currents). However, a previous study shows that such choice does not qualitatively change our results. The parameters in the model are set as follows, unless otherwise stated in the figures: C = 2, θ_{T} = 3; κ = 8; τ = 400. The strength of adaptation g_{z} is varied in the simulations.

### The phase-response curve as a measure of neural response

The phase-response curve (PRC) provides a useful description of the neural response in particular and of nonlinear oscillators in general (Ermentrout 1996; Hansel et al. 1995). Suppose that a neuron is firing regularly with an interspike interval (ISI) T. The PRC measures the degree to which the time of the next spike of a neuron is changed when it is briefly stimulated at different times in the firing cycle. Let t denote the time since the last spike, let A parameterize the amplitude of the perturbing stimulus, and let T′(A, t) denote the time of the next spike after the stimulus. Note that T′(0, t) = T, the baseline firing rate. We can calculate the PRC in the following way Let T″(A, t) denote the time of the second spike after the stimulus. The usual assumption of PRC theory is that T″(A, t) − T′(A, t) = T, i.e., the stimulus is “forgotten” after one spike. In practice, this is not usually the case and we can define the second-order PRC (see also Oprisan et al. 2004) as Herein we use square current pulses of short-duration with amplitude A. Although this amplitude can have a qualitative effect on the PRC, we choose to focus on the timing of the perturbation and on the intrinsic firing properties of the neuron. We thus set the amplitude to be small enough not to cause too many “direct” crossings and leave the issue of the strong perturbations to future investigation. We must note that in the experimental work of Reyes and Fetz (1993a), the perturbations were quite strong in some of the experiments. Nevertheless, with our canonical model we can follow the dependency of PRC shape on the amplitude of the perturbation as seen in Reyes and Fetz (1993b) (simulations not shown). Our goal here is not necessarily to reproduce their results quantitatively, but to provide a general mechanism for the qualitative trends in the PRC shapes.

For sufficiently small A, one expects that PRC_{1}(A, t) will be nearly linear in A, so we can write PRC_{1}(A, t) = Aν*(t), where ν*(t) is independent of A. ν*(t) is called the infinitesimal PRC. The advantage of this approximation is that the function ν*(t) can be easily computed numerically and furthermore, in some instances, can be explicitly written down. In particular for the model *Eq. 2* (5) For *Eq. 4*, the function ν*(t) cannot be explicitly determined and we find it numerically using the numerical package XPPAUT (Ermentrout 2002). When the PRC was computed directly from simulations, the duration of the perturbation was set at 2 ms and amplitude A = 0.1.

### PRC and PSTH

The poststimulus time histogram (PSTH) is the standard way to measure the response of a neuron to a stimulus. A schematic PSTH is shown in Fig. 3. Typically the PSTH has a peak that rises above the base firing rate (1/T) at the onset of the stimulus, followed by a trough (Fetz and Gustafsson 1983). The PSTH is directly related to the PRC if the stimuli are sufficiently small. Suppose that a neuron is firing with a mean ISI of T. Then the PSTH before the stimulus will be flat with a magnitude of 1/T in spikes per millisecond. To relate the PSTH to the PRC, we need to make two assumptions about the strength of the stimulus:

*1)* it is weak enough so that no new spikes are added on average (Fetz and Gustafson 1983), and

*2)* the slope of the PRC is > −1/T.

Both of these assumptions imply that the stimulus is small. The first assumption means that the transient stimulus advances the occurrence of the second spike in the ISI without directly evoking a third spike. The second assumption means that the falling edge of the PRC is shallower than 1. Moreover, T′(A, t_{1}) > T′(A, t_{2}) if t_{1} > t_{2}, i.e., a late stimulus can never make a neuron fire sooner than an early stimulus. To relate the PRC and PSTH, consider Fig. 3, which shows the PSTH during the interval of length T before the stimulus and the interval right after the stimulus. *Assumption 1* implies that the area under the PSTH over the second interval is 1. Suppose a spike occurs at time s in interval 1. Without the stimulus, we expect the spike to occur at the same time in the second interval because the average ISI is T. However, the stimulus causes the spike to occur at a time s′ instead. We can determine the time s′ from the PRC. The stimulus arrives at a time T − s after the spike at s in the first interval. Thus the phase will be changed by an amount PRC(T − s). Thus the time of the spike will be changed by TPRC(T − s) because the PRC measures phase rather than time. Thus because the PRC measures the spike advance. The function G(s) maps the prestimulus interval into the poststimulus interval; from *assumption 2*, the function G is invertible. [G′(s) = 1 + TPRC′(T − s) is positive from *assumption 2*.] The PSTH tells us the probability that a spike occurs at time s_{0} in the interval after the stimulus. Thus where the last expression arises from the fact that the spikes are uniformly distributed before the stimulus. This shows that, given the PSTH, we can obtain G^{−1}(t) and thus the PRC. A similar procedure was used to reconstruct the S-D curve (the experimental correlate of the PRC) from the cross-correlation histogram (Reyes 1990).

### Numerical methods for noise, variable inputs, and input–output cross-covariance

For the simulations that examined the effects of noise on the PRC (see results), the model was stimulated with DC input + noise (variance = 0.01 in units of θ) to yield randomized spike times (1,000 sample paths lasting 1 s each). The mean firing rate was computed to determine the control ISI duration to be used in computing the PRC. We again ran the simulation but with transient inputs. We then compiled and plotted the changes in the spike times relative to control for each run. The spike times from 1,000 sweeps were averaged to construct the PRC.

For the slow-current–modulated PRC calculations (see results), we computed the times of each of the action potentials arising from the slowly varying current. We then picked several ISIs of various lengths and applied the perturbing stimulus at different times within each of these intervals. From this, we compute the PRC for this interval (see above) and compare it to that computed from the stationary firing rate. Note that during the interval, the bias is changing. The applied slowly varying current is which is quasi-periodic. We apply it over 1 s with *I*_{0} = 0.65, *I*_{1} = 0.2, and find a range of ISIs from about 18 to 150 ms. From this data set, we choose two intervals with an ISI of about 18 ms and one with an ISI of about 48 ms.

For the relationship between transient input timing, PRC shape, and output spike timing, we generated random input times by sampling a Gaussian distribution centered on a predetermined mean time and with a given variance (=10 ms for all simulations shown). The mean times were set so as to correspond with different levels of the PRC (e.g., the peak and the minimum). The inputs were square injections of amplitude 0.01 and duration 1 ms.

We first computed the noiseless PRC for a variety of ISI durations and adaptation levels. For a given base ISI duration and SFA level, the mean of the input times was set to the firing phase (time relative to the first spike) that corresponded to either the minimum or the maximum of the PRC. For the strong adaptation and firing rate of 40 Hz this corresponded to input times centered on 5 and 20 ms, relative to the first spike. Then the model was integrated for one firing cycle (see methods on constructing the PRC) with the randomly timed input. This generated a distribution of the output spike times. We then computed the covariances between the output spike times and the input times using standard methods. This measure allowed us to see whether the maxima/minima of the PRC can predict the timing correlations between the inputs and the output spikes.

### PRCs for conductance-based model

We used Traub’s model (Traub et al. 1999) to compute the infinitesimal PRCs at various firing rates; the parameters of the model were set as in Ermentrout et al. (2001). The slow potassium current conductance was set to zero to ensure that the model showed no SFA. The amplitude of the input current was chosen to yield various firing rates. The PRC was then computed using the standard method from XPPAUT (Ermentrout 2002) and normalized to give the maximum amplitude of 1.

## RESULTS

### Properties of the θ-neuron

There has been accumulating evidence that pyramidal neurons obey type I spike-generating dynamics (e.g., Fricker and Miles 2000; Stiefel et al. 2003; Stoop et al. 2000; Tateno et al. 2004). We previously showed that the θ-neuron—and thus the type I theory—successfully reproduces key input–output characteristics of cortical neurons: the frequency–current (*F*/*I*) relation, the slow approach to the spike threshold, the all-or-none action potentials, refractory periods, changes in the input resistance during firing, and ISI variability (e.g., see Gutkin and Ermentrout 1998; Gutkin et al. 2003). This model has also been used to study the effect of SFA on emergent network synchronization (Ermentrout et al. 2001) and sustained activity in models of cultured neurons (Latham et al. 2000). We will use this model to describe the complex responses of neurons when perturbed by small EPSP-like inputs. The model is derived from conductance-based models with type I excitability (see methods for formal arguments) and can be generally applied to all neurons that exhibit type I excitability.

### Theoretical PRCs reproduce experimental data

PRCs were measured experimentally in cortical neurons (Reyes and Fetz 1993a,b; PRCs were termed shortening-delay or S-D curves). In these studies, the PRCs were constructed from repetitively firing neurons (10–30 Hz) that were probed with brief current pulses. The current transients shortened the ISI by an amount that depended on which part of the ISI the input was delivered. Early in the ISI, when membrane potential was well below threshold, the input shortened the ISI by an indirect albeit regenerative mechanism (Reyes and Fetz 1993a). When delivered late in the ISI, the input shortened the ISI by directly crossing threshold and evoking a short latency spike. Figure 4*A* (rectangles) shows a representative example from Reyes and Fetz (1993a). Although most of the PRCs had positive values, some showed small negative regions for inputs arriving at the beginning of the ISI. Another notable feature was that the peaks of most PRCs were considerably skewed to the right. Note that for the PRC in Fig. 4*A*, the slope of the falling edge (at phase delay >0.7) is −1 toward the end of the interval, indicating that the neuron fired with a very short delay after stimulation. Weaker stimuli do not result in instantaneous firing.

As shown by *Eq. 5*, the “pure” θ-neuron model has a symmetric PRC that is strictly positive. Our previous work (Ermentrout et al. 2001) suggested that the rightward skew in the PRC peak was a result of slowly activating potassium currents that caused the firing rate of neurons to adapt. When activated these currents differentially decrease the effectiveness of the perturbing stimulus at the start of the ISI where the adaptation current is maximally activated. As the current deactivates with time in the ISI, the effectiveness of the input recovers. As a result, the period of maximum ISI shortening (i.e., the PRC peak) occurs toward the end of the ISI.

To examine the effects of slowly activating outward currents in a general manner, we modified the θ-neuron to include a negative-feedback current that replicates adaptation. As described earlier (see methods) this slow feedback may be viewed as an amalgam of voltage-dependent potassium currents that underlie SFA in cortical pyramidal neurons.

Figure 4*B* shows the PRC for θ-neuron with different SFA levels. The bias current was adjusted to produce comparable firing rates. The PRC is computed in this and all remaining simulations by delivering a square pulse (duration: 0.01 ms; amplitude: 1). Here we set the neuron to fire at 20 Hz. The skew in the PRC peak increases with increasing levels of g_{z} (*Eq. 4*), and the maximal spike displacement (i.e., the maximal amplitude of the PRC) is reduced with increasing g_{z}. By adjusting the maximal strength of the “adaptation-current” (g_{z}), we were able to match the experimentally measured PRC with one computed using the θ-neuron (Fig. 4*A*). The steady-state firing rate of the model and cortical neuron were comparable. In Ermentrout et al. (2001), we matched the experimental results to a specific conductance-based model. Here we have generalized this result by showing that the θ-neuron with adaptation can fit the data. In fact, because the θ-neuron arises through a reduction of any type I neuronal model, the present results show that any type I conductance-based model with adaptation could be fit to the data, giving a general result.

### Factors that affect the PRC time course: effects of firing rate

We now show how the PRC can be used to elucidate how the effects of transient input on spike timing are related to the firing rate of the cell and the underlying intrinsic ionic conductances. For a neuron with no or weak adaptation (e.g., g_{z} < 0.5), the steady-state firing rate does not affect the overall shape of the PRC. Thus we see that an incoming EPSP will always advance a spike (see Fig. 5*A*). An EPSP will displace an individual spike (shorten the ISI) much more at the lower firing rates than at higher rates. For the neuron with a significant adaptation (e.g., g_{z}) the situation is different. The skew in the PRC is most prominent at low firing rates (≤20 Hz) and decreases substantially at higher firing rates (≥40 Hz). Thus the window of opportunity for an incoming fast input to affect the timing of the spike narrows relative to the total ISI, although within that window an incoming EPSP has a stronger effect (Fig. 5*B*).

A possible heuristic biophysical explanation for this effect is as follows. We suggested earlier that the skew in the PRC depends on the time course of activation of K^{+} currents (e.g., M-current) in the ISI. In effect, activation of the K-current at the beginning of the interval renders the neuron less responsive to the transient input in the earlier as compared with the later part of the ISI. The transient becomes more effective toward the end of the ISI as the current deactivates. Furthermore, at least in pyramidal cells, late inputs can activate voltage-gated inward currents that enhance their effectiveness. As shown by Reyes and Fetz (1993a), the amplitude of the voltage responses to transient input is minimal early in the ISI (less than the response measured at rest) and maximal late in the ISI (exceeds the response at rest). These along with the ramp increase in voltage trajectory in the ISI mean that input is more effective toward the end of the ISI, resulting in skewed PRCs.

### Long-lasting effects of the transient inputs: PRC_{2}

In cortical neurons, the transient pulse that shortened the first ISI, lengthened the following ISI (Fig. 6, negative-going trace; Reyes and Fetz 1993a). Similar results were obtained for the θ-neuron with adaptation (Fig. 6), where the PRC calculated for the following ISI (PRC_{2}) is predominantly negative. In general, when the first-interval PRCs of the model and experiments matched, so did those of the following ISI (analysis not shown).

The mechanism for this effect can be easily understood. The shortened first ISI leads to a transient activation of the slow potassium current, causing a negative bias that lasts into the following spike cycle. This bias lengthens the following ISI by effectively putting a forward brake on the cell excitability (an effect noted before in a more complex SFA model; Wang 1998). The ISIs recover to control levels at the timescale of the adaptation.

### Effects of noise and time-varying inputs on the computed PRC

In vivo, neurons are constantly bombarded by synaptic inputs, resulting in background noisy voltage fluctuations. It is therefore important to determine whether the PRC can be recovered in the presence of noise. In these simulations, noise was added to the DC bias so that the average firing rate was maintained at a specific level but with greater variability in the spike times. The PRCs can be readily recovered from noisy data by averaging (Fig. 7) for both high and low firing rates. In general, the recovered noisy PRC reproduces the deterministic PRC even when the variance of spike times for individual trials was quite high.

PRC theory requires that the underlying system be periodic. However, under in vivo conditions, the firing rate of an individual neuron varies over time. To determine the validity of our analyses under nonstationary conditions, we apply a slowly varying (rather than constant) current. We compute the times of each of the action potentials arising from this slowly varying current. We then pick several ISIs of various lengths and apply the perturbing stimulus at different times within each of these intervals. From this, we compute the PRC for a particular interval and compare it to that computed under the stationary condition. Note that during the interval, the bias is changing (see methods). We apply it over 1 s with *I*_{0} = 0.65, *I*_{1} = 0.2, and find a range of ISIs from about 18 to 150 ms.

From this data set, we choose two intervals with an ISI of about 18 ms and one with an ISI of about 48 ms. For the short ISI, the stationary PRC and the transient PRC are quite close both in the beginning and at the end of the spike train (Fig. 8, *A* and *B*). This suggests that slow nonstationarities, resulting say from the buildup of adaptation, do not affect short-ISI PRCs. For the longer ISIs, there is some discrepancy, but most of the qualitative features we described above for the steady-state PRC, such as the pronounced skew, remain. We note that this also confirms a point made by Wang (1998) that adaptation acts as a differentiator; the neuron is even more insensitive to stimuli after a period of excitation.

In Fig. 8*C* low-pass filtered white noise was added to the DC bias. As with the time-varying input (Fig. 8*B*), the resultant PRC is more skewed than that for constant firing rate, although the basic shapes remain very similar.

### Relationship between input and output spike-time distributions are predicted by PRC

The results above clearly show that under low firing rates and high SFA conditions the cell becomes differentially sensitive to inputs in a specific part of its firing cycle: toward the end of the ISI, when the membrane excitability has partially recovered from the slow negative-feedback current. Thus the efficacy of the presynaptic inputs, albeit weak in absolute strength, is dependent on their correlation with the postsynaptic spikes. This suggests that the strongly adapting neuron at low firing rates acts to detect “coincidence, ” in a loose sense of the word, between the presynaptic inputs and its own firing times.

We quantify this by computing the relationship between the distribution of spike times in our model (under different excitability conditions) and the distribution of synaptic input times. To do this, the model, biased to fire at a predetermined mean firing rate, was stimulated by inputs with a set amplitude (see methods) and times that were Gaussian distributed with a preset variance. The mean input time was centered either at the time corresponding to the maximum or the minimum of the deterministic PRC at the low firing/high adaptation condition. If the inputs at the time when the PRC peaks are indeed more efficient, the distribution of the postsynaptic spike times should depend strongly on the distribution of such inputs. For those inputs that have the mean at the minimum of the PRC, such dependency should be relatively weak. The model was biased to fire either at 10 or 40 Hz (Fig. 9 *left*/*right*, respectively) and with g_{z} values of 0 (Fig. 8*A*) or 4 (Fig. 9*B*). The input time variance was set at 10 ms and the mean set at either 20 or 80 ms, for the 10-Hz condition; 5 or 20 ms for the 40-Hz condition. To quantify the input–output timing dependency, we computed covariances in the input arrival times and the output spike times (Fig. 9).

For the strongly adapting cell firing at 10 Hz the peak of the PRC is about 80 ms (80% of ISI-Control) and the minimum is near 20 ms (20% of ISI-Control). Clearly, the covariance of spike times and input times centered at 80% of ISI-Control is higher (tenfold) than the same centered at 20% of ISI-Control. For higher firing rates and for the nonadapting neuron the differences between the respective covariances are much lower. Furthermore for the nonadapting neuron the change in the two respective covariances is largely insensitive to the mean firing rate. We also calculated the mutual information measure between the output spike times and the input times (results not shown) and found it to agree with the covariance calculations: inputs centered at the peak of the PRC yield higher mutual information with the output spikes than those that fall in the minimum of the PRC. We suggest these simulations support our conjecture that currents that produce SFA can produce a “coincidence”-detector behavior at lower firing rates.

### Relationship between the θ-neuron PRC and PRC of conductance-based model

A valid question to pose is: What might be the caveats on the mechanism we had proposed and how does it relate to conductance-based models? The θ-neuron infinitesimal PRC with no adaptation is symmetric and bell shaped and its shape is invariant with firing rates. As indicated earlier the derivation presupposes low firing rates and implies that the θ-neuron PRC should theoretically be a valid description of neural PRCs at relatively low firing rates (near the onset of the repetitive firing or rheobase). However, the mathematics do not define quantitatively the region of validity. We thus compared the θ-neuron PRCs with those computed for a conductance-based model for a pyramidal neuron. We use the Traub model (see methods), set the adaptation to zero (g_{m} = 0), and compute its PRCs for a variety of firing rates. We see that for a relatively large range of firing rates the PRC has a relatively small skew (Fig. 10). At higher firing rates a skew develops in the PRC, which is attributed to the influence of the relative refractory period brought about by the fast potassium current (the delayed rectifier). Thus we see that the K-based mechanism we propose for the skew at low firing rates also carries over to high rates. The key of course is in the different timescales of the current dynamics: the slow currents act at low rates; the fast ones do so at high rates. At the higher firing rates, in addition to the mechanism we propose it is likely that a sodium-based mechanism helps to induce the skew: the partial inactivation of the sodium current after the spike (Reyes and Fetz 1993a).

### PRC shape at low firing rates in strongly adapting neurons: when excitation delays the spike

For the sufficiently strongly adapting neurons, a negative region may appear in the PRC at the start of the ISI. A negative region in the PRC can occur at very low firing rates (see Fig. 11). Thus the EPSP can actually advance or delay an action potential depending on its timing with respect to the previous (“conditioning”) spike. Note that the delay in the spike time would correspond to a decrease in firing probability (Reyes and Fetz 1993a,b). The presence of this negative region depends on the steepness of the θ dependency of the z activation: the more partial activation is present at resting θ levels, the more likely is the negative part to appear. In other words adaptation with lower activation thresholds can yield negative regions in the PRC. This is in contrast to, for example, the adaptation currents that activate only when action potentials are present (those with high activation thresholds, such as Ca^{2+}-dependent K^{+}-currents), which do not cause a negative region. In terms of the modified θ-neuron, the spike-dependent adaptation occurs when C is large and θ_{T} is close to π in *Eq. 4* so that the dynamics of z approach those of *Eq. 3*. Previously we presented a dynamical mechanism for this apparently puzzling effect of low-threshold adaptation in a conductance-based model of cortical neurons: a change in the bifurcation structure of excitability (Ermentrout et al. 2001). Although we expect this mechanism also holds for our modified θ-neuron, such analysis is well beyond the scope of this report.

In biophysical terms we can offer the following a posteriori scenario. We need the perturbing stimulus to cause enough additional transient activation of this “negative” current to delay the subsequent spike. We suggest that this effect is explained by the relative size of the transients in the level of *I*_{z} and at higher firing rates this negative region does not appear in the PRC because *I*_{z} is saturated and shows relatively small spike and stimulus-evoked transients. As mentioned earlier, the skew in the PRC is produced by the state-dependent time dynamics of the adaptation feedback, after the neuron has reached its steady-state firing.

## DISCUSSION

In this report we have explored the effects of small EPSP-like perturbations on the firing characteristics of the θ-neuron. The θ-neuron is a canonical model for the type I spike-generation class of neural oscillators. We have augmented this basic model with a slow negative feedback to model low-threshold spike-frequency adaptation (SFA), thereby arriving at a two-dimensional canonical model for adapting regular-firing neurons. In principle such slow feedback can model any number of currents that change the excitability of the neuron as a function of its firing rate. We then used this canonical model as a tool to explore the behavior of the whole dynamical class. We can view the small input perturbations as short-term temporal correlations in the inputs that separate themselves from the background barrage of PSPs. To characterize the responses of a periodically firing neuron to transient inputs we computed or simulated phase-response functions. Previously, Ermentrout (1996) analyzed the neurons’ phase-response curve (PRC) and showed that the PRC was strictly positive and symmetric for type I neurons without SFA. Here we show that the SFA controls the shape of the PRC in a firing-rate–dependent manner.

Using the θ-neuron, we have studied the critical role that the slow hyperpolarizing currents play in shaping the changes in spike timing in response to transient inputs and the neurons’ basal firing rate. Given a neuron driven to firing by some slow input (e.g., a large number of asynchronous synaptic inputs), the slow hyperpolarizing currents reduce the portion of the ISI where the EPSP would have a significant effect, thus skewing the response function toward the end of the ISI. This implies that neurons with strong SFA are most sensitive to inputs that co-occur just before the postsynaptic spike. Such neurons would tend to act as pre-/postsynaptic coincidence detectors. Nonadapting neurons, on the other hand, are more like integrators; their response is less dependent on the time of the inputs. We expect that spike frequency adaptation is only one example of a number of processes that can produce a skew in the PRC, and that neurons with such skewed PRCs are more likely to show coincidence behavior than neurons with symmetric PRCs.

Furthermore, a negative region in the PRCs of the strongly adapting neurons at low firing rates shows that EPSPs, arriving right after the cell fires, have in fact a delaying effect on timing of the subsequent spike. This happens because the spike-evoked transients in the slow currents are large relative to the average amount of that current. During these transients the slow K-currents are more liable to decrease differentially the excitability of the membrane during the first part of the ISI; in other words the membrane spends more time in a state where it is not responsive to the transient input because the K-currents represent a relatively larger proportion of the total cross-membrane current. These transients can be quite fast because they come about as a result of the changes in the driving force of the K-currents. As the basal firing frequency of the neuron increases, these transients become smaller, relative to the total amount of the slow K-current. So the K-current acts like another DC input, PRC becomes less skewed, and the neuron becomes less of a coincidence detector. This point is further supported by analyzing spike-time responses with randomly timed inputs (see above).

Previously it has been shown that repeated “frozen-noise ” inputs can lead to precise timing of individual spikes across multiple trials (e.g., Calvin and Stevens 1968; Mainen and Sejnowski 1995). We examine the effects of a single EPSP and show that even a small input transient (overlaid on DC background) can effectively shift spike times. When such input is repeated at the same time across multiple trials, it will shift the spike times reliably and lead to repeated spike times. The results on the strongly adapting case further suggest that transients that come at the end of the neurons’ firing cycle should be most efficient at setting the spike times, and such effect should increase with decreasing firing rates. Thus we can view our results as an extension of the spike-timing hypothesis laid out in Mainen and Sejnowski (1995). Our theory shows us how the mean firing rate of a neuron interacts with the effects of transient inputs on the timing of individual spikes. For a neuron with a given amount of spike-frequency adaptation, at lower firing rates, the PRC is strongly skewed. The maximal spike displacement stemming from a transient input is greater than that for the same neuron firing at higher rates. Thus as expected, at lower firing rates, inputs have a much greater effect on the spike timing than for high firing rate. At the lower firing rates, however, these inputs must be more temporally correlated with each other and, most important, with the postsynaptic neuron, than at the higher firing rates. As the relative contribution of the slow potassium currents to the overall firing pattern of the neurons decreases, the need for temporal correlations in activity decreases, resulting in a weakly adapting cell, and the firing rate has little effect on the skew of the PRC.

The reason that a neuron is less responsive to transient excitatory inputs at a higher firing rate is quite intuitive. At the lower firing rates, the neuron spends most of its time in the subthreshold voltage range where the cell membrane is at maximal input resistance levels (as compared with that during the spike or right after the spike). That is, in the ISI the membrane has a chance to recover from the influence of the fast currents that underlie the spike and also the slow currents, such as AHP. In fact at lower firing rates, the average level of such currents is relatively low. When the neuron is firing at a higher firing rate, it spends much less time near the rest and much more time in the hyperpolarized or depolarized state, where voltage-dependent currents are active to various degrees. Such currents dominate the dynamics of the membrane and the inputs are decreased in their ability to perturb the phase of the oscillation (e.g., shift the spike times). The voltage-dependent conductances that underlie such active currents decrease the input resistance of the membrane, and the transient inputs are shunted. Thus we are inclined to suggest that the neuron’s overall sensitivity to inputs depends on the overall level of the currents and conductance set by the overall firing rate. Note that this argument also implies that at higher firing rates neurons should be also less sensitive to transient inhibitory inputs.

We have used the PRC theory and a canonical neural model to study the interaction between mean firing rate and the timing of spikes. When a neuron is undergoing stimulus-evoked repetitive firing (or oscillations), one can argue that it is in the rate-coding regime. It responds by modulating its overall firing frequency to the slow input current resulting, for example, from a large number of asynchronous synaptic inputs (e.g., Destexhe and Par é 1999; Harsch and Robinson 2000), whereas when the neuron is responding to repeated presentations of specific input patterns, it produces precisely timed action potentials (Bair and Koch 1996; Calvin and Stevens 1968; Fellous et al. 2000; Mainen and Sejnowski 1995; Nowak et al. 1997), and is in the spike-time–response regime. In this mode the membrane characteristics of the neuron undergo slow modulation as a result of the activation of the voltage-dependent conductances. These fluctuations in turn affect the way in which rapid synaptic inputs set the timing of the individual action potentials. We have quantified such spike-time setting interaction by the shape of the PRC. We showed that this interaction is firing rate dependent. Thus the PRC formalism allowed us to quantify the interaction between the rate and spike-time response modes. Finally, the effects of SFA and the firing rate on the shape of the PRC also have implications for the structure of neural activity at the circuit level.

Previously, Hansel et al. (1995) and Ermentrout (1996) gave general results, stating that all nonadapting neurons with type I excitability will antisynchronize (phase-lock half a cycle out of phase) with positive coupling. To arrive at this conclusion Ermentrout (1996) analyzed the neurons’ PRC and showed that antisynchrony is a result of the PRC’s being strictly positive and symmetric. Subsequently, Crook et al. (1998) showed that adapting type I neurons readily synchronize with excitation. At the time this seemed to present a contradiction. However, in Ermentrout et al. (2001) we analyzed the PRC shape for a specific type I biophysical model and found that increasing the skew in the PRC by slow hyperpolarizing currents can lead to stabilization of synchrony. A similar role has been proposed for AHP in stabilizing the alpha-range activity in the hippocampus (Kopell et al. 2000). This stabilizing effect of the SFA was linked with the shape of the PRC in conductance-based models of cortical and hippocampal neurons. Our results clearly show that the skew in the PRC for the canonical model depends on the SFA level and thus the results presented in Ermentrout et al. (2001) are generalizable to the whole type I class of models. Furthermore, this work should be viewed as complementary to the previous studies because we focus on the meaning of the PRC with respect to the spike timing of neural responses to synaptic stimuli as opposed to a technical tool for analysis of synchrony.

Herrmann and Gerstner (2001) developed a theory for the dependency of the PSTH on both the shape of the PSP and the amount of background noise. They obtained analytical results by linearizing the full nonlinear system and were thus able to use an impulse response function to obtain the response to any input. Their methods are applicable to spike-response models and, in particular, the leaky integrate-and-fire model. They are mainly interested in the effects of noise amplitude on the PSTH, which in our case could be translated to background firing rate. We have shown a strong dependency of the PRC on the background rate, which implies a similar strong dependency of the PSTH.

Overall, the PRC analysis has shown that SFA has three main effects on how small transient (e.g., single synaptic) inputs determine the timing of the subsequent spike. First, in a neuron with strong adaptation, the inputs that come later in the ISI have a stronger effect on the timing of the subsequent spike. Second, in neurons where adaptation results from voltage-dependent currents (e.g., M-type potassium current), an input early in the ISI can actually delay the next spike. Third, the total effect an input may have on the timing of a spike decreases with the increasing strength of SFA. Thus for stronger adapting cells, inputs must be more correlated and stronger to affect the timing of a spike.

For a given level of SFA, increasing the firing rate decreases the skew in the PRC, implying that the stability (or robustness) of synchrony is also dependent on the average activity in the circuit. In fact we speculate that at higher frequencies SFA-dependent synchrony may be lost, and another mechanism, such as local inhibition, may play the dominant role at higher frequencies. Finally, the level of SFA is affected by a number of neuromodulators. For example, acetylcholine tends to decrease muscarinic-sensitive currents as well as other slow K-currents in neural cells (e.g., see Schobesberger etal. 2000; Tang et al. 1997) as well as other hyperpolarizing currents (Haj-Dahmane and Andrade 1996; Klink and Alonso 1997). Dopaminergic agonists also appear to decrease SFA by turning off slow potassium currents (e.g., Henze et al. 2000; Yang and Seamans 1996). We speculate that such modulatory effects may lead to a decrease in the strength of synchrony at the lower-frequency range. This is indeed hinted at by experiments where the activation of cholinergic subcortical structures reduces global synchrony at the lower-frequency band of the cortical EEG (e.g., Metherate et al. 1992; Munk et al. 1996) and increases the mean firing rate and induces local synchronization at the higher frequency bands (e.g., Herculano-Houzel 1999; Munk et al. 1996; Steriade et al. 1996). This is again consistent with our hypothesis that the action of the modulators (e.g., acetylcholine or dopamine) would lead to decreasing skew in the PRC and decrease stability of global synchrony in the lower-frequency band.

It should be pointed out that recently PRCs were measured for human motor neurons (Mattei and Schmied 2002), with results essentially similar to those found by Reyes and Fetz (1993a,b). The same study looked at long-duration effects and showed that the shortening and lengthening alternates for successive ISIs. Also, Stoop et al. (2000) measured the response functions attributed to inhibitory inputs in rat visual cortex slices and found that these are inconsistent with integrate-and-fire neurons. Oprisan et al. (2004) focused on the second-order phase resetting and showed how the shape of the second-order PRC affects synchronization of neural oscillators. We note here that, on visual inspection, the experimental PRCs in those studies are consistent with type I models.

In conclusion our theory allows us to predict that spike-time setting arising from transient inputs should be the strongest in neurons at lower firing rates, and that in neurons with strong spike-frequency adaptation, inputs should be temporally correlated to the postsynaptic neuron to have significant effects on spike timing. Once again we point out that SFA is not the only way to produce a skewed PRC. Such skew may develop as the result of any process that changes a neuron’s excitability during specific portions of its firing cycle (e.g., inactivation of fast sodium currents, persistent sodium currents, etc.), and such a neuron would also be sensitive to temporal correlations between the inputs and its own firing times. However, our theory makes a strong experimentally testable prediction: the skew of the PRC can be changed by blocking or activating the slow potassium currents that underlie SFA. In fact, preliminary experimental data confirm the prediction of our theory that blocking of the M-current by muscarinic agonists can decrease the skew of the PRC (Stiefel et al. 2003).

We show that type I phase-response curves, produced by the θ-neuron, match the measured PRCs of cortical pyramidal neurons, and reproduce a number of effects of intrinsic currents, stimulus size, and timing observed in experiment. Phase-response curves are also a clear way to summarize the interactions between the intrinsic membrane properties, the mean firing rate, and the effects of transient inputs on spike timing. The transient PRCs are closely matched by steady-state PRCs if the stimulus changes sufficiently slowly. We maintain that this also provides further evidence that membrane excitability and spike-generating dynamics in cortical neurons are of type I as defined by Hodgkin (1948) and Rinzel and Ermentrout (1998). In other words, cortical neurons with weak spike-frequency adaptation give a biological instantiation of dynamical systems with local saddle-node bifurcations.

## GRANTS

This work was supported by National Science Foundation (NSF) grants to G. B. Ermentrout and B. S. Gutkin, NSF Bioinformatics Postdoctoral Fellowship to B. S. Gutkin, National Institute on Deafness and Other Communication Disorders Grant DC-005787-01A1 and NSF Grant IBN-0079619 to A. D. Reyes, and a Gatsby Foundation grant to B. S. Gutkin.

## Acknowledgments

We are grateful to P. Dayan and P. Latham for comments on earlier versions of the manuscript.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2005 by the American Physiological Society