## Abstract

Although many studies in neuroscience are based on comparing neuronal responses to single, isolated sensory or motor events, multiple events frequently occur in close temporal proximity in freely moving animals. This often obscures the precise temporal correlation between each event and the relevant brain activity. By simulating neuronal responses in multi-event tasks, we show that perievent time histograms (PETHs) greatly distort the underlying true responses. We propose a multi-event deconvolution method that can separate the contribution of each event to the overall neuronal activity. The improvements over PETH in analyzing real data are demonstrated using simulated data and a sample electrophysiological recording obtained from rats in a task involving responses to a reward predictive cue.

## INTRODUCTION

An important goal in neuroscience is to understand how brain activity generates sensation and behavior. A major step toward achieving this goal is to determine how changes in neuronal activity patterns correlate with transient stimuli or motor actions. In many situations, several events occur closely in time, making it difficult to determine to which event an observed change in neuronal activity is correlated. In some situations, this problem can be addressed by experimentally separating the events in time or by using experimental conditions where only one event is present. However, in many other situations, the behavioral response to one event as well as its related activity is contingent on other events and therefore cannot be separated experimentally. (Hollerman et al. 1998; Lau and Glimcher 2008; Lauwereyns et al. 2002; Nicola et al. 2004b; Samejima et al. 2005). For instance, to study how the brain translates sensory information into immediate actions, one has to determine (in each relevant brain region) whether the neuronal activity between the stimulus and the action codes for the sensory stimulus or the action that follows it. It is also possible that in some regions the observed activity specifically encodes some unique function of the relation of the two events.

The ambiguities that arise from response overlap are known in the field of electrophysiology and some methods have been put forward to address it (Jin et al. 2009; Lau and Glimcher 2008; Nicola et al. 2004a; Taha et al. 2007). Most of these methods assume that if the average neural activity is more correlated to one event than another, the response to that event should be of shorter duration and larger magnitude compared with the response to the other event. Some have extended this concept by comparing the rise time of the responses to each event to determine the most correlated event (Taha et al. 2007). Another group has used a shuffling of event timings between trials for each event to separate significant changes in the firing patterns that indicate response correlation with one event compared with others (Jin et al. 2009). Although intuitively accurate, in many cases, these methods can only report which of the events correlates better with the neural response. They cannot rule out the dependency of the response to the event with the lower correlation nor can they inform the precise relation of each event to the neuronal response. Here we develop a method to dissociate these activities and test it on extracellular recording of single neurons obtained in behaving rats. Importantly, this method could in theory apply to any type of time varying brain activity such as voltammetry, event-related potentials or functional-MRI.

## METHODS

### Deconvolution model

The deconvolution model assumes that the total firing rate of a neuron in each trial is equal to the linear sum of the contributions of each event-related firing (ERF) delayed by the event latencies in that trial. It further assumes that the shapes of the ERFs are independent of inter-event intervals and of time. This means that properties of ERFs such as duration or amplitude do not change by changing the spacing between events and are also stationary over trials. Given these assumptions, the firing of a neuron in trial *m* with a total of **K** event types and **T t**ime samples can be written as
*f*_{m}(*t*) represents the firing of the neuron in trial *m* as a function of time and *h*_{k}(*t*) is the neuron's ERF for event **k.** The term *z*_{mk}(*t*) is set to one when event *k* occurs and is zero otherwise. Furthermore η_{m}(*t*) represents the noise in trial *m* and can model both the inherent randomness of neuronal firing as well as the contributions of all other potential events in the task that are not explicitly included in the model.

If the firing rates and noise terms are collected over time and over all the **M** trials in column vectors *F* = [*f*_{1}(1)…*f*_{1}(*T*)…*f*_{M}(1)…*f*_{M}(*T*)]^{T} and *N* = [η_{1}(1)…η_{1}(*T*)…η_{M}(1)…η_{M}(*T*)]^{T} and the ERFs to be estimated are also collected into another column vector *H* = [*h*_{1}(1)…*h*_{1}(*T*)…*h*_{K}(1)…*h*_{K}(*T*)]^{T} and correspondingly the event timings *z*_{mk} are turned into a **MTxKT** matrix referred to as *Z*, the convolution equation can be written as linear regression
*Equation 2* shows the full model for the convolution of ERFs over all the trials for a given neuron. The goal of deconvolution is to find an unbiased estimator *Ĥ* to approximate the true ERFs in *H* given that the firing of a neuron in all trials (a.k.a *F*) as well as the timing of all events (a.k.a *Z*) in each trial are known.

Raw PETHs aligned to each event can be derived from *Eq. 2* by multiplying both sides of the equation with *Z*^{T} and normalizing (dividing) by the number of trials.
*H̄* = [*h̄*_{1}(1)…*h̄*_{1}(*T*)…*h̄*_{K}(1)…*h̄*_{K}(*T*)] is a **KTx1** column vector containing the raw PETHs for all events. *V* called hereafter the convolution matrix is defined as *Z*^{T}*Z* and is a **KTxKT** square matrix that is mapping the **K** ERFs to the corresponding **K** PETHs.

