Microsoft Word – Problem Set #4.docx
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.
2. The cross-covariance function between two point process is given by
( )( )(1) (1) (2) (2)0 0( ) t t
t
CCF N t N tττ λ λ−= Δ − Δ Δ − Δ∑ , where ( )itNΔ indicates whether the ith
neuron fired a spike in ( , ]t t t+Δ and ( )0
iλ is the average spiking rate for the ith neuron over
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
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?