Spike Train Rasters and Firing Irregularity

Import Libraries

import glob
import scipy.io
import numpy as np
import matplotlib.pyplot as plt
import elephant.statistics
import elephant.spike_train_generation
%matplotlib inline

Utility Functions

Run the cell below to load the utility functions used in the exercises.

class utils:
    @staticmethod
    def load_neuron(mo_path, ts_path=None):
        mat = scipy.io.loadmat(mo_path)
        sf  = mat['SparseFormat'][0, 0]
        dt  = float(sf['TimeResolutionMS'][0, 0])
        def parse(sf_field):
            interval = sf_field['CutIntervalMS'][0]
            t_start  = float(interval[0])
            t_stop   = float(interval[1])
            result   = []
            for d in range(sf_field['Data'].shape[1]):
                matrix = sf_field['Data'][0, d].toarray()
                sts    = []
                for trial in range(matrix.shape[1]):
                    idx = np.where(matrix[:, trial])[0]
                    sts.append(t_start + idx * dt)   # numpy array of spike times in ms
                result.append(sts)
            return result, t_start, t_stop
        parsed, t_start, t_stop = parse(sf)
        out = {'mo': parsed, 't_start_ms': t_start, 't_stop_ms': t_stop}
        if ts_path is not None:
            mat2 = scipy.io.loadmat(ts_path)
            out['ts'] = parse(mat2['SparseFormat'][0, 0])[0]
        return out

    @staticmethod
    def pool_isis_in_window(sts, t0_ms, t1_ms):
        all_isi = []
        for st in sts:
            epoch = st[(st >= t0_ms) & (st < t1_ms)]
            if len(epoch) > 1:
                all_isi.extend(np.diff(epoch))
        return np.array(all_isi)

    @staticmethod
    def psth(sts, bin_ms=50, t_start=-1000, t_stop=1000):
        bins      = np.arange(t_start, t_stop + bin_ms, bin_ms)
        all_times = np.concatenate(sts)
        counts, _ = np.histogram(all_times, bins=bins)
        rate_hz   = counts / (len(sts) * bin_ms / 1000)
        centres   = (bins[:-1] + bins[1:]) / 2
        return centres, rate_hz

    @staticmethod
    def simulate_rate_step_trains(base_hz, step_hz, n_trains=200, t_half_ms=1000):
        """Poisson trains stepping from base_hz to step_hz at t_half_ms. Returns list of numpy arrays (ms)."""
        import quantities as pq
        import neo
        rate_arr = np.concatenate([np.full(t_half_ms, base_hz), np.full(t_half_ms, step_hz)])
        rate_sig = neo.AnalogSignal(rate_arr * pq.Hz, sampling_period=1 * pq.ms, t_start=0 * pq.ms)
        return [elephant.spike_train_generation.NonStationaryPoissonProcess(rate_sig)
                .generate_spiketrain().rescale(pq.ms).magnitude
                for _ in range(n_trains)]

    # --- plotting utilities ---------------------------------------------------

    @staticmethod
    def plot_raster(sts, ax=None, alignment_label='t = 0', title=None, ms=4, alignment_ms=0):
        """Plot a spike raster for a list of spike trains. Returns the axis."""
        if ax is None:
            _, ax = plt.subplots(figsize=(8, 4))
        for trial, st in enumerate(sts):
            ax.plot(st, np.full(len(st), trial), '|k', ms=ms)
        ax.axvline(alignment_ms, color='r', ls='--', label=alignment_label)
        ax.set_xlabel('Time (ms)')
        ax.set_ylabel('Trial')
        if title:
            ax.set_title(title)
        ax.legend()
        return ax

    @staticmethod
    def plot_cv_hist(cvs, title=None):
        fig, ax = plt.subplots(figsize=(5, 3))
        ax.hist(cvs, bins=10, color='steelblue', edgecolor='white')
        ax.axvline(1.0, color='r', ls='--', label='Poisson (CV = 1)')
        ax.set_xlabel('CV'); ax.set_ylabel('Count')
        if title: ax.set_title(title)
        ax.legend()

    @staticmethod
    def plot_cv_barchart(mean_cvs, xlabel='Direction'):
        fig, ax = plt.subplots(figsize=(6, 3))
        ax.bar(range(len(mean_cvs)), mean_cvs, color='steelblue')
        ax.axhline(1.0, color='r', ls='--', label='Poisson (CV = 1)')
        ax.set_xlabel(xlabel); ax.set_ylabel('Mean CV')
        ax.set_xticks(range(len(mean_cvs))); ax.legend()

    @staticmethod
    def plot_isi_histograms(isis_list, labels, n_bins=40):
        n = len(isis_list)
        fig, axes = plt.subplots(1, n, figsize=(4 * n, 3), sharey=True)
        if n == 1: axes = [axes]
        for ax, isis, label in zip(axes, isis_list, labels):
            cv = isis.std() / isis.mean()
            ax.hist(isis, bins=n_bins, density=True, color='steelblue', edgecolor='white')
            ax.set_title(f'{label}\nCV = {cv:.3f}')
            ax.set_xlabel('ISI (ms)')
        axes[0].set_ylabel('Density')
        fig.tight_layout()

    @staticmethod
    def plot_psth(sts, bin_ms=50, t_start=-1000, t_stop=1000):
        centres, rate_hz = utils.psth(sts, bin_ms, t_start, t_stop)
        fig, ax = plt.subplots(figsize=(8, 3))
        ax.bar(centres, rate_hz, width=bin_ms, color='steelblue', edgecolor='white')
        ax.axvline(0, color='r', ls='--', label='t = 0')
        ax.set_xlabel('Time (ms)'); ax.set_ylabel('Firing rate (Hz)')
        ax.legend()
        return centres, rate_hz

    @staticmethod
    def plot_measures_vs_rate(jump_rates, measures):
        _colors  = {'CV': 'steelblue', 'CV2': 'coral', 'Lv': 'seagreen'}
        _markers = {'CV': 'o', 'CV2': 's', 'Lv': '^'}
        fig, ax = plt.subplots(figsize=(6, 3))
        for label, values in measures.items():
            ax.plot(jump_rates, values,
                    _markers.get(label, 'o') + '-',
                    color=_colors.get(label, 'k'),
                    label=label)
        ax.axhline(1.0, color='k', ls='--', lw=0.8, label='Poisson expectation')
        ax.set_xlabel('Second-half rate (Hz)'); ax.set_ylabel('Mean value')
        ax.legend()

    @staticmethod
    def plot_cv_cv2_scatter(cv_all, cv2_all):
        lim = max(max(cv_all), max(cv2_all)) * 1.05
        fig, ax = plt.subplots(figsize=(5, 5))
        ax.scatter(cv_all, cv2_all, s=15, alpha=0.5, color='steelblue')
        ax.plot([0, lim], [0, lim], 'k--', lw=1, label='CV = CV2')
        ax.axhline(1.0, color='r', ls=':', lw=1, label='Poisson (CV2 = 1)')
        ax.set_xlabel('CV'); ax.set_ylabel('CV2')
        ax.set_xlim(0, lim); ax.set_ylim(0, lim)
        ax.legend()

