## Abstract

Phase response curves (PRCs) are a simple model of how a neuron's spike time is affected by synaptic inputs. PRCs are useful in predicting how networks of neurons behave when connected. One challenge in estimating a neuron's PRCs experimentally is that many neurons do not have stationary firing rates. In this article we introduce a new method to estimate PRCs as a function of firing rate of the neuron. We call the resulting model a parameterized PRC (pPRC). Experimentally, we perturb the neuron applying a current with two parts: *1*) a current held constant between spikes but changed at the onset of a spike, used to make the neuron fire at different rates, and *2*) a pulse to emulate a synaptic input. A model of the applied constant current and the history is made to predict the interspike interval (ISI). A second model is then made to fit the modulation of the spike time from the expected ISI by the pulsatile stimulus. A polynomial with two independent variables, the stimulus phase and the expected ISI, is used to model the pPRC. The pPRC is validated in a computational model and applied to pyramidal neurons from the CA1 region of the hippocampal slices from rat. The pPRC can be used to model the effect of changing firing rates on network synchrony. It can also be used to characterize the effects of neuromodulators and genetic mutations (among other manipulations) on network synchrony. It can also easily be extended to account for more variables.

- phase response curves
- ARX models
- ARMAX models
- patch clamp

many normal brain rhythms and neuropathologies are thought to be emergent behaviors of populations of neurons (Buzsáki 2006; Llinás et al. 1999). Epilepsy, Parkinson's disease, and Alzheimer's disease are pathologies in which synchrony between different cortical areas is affected (Uhlhaas and Singer 2006). Changes in these rhythms may be caused by changes in the dynamics of the neurons in the networks producing these rhythms or in the network coupling. However, relating changes in neuronal activity at the cellular level to network level changes remains a challenge central to understanding many diseases and developing new treatments.

One tool for linking the dynamics of individuals to population dynamics is the phase response curve (Glass and Mackey 1988; Tass 1999; Winfree 1980). The phase response curve (PRC) is a simple measure of an oscillatory system's response to an input. For a periodically firing neuron, the PRC is the average advance of the next spike time given the phase in the neurons cycle when a stimulus is applied. Given the PRC, it is possible to predict how a network of neurons will behave when coupled together (Achuthan and Canavier 2009; Canavier et al. 1997; Crook et al. 1998; Netoff et al. 2005; Van Vreeswijk et al. 1994).

The power of PRC model is its generalizability and flexibility. Any change in the neuron's dynamics, or stimulus waveform, that affects a neuron's spike timing can be characterized with a PRC. It can be used to characterize a neuron's response to pulsatile inputs, such as synaptic inputs (Netoff et al. 2005; Reyes and Fetz 1993a, 1993b), electric fields from deep brain stimulation (Wilson et al. 2011), or continuous waveforms, such as white noise (Ermentrout et al. 2008; Torben-Nielsen et al. 2010). Many different perturbations to a neuron can affect its PRC, such as changes in firing rate (Fink et al. 2011), position of synapse on the dendrite (Crook et al. 1998; Goldberg et al. 2007; Lewis and Rinzel 2004), neuromodulators (Stiefel et al. 2008a, 2008b), or synaptic plasticity (Lengyel et al. 2005), but there is not yet formalism for modeling how the PRC may change with a change of these parameters. Instead, a PRC model is usually generated for each condition.

Animal models of seizure-like activity show that both the neuron's firing rate and synchrony change over the seizure (Netoff and Schiff 2002). In computational models, it has been shown that the PRCs of neurons depend on the firing rate (Gutkin et al. 2005). Network synchrony has been shown to change dramatically, depending on the neuron's dynamics (Fink et al. 2011), and changes in synchrony have been used in models to explain shifts in population synchrony (Lewis and Rinzel 2003). Recently, we have used PRCs to understand the transition between the tonic phase of a seizure, where firing rate of neurons is high but synchrony is low, and the clonic phase, where the firing rate is lower but synchrony is high. We hypothesize that at high firing rates, a neuron's spike time is determined by the kinetics of the voltage-gated channels and is insensitive to synaptic inputs. As the firing rate slows, the dynamics of the neuron changes, increasing its sensitivity to synaptic inputs, which effectively increases coupling, until the network spontaneously synchronizes (Wilson et al. 2011).