*Equation 3* shows that PETHs are in general themselves only a convoluted version of true ERFs. The further the convolution matrix (1/*M*)*V* is from identity matrix the more distorted *H̄* will be from *H*. This in turn will depend on the distribution of event timings in each trial. *Equation 3* not only confirms our main argument that PETHs can be very different from the actual neuronal ERFs but also provides a way to formulate this difference simply as a multiplication between a matrix containing all the event timings and the ERFs. To find the best estimate for the ERFs, one needs to undo all the time shifts that are applied to *H* by the convolution matrix (1/*M*)*V*. This can be done either iteratively (see *Eq. 6A*) or in a single pass using the least squares technique (see *Eqs. 8–10*). For computational reasons, the least squares technique is more efficiently implemented in the frequency domain, whereas for the iterative approach, we find it more intuitive to stay in the time domain. However, the iterative approach and the least square approach are theoretically applicable in either domain.

### Iterative deconvolution in time

To undo the effect of convolution matrix *V* on the ERFs, an inverse matrix *V*^{−1} needs to found. It can be shown that for an invertible and properly scaled matrix **B**, **B**^{−1} can be calculated using the Neumann series (Stewart 1998)
*V,* which is equal to *Z*^{T}*Z* is a symmetric positive semi-definite matrix and when *V*^{−1} exists, it can almost always be properly scaled for the Neumann series to converge. To properly scale the convolution matrix(1/*M*)*V* in *Eq. 3* one can divide it by any real positive number *S*> (λ_{max}/2) where λ_{max} is the largest eigenvalue of (1**/***M*)*V*.