import os
import owncloud
from tqdm import tqdm

oc = owncloud.Client.from_public_link('https://uni-bonn.sciebo.de/s/m2pCDzgb2wq9d9a', folder_password='ibots')

os.makedirs('data', exist_ok=True)
for f in tqdm(oc.list('/'), desc='Downloading Files to data/'):
    if not os.path.exists(f'data{f.path}'):
        oc.get_file(f.path, f'data{f.path}')

Load Recording Data

data = utils.load_neuron('data/joe096-5-C1-MO.mat', 'data/joe096-5-C1-TS.mat')

In the present exercises, you will analyze single unit recordings from the primary motor cortex (M1) of one monkey (monkey 1 in Rickert et al., 2009 and monkey M1 Rostami et al. 2024), which performed a delayed center-out task. The experiments were carried out in the lab of Alexa Riehle. The task and the data are described in detail in Rickert et al., (2009). For an introduction to directional tuning in the motor cortex an in-depth coverage is found in Riehle & Vaadia (2005).

Each trial of the experiment the monkey moved a cursor to one of 6 target directions. Spike trains are cut around two alignment events: movement onset (MO) and trial start (TS). The cell above loads the data of one neuron (neuron 5) of one recording session: data['mo'][d] is the list of spike trains of the individual trials of that neuron (numpy arrays of spike times in ms) for direction d (0-5) aligned to MO; data['ts'][d] provides the same trials aligned to TS.

Code Reference

Exercises

Example: Plot a spike raster for direction 0, aligned to MO.

utils.plot_raster(data['mo'][0], alignment_label='MO', title='Direction 0, MO-aligned')

Exercise: Plot MO-aligned rasters for all 6 movement directions in a 2×3 grid with shared axes.

Solution