However, a limitation in testing these theories experimentally has been in measuring PRCs from the same neuron across frequencies. Generally, PRCs are measured experimentally by inducing the neuron to fire periodically, by injecting constant current to the soma through a patch-clamp electrode, and then adding a stimulus to measure how the input may affect its next spike time (Reyes and Fetz 1993a, 1993b). To keep the neuron firing at the same period over the duration of the experiment, we have implemented a closed-loop spike rate controller that adjusts the applied current to compensate for drift in the neuron's interspike interval (ISI) (Miranda-Domínguez et al. 2010; Netoff et al. 2005). The use of a periodically firing neuron is not essential, but it is just the easiest condition to identify changes in the expected spike time. An alternative to modulating the applied current used to hold the neuron's firing rate constant over the experiment is to generate an adaptive model of the neuron's ISI that can predict the next spike time with high accuracy.

In this article we present a new method for measuring PRCs as the current applied to the neuron is varied. We excite the neuron using a current that has two components. The first component is a current that is changed at each spike time but held constant between spikes. The second component is a current that emulates a synaptic current. This transient current is applied at each ISI randomly. To analyze the resulting data, first, we generate an autoregressive model with exogenous input (ARX model) to predict the next ISI as the current applied to the neuron is varied. Given an accurate prediction of the unperturbed spike times, we then relate the mismatch between the predicted and the measured ISI to the timing of the transient input. This approach allows us to generate parameterized PRCs (pPRCs), a two-dimensional PRC where phase advance is predicted as a function of the stimulus phase and the firing rate of the neuron. We demonstrate our approach for estimating pPRCs in a mathematical neuron model proposed by Golomb and Amitai (1997) and in pyramidal neurons in the CA1 region of the hippocampal slices from rat.

## METHODS

The experiment is to vary the neuron's ISI by perturbing the applied current at the onset of each spike, and then a synaptic conductance is added at a random phase, as illustrated in Fig. 1. From these data we generate a model that predicts the ISI given the applied current and timing of the applied synaptic conductance. This model is created in two stages: In the first stage, a model is generated to predict the ISI given the applied current on this ISI and the history of applied currents and ISIs. This is done by recursively fitting an ARX model to the ISI as a function of the applied current. The ARX model can then be used to predict the firing rate as a function of the applied current. Because the ARX model is fit recursively, its parameters are able to adapt to nonstationary behaviors of the neuron over the duration of the experiment. The second stage is to fit a model to the phase advance from the expected ISI given the phase of the synaptic input. Because the PRC is dependent on the expected ISI, we generate a two-dimensional model that fits the advance given the phase and the expected ISI. The following sections describe the design of the experiment and the two steps of the analysis.

### Computational Model

To test this approach for generating models, we first measured pPRCs in a computational model of a CA1 pyramidal neuron developed by Golomb and Amitai (1997), which we refer to as the GACell. The full model is described in their report in detail, but we briefly describe the model here. The GACell is a single-compartment Hodgkin-Huxley-type model of coupled differential equations.

#### Current balance equation.

The total current in the GACell can be summarized in the following equation:
*C* = 1 μF/cm^{2} is the membrane capacitance and *V*(*t*) is the membrane potential of a neuron at time *t*. The right side consists of several time-varying currents, such as applied current, *I*_{app}, and voltage-gated intrinsic currents: the fast sodium current, *I*_{Na}; the persistent sodium current, *I*_{NaP}; the delayed-rectifier potassium current, *I*_{Kdr}; the A-type potassium current, *I*_{KA}; the slow potassium current, *I*_{K,slow}; and the leak current, *I*_{L}. For brevity, the details of the voltage dependence of these currents is left to the report by Golomb and Amitai (1997). The model was integrated using a fourth-order Runge-Kutta algorithm in the Real-Time eXperiment Interface (RTXI) environment (described in *Brain Slice Preparation*) at 100 kHz and sampled at 5 kHz. Current was applied so that the neuron fired at a steady-state firing rate of 10 Hz (i.e., a period of 100 ms). Experimental protocol for the GACell was exactly the same as for the real neuron with the exception that all stimuli were applied to the model rather than to the patch-clamp amplifier (described in *Experimental Design*).

### Brain Slice Preparation

PRCs were also measured in vitro from pyramidal neurons in brain slices prepared from the hippocampal formation (CA1) of rats. All experiments were performed following approved guidelines from Research Animal Resources of the University of Minnesota using an Institutional Animal and Use Care Committee-approved protocol. Long-Evans rats 14–21 days old of either sex were deeply anesthetized with isoflurane and decapitated. The brains were removed and placed into chilled artificial cerebral spinal fluid (ACSF) containing (in mM) 126 NaCl, 25 NaHCO_{3}, 25 glucose, 1.25 NaH_{2}PO_{4}·H_{2}O, 2.5 KCl, 1 MgCl_{2}, and 2 CaCl_{2}. Slices of 400-μm thicknesses were cut using a Vibratome and perfused with ACSF. The slices were aerated with O_{2} and incubated at 37°C for at least 1 h or until use.

