AT2 – Neuromodeling: Problem set #3 SPIKE TRAINS

PROBLEM 1: Poisson spike trains.

Our tour of the brain starts by looking more closely at the brain's basic signaling elements, action potentials or "spikes", and by constructing some simple tools that will allow us to deal with them.

In [1]:
# /!\ I'm using Python 3

import plotly
import numpy as np
import plotly.plotly as py
import plotly.graph_objs as go

import matplotlib.pyplot as plt
import scipy.io as sio
from scipy import stats
%matplotlib inline

plotly.offline.init_notebook_mode(connected=True)
from IPython.core.display import display, HTML, Markdown
# The polling here is to ensure that plotly.js has already been loaded before
# setting display alignment in order to avoid a race condition.
display(HTML(
    '<script>'
        'var waitForPlotly = setInterval( function() {'
            'if( typeof(window.Plotly) !== "undefined" ){'
                'MathJax.Hub.Config({ SVG: { font: "STIX-Web" }, displayAlign: "center" });'
                'MathJax.Hub.Queue(["setRenderer", MathJax.Hub, "SVG"]);'
                'clearInterval(waitForPlotly);'
            '}}, 250 );'
    '</script>'
))

# Colorscales
def colorscale_list(cmap, number_colors, return_rgb_only=False):
    cm = plt.get_cmap(cmap)
    colors = [np.array(cm(i/number_colors)) for i in range(1, number_colors+1)]
    rgb_colors_plotly = []
    rgb_colors_only = []
    for i, c in enumerate(colors):
        col = 'rgb{}'.format(tuple(255*c[:-1]))
        rgb_colors_only.append(col)
        rgb_colors_plotly.append([i/number_colors, col])
        rgb_colors_plotly.append([(i+1)/number_colors, col])
    return rgb_colors_only if return_rgb_only else rgb_colors_plotly

In many neurons of the brain, the omission of individual spikes appears to occur almost at random, similar to the clicks in a Geiger counter. We will therefore start by creating some artificial, random spike trains. We will quantify such a spike train as a series of $0$s (no spikes) and $1$s (spike).

(a) To begin, create a vector of $1000$ elements so that, on average but in the most irregular manner, every fourth element in the vector is a spike. Plot the vector, by using one dot for each spike, for example, or the matplotlib function eventplot.

In [2]:
N = 1000
spikes = np.zeros(N)


spikes[np.random.rand(N) < 1/4] = 1

trace_sp = go.Bar(
    y = spikes,
    opacity=0.75,
    name = 'Spikes'
)

layout_bar = go.Layout(
    title='Spike train',
    xaxis=dict(
        title='Time step'
    ),
    yaxis=dict(
        title='Spike',
        tick0=0,
        dtick=1
    )
)

plotly.offline.iplot(
    go.Figure(
        data=[trace_sp], 
        layout=layout_bar)
)

(b) In the next step, we want to introduce time units. We will associate every $0$ or $1$ with a time bin of length $∆t$ msec. Choose $∆t = 2$ msec and create a spike train of length $1$ sec with an average rate of $25$ spikes/sec. (Careful with the units!) Plot the spike train similar to above, but now use the correct time axis.

In [3]:
delta_t = 2e-3
t_max = 1
rate = 25 # number of spikes/sec

spikes = np.zeros(int(t_max/delta_t))

spikes[np.random.rand(int(t_max/delta_t)) < rate*delta_t] = 1




trace_sp = go.Bar(
    x = 2*np.arange(len(spikes)),
    y = spikes,
    opacity=0.75,
    name = 'Spikes'
)

layout_bar = go.Layout(
    title='Spike train',
    xaxis=dict(
        range = [0, 2*int(t_max/delta_t)],
        title='Time (ms)',
        tick0=0,
        dtick=200
    ),
    yaxis=dict(
        title='Spike',
        tick0=0,
        dtick=1
    )
)

plotly.offline.iplot(
    go.Figure(
        data=[trace_sp], 
        layout=layout_bar)
)

(c) Generate $N = 200$ spike trains with firing rate $25$ Hz and count the total number of spikes in each of them. Plot $50$ of these trials as a rastergram. Plot a histogram of the spike counts.

In [27]:
N = 100000
delta_t = 2e-3
t_max = 1
rate = 25 # number of spikes/sec

spikes = np.zeros((N, int(t_max/delta_t)))
spikes[np.random.rand(N,int(t_max/delta_t)) < rate*delta_t] = 1


# Histogram of spike counts

trace_hist = go.Histogram(
    x=np.sum(spikes,axis=1),
    opacity=0.75,
)