fig, axes = plt.subplots(2, 3, figsize=(12, 6), sharex=True, sharey=True)
for d, ax in enumerate(axes.flat):
    utils.plot_raster(data['mo'][d], ax=ax, alignment_label='MO', title=f'Direction {d}')
fig.tight_layout()

Exercise: Compare MO-aligned and TS-aligned rasters for direction 0 side by side using utils.plot_raster. How does the choice of alignment event affect the apparent structure of the spike trains across trials?

Solution

fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True)
for ax, key, label in zip(axes, ['mo', 'ts'], ['MO', 'TS']):
    utils.plot_raster(data[key][0], ax=ax, alignment_label=label, title=f'Aligned to {label}')
fig.tight_layout()
# MO alignment concentrates peri-movement spiking near t=0.
# TS alignment spreads the same spikes because reaction time varies across trials.

The rasters above already hint that firing rate is not constant within a trial: spike density visibly increases around movement onset. Before measuring irregularity, it is useful to see what a rate step looks like in a controlled, simulated setting (where the ground truth is known) and then to confirm the same structure in the real data.

Code Reference

Exercises

Exercise: Simulate 30 Poisson spike trains stepping from 20 Hz to 80 Hz at t = 1000 ms using utils.simulate_rate_step_trains(20, 80, n_trains=30). Plot them as a raster with utils.plot_raster, placing the alignment marker at the rate step. What do you observe?

Solution

sts_step = utils.simulate_rate_step_trains(base_hz=20, step_hz=80, n_trains=30)
ax = utils.plot_raster(sts_step, alignment_label='Rate step (80 Hz)', alignment_ms=1000)
ax.set_xlim(0, 2000)
# Spike density visibly increases after t = 1000 ms.
# The raster reveals the rate jump without any computation.

(0.0, 2000.0)

Exercise: Compare rasters for three conditions side by side: 20→20 Hz (stationary control), 20→40 Hz, and 20→80 Hz, each with 30 trains. Which conditions show a visible rate change at t = 1000 ms?

Solution

conditions = [(20, 20), (20, 40), (20, 80)]
fig, axes = plt.subplots(1, 3, figsize=(14, 4), sharey=True)
for ax, (base, step) in zip(axes, conditions):
    sts = utils.simulate_rate_step_trains(base, step, n_trains=30)
    utils.plot_raster(sts, ax=ax, alignment_label='Rate step', title=f'{base}→{step} Hz', alignment_ms=1000)
    ax.set_xlim(0, 2000)
fig.tight_layout()
# 20→20 Hz: uniform spike density throughout (stationary Poisson).
# Rate increase becomes progressively more visible for 20→40 and 20→80 Hz.

The inter-spike interval (ISI) is the time between consecutive spikes. For a stationary spike train, the ISI distribution is interpretable in terms of point process theory, but this requires that the firing rate is constant within the analysis window. We begin with this assumption and will see shortly where it breaks down.

The coefficient of variation CV = σ_ISI / μ_ISI summarises the ISI distribution in a dimensionless number independent of the mean rate. For a homogeneous Poisson process (exponential ISIs), CV = 1. CV < 1 indicates more regular (sub-Poisson) firing; CV > 1 indicates bursty or clustered activity. Cortical neurons in the awake animal typically show CV ≈ 0.5-1.2 (Softky & Koch 1993; Shadlen & Newsome 1998).

For any stationary renewal process (each ISI drawn independently from the same distribution), the Fano factor and the squared CV are asymptotically linked (Cox & Lewis 1966; Nawrot 2010):

For a Gamma(κ) renewal process, CV = 1/√κ and FF ≈ 1/κ. The Poisson process is the special case κ = 1.

A fundamental limitation: ISI analysis is valid only when the firing rate is stationary within the observation window. In vivo, the rate often changes systematically within a trial (visible in the rasters: firing ramps up around movement onset). Pooling ISIs across such a window mixes intervals from different rate epochs and inflates CV beyond its intrinsic value. Two remedies are possible: (1) restrict analysis to a stationary sub-window identified from a PSTH, or (2) use local irregularity measures that compare only consecutive ISI pairs and are insensitive to slow rate drifts (next section).

You explored this rate-inflation effect in the homework. Here we quantify it directly by computing CV in short windows placed at and around the rate step, then confirm the pattern with a controlled rate-step simulation.

Code Reference

Exercises

Example: Compute and print mean ISI and CV for trial 0, direction 0, MO-aligned, ignoring that the data is highly non-stationary and therefore will yield an inflated CV.