Patch-clamp electrodes (6–10 MΩ) were pulled from borosilicate glass (P-97 micropipette puller; Sutter Instrument) and filled with intracellular recording fluid consisting of (in mM) 120 K-gluconate, 20 KCl, 10 HEPES, 0.2 EGTA, 2 MgCl_{2}, 4 Na_{2}-ATP, 0.3 Tris-GTP, and 7 Na_{2}-phosphocreatine. Pyramidal cells in the CA1 region of the hippocampal formation were whole cell patched and recorded using a current-clamp amplifier (MultiClamp 700B; Axon Instruments, Molecular Devices, Sunnyvale, CA). Data were acquired and stored using a real-time Linux-based system (Dorval et al. 2001) called RTXI, which allows implementation of complex closed-loop protocols running in hard real time and is publicly available (http://www.rtxi.org). Data were sampled with 16-bit resolution at 8 kHz.

### Experimental Design

Traditionally, we have measured PRCs by holding the neuron at a single firing rate over the duration of the experiment. However, in this method the current is held constant over the duration of the ISI but is changed at the onset of each action potential. This provides a data set in which the PRC can be estimated as a function of both the stimulus phase and the neuron's ISI. First, a baseline current (*I*_{mean}) is selected manually, which makes the neuron fire approximately at the desired ISI. Next, at each subsequent spike *i*, the current applied is selected randomly from a set of possible perturbations to *I*_{mean}, denoted as *I*_{DC}(*i*), and is held constant over the next ISI. The possible values of *I*_{DC} are selected randomly at each spike time from the set [0.9*I*_{mean}, 0.95*I*_{mean}, *I*_{mean}, 1.05*I*_{mean}, 1.1*I*_{mean}]. In addition to *I*_{DC}(*i*), a synaptic conductance input (*I*_{syn}) is added on each cycle. The timing of *I*_{syn} is selected using a uniform distribution over the expected ISI. The entire experiment is depicted in Fig. 1 (results presented are from the GACell model of Golomb and Amitai 1997).

*I*_{DC}(*i*) can adopt different discrete values. We chose five applied current values equally spaced around a mean current value *I*_{mean}. The difference in current values was chosen to be large enough to produce changes in ISI greater than the random variability in the neuron's ISI. The value of *I*_{DC}(*i*) is randomly selected at each spike time using a Markov process. A Markov process was used to choose the current values so that sequences of the same ISI may occur with reasonable probability so that steady-state behaviors could be estimated accurately. The Markov process also was designed to limit large jumps in the applied current from one cycle to the next. To achieve this, the probability of the current staying at the same value on the subsequent ISI was set equal 0.91, and the probability of jumping up or down one current value was set at 0.08, and two current values at 0.01. These transition probabilities were set so that the expected distribution across all current values was uniform.

The random sequence for *I*_{DC}(*i*) was generated in Matlab, stored to a file, and read by the RTXI module. Figure 1*C* shows the time course evolution of the applied current as a function of spike time along with a histogram of the amplitudes. A synaptic input is then applied at a randomly chosen phase on each ISI, as shown in Fig. 1, *A* and *D*. The applied synaptic current is calculated using a synaptic conductance waveform and membrane potential of the neuron as follows:
*t* is the time from the onset of the start of the synaptic input, *I*_{syn} = 0 for *t* < 0, τ_{fall} is the falling time constant (6.23 ms), τ_{rise} is the rising time constant (2.16 ms), *R*_{p} is the synaptic reversal potential, which we set at 0 mV to simulate excitatory synapses, and *G* scales the conductance to current. The value of *G* was empirically determined for each neuron to induce a measurable spike time advance without evoking action potentials, generally on the order of 400 pS. If the cell fires before *I*_{syn} is applied, no synaptic input is applied to the neuron and the interval is not used in estimating the shape of the PRC.

### Predicting the Neuron's Unperturbed ISI

To measure how a synaptic input perturbs a neuron from its normal period, an accurate prediction of its unperturbed ISI is required. The easiest model is to stimulate the neuron with a constant input and wait until it has settled to a steady-state ISI and then apply the stimulus. However, we are interested in PRCs at different firing rates. This could be done serially by stepping the current and measuring the PRC, then stepping the current again. Alternatively, this can be done by applying a different current at each spike over the duration of the ISI and fitting a model to predict the ISI based on the history of the applied current and the neuron's past ISIs. Here we predict the neuron's ISI as a function of the applied current and past ISIs using a quadratic ARX model:
*i* is the index of the current ISI and *i* − α the preceding α ISIs, and *I*_{DC}(*i*) is the applied current on the cycle to be predicted. The number of historical ISIs and applied currents are determined by the variables *m*_{k} and *n*_{k}, respectively.

Once the parameters (*m*_{k} and *n*_{k}) are selected, the *a* and *b* values can be determined by fitting the model to the data by solving for the least-squared solution. However, as the ISI drifts, it is prudent to estimate the parameters recursively, updating them at each spike to minimize the mismatch between the predictions and the data. Through recursive estimation of the model, drift in the neuron's dynamics is compensated.

The recursive estimation can be repeated multiple times to start with a realistic set of initial conditions each time the estimation is repeated. The recursive least-squares estimation of the parameters (i.e., the *a*′s and *b*′s) of *Eq. 3* is performed following the formulation presented by Simon (2006):
*y*_{i} represents ISI. *H* is a matrix that contains the values of the ISIs and *I*_{DC}(*i*) required to predict the ISI, as indicated by *Eq. 3*.

To select the number of a's and *b*'s, an initial guess for *a* and *b* values are stored in *P*_{0}:
*Eq. 7*: predicted ISI =

### Calculating the Neuron's Parameterized PRC

Once the next spike time can be accurately predicted, we can then add a synaptic input at a phase and measure how it perturbs the spike time. The “spike time advance” (STA) can be calculated from the residuals of the predicted ISI from the measured perturbed ISI given the phase *I*_{syn} was applied. The stimulus “phase” (*P*) is defined here as the stimulus time since the last spike normalized to the unperturbed ISI. Hence, *P* is bounded between 0 and 1.

The STA is calculated as a function of both the predicted ISI and the phase of the applied stimulus. We will use a polynomial with two independent variables to fit the STA as a function of the predicted ISI and the stimulus phase. The variables used to calculate the STA are the *na* recursive terms of the STA, i.e., the historical values of STA, the phase of the stimulus (*P*) and the predicted ISI given the applied current (ISI_{pred}). To solve for the PRC, these variables are concatenated into a matrix B. B has *n* rows, one row for each observation, and *m* columns, one for each independent measure and the cross terms.

From this matrix B and the set of appropriate weights, the STA can be calculated as follows:
*n* × 1 vector containing the *n* observations of spike time advances, and *m* × 1 vector containing the weights for each of the independent variables. These weights are the coefficients to the polynomial function from which the PRC can be calculated given ISI_{pred} and the stimulus phase, *P*.

The matrix B is determined as follows: the first *na* columns of *B*_{i×m} consist of the history of the STAs:
*np* historical values of the predicted ISIs, *Y* = ISI_{pred}, and the *nb* historical values of the stimulus phase, *P*. To obtain a polynomial, there are columns containing these values squared and cubed, and all the cross products of these terms. For *np* = 2 and *nb* = 2, the remaining columns of B are then
*m* = *na* + *nb*[(*np* + 1)^{2} − 1]. The unknown set of weights *Eq. 8* using linear algebra. However, simply solving for *Eq. 8* is ill posed. To correct for this we used the truncated singular value decomposition (TSVD) in B to solve for

The TSVD approach to solving ill-posed problems is to solve for the minimum norm solution for _{1}, …, _{m} and V = _{1}, …, _{m} have orthonormal columns, known as left and right singular vectors, respectively, and Σ = diag(σ_{1}, …, σ_{m}) contains the singular values of B on the diagonal. The singular values of B are ordered such that
*Eq. 11*, given the singular value decomposition of B (*Eq. 12*), is

The goal is to eliminate vectors that have small singular values, which generally have more sign changes in the elements of the singular vectors, thereby reducing the noise in the PRC estimate. If we choose *l* = *m*, then this solution is the same as the pseudo-inverse method. The singular values of B are ordered from major to minor, and truncating *Eq. 14* to an appropriate value of *l* eliminates the noisier components. The trade-off is between in-sample predictive power and model robustness to noise. If successful, a small decrease in in-sample prediction will actually increase out-of-sample predictive power.

## RESULTS

We first present results predicting the spike times as current applied to the neuron is changed at each spike time given the history of applied currents and ISIs. Given an accurate prediction of the unperturbed spike time, we then measure how the phase of a synaptic input perturbs the neuron from its period. We present results from a computational model of a pyramidal neuron, the GACell model, and then in pyramidal neurons recorded from hippocampal brain slices using dynamic clamp.

### Unperturbed Spike Time Prediction

The first step is to establish a model of the ISIs provided the history of the neuron's ISIs and the applied currents. Figure 2 illustrates the accuracy of our unperturbed spike estimation in the GACell (Golomb and Amitai 1997) model using the recursively fit ARX model (*Eq. 7*). The voltage trace (Fig. 2*A*) and applied current (Fig. 2*B*) during the first nine spikes are shown. The variables *m*_{k} and *n*_{k}, determine the number of historical ISIs and applied current. We set them to 5. They indicate the length of the memory of the neuron to changes in current and history of ISIs. For example, if there is a transient, it may take five spikes to settle; therefore, a history of five spikes may be necessary to make an accurate prediction. Predicted and measured ISIs are plotted in Fig. 2*C*, and Fig. 2*D* shows the predicted ISI against the measured ISI.

To characterize how well the model predicts the data, we measure the Pearson *R* correlation between the predicted and measured ISIs. The square of the *R* value is a measure of the percentage of variance in the data explained by the model. In the GACell, the *R* value between the predicted ISI from the ARX model and the actual spike times is 0.9942, indicating that the model explains 98.8% of the variance (*R*^{2}), leaving only 1.2% of the variance unexplained. The residuals from the ARX model can be attributed to higher order dynamics and the effect of the applied synaptic inputs.

It is interesting to note that although neurons are highly nonlinear, a simple linear ARX model is able to describe more than 99% of the ISI variability as a function of the applied current. Figure 3 shows the actual and predicted ISIs and their correlation from a pyramidal neuron in CA1. In this case, the correlation coefficient between the predicted and measured ISIs is *R* = 0.8346, explaining 69.6% of the variance, leaving only 30.4% unexplained. The decreased *R* value reflects nonlinear dynamics in the neuron that are not considered in the model. The ISI prediction was performed on 45 different cortical neurons. In more than 40 cells, the total variance explained was greater than 50% (see Fig. 7). The remaining five cells did not fire periodically; therefore, a linear model could not provide highly accurate predictions.

In this section we have presented the results of the prediction of the unperturbed spike times. It is important to note that this model is being fit and predicted concurrently with the application of synaptic inputs. Therefore, the same data set is used in the estimation of the unperturbed spike times as well as the measurement of the STAs with respect to the stimulus phase.

### Spike Time Advance Measurement

Phase advance was then determined from the perturbation by the stimulus input of the neuron from the expected period, as calculated from the residuals of the ARX model. Before fitting the PRC function to the residuals, we removed outliers that were greater than 3σ from the mean and replaced them with their expected values. These outliers usually were from missed action potential, where the period was much longer than expected.

Ideally, the residuals from the ARX model represent the perturbation in the spike times caused by the applied synaptic input. To make a model of this phase advance, we fit a polynomial with two independent variables to predict the spike advance as a function of the stimulus phase and the expected unperturbed ISI. Figure 4 shows the estimated PRC for the GACell model used in Fig. 2. Here we set *np* = 4 (the polynomial order), *na* = 0 (the number of synaptic inputs from the previous cycles), and *nb* = 1 (the number of currents applied over previous cycles). For the TSVD solution, we used the first seven singular values. The correlation coefficient between the model and the phase advance is 0.8477. This correlation explains 72% of the residual 1.2% not explained by the ARX model. The ARX model and the PRC together explain a total of 99.1% of the variance. In the resulting pPRC, it can be seen that the shape changes as a function of the firing rate (Fig. 4, *C* and *D*).

Figure 5 shows the corresponding PRC for the pyramidal neuron used in Fig. 3. For this model, the parameters were set at *np* = 4, *na* = 0, and *nb* = 1, as was done with the neuron model. Again, only the first seven singular values of B were used. In this experiment, the correlation coefficient of the predicted and observed STA was 0.63, explaining 40% of the remaining 72% from the ARX model, resulting in a total of 88.8% of the spike time variance explained. Like the GACell, the real neuron also shows a dramatic change in the estimated PRC with firing rate.

Figure 6 shows the PRC estimations that have the highest total *R* value. Figure 6 shows the ISIs over time, the data and fit PRCs in three dimensions, and a projection across the predicted ISIs in two dimensions, as well as the estimated PRCs for the fast, medium and slow spiking rates. Neurons were selected that were stable over the duration of the experiment, that had a low coefficient of variation of the ISIs, indicating that the neuron was a regular spiker, and that were recorded long enough to have at least a few hundred spikes. Notice that even some of the best cells have nonstationary firing rates (*cell 6*), and some of them skipped spikes (*cell 3*). The recursive ARX model was able to compensate for the nonstationarities, but as mentioned, the skipped spikes were removed as outliers before the models were fit. Across all cells, the estimated PRCs peaked to the right of center. Furthermore, the peaks of the PRCs were further to the right for the longer ISIs (slower firing rate) than for shorter ISIs (faster firing rate).

Figure 7 summarizes the results obtained from 45 pyramidal neurons. For each neuron, an adaptive ARX model was first fit to the data to infer the unperturbed ISI. The correlation coefficient (*R*_{ARX}) between the predicted and observed ISI is plotted using a solid gray line. The correlation coefficient between the predicted and the observed STA is shown using a dashed gray line (*R*_{REG}). The variances explained by the ARX model and the PRC model on the residuals can be added to account for the total variance explained. The cumulative variance is calculated as follows: *R*_{t}^{2} = *R*_{ARX}^{2} + *R*_{REG}^{2}(1 − *R*_{ARX}^{2}). The total variance is shown as a solid black line. In 30 cells, the pPRC model can account for more than 50% of the neuron's variance. The example from the neuron shown in Figs. 3 and 5 was the best case, corresponding to *neuron 1* in Fig. 6. The cases where the model could only explain a small fraction of the total variance were neurons that missed spikes regularly. These were neurons where 10% of the recorded ISIs had a value at least twice that expected. Although this is a common behavior of neurons, the model assumes the neuron fires periodically in response to a constant input and cannot generate an accurate model to predict the spike times.

### Comparison Between the Parameterized and Infinitesimal PRCs

Finally, to test the accuracy of the estimated pPRCs for the GACell, we compared them with the infinitesimal PRC (iPRC) calculated using the adjoint method. The adjoint method (Ermentrout 2002) estimates the iPRC from the GACell equations. The iPRC is an estimate of the neuron's response to an infinitesimally short stimulus pulse, a Dirac-delta function. The iPRC was measured using the XPP software package (http://www.math.pitt.edu/∼bard/xpp/xpp.html) at three ISIs (85, 100, and 125 ms) and solved at 0.1-ms time step resolution. The resulting iPRCs are plotted as dashed black lines in Fig. 8. Some discrepancies between our measured pPRCs (solid black line) and the iPRC (dashed black line) can be seen. This might be attributed to three different causes: the first is that the estimate of the pPRC is a polynomial fit to the data, which may smooth the curve a little; the second is that the pPRC was estimated using a finite stimulus waveform; and the third is that the pPRC was estimated with a finite data set, and where the response of the neuron is noisy, the error in the fit may be high. To address the effect of using a polynomial to fit the data, we fit a fourth-order polynomial, the same order used in the pPRC, to the iPRCs and plotted the resulting curves as dashed gray lines in Fig. 8. To address the effect of the finite stimulus waveform, we convolved the polynomial fit with the synaptic current used as perturbation in our experiments, plotted as solid gray lines in Fig. 8. Notice that the adjoint method is able to capture the biphasic PRC of the GACell. However, when the iPRC is convolved with the synaptic current, the negative component of the STA is blurred, as happened with our pPRC method.

## DISCUSSION

The goal of this article is to present a new approach to estimate PRCs from neurons as a function of the neuron's firing rate. The method used to predict the ISI has two stages: *1*) a recursive linear model is fit to predict the ISIs as a function of the applied current and its history, and *2*) a polynomial function is fit to predict how a synaptic input will perturb the neuron from the expected ISI given the phase of the synaptic input and the expected ISI. We have shown that this method can explain nearly 99% of the spike time variance over a range of firing rates in a computational model. We have also shown this method can be applied to real neuronal data obtained from pyramidal neurons in the hippocampal formation (CA1), explaining over 60% of the variance in over half the neurons we recorded. Coefficients estimated from the best seven neurons we recorded, presented in Fig. 6, are available for download from http://neuralnetoff.umn.edu/pPRC.html. A limitation of the method is that accurate models could not be generated in irregularly spiking neurons. This limitation will be addressed in future work.

Many different methods for measuring the PRC already exist, from the theoretical approaches by solving the adjoint equation (Ermentrout and Kopell 1991; Ermentrout 1996, 2002), new numerical methods for estimating the PRC (Govaerts and Sautois 2006), and ones that can be implemented experimentally, such as the direct method, or through a white noise stimulation (Ermentrout et al. 2007; Ota et al. 2009a, 2009b; Torben-Nielsen et al. 2010). Our approach is an extension of the direct method. A major advance of our approach is that instead of holding the neuron at a fixed firing interval over the duration of the experiment, we use an adaptive prediction algorithm to estimate the next spike time given the history of the neuron.

There are also many different functions that can be fit to the data to estimate the PRC. In the past we have used a polynomial function (Netoff et al. 2005). Again, here we chose a polynomial function because it could easily be generalized to two dimensions. Truncated Fourier series (Galán et al. 2005) also have been used for fitting PRC data. It is possible to extend the Truncated Fourier series to higher dimensions, but we found that the results were not significantly improved over the use of a polynomial fit (results not shown).

This work shows that spike times of neurons can be predicted with sufficient fidelity while input current is varied that it is possible to detect perturbations by synaptic-like inputs. This is the first time we have seen a model that takes into account the prior stimuli and ISIs to predict the next spike. Furthermore, using recursive estimations, we were able to make this model adaptive, which further increased the prediction accuracy as the dynamics of the neuron changed over the duration of the experiment.

These results are also the first experimental evidence measuring the PRCs from neurons as a function of firing rate and stimulus phase. The PRC model is extremely useful in predicting how networks of neurons will synchronize, but it is limited by the fact that it is usually estimated around a narrow set of parameters. This limits its applicability to predict network behaviors as parameters, such as firing rate changes. Several computational models have been used to study how synchrony changes as a function of firing rate. In addition, other groups have developed models showing the influence of firing rate and cell type on synchrony (Fink et al. 2011). This article supports modeling results predicting that PRCs will decrease in amplitude as firing rate increases when neurons are oscillating around 10 Hz, or close to the theta frequency. These predictions suggest that if coupled through excitatory connections, neurons will not synchronize as strongly at higher firing rates compared with slower firing rates. We have used this dependency of synchrony on firing rate to explain changes from asynchronous tonic phase to synchronous clonic phases in seizures (Wilson et al. 2011). Changes in synchrony as a function of firing rate also have been used to explain changes in gamma oscillations (Lewis and Rinzel 2003). The PRC estimated in these experiments can be directly incorporated into these kinds of models to explain changes in synchrony as dynamics of neurons change.

The nonlinear nature of neuronal responses to stimulus input led to some prediction errors by the linear ARX model. A quadratic ARX model decreases this error significantly for the GACell; however, we found no improvement when this was applied to the real neuron where the nonlinearities were obscured by the noise. Presumably, the noise in the real neuron is due to synaptic inputs and thermal noise. To fit the greater number of terms in the quadratic ARX model, larger data sets are required, and these may not be acquired in a reasonable amount of time. Because of this, we decided to only present the linear ARX model results in this article. However, if larger perturbations to the ISI were induced by using larger current changes, we may expect that a quadratic or even cubic polynomial might significantly improve prediction accuracy.

Our method works best where the change in frequency can be best predicted by using a quadratic model, but as nonlinearities accumulate, the prediction accuracy will degrade. The quality of the prediction depends also on the length of the history used. If the history is too short, effects from changes in current prior to the window of data used can effect the spike times but are not accounted for. If the history is too long, the model may try to utilize data that just add noise to the prediction. The optimal history length should be long enough to minimize the effect unmeasured variables, such as real synaptic inputs. At the other extreme, the cell's behavior is not stationary and has outliers; a long history can make the predictions less accurate by using data that are no longer predictive. Although this is a two-dimensional model, we found that we could get a good model by recording around 1,000 stimulations, about 10 min of recording. However, we also have obtained models where the number of spikes recorded was fewer than 400 stimuli.

The estimated PRCs we have observed in this study have the typical shape of the curves seen using traditional methods. The PRCs calculated from the computational model and the cortical neuron (Figs. 4*D* and 5*D*) show that the PRC skews toward the right and the peak goes higher as the mean ISI increases (firing rate decreases). This qualitative behavior is also reported by Gutkin et al. (2005) in a computational model of excitatory neurons. Gutkin et al. point out in their model that the displacement of the peak to the right as the ISI decreases is a consequence of slowly adapting potassium currents that cause the firing rate of the cell to change. The adaptation current is at its maximum shortly after an action potential and suppresses the effect of synaptic inputs and deactivates later in the ISI. As this current shunts the effect of the inputs, the maximum effect of a synaptic input will occur at the end of the ISI. As the ISI decreases, the adaptation current shunts inputs for a larger portion of the phase.

Some of the PRCs presented in Fig. 6 do not return to zero and therefore are not periodic, but in theory they should be. The polynomials fit were not constrained to be periodic. The differences can be due to errors in fitting the data at the beginning of the phase, where the phase advance is noisy. They also can be aperiodic because the effect of the stimulus can last beyond the period to which it is applied. In these cases, the first-order PRC, the effects seen in the cycle in which the stimulus is applied, is not sufficient to explain the full dynamics of the neuron, and the second (the effects on the next cycle) may need to be included to get a periodic solution.

In this study, we estimated the PRC by fitting a polynomial to the data by a least-squares prediction error. We found that using a TSVD in finding the least-squared solution dramatically reduced the noise in the fit. Alternative methods have been used to fit models to data that may be applicable, such as a Bayesian approach (Ota et al. 2009b).

In fitting PRCs, higher order effects, i.e., influences from inputs applied on previous cycles, can influence the measure. Although the largest effect is generally in the period that the stimulus is applied, prediction accuracy of network activity can be improved by including the effects of the higher order PRCs (Achuthan and Canavier 2009; Oprisan and Canavier 2001). To minimize these higher order influences while measuring PRCs in the past, we have applied a synaptic current every sixth spike. Another approach is to saturate the second-order effects by stimulating the neuron on every cycle at the same time lag and measure the resulting steady-state ISI; the resulting PRC has been called a functional PRC (fPRC) (Cui et al. 2009). An advantage of our approach over use of the fPRC is that the polynomial used in our model can easily be extended to account for the higher order effects in a single model by adding more terms including the historical phases and ISIs.

There are differences between the PRCs resulting from the use of our method and the iPRCs measured using the adjoint method from the GACell. To make the iPRC comparable to the functions resulting in pPRC, we fit the iPRC with a polynomial and then convolved the resulting function with the synaptic current waveform. The convolution removes the phase delay (negative section of the PRC) observed at early phases using the iPRC. Because of the relatively low numbers of samples and the noisy measures, our method did not have the fidelity to capture the small negative phase at the action potential seen in the iPRC. Furthermore, the peaks of the pPRCs are shifted to the left with respect to those of the convolved iPRCs; we attribute the differences to the noise in the measurements in the early phase and the finite data lengths used to estimate the pPRCs. However, the change in the PRC shape predicted by the adjoint method is followed by our proposed method.

This method can also be extended to relate the STA to other variables besides the firing rate, such as the concentration of a neuromodulator. It is also possible to extend this formalism to relate the STA to three or more variables. For example, one option for a third variable could be measuring the effect of stimulus intensity in addition to the phase and expected unperturbed ISI of the neuron.

This method enables us to more accurately model the effect of antiepileptic drugs over a wide range of neuronal behaviors.

### Conclusion

This article presents new experimental methods for estimating PRCs as a function of the stimulus phase and firing rate of the neuron, which we term a parameterized PRC (pPRC). The advantage of this pPRC is that it is not limited to modeling the neuron at one fixed condition, such as ISI, but can be used to model a neuron's dynamics over a range of parameters. This model can be used for large-scale and efficient models where ISIs may vary across neurons and/or across time. These models may be used to generate efficient network simulations. The experiments demonstrate that neuron spike times can be varied considerably by applying current steps but still can be predictable given the history of the inputs and neuron's spikes. The experimentally measured pPRCs in this study support the hypothesis generated in many modeling studies that the phase response curve is smaller in amplitude at higher firing rates, predicting that neurons coupled through excitatory synapses will not synchronize as well when firing rapidly as they will when they are firing slowly. One advantage of this approach of the method used to estimate the pPRCs is that it can easily be extended to account for other perturbations to the neuron or to account for more history and higher order nonlinear responses. The use of TSVD in solving the coefficients for the polynomials increases the robustness of the estimates under noisy and sparse estimates.

## GRANTS

We acknowledge funding from Tecnológico de Monterrey, a Fellowship of Design of Medical Devices, University of Minnesota; a Grant-in-Aid from the University of Minnesota; the laboratory of Dr. Damien Fair (Oregon Health & Science University), and National Science Foundation Career Award 0954797 (to T. Netoff).

## DISCLOSURES

No conflicts of interest, financial or otherwise, are declared by the authors.

## ENDNOTE

At the request of the authors, readers are herein alerted to the fact that additional materials related to this manuscript may be found at the institutional website of one of the authors, which at the time of publication they indicate is http://neuralnetoff.umn.edu/pPRC.html. These materials are not a part of this manuscript and have not undergone peer review by the American Physiological Society (APS). APS and the journal editors take no responsibility for these materials, for the website address, or for any links to or from it.

## AUTHOR CONTRIBUTIONS

O.M.-D. and T.I.N. conception and design of research; O.M.-D. performed experiments; O.M.-D. and T.I.N. analyzed data; O.M.-D. and T.I.N. interpreted results of experiments; O.M.-D. prepared figures; O.M.-D. and T.I.N. drafted manuscript; O.M.-D. and T.I.N. edited and revised manuscript; O.M.-D. and T.I.N. approved final version of manuscript.

- Copyright © 2013 the American Physiological Society