The normalized PETH formed can be represented as
*H̄*_{/S}(*t*) to represent the normalized vector of raw PETHs by *S* . The inverse of convolution matrix (1/*MS*)*V* can then be found using *Eq. 4*. The result can be convolved with both sides of *Eq. 5* after which we have

*Equation 6A* shows that *Ĥ* can be constructed by adding the terms of the Neumann series in each iteration with superscripts on *Ĥ* representing the iteration number and the zero iteration being the normalized PETHs themselves. *Equation 6B* shows an estimate of the deconvolution error in predicting the true ERFs and its dependence on the trial noise terms.

### Complete deconvolution in frequency

Given the fact that convolution in time domain is equivalent to multiplication in frequency domain (Oppenheim et al. 1997), we can rewrite *Eq. 7* this time in the frequency domain, we have

One of the benefits of representing the deconvolution problem in the frequency domain is that as can be seen in *Eq. 7* the value of *f*_{m}(ω) at ω = ω_{0} depends on the values of *z*_{mk}(ω), *h*_{k}(ω), and η_{m}(ω) only at ω = ω_{0} and as such is decoupled from all other frequencies and can be studied independently from them. Therefore if the firing rates and noise terms over all the **M** trials are collected for each frequency in column vectors *F*(ω) = [*f*_{1}(ω)…*f*_{M}(ω)]^{T} and *N*(ω) = [η_{1}(ω)…η_{M}(ω)]^{T} and the ERFs to be estimated are also collected into another column vector for each frequency *H*(ω) = [*h*_{1}(ω)…*h*_{K}(ω)]^{T} and correspondingly the event timings *z*_{mk} are turned into a **MxK** matrix referred to as *Z*(ω), the convolution equation can again be written as linear regression this time in the frequency domain
*Eq. 8* by
*V*(ω) is a **KxK** the convolution matrix and is equal to *Z*(−ω)^{T}*Z*(ω). To deconvolve, we need to find the inverse of (1/*M*)*V*(ω) at each frequency. Multiplying both sides of *Eq. 9* with this inverse results in the deconvolved ERF vector in the frequency domain

^{2} is the noise variance.

*Equations 6* and *10* should yield equivalent results for deconvolving the ERFs. In theory the frequency domain approach is more efficient for complete deconvolution due to decoupling the convolution effect on each frequency. However, in practice due to singularities of the convolution matrix at some frequencies or violations of the model assumptions, the partial deconvolution in time domain might provide more accurate results.

The deconvolution requires the convolution matrix to be invertible whether in time or in frequency domain. This means that the data should have enough number of trials with some minimum amount of jitter in the event timings so as to provide the deconvolution model with enough power to parse out ERFs over all frequencies. This is especially true for the lower frequencies where larger jitters are needed to reveal their amplitude changes over time. In fact in the case of single occurrence of each event in a given trial it can be shown that the DC mapping of the ERFs to the PETHs is not invertible because *V*(ω = 0) is a **KxK** matrix with all its elements set to one. In such cases where deconvolution is not well defined for low frequencies, the deconvolved ERFs may show spurious up- or downward DC shifts that can be corrected back to zero. Having a variable repetition of at least some events during each trial can resolve the ambiguity of the ERFs DC content because the changing number of events in each trial can provide another dimension of information beside the temporal jitters for the deconvolution of the ERFs.

### Cross-validation

To estimate the ERF estimation error in the absence of the ERFs across-validation technique can be used. In cross-validation, we divide the trials into two groups, one that usually includes the majority of trials is used for training and one with the remaining trials for validation of the fit. We have
*t* represents the training trials and subscript *v* represents the validation trials. The cross-validation MSE for iteration *n* is defined as
*Eqs. 12* and *13* the cross-validation error is similar to ERF error except for an added noise term and scaling by *V*_{v}. It can be shown that if ERF estimation error goes down to zero then the cross-validation error goes down to its lower bound set by noise.

Alternatively we can divide the data randomly into two equal halves and use the ERFs estimated by one part and the plug in principle to try to directly estimate the ERF estimation error as the following
*Eq. 13.*

### Behavior and electrophysiology

Rats were trained with reward predictive auditory cues (≤10 s) randomly presented on a variable interval schedule with an average interval of 60 s. A cartoon of the experimental apparatus is shown in Fig. 5*A* (For a detailed description of the ex perimental paradigm, see Ambroggi et al. 2008; Nicola et al. 2004b). Active lever presses only during presentation of the cue resulted in the termination of the cue and delivery of a 10% sucrose reward into an adjacent reward receptacle. The timing of cue presentations, lever presses, and the entries into and exits out of the reward receptacle were all simultaneously recorded together with the firing of ventral striatal neurons using chronically implanted micro-array electrodes. A subset of randomly chosen neurons (*n* = 106) from a previous study (Ambroggi et al. 2008) is used here.

### Simulation and data analysis

All simulations and data analyses were done in MATLAB 7.9.0. The algorithms used for simulation and analysis are available on request.

## RESULTS

### Deconvolution intuition

In electrophysiology, the conventional method to estimate a neuron's response to a stimulus is to shift and time-lock its spike train to that stimulus and then make a PETH by averaging the time-locked activity over many trials. An example PETH for an actual neuron is shown in Fig. 1. The firing of this neuron is time-locked to both a sensory cue and the motor action that follows it (Fig. 1*A*). The PETHs in Fig. 1*B* show a broad excitation after the cue and a sharp peak of activity around the motor action. However, the firing rasters show that the activity of this neuron is more tightly locked to the motor action compared with the cue, but because the motor actions happen in temporally close yet variable intervals after the cue a broad excitation appears shortly in the cue PETH following cue onset. Therefore based on the PETH results, one could hypothesize that this neuron is excited to the cue presentation as well as the motor action that can lead to an error in response classification.

Our proposed deconvolution method can be used to address this issue by teasing apart the motor activity from the cue activity as the following. If the neuron is only responsive to the action and not the cue, one should be able to account for the cue PETH by using the action PETH in each trial shifted to the time of the action in that trial and then adding the shifted versions of the action PETH over all the trials as shown in Fig. 1*C*. A simple overlay of the predicted cue PETH constructed by this technique on the raw cue PETH shows a good degree of overlap between the two responses (Fig. 1*D*, *left*), indicating that the broad excitation appearing on the cue PETH is likely the result of the neuron's response to the action shifted over many trials. Subtracting the predicted cue PETH from the raw cue PETH then successfully eliminates the influence of motor action on cue related firing (Fig. 1*E*, *left*, baseline is added back for comparison). On the other hand because most of the neuron's response is due to the action and not the cue, the same procedure performed on the raw motor PETH (subtracting the predicted PETH made by shifting the cue PETH in each trial from the motor PETH) does not account for the activity of the neuron around the motor action (Fig. 1, *C–E, right*; Fig. 1*E* baseline is added back for comparison). Therefore this neuron is more accurately categorized as action-encoding but not cue-encoding after we account for the response to the motor action in the raw cue PETH.

In general if we define an event-related firing (ERF) to be the true firing of a neuron to one event independent of other event timings, then it follows that PETHs themselves are a superimposition of the ERFs that are time shifted and mixed based on the relative timings of the events. In a linear and time independent regime, PETHs represent a convolution of the ERFs with the event timings. To recover the ERFs, one needs to reverse the convolution or in other words deconvolve the PETHs. This can be done either by trying to solve the inverse problem of convolution in one step (complete deconvolution) or iteratively where the current ERF estimates in each step are used together with the event timings to improve the ERF estimates in the next step (partial deconvolution). The general formulation for the deconvolution is laid out in methods where the procedure presented in the preceding text in Fig. 1 is shown to constitute only the first iteration in using raw PETHs and timing of all events repeatedly for partial deconvolution of the ERFs. This iterative approach is carried out in the time domain and will be referred to as deconvolution in time or temporal deconvolution (see methods, *Eq. 6*). The complete deconvolution that aims to solve the deconvolution in one step through a direct matrix inversion is carried out in the frequency domain and will be referred to as deconvolution in frequency or frequency deconvolution (see methods, *Eq.10*). It is important to note that due to the correspondence of time and frequency domain the choice of performing partial deconvolution in time and complete deconvolution in frequency is rather arbitrary and therefore the true comparison in the following sections is between partial and complete deconvolutions not between time and frequency domains per se. The following sections will examine the deconvolution results first for simulated data and then for actual single unit recordings.

### Simulation results

In this section, we first show how response overlaps can distort ERFs obtained for an event by simulating two hypothetical experiments. One of the examples involves multiple events happening closely in time and another involves high-frequency repetitions of a single event. We show that in both cases deconvolution in time or in frequency can recover the ERFs successfully.

##### MULTI-EVENT DECONVOLUTION.

For the first experiment, we assume to have four events A–D. For simplicity we assume that each event happens only once during each trial and always in the order ABCD with a total of 40 trials. The inter-event intervals between events in a given trial are drawn randomly from Poisson distributions. We assume to have a neuron with ERFs that show excitation to event A, excitation then inhibition to event B, inhibition to event C, and inhibition then excitation to event D (Fig. 2*C*, dotted lines). The firing of this neuron in two different trials in the presence of all four events is shown in Fig. 2*A*. The PETHs are obtained by convolving the ERFs with the event timings (*Eq. 1* with Gaussian noise, SD = 1). Fig. 2*C* shows the event-locked PETHs overlaid on their corresponding ERFs. This figure clearly demonstrates that the raw PETHs built by simply averaging the firing over trials can deviate considerably from the true ERFs and as such provide poor estimates of the true ERFs.

It is important to note that in this example and in other multi-event experiments where the number of each event per trial does not change between trials, the DC contents of the ERFs are not deconvolvable. This can be verified in the current example by looking at the determinant of the convolution matrix and noticing that it is zero (and therefore not invertible) at or near DC frequencies (see Fig. 2*B*). This is because the DC contents of the ERFs are not affected by the timing jitters that are used to deconvolve the higher frequency contents and as such different trials here do not provide additional information on the value of the ERFs DC content. This leaves us with an ill posed problem of one equation for the total DC content and four unknowns corresponding to each ERF. Nevertheless after correcting for DC shifts in the results, the application of deconvolution method successfully recovers the ERFs with high fidelity in time or in frequency as shown in Fig. 2, *D* and in *E*, respectively.

The ERF estimation error can be measured by the average mean squared error (hereafter MSE) for all four events integrated over time and as a function of deconvolution iterations (Fig. 2*F*). The ERF estimation error shows >70% improvement just after the first deconvolution iteration in time compared with the raw PETHs (MSE = 0.7 vs. 2.6). The ERF estimate continues to converge to the true ERFs, and the error goes down by another order of magnitude after 1,000 iterations comparable to average PETH firing noise (MSE = 0.03 vs. 0.025 for noise variance in 40 trials). Our ability to fit the raw firing PETHs in the presence of all four events based on the deconvolved ERFs also greatly improves over deconvolution iterations with errors as low as 10^{−7} after 1,000 deconvolution iterations. (see Fig. 2*G*). A final deconvolved signal (equivalent to infinite iterations in time domain) can be obtained from frequency deconvolution with errors that are not distinguishable from the average firing noise (MSE = 0.027 vs. 0.025 for noise variance in 40 trials; see Fig. 2*F*).

##### FAST OCCURRING SINGLE EVENT DECONVOLUTION.

For the second experiment, let us assume an event A happening 10 times in each trial with randomly spaced inter-event intervals from a Poisson distribution. We assume to have a neuron with an ERF that shows excitation then inhibition to this event and inter-event intervals that are shorter than the ERF duration. The firing of this neuron in two different trials in the presence of 10 event repetitions is shown in Fig. 3*A*. The raw PETH is again obtained by convolving the ERF with the event timings (*Eq. 1* with Gaussian noise, SD = 1).The close temporal spacing between events in this case results in a broad distortion of ERF estimate based on the raw PETH. Figure 3*C* demonstrates a gross overestimation of both amplitude and duration of the ERF. Applying the deconvolution method in time or in frequency successfully recovers the actual neuronal response to a single occurrence of event A (Fig. 3, *D* and *E*). Although small jitters in the event timings tend to reduce the determinant of the convolution matrix in low frequencies, in this case, the convolution matrix is invertible over all frequencies making the deconvolution in the frequency domain feasible for all frequencies (Fig. 3*B*).

ERF estimation error measured by MSE again shows considerable improvement close to an order of magnitude only after the first deconvolution iteration compared with the raw PETH (MSE = 1.5 vs. 9.8). This error continues to go down after some 500 iterations to levels undistinguishable from average firing noise. The MSE resulting from the frequency domain deconvolution (equivalent to infinite iterations in time domain) is also comparable to the average PETH firing noise (MSE = 0.0022 for both time and frequency deconvolution vs. 0.0025 for noise variance in 40 trials containing 10 event repetitions each; Fig. 3*F*). Also the raw PETH fitting accuracy using the deconvolved ERF improves rapidly over deconvolution iterations where after 1,000 iterations, the MSE for predicting PETHs goes down to a negligible 10^{−12} (Fig. 3*G*).

### Cross-validation of deconvolution

The simulation results in the previous section demonstrate the feasibility of complete deconvolution of the ERFs when the data are generated based on the deconvolution model assumptions with the ERF estimation error bounded only by a linear combination of the raw PETH noise variances (see methods, *Eqs. 6B* and *10B*).

In dealing with the actual firing data, the assumptions of the model can be violated to varying degrees. For instance, the ERFs might not be stationary over time or might show dependence on the inter-event intervals (a.k.a. ERF interactions). Furthermore in dealing with the actual data, we might be limited by the number of available trials to average out the trial noise that would raise the lower bound of deconvolution error. These factors can adversely affect our ability to fully recover the true ERFs. If the ERFs were known, one could simply measure the MSE of the deconvolved ERFs in predicting the true ERFs and obtain a quantitative measure of how the assumption violations or trial noise are affecting the deconvolution performance over the iterations. Such a quantitative performance measure is in fact shown for the simulation results in Figs. 2*F* and 3*F* where the ERF estimation error was calculated and plotted over deconvolution iterations. In dealing with real data, obviously the true ERFs are not known and therefore another approach is needed to monitor the model performance over the deconvolution iterations.

One way of characterizing the effectiveness of the deconvolution for real firing data are by cross-validation of deconvolved ERFs. Cross-validation is a simple but powerful statistical technique that can be used to choose between models (Stone 1974), which in our case would be between performing or not performing the deconvolution and also between different levels of deconvolution for a given neuron. Cross-validation can be done by randomly dividing the trials into two groups, one used for training the model and finding the ERFs using a certain deconvolution iteration and the other one for validating the estimated ERFs to see if they can convolve back to predict the PETHs made from the test trials (see methods, *Eq. 11*).

We first tested the sensitivity of the deconvolution to noise by decreasing the signal-to-noise ratio in the four event deconvolution experiment shown in Fig. 2. The signals to noise values tested were 22, 12, 6, and 4 dB. Increasing the noise level resulted in increasing levels of cross-validation error as expected. Although in all cases, some levels of partial deconvolution outperformed the raw PETH in cross-validation, the deconvolution iteration resulting in the lowest cross-validation error decreased in conjunction with the signal-to-noise ratio. (Fig. 4*A*, *left*). Because in the case of simulation we a priori know the true ERFs, MSE for ERF prediction error at each noise level can be calculated as well. One can see that the trend of changes for the ERF prediction error is qualitatively the same as that of cross-validation (corrcoef = 0.85, *P* < 0.001) though the iteration with minimum error for each noise level is not necessarily the same for the two error measures (see Fig. 4*A*, *right*). This confirms that in the absence of the ERF prediction error, cross-validation error can provide a comparable criterion for finding the appropriate level of deconvolution in a given neuron (for a derivation of the relationship between ERF and cross-validation MSEs look at *Eqs. 12* and *13* in the methods).

It is important to note that the changes in the cross-validation and ERF MSE reflect changes in both the estimates bias and variance. In cases of low signal-to-noise ratio, the amount of noise that is added after a few iterations integrated over time more than cancels out the improvement in the ERF estimate bias that results in the rising of the cross-validation MSE after only a few trials as can be seen in Fig. 4*A*, both *right* and *left columns* (also see methods *Eq. 6B* for added noise terms in each trial).

We also tested the sensitivity of the deconvolution to increasing levels of ERF interactions in one of the many possible ways that such interactions may occur. In this case, we assumed that the amplitude of the ERFs would decrease and their duration would increase with the interval between the events. The performance of the deconvolution model was compared when there were no interactions versus when there was interactive modulation of the first ERF, the first and second ERFs, up to the point that all four ERFs were modulated by the spacing between their events and the next event except for the last event's ERF that was modulated by the interval between itself and the previous event. Again increasing levels of ERF interactions resulted in increased cross-validation error as expected. And although in all cases deconvolution outperformed the raw PETH in cross-validation, the deconvolution iteration resulting in the lowest cross-validation error decreased by increasing ERF interactions (Fig. 4*B*, *left*). Here too the ERF prediction error showed a similar trend of changes as a function of ERF interaction and deconvolution level to the cross-validation (corrcoef = 0.6, *P* < 0.001), but the iteration with minimum ERF prediction error was not necessarily the same as that of the cross-validation error (Fig. 4*B*, *right*).

It is also worth pointing out that in both cases of low signal-to-noise ratio or high degree of ERF interactions, the estimate provided by frequency deconvolution (equivalent to complete deconvolution) can be inferior to a partial deconvolution (obtained after some deconvolution iterations). A fact that is demonstrated both by the ERF estimation error and the cross-validation error for varying levels of noise or ERF interactions. This can be partially due to the singularity of the convolution matrix (Fig. 3*B*) in this experiment for the lowest frequencies as well as the discussed effect of noise in low signal-to-noise ratio signals in raising the estimate errors.

In using the cross-validation as a tool for selecting the appropriate deconvolution iteration we define the “optimal iteration” to be the iteration resulting in the lowest significant cross-validation error compared with the raw PETH error and the “maximal iteration” to be the highest iteration having a cross-validation error significantly lower than that of the raw PETH (paired *t*-test *P* < 0.01 for significance). Optimal and maximal iterations therefore define the most and least conservative criteria respectively for the degree of deconvolution that can be safely applied to the data.

Finally it is worth mentioning that we could have chosen an alternative way of characterizing the effectiveness of the deconvolution by dividing the data randomly into two equal parts and using the complete deconvolution from one half to estimate the ERF MSE directly based on the plug in principle. This direct method, which normally yields similar results to cross-validation (see methods, *Eq. 14*), requires double the amount of deconvolution compared with cross-validation and therefore is more computationally intensive. Furthermore in many cases, complete deconvolution of ERFs might not be available or might show large fitting errors (see next section, Figs. 6*B* and 7*A*). Therefore we have used cross-validation as the faster and more robust approach in this paper thought when applicable all of our results with regard to cross-validation should be replicable using this direct method as well.

### Experimental results

In this section we demonstrate the use of the deconvolution method to parse out actual neuronal responses and will use the cross-validation technique introduced in the previous section to calibrate the number of deconvolution iterations applied on each neuron given its firing pattern and the degree of noise.

We used a subset of 106 neurons recorded in the nucleus accumbens (NAc) in a discriminative stimulus task (Ambroggi et al. 2008). Briefly, two auditory cues, discriminative stimuli (DS) and neutral stimuli (NS) were randomly presented (with durations of ≤10 s for the DS and 10 s for the NS) on a variable interval schedule with an average interval of 30 s. A lever press was required during DS presentation to cause the delivery of a reward into an adjacent reward receptacle and responding during NS was never rewarded. For simplicity, only the DS is displayed in Fig. 5 with intervals between the two DS's shown on average to be 60 s.

Firing of neurons on trials when the animal successfully responded to the cue can be modulated by four events that happen shortly one after the other. These events include presentation of the DS, lever press (LP), receptacle entry (Ent), and receptacle exit (Ext). Many of the neurons recorded in this task show changes that are time-locked to more than one event. In fact complex responses to multiple events are quite common in the rat ventral striatum (Humphries and Prescott 2009; Nicola 2007; Nicola et al. 2004b; Schultz et al. 1992). Figure 6*A* shows the firing raster of a putative medium spiny neuron recorded in rat NAc during the rewarded trials. This neuron is excited after the DS, shows inhibition during the consummatory period between receptacle entry and exit and then has a brief excitation following receptacle exit. Although the PETH fitting error decreases rapidly over deconvolution iterations, cross-validation for this neuron reaches a minimum after about three deconvolution iterations (optimal deconvolution) but then increases back with increasing degree of deconvolution with the maximal deconvolution happening at iteration 22 (Fig, 6, *B* and *C*). The result of partial deconvolution with optimal and maximal iterations is shown together with the raw PETHs in Fig. 6*A*. One can see that the excitation to DS based on the deconvolved ERF has a shorter duration compared with what was previously estimated by the raw PETH. This can be verified by looking at the firing raster time locked to the DS where part of the excitation can be seen to be due to lever press preparation. Deconvolution also removes the broad consumption-related inhibition from the DS PETH as well as the lever press and cue related excitations from the Ent and Ext PETHs. The presence of the excitation following the receptacle exit in the Ext ERF shows that in fact this activity is correlated with exiting the receptacle and not a contamination by other events.

The fact that the partial deconvolution outperforms the raw PETHs in cross-validation shows that the neuron's activity obeys at least partially the deconvolution model (see methods, *Eq. 1*). But possibly due to interactions, nonlinearities and moderate signal-to-noise ratio (∼10 dB, see Fig. 4*A*), complete deconvolution of responses is not feasible. Another fact that further complicates the complete deconvolution of this neuron's firing is the fact that the determinant of the convolution matrix made from event timings shows a rather wide notch for low frequencies that will make complete deconvolution at those frequencies erroneous. This results in high PETH fitting errors for frequency deconvolution in Fig. 6*B*. (see Fig. 6*D* for the convolution determinant).

Among the 106 neurons recorded, 85 were responsive to at least one event in rewarded trials. As shown in Fig. 7*A*, the average PETH fitting error for responsive neurons decreases rapidly over deconvolution iterations while on average the frequency deconvolution shows large errors in deconvolving the ERFs. Once more the poor performance of the frequency deconvolution can partially be explained by average determinant of the convolution matrix which shows a wide notch over low frequencies (see Fig. 7*C* for the average convolution determinant). As discussed previously, this is a direct consequence of multi-event experimental paradigms with unchanging number of occurrences for each event per trial in which case a partial temporal deconvolution usually outperforms frequency deconvolution.

The optimal deconvolution for each of the responsive neurons can be determined by finding the minimum of the cross-validation error term. The population average of the cross-validation error for the responsive neurons indicates that on average partially deconvolved ERFs outperform raw PETHs in predicting the validation trial firing patterns with the average cross-validation error being significantly lower than the average raw PETH errors for up to six iterations (multiple comparison *P* < 0.01 Bonferroni) and a maximum improvement of ∼8% in the average cross-validation MSE (Fig. 7*B*).

For 25% of the neurons, the minimum of the cross-validation error happens for one or more iterations. For the other neurons, this minimum happens at the zero level deconvolution. The deconvolved ERFs at zero iteration are simply the raw PETHs from the training trials normalized by a factor, which, in this case, is the number of event types (see methods, *Eq. 5*). The fact that in all the responsive neurons partial deconvolution outperforms raw PETHs in predicting the validation trial firings means that the raw PETHs time-locked to each event are influenced by the timing of the surrounding events. The optimal deconvolution found for each neuron is correlated with the signal-to-noise ratio (corrcoef = 0.71 *P* < 0.01) with higher degrees of deconvolution corresponding to higher signal–to-noise ratios (Fig. 7*D*). The maximal deconvolution iteration for each neuron shows a wider distribution over iterations with an average of about six iterations (mode equal to 4 iterations). The maximal iteration also shows dependence on the level of signal to noise with higher iterations corresponding to higher levels of signal-to-noise ratio. (corrcoef = 0.40 *P* < 0.01; Fig. 7*E*).

The neuronal responses shown in FIG. 8. a demonstrate that in fact in comparison with the raw PETHs the degree of contamination by non-time-locked events is reduced for all of the deconvolved ERFs after partial deconvolution (results of maximal iteration shown). For instance the duration of excitations following the DS is reduced significantly based on the deconvolved DS ERF compared with the raw DS PETH (average 1.51 vs. 2.18 s, paired *t*-test *P* < 0.05), pointing to the fact that the post-DS excitations in these cells were probably contaminated by excitations preceding the lever presses (Fig. 8*B*). Interestingly some responses were not affected by deconvolution. Indeed receptacle exit ERFs maintained postexit excitations after deconvolution (*P* > 0.05 paired *t*-test for average response from −1 to 1 s around receptacle exit) while the cue and lever press contaminations on the exit ERFs are attenuated after deconvolution (*P* < 0.01 paired *t*-test for average response from −7 to −5 s around receptacle exit corresponding to average time of cue presentation and lever press, respectively). Also inhibitions during the consumption period are shown to be mainly correlated with the receptacle entry and not with the receptacle exit since the sustained inhibition prior to the exit is removed after deconvolution from Ext ERFs. On the other hand, Ent ERFs still show inhibition during the consumption even after the deconvolution (Fig. 8*B*). However, it is important to be cautious when interpreting the appearance or disappearance of broad sustained changes compared with the faster phasic changes in the ERFs following the deconvolution due to the problematic nature of low frequency deconvolution in this experiment (Fig. 7*C*).

Besides changes in response dynamics, there is a considerable change in the percentage of responsive neurons to each event after deconvolution. The total number of neurons responsive to each event was reduced by 36% on average (averaged over all 4 events; see Fig. 8*A*, *right* vs. *left*). This reduction happens while the total number of neurons that are responsive to at least one event stays roughly the same before and after deconvolution (85 neurons before vs. 79 neurons after). This indicates that the deconvolution method is effective in reducing the neurons that are falsely categorized as responsive to multiple events (Fig. 8*C*).

## DISCUSSION

The ambiguity that typically occurs in the analysis of neuronal coding of temporally close events can be viewed as a problem with the PETHs being the convolved readout of several underlying ERFs that are time-shifted with respect to each other over trials. We first demonstrated that the raw PETHs in some cases provide a poor estimate of the underlying ERFs due to their overlap. We then presented a deconvolution formulation that is suited to the specific problem of finding ERFs in a multi-event task based on raw PETHs. We derived a pair of theoretically equivalent formulations for deconvolution in the time and in the frequency domains and showed through simulation that the formulation proposed can successfully parse out individual ERFs in multi-event and in fast occurring single event tasks. Because convolution in time is equivalent to multiplication in frequency (Oppenheim et al. 1997), deconvolution in the frequency domain theoretically simplifies the inversion of the convolution matrix and so is used for complete deconvolution. On the other hand, the temporal approach involves a step-by-step deconvolution to approximate this inverse. Although the two methods are theoretically equivalent, the iterative approach in the temporal deconvolution is more robust when applied to actual data. This is due to the fact that the deconvolution model entails validity of certain assumptions such as linearity and stationarity of responses and nonsingularity of the convolution matrix that are not necessarily true when dealing with actual event timings and neuronal responses. In dealing with neuronal recordings, the temporal approach has the advantage of permitting a partial deconvolution given the level of the noise and the degree to which the data deviates from model assumptions. In this sense, the iterative temporal approach can be more flexible than the complete deconvolution approach in the frequency domain.

Noise can also play a limiting role on the degree of deconvolution that can be successfully applied to the actual PETHs (Figs. 4 and 7). It is therefore recommended to use a low-pass filter on the data by applying Gaussian smoothing or Wiener filters (Papoulis 1977) on the data or have bin sizes that are large enough to limit the noise power before applying the deconvolution. Noise levels can also be lowered proportional to the square root of the number of trials. As shown in the result section, neurons with higher signal–to-noise ratios tend to have higher optimal and maximal deconvolution iterations.

It is also important to note that the experimental design and the amount of jitter in the event timings can influence the effectiveness of the deconvolution procedure, especially in deconvolving the broad tonic changes in the ERFs (which correspond to low frequency content of the signals). In fact there are two channels of information that are used for deconvolution, one is the changes in the number of occurrences of each event per trial and the other is the changes in the inter-event intervals between different events over trials. Of these two factors, the former is more important when deconvolving the tonic changes as they are less affected by the event timing jitters due to their slow temporal dynamics. For example in the discussed experimental paradigm where each event occurred only once per trial, the deconvolution of the lowest frequencies was not well defined. In such cases, caution is recommended when drawing conclusions about changes in tonic components of the ERFs. The best case scenario for successful deconvolution ERFs therefore includes variable number of occurrences of events per trial with extensive event timing jitter between trials.

Cross-validation can be used as an effective criterion to choose between using raw PETHs without deconvolution or a partial deconvolution to obtain the best possible estimates for the ERFs. Depending on the error tolerance, a range of partial deconvolutions from optimal to maximal can be chosen for each neuron. The result of the partial deconvolution should provide ERF estimates that are closer to the true ERFs. This in turn can inform the response detection algorithms to provide a better estimate of the responsive neurons to each event.

Deconvolution on real data obtained from extracellular recordings provides the proof of concept that this method can modify the conclusions obtained from such data. First it allows better categorization of the response types. Indeed as shown in Figs. 1 and 6, based on raw PETHs, neurons can be erroneously categorized as responsive to a given event if they respond strongly to another close event. The fact that cross-validation finds the iteration zero on average to fit the validation data better than the raw PETH (Fig. 7*B*) supports the idea that many of the raw PETHs are dependent on the responses to the surrounding events. Categorizing the neurons based on their deconvolved ERFs addresses this issue. As expected using the deconvolved ERFs instead of the raw PETHs reduces the number of events the neurons are responsive to (Fig. 8*C*). Furthermore, the deconvolution allows a better estimation of the response dynamics, both in duration and amplitude. Notably if a neuron is responsive to two temporally close events, the responses add up. For example we find that the duration of the neuronal responses to reward predictive cues (DS) is reduced after deconvolution. This results from correcting for the preparatory responses to the lever press following the cue onset. Therefore deconvolution provides a better estimate of duration in the presence of temporally proximal events.

The ambiguities that arise from response overlap are not only relevant in electrophysiology but in neuroscience as a whole. In fact there are studies in functional MRI and EEG literature that have identified this issue and provided methods for deconvolution of responses (Bohórquez and Ozdamar 2006; Dale 1999; Glover 1999; Hinrichs et al. 2000) However, all of these methods only allow for a complete deconvolution (like the one presented here in frequency domain) and as such are prone to large errors in the presence of interactions and convolution matrix singularities as discussed and demonstrated in this paper. Furthermore some of these studies are only appropriate for a single event deconvolution and are not general enough to allow for multiple event types (Bohórquez and Ozdamar 2006; Glover 1999). To our knowledge, the method presented here is the most general formulation of deconvolution for trial based experimental paradigms. It is also the first to allow for a partial deconvolution with a criterion based on cross-validation to address the success of varying degrees of deconvolution in dealing with the real data.

Finally our formulation of deconvolution assumes the timings of all the events in a task to be available. It may be the case in some experiments that the existence and/or the timings of some events are not known or recorded. These unknown event occurrences, especially when correlated in time with the recorded events, can introduce significant errors in the ERF estimations. In these cases, other techniques such as blind deconvolution or source separation (Comon 1994; Haykin 1994) can be used. These techniques not only estimate the original signals (in our case the ERFs) but also try to approximate the mixing or convolution matrix. Methods such as independent component analysis (ICA) that assume statistical independence between the original signals are successfully applied in many source separation problems. The number of unknown events that needs to be estimated in this way might in turn be decided by a model selection technique such as cross-validation. Such an extensive model that includes known events and optimizes for unknown events and their timings, however, would be heavily parameterized and computationally intensive. The usefulness of such models in deconvolving the ERFs remains to be tested in the future.

In conclusion, we have described and tested a novel technique that allows separation of the components of the neuronal activity that are associated with temporally close but distinct behavioral events. We have demonstrated the value of the deconvolution method using simulated data showing that it can significantly improve estimates of the event-related neuronal responses compared with the raw PETHs. We then applied our method to single unit recordings of nucleus accumbens neurons obtained in behaving rats showing that the deconvolution of the neuronal responses disambiguates the classification of responsive neurons to each event and adjusts each event-related response by reducing the response component that is due to the overlap of temporally adjacent events. We believe that this novel analytical tool can greatly improve the interpretability of data and as such is essential for understanding the temporal relationship of neuronal activity to behavioral events.

## GRANTS

This work was supported by funds provided by the state of California for medical research on alcohol and substance abuse through the University of California, San Francisco.

## DISCLOSURES

No conflicts of interest are declared by the authors.

## ACKNOWLEDGMENTS

The authors are grateful to G. Hjelmstad and A. Amini for helpful discussions and comments on this work.

- Copyright © 2010 the American Physiological Society