st  = data['mo'][0][0]
isi = elephant.statistics.isi(st)            # numpy array of ISIs in ms
print(f'Mean ISI: {isi.mean():.1f} ms')
print(f'CV:       {elephant.statistics.cv(isi):.3f}')

Mean ISI: 44.5 ms
CV:       1.215

Exercise: Compute the per-trial CV for direction 0 (skip trials with fewer than 3 spikes). Store the values in cvs_d0 and plot the CVs across trials with utils.plot_cv_hist(). Print the mean CV.

Solution

cvs_d0 = [elephant.statistics.cv(elephant.statistics.isi(st))
          for st in data['mo'][0] if len(st) > 2]
utils.plot_cv_hist(cvs_d0)
print(f'Mean CV direction 0: {np.mean(cvs_d0):.3f}')

Mean CV direction 0: 1.105

Exercise: Let’s explore the non-stationarity of the data in more detail. Build a PSTH for direction 0 (all MO-aligned trials pooled, bin size 50 ms) using utils.plot_psth(). When does the firing rate begin to rise? When does it peak, and when does it start to stabilise? Use this as a guide to understand which post-movement windows are most non-stationary.

Solution

centres, rate_hz = utils.plot_psth(data['mo'][0], bin_ms=50)
# Rate is roughly flat before -200 ms and rises sharply around MO.
# The pre-movement window -900 to -200 ms looks approximately stationary.

Exercise: The PSTH above shows that firing rate shows both variable and stable periods. Let’s see how the CV behaves in smaller windows, where activity is more stationary, thus likely to lead to less inflation. Pool ISIs from direction 0 (MO-aligned, all trials) in 400 ms windows that slide across the post-movement period. Use utils.pool_isis_in_window and display the five ISI distributions side by side with utils.plot_isi_histograms. How does the ISI shape and CV change as you move away from the rate step? How does the CV differ from the one calculated across the whole data segment in the first exercises of this section?

Solution

window_starts = range(-1000,600,300)
window_ms     = 400

isis_list = [utils.pool_isis_in_window(data['mo'][0], t0, t0 + window_ms)
             for t0 in window_starts]
labels    = [f'{t0}–{t0 + window_ms} ms' for t0 in window_starts]

utils.plot_isi_histograms(isis_list, labels)
# Much lower CVs observed

Exercise: Let’s investigate in detail how a rate step can affect the CV measurement. For each step-to rate in jump_rates = [20, 30, 40, 60, 80] Hz, simulate 200 Poisson trains with utils.simulate_rate_step_trains(base_hz=20, step_hz=jump), similar to the rasters in the first section of this exercise. Compute the mean CV per condition, store in jump_cvs, and plot with utils.plot_measures_vs_rate(jump_rates, {'CV': jump_cvs}).

Solution

base_rate  = 20                     # Hz, first half of each train
jump_rates = [20, 30, 40, 60, 80]  # Hz, second half (20 = stationary control)

sts_by_jump = {}  # stored for reuse in the next exercise and section 3
jump_cvs    = []
for jump in jump_rates:
    sts = utils.simulate_rate_step_trains(base_rate, jump, n_trains=200)
    sts_by_jump[jump] = sts
    cvs = [elephant.statistics.cv(elephant.statistics.isi(st)) for st in sts if len(st) > 2]
    jump_cvs.append(np.mean(cvs))

utils.plot_measures_vs_rate(jump_rates, {'CV': jump_cvs})
# CV grows monotonically with the rate jump.
# Even though spiking is locally Poisson, mixing slow and fast ISIs inflates CV above 1.

Exercise: Using sts_by_jump, pool ISIs across all 200 trains separately for three conditions: 20->20 Hz, 20->40 Hz, and 20->80 Hz. Pass them as isis_list to utils.plot_isi_histograms().

Solution

isis_list = [
    np.concatenate([elephant.statistics.isi(st) for st in sts_by_jump[j] if len(st) > 1])
    for j in [20, 40, 80]
]
utils.plot_isi_histograms(isis_list, ['20->20 Hz', '20->40 Hz', '20->80 Hz'])
# Stationary (20->20 Hz): near-exponential ISI distribution, CV close to 1.
# Larger jumps: bimodal distribution mixing pre-jump and post-jump ISIs.
# The histogram shape reveals non-stationarity; the high CV reflects it numerically.

When the firing rate changes within a trial, CV is inflated by rate modulation, as demonstrated in the previous section. Two closely related local measures address this by comparing only consecutive ISI pairs normalised by their local sum. This makes them insensitive to slow rate drifts: if the rate changes slowly, consecutive ISIs change together and their normalised difference remains close to its stationary value.

