Analyzing Rhythms Part 3 (Spectra of spike trains)

# Load modules we'll need.
import numpy as np
import matplotlib.pyplot as plt

Example: a randomly spiking neuron

To start, let’s create a fake spike train for a randomly spiking neuron, and compute the autocovariance and spectrum.

N  = 5000;                          # Number of bins.                   
dt = 0.001;                         # Duration of each bin [s].
T  = N*dt;                          # Total time of observation [s].
tm = np.arange(0,N)*dt;             # Time axis for plotting

lambda0 = 5                         # Average firing rate [Hz]
p0      = lambda0*dt;               # Probability of a spike in a time bin
dn      = np.random.binomial(1,p0,N)# Create the spike train as "coin flips"

plt.plot(tm, dn)                    # Plot it.
plt.xlabel('Time [s]');

# Compute the autocovariance.

ac_xx = "SOMETHING

                                    # Plot it.
lags  = np.arange(-N + 1, N)        # Create a lag axis,
plt.plot(lags * dt, ac_xx)          # ... and plot the result.
plt.xlabel('Lag [s]')
plt.ylabel('Autocovariance');

print('lambda0*dt = ',lambda0*dt)  # Compare expected r_{nn}[0] to computed value.
print('r_{nn}[0]  = ',ac_xx[N-1])

# Compute the spectrum.

Dj = "SOMETHING"                            # Compute the FT,
Pj = "SOMETHING"                            # ... and the spectrum.

f = np.fft.fftfreq(N, dt)                   # Create frequency axis.

plt.plot(f, np.real(Pj))                    # Plot the spectrum.
in_class_guess = "SOMETHING"                # And our guess from in-class analysis.
plt.plot(f,in_class_guess*np.ones(N), 'b')
plt.xlabel('Frequency [Hz]');

# Repeat the entire simulation many times, and plot the average autocovariance and spectrum

K  = 1000                                       # Number of times to repeat the simulation.
lambda_est = np.zeros(K)                        # Vector to store estimate of lambda.
P  = np.zeros([K,np.size(Pj)])                  # Matrix to store spectra,
AC = np.zeros([K,2*N-1])                        # ... and AC.

for k in np.arange(K):                          # For each repeat,
    dn = np.random.binomial(1,p0,N)             # ... create a new spike train.

    lambda_est[k] = "SOMETHING"
                                                # Compute the AC,
    ac_xx = "SOMETHING"
    AC[k,:] = ac_xx                             # ... and save the result.
    
    Dj = "SOMETHING"                            # Compute the FT,
    Pj = "SOMETHING"                            # ... and the spectrum,
    P[k,:] = Pj                                 # ... and save the result.

plt.figure()                                    # Plot it.
plt.plot(lags*dt, np.mean(AC,0))
plt.xlabel('Lag [s]')
plt.ylabel('Autocovariance');

plt.figure()
plt.plot(f, np.mean(P,0))                       # Plot the spectrum, averaged over repeats
in_class_guess = "SOMETHING"                    # Use lambda_est to compute our guess from in-class analysis.
plt.plot(f,in_class_guess*np.ones(N), 'b')
plt.xlabel('Frequency [Hz]');

Example: a randomly spiking neuron + refractory period

Now, let’s create a fake spike train for a randomly spiking neuron with a refractory period, and compute the autocovariance and spectrum.

N  = 5000;                          # Number of bins.                   
dt = 0.001;                         # Duration of each bin [s].
T  = N*dt;                          # Total time of observation [s].
tm = np.arange(0,N)*dt;             # Time axis for plotting

lambda0 = 5                         # Average firing rate [Hz]
p0      = lambda0*dt;               # Probability of a spike in a time bin
dn      = np.random.binomial(1,p0,N)# Create the spike train as "coin flips"
refractory_period = 10              # Add a refractory period
for i in np.arange(N):
    if dn[i]==1:
        dn[i+1:i+refractory_period] = 0

plt.plot(tm, dn)                    # Plot it.
plt.xlabel('Time [s]');

# Compute the autocovariance.

# Compute the spectrum.

# Repeat the entire simulation many times, and plot the average autocovariance and spectrum