|Home | About | Journals | Submit | Contact Us | Français|
We demonstrate that two key theoretical objects used widely in Computational Neuroscience, the phase-resetting curve (PRC) from dynamics and the spike triggered average (STA) from statistical analysis, are closely related when neurons fire in a nearly regular manner and the stimulus is sufficiently small. We prove that the STA due to injected noisy current is proportional to the derivative of the PRC. We compare these analytic results to numerical calculations for the Hodgkin-Huxley neuron and we apply the method to neurons in the olfactory bulb of mice. This observation allows us to relate the stimulus-response properties of a neuron to its dynamics, bridging the gap between dynamical and information theoretic approaches to understanding brain computations and facilitating the interpretation of changes in channels and other cellular properties as influencing the representation of stimuli.
Dynamical systems models like the Hodgkin Huxley model of the squid giant axon are very effective at replicating the firing properties of individual neurons. Such models have been extremely useful tools for understanding the mechanisms of neural excitability and for simulating neural circuits. However, these models have been less successful in aiding the understanding of what neurons compute. Given a model that perfectly describes the behavior of a neuron or a network, we are in most cases still at a loss to say what computation this neuron or circuit is performing. Rather, the natural concepts and objects of the theory of computation (e.g. stimuli, features, coding) seem only distantly related to those of the theory of dynamical systems (e.g. differential equations and attractors). Here we relate two mathematical objects central to these two approaches: the phase resetting curve (PRC) and the spike triggered average (STA). We start by introducing these objects and then prove that the STA is proportional to the derivative of the PRC in the weak stimulus limit. We show that this approach allows us to efficiently and accurately compute the STA from the PRC and vice versa, in the case of numerical simulations as well as in the case of real neurons. The ability to compute these functions from each other allows us to make some progress in relating dynamics of neurons to their ability to code features.
The PRC describes how the spiking of a regularly firing neuron is altered by incoming input, that is, how the time of the next spike is shifted as a function of the stimulus time relative to the previous spike: Δ(t) (T - (t))/K where T is the natural period and (t) is the time of the spike given a stimulus at time t after the last spike. The constant K is the magnitude of the stimulus, e.g. 1 mV since the PRC depends on the amplitude of the perturbation. [1, 2] show that for small stimuli, x(t), any stable limit cycle oscillator can be reduced to a scalar model for its phase, θ:
where θ [0, T) and Δ(θ) is the PRC. For neurons, the stimulus has dimensions of millivolts per millisecond. We take θ = 0 to be the time of spiking. The PRC is valid for neurons which are repetively firing (that is, on a stable limit cycle), but the concept of the PRC can be applied even when the neurons are quite noisy, often permitting a PRC to be measured experimentally from neurons in vitro [3, 4]. Δ(θ) is readily computed for differential equation models, and can also be computed for neural models either via direct perturbation- or indirectly [3, 4]. Equation (1) shows that for a regularly firing neuron, the PRC provides, an answer to the question of “how will a stimulus influence when the next spike will come?”
A contrasting approach to understanding neural computation has focused on neural coding, by which we mean determining what features of stimuli are represented by single spikes and spike trains . Such analysis of the stimulus dependence of spike trains often includes the calculation of the STA, which is related to the reverse correlation [8, 9]. The STA is defined as the average stimulus with a given prior statistics preceding an action potential in a neuron. For our purposes, by stimulus we refer to the current injected into a neuron (divided by the capacitance). In other experimental protocols, the stimulus may be a sensory stimulus presented to an animal while a neuron is recorded. If x(t) is the stimulus:
where the average is taken over all spike times, τj. In the present context, the STA is an answer to the question of “on average, what temporal features of the stimulus lead to spiking?”
To derive the main result which applies when the stimulus is small, we write x(t) = σξ(t) where σ is the magnitude of the signal and ξ(t) is a zero-mean stationary stochastic process with unit variance. (Note that for white noise applied to the voltage equation for a neuron model, σ2 has dimensions of mV2/ms.) We will use the smallness of σ to estimate the time of a spike in equation (1) and from this obtain the STA from (2). With θ(0) = 0 as the initial condition (we assume the neuron has just spiked) we write
which upon integration yields
We want to determine the time τ at which the oscillator spikes again, that is θ(τ) = T. As with θ, τ depends on σ, so we write τ = τ0 + στ1 + . . ., and find τ0 = T and
Thus, we find that to order σ,
where C(t) = ξ(t)ξ(0) is the auto correlation of the stimulus. To get the fourth line, we note that the expected value of ξ(t) is zero; to get the last line, we have integrated by parts. Dropping the higher order terms and assuming the stimulus is white noise, we obtain the main result:
where we use the T-periodicity of Δ(t) to drop the T. This result shows that the dynamics of a neuron, as captured by the PRC, can be used to predict the STA and conversely, given the STA, we can estimate the PRC of a repetitively firing neuron. Furthermore, (4) is valid for any small, zero-mean, stationary, stochastic stimulus.
Two key issues arise in this analysis. First, we must translate the periodic PRC to the generally aperiodic STA. This conversion is most straightforward when neuron is firing at a nearly constant rate and the PRC is well-defined. While the STA is in principle aperiodic, in reality it is only sensible to define the STA over the time interval prior to a spike in which there are no other spikes. Thus the time over which both the STA and the PRC can be clearly defined is the interval between spikes, i.e. the average period. Second, we must note that because the PRC is the integral of the STA, it is defined only up to an additive constant term. However, if we assume that the PRC vanishes at t = 0, T (as is common in neurons), then we can determine the integration constant.
To test the theory, we compute the STA for the Hodgkin-Huxley (HH) equations and then use equation (5) to compute the PRC subject to the constraints that it vanishes at t = 0 and at the average period, t = T. We drive the four variable biophysical Hodgkin-Huxley model with a constant bias current (I = 10μA/cm2)to make it fire at 70 Hz, inject noise and compute the spike triggered average. (7000 spikes were used.) We then numerically integrate the STA and time-reverse it to reconstruct an approximate PRC. We compute the exact PRC using the method in  and compare the two methods for two different amplitudes of noise. Figure 1 shows that in both the low and high noise case, the PRCs calculated from the STAs are almost identical to the actual PRC. Later (see figure 3), we will systematically quantify the dependence of the reconstruction on the statistics of the spike trains for both the HH and phase models.
We next tested this transformation on real neurons. We performed whole cell recordings from olfactory bulb mitral cells. We injected these cells with DC current to cause them to fire repetitively at 50 ± 6 Hz, added noisy current (with an amplitude 10% of the DC) and recorded spike times over intervals of 2-2.5 seconds, repeated 100-120 times. From the recorded spike trains we first eliminated the initial period (250-600 ms) of spike frequency accommodation and then calculated the STA. We calculated the PRC from the STA as per eq. (5). For comparison, we also calculated the PRC using a method based on injection of aperiodic perturbing pulses . Figure 2 shows the estimated PRC obtained from these two methods from the an olfactory bulb mitral cell. In the PRCs estimated by both methods, there is a clear negative region after the spike followed by a larger positive region, consistent with our earlier estimation of the PRC for olfactory bulb mitral cells. We also were interested in the fact that the STA-based method was possible despite the fact that cells were not firing in a precisely oscillatory manner. We measured the standard deviation of the interspike intervals in our recordings and found it to be approximately 10% of the firing rate. Thus our method is robust for at least this level of variation in the firing rate. We also applied these methods to recordings of neurons in the mouse somatosensory cortex (data not shown) with similar results.
To explore further how this relationship between the PRC and the STA depends on the regularity of the periodic firing, we drive both the HH model and simple phase models of the form of (1) with larger and larger amplitude noise to make their firing more irregular. We examined the effects of noise on phase models with two commonly used PRCs (Δ(t) = sin t, Figure 3A1; Δ(t) = 1 - cos t, Figure 3A2). Calculating the STA and the PRC for these higher noise simulations (after correcting for the change in average firing rate) resulted in PRCs that less and less closely resembled the actual PRC for these models although the general shape of the PRC was somewhat maintained (Figure 3A1 and A2). We injected increasingly larger noise (σ = 0.25 - 16) into the HH model. While the shape of the PRC is degraded, the zero crossing is preserved remarkably well. We quantified this degradation in the quality of the estimated PRC by calculating the correlation coefficient (R) between the actual PRC and the PRC estimated from the STA and plotted this against the CV of the ISI distribution. We observed that for both phase oscillator models the estimate provided a good approximation (R>0.75) of the PRC up to CV = 0.4. For CVs above 0.4, the correlation between the actual and estimated PRC declined rapidly for both models. For the HH model, it was difficult (even with the large noise applied) to increase the CV beyond about 0.3, when in the oscillatory mode. Trying larger noise values, led to numerical difficulties. We remark that in figure 3, we have normalized the magnitude of the PRC to a maximum of 1 and that the un-normalized reconstructed PRCs for the HH model had amplitudes that were smaller at large noise values.
Relating neural dynamics to neural coding has been termed a grand challenge for neuroscience and we believe that our work describes an important step towards relating these two subjects, albeit for a restricted set of stimuli and neural operating range. Strengthening this connection will provide valuable means of relating the wealth of data on biophysical properties of neurons (as captured in dynamical systems models of neuronal properties) to questions about the properties of stimuli that are being computed by neurons. There are two limitations to our approach. However, our methods can still be applied to many situations. First, the PRC, as useful as it is, has limits to its applicability; specifically, in this context, it is defined only for nearly periodically firing neurons where only the timing of spikes (and not the rate) are altered. Nonetheless, neurons in many brain networks are active spontaneously and the strength of any one input is weak. Thus, the assumption that inputs modulate the timing of spikes rather than adding more spikes may hold. For example,  demonstrated that realistic synaptic conductances in the aplysia satisfy the mathematical criteria of “weak coupling” in the sense that the notion of phase still makes sense.
Secondly, traditional STA analysis refers to sensory stimuli rather than directly injected neurons. However, a number of authors have used current injections to understant the “neural code.”  were among the first to try to extract statistical information such as the STA and the spike triggered covariance (STC) from a biophysical model for a neuron. Using simulations of the Hodgkin-Huxley equations, they compute both the STA and the STC in a situation where the stimulus is sufficiently strong to elicit spikes. They show that there are important features of the stimulus encoded in the STC which are not evident in the STA. The calculations in section II could be extended to incorporate higher-order effects and thus relate the STC to other aspects of the PRC .  also consider the effects of the amplitude of the stimulus on the STA; such effects would be reflected in changing the shape of the reconstructed PRC (as in figure 3A3). While our theory is linear (with respect to the effects of the stimulus), it is known that the shape of the PRC is also affected by stimuli which are sufficiently large , so that some of the shape effects in figure 3 may be due to pushing the simulations beyond the linear range. More recently,  derived the STA for the integrate-and-fire model and produced formulae in the limit of small noise.  characterized the spike triggered average voltage (a related quantity which requires intracellular recording of the membrane voltage fluctuations) for several different models and related it to the underlying dynamics while  studied the relationship between the STA and STC to subthreshold dynamics.  compute the PSTH and STA of a general class of spike-response models, thus providing some insight into the relationships between dynamics of spiking and the neural code. This work is the closest to ours in its generality. They apply their results to the data in . Our methods, while in a more restricted situation (regularly firing) are very general in that they apply to any model of an experimental system for which a PRC is defined. In particular, PRCs encode dynamic information about subthreshold behavior [3, 20], which is not possible in spike-response models. The PRC can be directly computed from any biophysical model for periodic neuronal firing.
The authors were supported by NSF grants, DC005798, MH079504 and DMS0513500 and they would like to thank Brent Doiron and Nicolas Fourcaud-Trocme for their comments on this manuscript.
PACS numbers: 87.18.Sn, 05.45.Xt, 84.35.+i, 89.75.Hc