CV2 (Holt et al. 1996; Ponce-Alvarez et al. 2010):

Local Variation Lv (Shinomoto et al. 2003):

For a stationary Poisson process, CV2 = Lv = 1. For a stationary renewal process, CV2 ≈ Lv ≈ CV (all three converge). The key diagnostic: if CV » CV2 (or Lv), the spike train is non-stationary within the observation window.

The two measures penalise ISI pair differences differently: CV2 uses the absolute difference (L1 penalty), so all deviations count equally regardless of size. Lv uses the squared difference (L2 penalty), giving disproportionately larger weight to pairs with very unequal consecutive ISIs. As a result, Lv is slightly more sensitive to extreme burstiness or pausing, while CV2 is more robust to individual outlier pairs. In practice, for typical cortical neurons both measures yield nearly identical results, and Lv = CV2 for a Poisson process (Ponce-Alvarez et al. 2010).

Code Reference

Exercises

Example: For direction 0, trial 0, MO-aligned, compute and compare CV, CV2, and Lv.

st  = data['mo'][0][0]
isi = elephant.statistics.isi(st)
print(f'CV  = {elephant.statistics.cv(isi):.3f}')
print(f'CV2 = {elephant.statistics.cv2(isi):.3f}')
print(f'Lv  = {elephant.statistics.lv(isi):.3f}')

CV  = 1.215
CV2 = 0.730
Lv  = 0.600

Exercise: Using sts_by_jump and jump_cvs from the previous section, compute mean CV2 and mean Lv for each condition. Store in jump_cv2s and jump_lvs. Plot CV, CV2, and Lv together with utils.plot_measures_vs_rate(jump_rates, {'CV': jump_cvs, 'CV2': jump_cv2s, 'Lv': jump_lvs}).

Restart point: re-run the simulation cell in the Firing Irregularity section to regenerate sts_by_jump and jump_cvs.

Solution

jump_cv2s, jump_lvs = [], []
for jump in jump_rates:
    cv2s, lvs = [], []
    for st in sts_by_jump[jump]:
        if len(st) > 3:
            isi = elephant.statistics.isi(st)
            cv2s.append(elephant.statistics.cv2(isi))
            lvs.append(elephant.statistics.lv(isi))
    jump_cv2s.append(np.mean(cv2s))
    jump_lvs.append(np.mean(lvs))

utils.plot_measures_vs_rate(jump_rates, {'CV': jump_cvs, 'CV2': jump_cv2s, 'Lv': jump_lvs})
# CV climbs with the rate step; CV2 and Lv stay flat at 1 - the true Poisson irregularity.
# Local measures are blind to slow rate changes because they compare only neighbouring ISI pairs.

Exercise: For each of the 6 directions, compute mean CV and mean CV2 per trial (MO-aligned, skip trials with fewer than 4 spikes). Print a summary table. Is CV consistently higher than CV2?

print(f"{'Dir':>4}  {'mean CV':>8}  {'mean CV2':>9}")
for d in range(6):
    cvs  = [elephant.statistics.cv(elephant.statistics.isi(st))
            for st in data['mo'][d] if len(st) > 2]
    cv2s = [elephant.statistics.cv2(elephant.statistics.isi(st))
            for st in data['mo'][d] if len(st) > 3]
    print(f"{d:>4}  {np.mean(cvs):>8.3f}  {np.mean(cv2s):>9.3f}")

 Dir   mean CV   mean CV2
   0     1.105      0.580
   1     1.358      0.672
   2     1.020      0.790
   3     1.326      0.565
   4     1.230      0.589
   5     1.333      0.611

Exercise: Scatter CV vs. CV2 per trial for all 6 directions pooled (MO-aligned, ≥4 spikes). Add the diagonal CV = CV2 and a reference line at CV2 = 1. Most trials should fall above the diagonal, confirming that CV is inflated by within-trial rate modulation.

Solution

cv_all, cv2_all = [], []
for d in range(6):
    for st in data['mo'][d]:
        if len(st) > 3:
            isi = elephant.statistics.isi(st)
            cv_all.append(elephant.statistics.cv(isi))
            cv2_all.append(elephant.statistics.cv2(isi))

utils.plot_cv_cv2_scatter(cv_all, cv2_all)
# Most trials fall above the diagonal (CV > CV2): CV is inflated by within-trial rate modulation.
# CV2 < 1 for most trials: local spiking is sub-Poisson, consistent with Gamma(kappa > 1).

