Relating Neural Dynamics to Neural Coding

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 under a wide range of stimulus conditions. We prove that the STA 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.


I. INTRODUCTION
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, such models, though they are themselves computational, 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.

II. PRC AND STA.
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: is the time of the spike given a stimulus at time t after the last spike. [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. 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 [3,4]. ∆(θ) is readily computed for differential equation models, and can also be computed for neural models either via direct perturbation[5]- [7] 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 [8]. Such analysis of the stimulus dependence of spike trains often includes the calculation of the STA, which is related to the reverse correlation [9]. The STA is defined as the average stimulus preceding an action potential in a neuron. For our purposes, by stimulus we refer to the current injected into a neuron. In other cases, the stimulus may be a sensory stimulus presented to an anaesthetized animal. Thus, if x(t) is the stimulus: where the average is taken over all spike times, T j . The STA is an answer to the question of "what temporal features of the stimulus lead to spiking?"

III. DERIVATION OF THE MAIN RESULT.
To derive the main result, we use two previous results. A standard result in the theory of neural decoding [9,10,11] shows that the optimal linear approximation to the firing rate, r(t) of a neuron to a stimulus, I(t) is given by: where r 0 is the background firing rate. For stimuli which are brief, the response, r(t) is called the post-stimulus time histogram (PSTH) by experimentalists. For a δ(t) function stimulus, the PSTH is exactly the STA once the background rate, r 0 is removed: [3] relate the response of an oscillating neuron to a brief stimulus and the PRC as: 3 where G(s) = s − T ∆(T − s). If the signal is small then we can approximate the inverse of Differentiating (4), using approximation (5), and noting that r 0 = 1/T , we get STA(t) = −∆ ′ (T −t) for a unit impulse stimulus, I(t) = δ(t). For more general stimuli with C(t, s) =< I(t)I(s) >, we obtain: which in the case of white noise with strength σ 2 becomes 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.
Two key issues arise in this analysis. First, we must translate the PRC, which is a function of the phase of the oscillating neuron and thus periodic, to the STA, which is a function of the time since the last spike and thus aperiodic. This conversion is most straightforward in the case in which the neuron is firing at a nearly constant rate which is the condition in which the PRC was originally defined. In this case the phase of the oscillator can be translated to the time since the last spike simply by multiplying by the period of the oscillation. 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 a 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.

IV. EXAMPLES
To test the theory, we compute the STA for the Hodgkin-Huxley equations and then use equation (6)  and compute the spike triggered average. We then numerically integrate the STA and timereverse it to reconstruct an approximate PRC. We compute the exact PRC using the adjoint method [12]- [14] 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. Thus, even when the amplitude of the noise was sufficient to cause some spikes to drop out, the approximation (6) still results in a very accurate estimate of the PRC. Because the STA is robust to independent sources of noise, so is the reconstruction of the PRC from the STA.
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 (eq. (2)). We calculated the PRC from the STA as per eq. (6). For comparison, we also calculated the PRC using a method based on injection of aperiodic perturbing pulses [4]. Figure 2 shows the estimated PRC obtained from these two methods from the an olfactory bulb mitral cells. In the PRCs estimated by both using two different methods from recordings of an olfactory bulb mitral cell. The gray line shows estimate from method described previously [4]. Black line shows average PRC for the same cell calculated from the STA, as described above methods, there is a substantial 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. Thus, 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 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 3a; Figure 3b). Increasing the amplitude of the injected noise resulted in, first, a broadening and then a leftward shift in the shape of the interspike-interval (ISI) histograms (Figure 3a1 and a2). 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 3b1 and b2). We quantified this degradation in the quality of the estimated PRC by calculating the correlation coefficient (R) between the 6 actual PRC and the PRC estimated from the STA and plotted this against the CV of the ISI distribution. We observed that for both 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.

V. DISCUSSION.
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 cases. 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. The main limitation to our approach, which still allows it to be applied to many situations, arises from several assumptions that we have made and from the fact that the two mathematical objects that we have linked, the STA and the PRC, as useful as they are, have limits to their applicability. First, the STA provides only a first order approximation of the response of a neuron to stimulus. That is, it does not account for the fact that spikes interact with each other. Thus any analysis of neuronal dynamics that allows for estimation only of the STA is incomplete. Calculations of at least the spike triggered covariance are required to provide a more accurate description of the coding properties of many real neurons. Similarly, the definition of the PRC is based on the case of a regularly firing neuron receiving a weak input i.e. an input that does not change the neurons average firing rate but only the timing of spikes, although extensions to examples involving stronger stimuli have often proved fruitful in many contexts. [15] were among the first to try to extract statistical information such as the STA and the STC from a biophysical model for a neuron. Using simulations of the Hodgkin-Huxley equations, they compute both the STA and the spike triggered covariance (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. More recently, [16]derived the STA for the integrate-and-fire model and produced formulae in the limit of small noise. [17] 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 show PRCs estimated from the STA as described in the main text for different noise levels (as above). In both cases low noise levels (orange traces) produce most precise estimates of the PRCs.
(a3, b3) Plot of correlation between actual and estimated PRCs vs CV of the ISI distribution for the same date as plotted above. Note that correlation is good up to a CV of 0.4.
to the underlying dynamics. [18] compute the PSTH and STA of a general class of spikeresponse 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 [11]. 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. Further, there is a direct connections between the underlying dynamics of regularly firing neurons and their PRC [3,6,19].
In the context of direct sensory inputs to a resting neuron, it is clear that strong stimuli are necessary to evoke responses. However, 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. In a recent paper Preyer et al demonstrated that realistic synaptic conductances in the aplysia satisfy the mathematical criteria of "weak coupling." Mathematically, "weak" means that the inputs are sufficiently small so that the notion of phase persists. In [3], we showed that the idea of a PRC still works very well even when the system is nonstationary and the ISIs are slowly varying. In this paper, we have described a relationship between the first order statistics (spike triggered average) and the PRC. Because dynamical systems and codingbased approaches are both concerned with the fundamental question of what makes a neuron spike we anticipate that further analysis will allow this approach to be extended to cover situations in which the average firing rate of a neuron is rapidly varying.