Auto- and crosscorrelation functions for spike trains


Run this example as a Jupyter notebook:

See our guide for more information and troubleshooting.


A time bin of size tbin is centered around the time difference it represents. If the correlation function is calculated for tau in [-tau_max, tau_max], the pair events contributing to the left-most bin are those for which tau in [-tau_max-tbin/2, tau_max+tbin/2) and so on.

Correlate two spike trains with each other assumes spike times to be ordered in time. tau > 0 means spike2 is later than spike1

  • tau_max: maximum time lag in ms correlation function

  • tbin: bin size

  • spike1: first spike train [tspike…]

  • spike2: second spike train [tspike…]

import nest
import numpy as np


def corr_spikes_sorted(spike1, spike2, tbin, tau_max, resolution):
    tau_max_i = int(tau_max / resolution)
    tbin_i = int(tbin / resolution)

    cross = np.zeros(int(2 * tau_max_i / tbin_i + 1), "d")

    j0 = 0

    for spki in spike1:
        j = j0
        while j < len(spike2) and spike2[j] - spki < -tau_max_i - tbin_i / 2.0:
            j += 1
        j0 = j

        while j < len(spike2) and spike2[j] - spki < tau_max_i + tbin_i / 2.0:
            cross[int((spike2[j] - spki + tau_max_i + 0.5 * tbin_i) / tbin_i)] += 1.0
            j += 1

    return cross


nest.ResetKernel()

resolution = 0.1  # Computation step size in ms
T = 100000.0  # Total duration
delta_tau = 10.0
tau_max = 100.0  # ms correlation window
t_bin = 10.0  # ms bin size
pc = 0.5
nu = 100.0

nest.local_num_threads = 1
nest.resolution = resolution
nest.overwrite_files = True
nest.rng_seed = 12345

# Set up network, connect and simulate
mg = nest.Create("mip_generator")
mg.set(rate=nu, p_copy=pc)

cd = nest.Create("correlation_detector")
cd.set(tau_max=tau_max, delta_tau=delta_tau)

sr = nest.Create("spike_recorder", params={"time_in_steps": True})

pn1 = nest.Create("parrot_neuron")
pn2 = nest.Create("parrot_neuron")

nest.Connect(mg, pn1)
nest.Connect(mg, pn2)
nest.Connect(pn1, sr)
nest.Connect(pn2, sr)

nest.Connect(pn1, cd, syn_spec={"weight": 1.0, "receptor_type": 0})
nest.Connect(pn2, cd, syn_spec={"weight": 1.0, "receptor_type": 1})

nest.Simulate(T)

n_events_1, n_events_2 = cd.n_events

lmbd1 = (n_events_1 / (T - tau_max)) * 1000.0
lmbd2 = (n_events_2 / (T - tau_max)) * 1000.0

spikes = sr.get("events", "senders")

sp1 = spikes[spikes == 4]
sp2 = spikes[spikes == 5]

# Find crosscorrelation
cross = corr_spikes_sorted(sp1, sp2, t_bin, tau_max, resolution)

print("Crosscorrelation:")
print(cross)
print("Sum of crosscorrelation:")
print(sum(cross))

Gallery generated by Sphinx-Gallery