Exercise: Compare three estimates of firing irregularity pooled across all 6 directions: (1) CV from the stable pre-movement epoch (-600 to -200 ms), (2) mean per-trial CV2, and (3) mean per-trial Lv. Print a summary table. Do the three estimates agree?

isis_stable_all = np.concatenate([
    utils.pool_isis_in_window(data['mo'][d], -600, -200) for d in range(6)
])
cv_stable_all = isis_stable_all.std() / isis_stable_all.mean()

cv2s_all, lvs_all = [], []
for d in range(6):
    for st in data['mo'][d]:
        if len(st) > 3:
            isi = elephant.statistics.isi(st)
            cv2s_all.append(elephant.statistics.cv2(isi))
            lvs_all.append(elephant.statistics.lv(isi))

print(f'{"Measure":<45}  {"Value":>6}')
print(f'{"CV  (stable epoch -600 to -200 ms, all dirs)":<45}  {cv_stable_all:>6.3f}')
print(f'{"CV2 (mean per trial, all dirs)":<45}  {np.mean(cv2s_all):>6.3f}')
print(f'{"Lv  (mean per trial, all dirs)":<45}  {np.mean(lvs_all):>6.3f}')
# All three estimates cluster around 0.6-0.8, consistent with a Gamma(kappa ~2-4) process.
# CV2 and Lv match the stationary-epoch CV without requiring an explicit window choice.

Measure                                         Value
CV  (stable epoch -600 to -200 ms, all dirs)    0.821
CV2 (mean per trial, all dirs)                  0.636
Lv  (mean per trial, all dirs)                  0.469

Avila-Akerberg O, Chacron MJ (2011) Nonrenewal spike train statistics: causes and functional consequences on neural coding. Exp Brain Res 210: 353-371.

Cox DR, Lewis PAW (1966) The statistical analysis of series of events. Chapman & Hall, London.

Farkhooi F, Strube-Bloss MF, Nawrot MP (2009) Serial correlation in neural spike trains. Phys Rev E 79: 021905.

Farkhooi F, Muller E, Nawrot MP (2011) Adaptation reduces variability of the neuronal population code. Phys Rev E 83: 050905.

Holt GR, Softky WR, Koch C, Douglas RJ (1996) Comparison of discharge variability in vitro and in vivo in cat visual cortex neurons. J Neurophysiol 75: 1806-1814.

Nawrot MP (2010) Analysis and interpretation of interval and count variability in neural spike trains. In: Grün S, Rotter S (eds) Analysis of parallel spike trains. Springer, New York, pp 34-58.

Ponce-Alvarez A, Kilavik BE, Riehle A (2010) Comparison of local measures of spike time irregularity and relating variability to firing rate in motor cortical neurons. J Comput Neurosci 29: 351-365.

Rickert J, Riehle A, Aertsen A, Rotter S, Nawrot MP (2009) Dynamic encoding of movement direction in motor cortical neurons. J Neurosci 29: 13870-13882.

Shadlen MN, Newsome WT (1998) The variable discharge of cortical neurons: implications for connectivity, computation, and information coding. J Neurosci 18: 3870-3896.

Softky WR, Koch C (1993) The highly irregular firing of cortical cells is inconsistent with temporal integration of random EPSPs. J Neurosci 13: 334-350.

Shinomoto S, Shima K, Tanji J (2003) Differences in spiking patterns among cortical neurons. Neural Comput 15: 2823-2842.

Further Reading

Nawrot MP et al. (2008) Measurements of variability dynamics in cortical spike trains. J Neurosci Meth 169: 374-390.

Riehle A, Vaadia E (eds, 2005) Motor Cortex in Voluntary Movements. CRC Press, Boca Raton.

Nawrot MP, Riehle A, Aertsen A, Rotter S (2003) Variability of motor cortical activity explained by network dynamics on multiple time scales. PhD thesis, Universität Freiburg.

Riehle A, Brochier T, Nawrot M, Grün S (2018) Behavioral context determines network state and variability dynamics in monkey motor cortex. Front Neural Circuits 12: 52.

Rostami V, Rost T, Schmitt FJ, van Albada SJ, Riehle A, Nawrot MP (2024) Spiking attractor model of motor cortex explains modulation of neural and behavioral variability by prior target information. Nat Commun 15: 6304.

Shadlen MN, Newsome WT (1998) The variable discharge of cortical neurons. J Neurosci 18: 3870-3896.