Analyzing neural spike data

Today, we’ll get to analyzing spike data. In particular, we’ll deal with estimating firing rates, which are typically the variables of interest. Once we have a good firing rate estimate, we often want to know how such rates are related to behavior, either evoked by it (as in sensory areas) or driving it (as in the motor system).

We’ll begin with where we left off in the last section:

%run exploring_spike_data.ipynb

Estimating mean firing rates

As we saw before, the simplest method for estimating a firing rate for a cell is to average its responses over multiple trials. Here, we’ll practice that on real data.

Exercise

1. Using the spike histogram data, compute the average spike count in each bin across trials for unit 1.

2. Convert this from a mean spike count per bin to a mean firing rate.

3. Plot the result as a function of time.

Solutions:

Show code cell content Hide code cell content

# 1. Average spike counts for unit 1

# Recall that spk_bin is of shape n_timepoints/bin_width X n_trials X n_units.
# Let's use 1ms bins, and average across axis 1 (remember, Python follows zero-indexing)
bin_width = 1
spk_bin = spikes_to_bin_varwidth(spikes, t_min=0, t_max=n_timepoints, bin_width=bin_width)
avg_spk_count = np.mean(spk_bin[:,:,0], axis=1)

Show code cell content Hide code cell content

# 2. Convert this from a mean spike count per bin to a mean firing rate.

avg_rate = avg_spk_count/(bin_width*1e-3) # bin_width is ms, need to convert to seconds

Show code cell content Hide code cell content

# 3. Plot the result as a function of time.

plt.figure()
t = np.arange(0, n_timepoints, bin_width)
plt.plot(t, avg_rate)
plt.ylabel("Firing rate (Hz)")
plt.xlabel("Time (ms)")
plt.title("Unit 1");

Smoothing

Typically, even with trial averaging, firing rates have a lot of noise. Under the assumption that rates are not changing drastically from one bin to the next, we can often get a more interpretable result by smoothing.

The most common method of smoothing is calculating a moving average, which we can do using NumPy’s filtering commands. This demo is more complicated than you might need, so just focus on the key line, which involves numpy.convolve.

Exercise

1. Smooth the trial-averaged firing rate over a 20ms window.

2. Plot this as a function of time on the same plot as the unsmoothed average.

Solutions:

Show code cell content Hide code cell content

# 1. Smooth the trial-averaged firing rate over a 20ms window.

window_len = 20
w = np.ones(window_len)
y = np.convolve(w/w.sum(), avg_rate, mode='same')

Show code cell content Hide code cell content

# 2. Plot smoothed and unsmoothed together.

plt.plot(t, avg_rate, alpha=0.25)
plt.plot(y, color='black')
plt.ylabel("Firing rate (Hz)")
plt.xlabel("Time (ms)")
plt.title("Unit 1")
plt.legend(["Raw", "Moving average, 20ms"]);

Estimating firing rates for single trials

In cases where trial-to-trial variability (neural or behavioral) is important, we might want to have an estimate of firing rate on a single trial to correlate with behavior. This is a tricky proposition, since, again, firing rate is not something we directly observe. In statistical language, we observe events from a Poisson process, which is equivalent, as we have seen, to observing counts within time bins that are governed by a Poisson distribution. In both cases, we assume there is an underlying spike rate or spike density \(\lambda(t)\) that varies in time.

As discussed previously, this is a hard problem, but as seen in our investigation of the PSTH, one approach is to trade variance for bias by smoothing the binned single-trial spikes.

So let’s estimate some single-trial firing rates. The simplest kerenel is a boxcar or sliding window of width \(w\). This simply replaces each spike count of \(n\) with a spike count of \(n/w\) in each bin of the window (i.e., a moving average).

Exercise

1. Create a time points x trials x units array of firing rate estimates by performing a moving average. The code above gives a straightforward way to do this.

2. Experiment with different window widths. What tradeoff is made by using these different widths?

3. Plot the collection of single-trial firing rates for a given unit as the first panel of a multi-panel plot of the same form as last class.

Solutions:

Show code cell content Hide code cell content

# 1. Create the array of firing rate estimates 

def firing_rate_moving_avg(spikes, window_len=20):
    w = np.ones(window_len)
    estimate = np.zeros_like(spikes) # Create an empty array of the same size as spikes
    for i in range(n_trials):
        for j in range(n_units):
            y = np.convolve(w/w.sum(), spikes[:,i,j], mode='same')
            estimate[:,i,j] = y #/ (window_len * 1e-3)
    return estimate

Show code cell content Hide code cell content

# 2. Different window widths

