MA 568 – Statistical Analysis of Point Process Data
Problem Set #4 Due November 20, 2018
Most of the analyses we have developed so far have focused on point process data in the time domain. It is often useful to visualize spiking data and construct models in the frequency domain, as well. In this problem, we will perform a spectral analysis on the spiking activity of a pair of neurons in rat barrel cortex.
The vibrissae (whiskers) of a rat are arranged in an ordered pattern of rows and columns along its face. The barrel cortex contains a topographic representation of this pattern of whiskers. Stimulation of a single whisker causes neurons in the corresponding cortical column, known as a “barrel”, to fire spikes.
Please download the file BarrelSpikes.mat (BarrelSpikes.csv) from the course website. This dataset contains the spiking activity of two neurons in separate
barrels, while a bar was swept across its whiskers in a rostral to caudal direction (from the nose toward its tail). The stimulation was regular and periodic, but the experimenter forgot to write down the frequency of stimulation. Each spike train was recorded at a sampling frequency of 1 kHz.
1. Plot the spiking activity of both neurons in a way that clearly shows the relative timing of their spikes. Describe the structure of the spiking activity for each neuron and the relation of the spiking activity between them.
t
0
the entire observation interval. The autocovariance function for a single point process is simply the cross-covariance of one point process with itself. Compute and plot the autocovariance function for the spiking activity of each neuron as a function of lag, and the cross-covariance function between the two spike trains. What do these functions suggest about the stimulus driving these neurons and the relation between the two point processes?
Useful MATLAB (R) function: xcorr (acf)
3. An estimate of the power spectral density (PSD) of a point process can be computed in two ways: 1) by multiplying the Fourier transform of the original data by its complex conjugate, or 2) by taking the Fourier transform of its autocovariance function. Compute and plot the power spectral density of the spiking activity for each neuron. Be carful to label the frequencies on the x-axis correctly and indicate their units. What is the stimulation frequency? Over what range of frequencies is bursting behavior observed?
Useful MATLAB (R) functions: fft, conj
2. The cross-covariance function between two point process is given by CCF(τ)=∑(ΔN(1) −λ(1)Δt)(ΔN(2) −λ(2)Δt), where ΔN(i) indicates whether the ith
t0t−τ0 t
neuron fired a spike in (t,t + Δt] and λ(i) is the average spiking rate for the ith neuron over
page 2: MA 568 – Problem Set 5
4. An estimate of the complex cross-spectrum between two point processes can be computed either 1) by multiplying the Fourier transform of one process by the complex conjugate of the Fourier transform of the other, or 2) by taking the Fourier transform of their cross-covariance function. Compute the complex cross-spectrum between the two point processes and plot its magnitude and phase. What is the phase difference at the frequency corresponding to the bar stimulus? Which neuron comes from the barrel that represents the whisker closest to the nose?
Useful MATLAB (R) functions: abs, angle (ccf, mod, arg)
5. Using the results of this spectral analysis as a starting point, write down an appropriate conditional intensity model for the firing of each of these two neurons. What variables should be part of this model? Using any of the methods we learned in this class, find maximum likelihood estimates of the parameters of your model and measure the goodness- of-fit between your model and the data. How do the estimated model parameters relate to the structure observed in the frequency domain analysis?
6. (Extra credit) Compute spectrograms of the spectral power as a function of time and frequency for both neurons. Do you detect any nonstationarity in the spectral representations of the data? How is this reflected in your spectral estimates from parts 3- 4? How might you account for this in your model from part 5?