trace_hist2 = go.Histogram(
    x=np.sum(spikes,axis=1),
    histnorm='probability',
    opacity=0.75,
    name='Empirical distribution'
)

xs = np.arange(stats.poisson.ppf(0.0001, rate*t_max), stats.poisson.ppf(0.9999, rate*t_max))

trace_poiss = go.Scatter(
    x=xs,
    y=stats.poisson.pmf(xs, rate*t_max),
    name='Poisson distribution with mean {}'.format(rate*t_max)
)

layout_hist = go.Layout(
    title='Histogram of spike counts',
    xaxis=dict(
        title='Number of spikes'
    ),
    yaxis=dict(
        title='Number of spike trains having this number of spikes'
    )
)

layout_comparison = go.Layout(
    title='Probability distribution of the number of spike counts',
    xaxis=dict(
        title='Number of spikes'
    ),
    yaxis=dict(
        title='Probability'
    )
)

fig = go.Figure(data=[trace_hist], layout=layout_hist)
plotly.offline.iplot(fig)

fig2 = go.Figure(data=[trace_hist2, trace_poiss], layout=layout_comparison)
plotly.offline.iplot(fig2)
In [33]:
N = 200
delta_t = 2e-3
t_max = 1
rate = 25 # number of spikes/sec

spikes = np.zeros((N, int(t_max/delta_t)))
spikes[np.random.rand(N,int(t_max/delta_t)) < rate*delta_t] = 1


# Histogram of spike counts

trace_hist = go.Histogram(
    x=np.sum(spikes,axis=1),
    opacity=0.75,
)

trace_hist2 = go.Histogram(
    x=np.sum(spikes,axis=1),
    histnorm='probability',
    opacity=0.75,
    name='Empirical distribution'
)

xs = np.arange(stats.poisson.ppf(0.001, rate*t_max), stats.poisson.ppf(0.999, rate*t_max))

trace_poiss = go.Scatter(
    x=xs,
    y=stats.poisson.pmf(xs, rate*t_max),
    name='Poisson distribution with mean {}'.format(rate*t_max)
)

layout_hist = go.Layout(
    title='Histogram of spike counts',
    xaxis=dict(
        title='Number of spikes'
    ),
    yaxis=dict(
        title='Number of spike trains having this number of spikes'
    )
)

layout_comparison = go.Layout(
    title='Probability distribution of the number of spike counts',
    xaxis=dict(
        title='Number of spikes'
    ),
    yaxis=dict(
        title='Probability'
    )
)

fig = go.Figure(data=[trace_hist], layout=layout_hist)
plotly.offline.iplot(fig)

fig2 = go.Figure(data=[trace_hist2, trace_poiss], layout=layout_comparison)
plotly.offline.iplot(fig2)
In [29]:
# Rastergram
traces_spikes = []

for i, spike in enumerate(spikes[:50]):
    traces_spikes.append(
        go.Scatter(
            x = 2*np.arange(len(spike)),
            y = (i+1)*spike,
            mode = 'markers',
            name = 'spike #{}'.format(i+1),
            marker = dict(
                symbol = 'square'
            ),
            showlegend = False
        )
    )

layout = go.Layout(
    title='Rastergram of 50 spike trains with {}Hz firing rate'.format(rate),
    xaxis=dict(
        title='Time (ms)'
    ),
    yaxis=dict(
        range=[0.5, i+2],
        title='Trial',
        ticklen= 5,
        gridwidth= 2,
        tick0=1
    )
)

plotly.offline.iplot(go.Figure(data=traces_spikes, layout=layout))

(d) Advanced: Compute the histogram of interspike intervals for the above spike trains. Verify that the histogram of spike counts follows a normal distribution (why?) and the histogram of interspike intervals an exponential distribution.

In [47]:
N = 200
delta_t = 2e-3
t_max = 1
rate = 25 # number of spikes/sec

spikes = np.zeros((N, int(t_max/delta_t)))
spikes[np.random.rand(N,int(t_max/delta_t)) < rate*delta_t] = 1


# Histogram of interspike intervals

trace_hist = go.Histogram(
    x=delta_t*1000*np.hstack([np.diff(np.where(s==1)) for s in spikes])[0],
    opacity=0.75,
)

layout_hist = go.Layout(
    title='Histogram of interspike intervals for {} spike trains with {}Hz firing rate'.format(N, rate),
    xaxis=dict(
        title='Interspike interval (ms)'
    ),
    yaxis=dict(
        title='Number of times this interspike interval appears in the spike trains'
    )
)


fig = go.Figure(data=[trace_hist], layout=layout_hist)

plotly.offline.iplot(fig)