window_lens = [20,50,100,200]
for w_len in window_lens:
    lambda_ = firing_rate_moving_avg(spk_bin, w_len)[:,10,0]
    
    # divide by bin size to get rate
    plt.plot(lambda_/1e-3)
plt.legend(window_lens, title="Window\nlength (ms)");
plt.ylabel("Rate (spks/s)")
plt.xlabel("Time (ms)");

Show code cell content Hide code cell content

# 3. Plot all the rates

lambda_ = firing_rate_moving_avg(spk_bin, window_len=50)

# Plot the behavior as before (see Week 1, part 1)
def plot_summary(unit, mat, H, V):
    fig, axes = plt.subplots(3, 1, figsize=(4,9));
    im1 = axes[0].pcolormesh(mat[:,:,unit].T);
    im2 = axes[1].pcolormesh(H,cmap="jet")
    im3 = axes[2].pcolormesh(V,cmap="jet")
    plt.colorbar(im1,ax=axes[0])
    plt.colorbar(im2,ax=axes[1])
    plt.colorbar(im3,ax=axes[2])
    axes[0].set_title("Unit 1")
    axes[1].set_title("Horizontal position")
    axes[2].set_title("Vertical position")
    fig.tight_layout()

plot_summary(0, lambda_, H, V)

Gaussian smoothing

In contrast to the moving average, we might like to use a filter that does not “share” the spike equally among many time bins but keeps the bulk of its effect concentrated near the real data. That is, we’d like something more bump-like. A natural candidate is the Gaussian, or Normal, distribution.

For smoothing by an arbitrary filter, NumPy has many methods available, but the most straightforward one is to generate a vector representing the filter and use the convolution command (convolve) to perform the smoothing. This is the same method we’ve seen now several times, but here we’ll dig a little more into the details.

Exercise

1. Create a Gaussian filter.

   1. Start by creating a range of times over which to define the filter. I suggest -1000:1000, which is 1 second on either side of 0.

   2. Pick a smoothing width. This will be the standard deviation of the distribution. I suggest 20ms to start.

   3. Create the filter. The scipy.stats.norm.pdf function implements the classic bell curve.

   4. Normalize your filter so that it sums to 1. This is important so that we maintain the total spike density in the time series (i.e., the sum of the smoothed signal is still the total spike count.)

   5. Plot your filter as a function of time to make sure it’s right.

2. Perform the smoothing. This will work just like the moving average, except you will use your Gaussian window to generate the single-trial estimates. A few caveats:

   • The first argument to convolve should be the spike count time series. The second argument is the filter. While convolution as a mathematical operation treats both arguments the same (it is commutative), extra arguments to the function do make assumptions.

   • To wit: you will need to specify the 'same' argument to convolve. This is to ensure that NumPy returns a time series of the same length as the one you put in. It’s also the reason we needed to specify the spikes and filter in the order we did, and why the middle of the filter is treated as time lag 0.

3. Generate the same plot as above using these new smoothed firing rates.

4. What happens as you change the width of the filter?

Solutions:

Show code cell content Hide code cell content

# 1. Create a Gaussian filter.

import scipy.stats as stats

window_size = np.arange(-1000,1001,1)
window_std = 20
window = stats.norm.pdf(window_size, 0, window_std)
window /= window.sum()
plt.plot(window_size[750:1250], window[750:1250])
plt.xlabel("Time (ms)");

Show code cell content Hide code cell content

# 2. Perform the smoothing.

# Create an empty array of the same size as spikes
estimate = np.zeros_like(spk_bin)

# With mode='same' we need to truncate the window to the middle 500 time points
for i in range(n_trials):
    for j in range(n_units):
        y = np.convolve(window[750:1250], spk_bin[:,i,j], mode='same')
        estimate[:,i,j] = y

Show code cell content Hide code cell content

# 3. Generate the same plot as above using these new smoothed firing rates.

plot_summary(0, estimate, H, V)

Show code cell content Hide code cell content

# 4. What happens as you change the width of the filter?

def firing_rate_estimate(spikes, w):
    estimate = np.zeros_like(spikes) # Create an empty array of the same size as spikes
    for i in range(n_trials):
        for j in range(n_units):
            y = np.convolve(w, spikes[:,i,j], mode='same')
            estimate[:,i,j] = y
    return estimate

window_stds = [20,50,100,200]
window_size = np.arange(-250,250)

for std in window_stds:
    window = stats.norm.pdf(window_size, 0, std)
    window /= window.sum()
    lambda_ = firing_rate_estimate(spk_bin, window)[:,10,0]
    
    # convert to rate by dividing by bin size
    plt.plot(lambda_/1e-3)   
    
plt.legend(window_stds, title="Standard\ndeviation");
plt.ylabel("Rate (spks/s)")
plt.xlabel("Time (ms)");