[pointprocesses] updated some figures

This commit is contained in:
Jan Benda 2021-01-09 23:23:47 +01:00
parent 10f71d5661
commit e5fd1ecb32
5 changed files with 101 additions and 85 deletions

View File

@ -55,11 +55,11 @@ def plotisih( ax, isis, binwidth=None ) :
binwidth = 5e-4 binwidth = 5e-4
h, b = np.histogram(isis, np.arange(0.0, np.max(isis)+binwidth, binwidth), density=True) h, b = np.histogram(isis, np.arange(0.0, np.max(isis)+binwidth, binwidth), density=True)
ax.text(0.9, 0.85, 'rate={:.0f}Hz'.format(1.0/np.mean(isis)), ha='right', transform=ax.transAxes) ax.text(0.9, 0.85, 'rate={:.0f}Hz'.format(1.0/np.mean(isis)), ha='right', transform=ax.transAxes)
ax.text(0.9, 0.75, 'mean={:.0f}ms'.format(1000.0*np.mean(isis)), ha='right', transform=ax.transAxes) ax.text(0.9, 0.7, 'mean={:.0f}ms'.format(1000.0*np.mean(isis)), ha='right', transform=ax.transAxes)
ax.text(0.9, 0.65, 'CV={:.2f}'.format(np.std(isis)/np.mean(isis)), ha='right', transform=ax.transAxes) ax.text(0.9, 0.55, 'CV={:.2f}'.format(np.std(isis)/np.mean(isis)), ha='right', transform=ax.transAxes)
ax.set_xlabel('ISI', 'ms') ax.set_xlabel('ISI', 'ms')
ax.set_ylabel('p(ISI)', '1/s') ax.set_ylabel('p(ISI)', '1/s')
ax.bar( 1000.0*b[:-1], h, bar_fac*1000.0*np.diff(b), facecolor=colors['blue']) ax.bar( 1000.0*b[:-1], h, bar_fac*1000.0*np.diff(b), **fsA)
# parameter: # parameter:
rate = 20.0 rate = 20.0
@ -85,16 +85,18 @@ x[x<0.0] = 0.0
inhspikes = pifspikes(x, trials, dt, D=0.3) inhspikes = pifspikes(x, trials, dt, D=0.3)
fig, (ax1, ax2) = plt.subplots(1, 2) fig, (ax1, ax2) = plt.subplots(1, 2)
fig.subplots_adjust(**adjust_fs(fig, top=1.5)) fig.subplots_adjust(**adjust_fs(fig, top=0.5, right=1.5))
ax1.set_title('stationary') ax1.set_xlim(0.0, 150.0)
ax1.set_xlim(0.0, 200.0) ax1.set_ylim(0.0, 31.0)
ax1.set_ylim(0.0, 40.0) ax1.set_xticks(np.arange(0.0, 151.0, 50.0))
plotisih(ax1, isis(homspikes)) ax1.set_yticks(np.arange(0.0, 31.0, 10.0))
plotisih(ax1, isis(homspikes), 0.005)
ax2.set_title('non-stationary') ax2.set_xlim(0.0, 150.0)
ax2.set_xlim(0.0, 200.0) ax2.set_ylim(0.0, 31.0)
ax2.set_ylim(0.0, 40.0) ax2.set_xticks(np.arange(0.0, 151.0, 50.0))
plotisih(ax2, isis(inhspikes)) ax2.set_yticks(np.arange(0.0, 31.0, 10.0))
plotisih(ax2, isis(inhspikes), 0.005)
plt.savefig('isihexamples.pdf') plt.savefig('isihexamples.pdf')
plt.close() plt.close()

View File

@ -1,36 +1,41 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\chapter{Spiketrain analysis} \chapter{Spike-train analysis}
\label{pointprocesseschapter} \label{pointprocesseschapter}
\exercisechapter{Spiketrain analysis} \exercisechapter{Spike-train analysis}
\entermde[action potential]{Aktionspotential}{Action potentials} \entermde[action potential]{Aktionspotential}{Action potentials}
(\enterm[spike|seealso{action potential}]{spikes}) are the carriers of (\enterm[spike|seealso{action potential}]{spikes}) carry information
information in the nervous system. Thereby it is the time at which the within neural systems. More precisely, the times at which action
spikes are generated that is of importance for information potentials are generated contain the information. The waveform of the
transmission. The waveform of the action potential is largely action potential is largely stereotyped and therefore conveys no
stereotyped and therefore does not carry information. information. Analyzing the statistics of spike times and their
relation to sensory stimuli or motor actions is central to
The result of the pre-processing of electrophysiological recordings are neuroscientific research. With multi-electrode arrays it is nowadays
series of spike times, which are termed \enterm{spiketrains}. If possible to record from hundreds or even thousands of neurons
measurements are repeated we get several \enterm{trials} of simultaneously. The open challenge is how to analyze such data sets in
spiketrains (\figref{rasterexamplesfig}). order to understand how neural systems work. Let's start with the
basics in this chapter.
Spiketrains are times of events, the action potentials. Analyzing
spike trains leads into the realm of the so called \entermde[point The result of the pre-processing of electrophysiological recordings
are series of spike times, which are termed \enterm[spike train]{spike
trains}. If measurements are repeated we get several \enterm{trials}
of spike trains (\figref{rasterexamplesfig}). Spike trains are lists
of times of events, the action potentials. Analyzing spike trains
leads into the realm of the statistics of so called \entermde[point
process]{Punktprozess}{point processes}. process]{Punktprozess}{point processes}.
\begin{figure}[ht] \begin{figure}[bt]
\includegraphics[width=1\textwidth]{rasterexamples} \includegraphics[width=1\textwidth]{rasterexamples}
\titlecaption{\label{rasterexamplesfig}Raster-plot.}{Raster-plots of \titlecaption{\label{rasterexamplesfig}Raster plot.}{Raster plots of
ten trials of data illustrating the times of the action ten trials of data illustrating the times of action
potentials. Each vertical dash illustrates the time at which an potentials. Each vertical stroke illustrates the time at which an
action potential was observed. Each line displays the events of action potential was observed. Each row displays the events of one
one trial. Shown is a stationary point process (homogeneous point trial. Shown is a stationary point process (homogeneous point
process with a rate $\lambda=20$\;Hz, left) and an non-stationary process with a rate $\lambda=20$\;Hz, left) and an non-stationary
point process (perfect integrate-and-fire neuron driven by point process with a rate that varies in time (noisy perfect
Ohrnstein-Uhlenbeck noise with a time-constant $\tau=100$\,ms, integrate-and-fire neuron driven by Ohrnstein-Uhlenbeck noise with
right).} a time-constant $\tau=100$\,ms, right).}
\end{figure} \end{figure}
@ -43,7 +48,7 @@ spike trains leads into the realm of the so called \entermde[point
crosses some threshold. For example:\vspace{-1ex} crosses some threshold. For example:\vspace{-1ex}
\begin{itemize} \begin{itemize}
\item Action potentials/heart beat: created by the dynamics of the \item Action potentials/heart beat: created by the dynamics of the
neuron/sinoatrial node neuron/sinoatrial node.
\item Earthquake: defined by the dynamics of the pressure between \item Earthquake: defined by the dynamics of the pressure between
tectonical plates. tectonical plates.
\item Communication calls in crickets/frogs/birds: shaped by \item Communication calls in crickets/frogs/birds: shaped by
@ -59,21 +64,21 @@ spike trains leads into the realm of the so called \entermde[point
$T_i=t_{i+1}-t_i$ and the number of events $n_i$. } $T_i=t_{i+1}-t_i$ and the number of events $n_i$. }
\end{figure} \end{figure}
\noindent
A temporal \entermde{Punktprozess}{point process} is a stochastic A temporal \entermde{Punktprozess}{point process} is a stochastic
process that generates a sequence of events at times $\{t_i\}$, $t_i process that generates a sequence of events at times $\{t_i\}$. In
\in \reZ$. In the neurosciences, the statistics of point processes is the neurosciences, the statistics of point processes is of importance
of importance since the timing of neuronal events (action potentials, since the timing of neuronal events (action potentials, post-synaptic
post-synaptic potentials, events in EEG or local-field recordings, potentials, events in EEG or local-field recordings, etc.) is crucial
etc.) is crucial for information transmission and can be treated as for information transmission and can be treated as such a process.
such a process.
The events of a point process can be illustrated by means of a raster The events of a point process can be illustrated by means of a raster
plot in which each vertical line indicates the time of an event. The plot in which each vertical line indicates the time of an event. The
event from two different point processes are shown in event from two different point processes are shown in
\figref{rasterexamplesfig}. Point processes can be described using \figref{rasterexamplesfig}. In addition to the event times, point
the intervals between successive events $T_i=t_{i+1}-t_i$ and the processes can be described using the intervals $T_i=t_{i+1}-t_i$
number of observed events within a certain time window $n_i$ between successive events or the number of observed events within a
(\figref{pointprocessscetchfig}). certain time window $n_i$ (\figref{pointprocessscetchfig}).
\begin{exercise}{rasterplot.m}{} \begin{exercise}{rasterplot.m}{}
Implement a function \varcode{rasterplot()} that displays the times of Implement a function \varcode{rasterplot()} that displays the times of
@ -94,7 +99,7 @@ real-valued variables:
\begin{figure}[t] \begin{figure}[t]
\includegraphics[width=0.96\textwidth]{isihexamples}\vspace{-2ex} \includegraphics[width=0.96\textwidth]{isihexamples}\vspace{-2ex}
\titlecaption{\label{isihexamplesfig}Interspike interval \titlecaption{\label{isihexamplesfig}Interspike-interval
histograms}{of the spike trains shown in \figref{rasterexamplesfig}.} histograms}{of the spike trains shown in \figref{rasterexamplesfig}.}
\end{figure} \end{figure}
@ -144,7 +149,7 @@ the interval $T_i$. The parameter $k$ is called the \enterm{lag}
(\determ{Verz\"ogerung}) $k$. Stationary and non-stationary return (\determ{Verz\"ogerung}) $k$. Stationary and non-stationary return
maps are distinctly different \figref{returnmapfig}. maps are distinctly different \figref{returnmapfig}.
\begin{figure}[t] \begin{figure}[tp]
\includegraphics[width=1\textwidth]{returnmapexamples} \includegraphics[width=1\textwidth]{returnmapexamples}
\includegraphics[width=1\textwidth]{serialcorrexamples} \includegraphics[width=1\textwidth]{serialcorrexamples}
\titlecaption{\label{returnmapfig}Interspike interval \titlecaption{\label{returnmapfig}Interspike interval

View File

@ -1,5 +1,6 @@
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from plotstyle import *
def hompoisson(rate, trials, duration) : def hompoisson(rate, trials, duration) :
spikes = [] spikes = []
@ -64,23 +65,22 @@ x[x<0.0] = 0.0
# pif spike trains: # pif spike trains:
inhspikes = pifspikes(x, trials, dt, D=0.3) inhspikes = pifspikes(x, trials, dt, D=0.3)
fig = plt.figure( figsize=(9,4) ) fig, (ax1, ax2) = plt.subplots(1, 2, figsize=cm_size(figure_width, 0.5*figure_width))
ax = fig.add_subplot(1, 2, 1) fig.subplots_adjust(**adjust_fs(fig, left=4.0, right=1.0, top=1.2))
ax.set_title('stationary')
ax.set_xlim(0.0, duration)
ax.set_ylim(-0.5, trials-0.5)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Trials')
ax.eventplot(homspikes, colors=[[0, 0, 0]], linelength=0.8)
ax = fig.add_subplot(1, 2, 2) ax1.set_title('stationary')
ax.set_title('non-stationary') ax1.set_xlim(0.0, duration)
ax.set_xlim(0.0, duration) ax1.set_ylim(-0.5, trials-0.5)
ax.set_ylim(-0.5, trials-0.5) ax1.set_xlabel('Time [s]')
ax.set_xlabel('Time [s]') ax1.set_ylabel('Trial')
ax.set_ylabel('Trials') ax1.eventplot(homspikes, colors=[lsA['color']], linelength=0.8, lw=1)
ax.eventplot(inhspikes, colors=[[0, 0, 0]], linelength=0.8)
ax2.set_title('non-stationary')
ax2.set_xlim(0.0, duration)
ax2.set_ylim(-0.5, trials-0.5)
ax2.set_xlabel('Time [s]')
ax2.set_ylabel('Trial')
ax2.eventplot(inhspikes, colors=[lsA['color']], linelength=0.8, lw=1)
plt.tight_layout()
plt.savefig('rasterexamples.pdf') plt.savefig('rasterexamples.pdf')
plt.close() plt.close()

View File

@ -40,19 +40,21 @@ def pifspikes(input, trials, dt, D=0.1) :
spikes.append( times ) spikes.append( times )
return spikes return spikes
def isis( spikes ) : def isis( spikes ) :
isi = [] isi = []
for k in range(len(spikes)) : for k in range(len(spikes)) :
isi.extend(np.diff(spikes[k])) isi.extend(np.diff(spikes[k]))
return np.array( isi ) return np.array( isi )
def plotreturnmap(ax, isis, lag=1, max=None) :
def plotreturnmap(ax, isis, lag=1, max=1.0) :
ax.set_xlabel(r'ISI$_i$', 'ms') ax.set_xlabel(r'ISI$_i$', 'ms')
ax.set_ylabel(r'ISI$_{i+1}$', 'ms') ax.set_ylabel(r'ISI$_{i+1}$', 'ms')
if max != None :
ax.set_xlim(0.0, 1000.0*max) ax.set_xlim(0.0, 1000.0*max)
ax.set_ylim(0.0, 1000.0*max) ax.set_ylim(0.0, 1000.0*max)
ax.scatter(1000.0*isis[:-lag], 1000.0*isis[lag:], c=colors['blue']) isiss = isis[isis<max]
ax.plot(1000.0*isiss[:-lag], 1000.0*isiss[lag:], clip_on=False, **psAm)
# parameter: # parameter:
rate = 20.0 rate = 20.0
@ -79,11 +81,14 @@ inhspikes = pifspikes(x, trials, dt, D=0.3)
fig, (ax1, ax2) = plt.subplots(1, 2) fig, (ax1, ax2) = plt.subplots(1, 2)
fig.subplots_adjust(**adjust_fs(fig, left=6.5, top=1.5)) fig.subplots_adjust(**adjust_fs(fig, left=6.5, top=1.5))
ax1.set_title('stationary')
plotreturnmap(ax1, isis(homspikes), 1, 0.3) plotreturnmap(ax1, isis(homspikes), 1, 0.3)
ax1.set_xticks(np.arange(0.0, 301.0, 100.0))
ax1.set_yticks(np.arange(0.0, 301.0, 100.0))
ax2.set_title('non-stationary')
plotreturnmap(ax2, isis(inhspikes), 1, 0.3) plotreturnmap(ax2, isis(inhspikes), 1, 0.3)
ax2.set_ylabel('')
ax2.set_xticks(np.arange(0.0, 301.0, 100.0))
ax2.set_yticks(np.arange(0.0, 301.0, 100.0))
plt.savefig('returnmapexamples.pdf') plt.savefig('returnmapexamples.pdf')
plt.close() plt.close()

View File

@ -13,6 +13,7 @@ def hompoisson(rate, trials, duration) :
spikes.append( times ) spikes.append( times )
return spikes return spikes
def inhompoisson(rate, trials, dt) : def inhompoisson(rate, trials, dt) :
spikes = [] spikes = []
p = rate*dt p = rate*dt
@ -55,7 +56,9 @@ def plotserialcorr(ax, isis, maxlag=10) :
ax.set_ylabel(r'ISI correlation $\rho_k$') ax.set_ylabel(r'ISI correlation $\rho_k$')
ax.set_xlim(0.0, maxlag) ax.set_xlim(0.0, maxlag)
ax.set_ylim(-1.0, 1.0) ax.set_ylim(-1.0, 1.0)
ax.plot(lags, corr, '.-', markersize=15, c=colors['blue']) ax.plot([0, 10], [0.0, 0.0], **lsGrid)
ax.plot(lags, corr, clip_on=False, zorder=100, **lpsAm)
# parameter: # parameter:
rate = 20.0 rate = 20.0
@ -87,6 +90,7 @@ plotserialcorr(ax1, isis(homspikes))
ax1.set_ylim(-0.2, 1.0) ax1.set_ylim(-0.2, 1.0)
plotserialcorr(ax2, isis(inhspikes)) plotserialcorr(ax2, isis(inhspikes))
ax2.set_ylabel('')
ax2.set_ylim(-0.2, 1.0) ax2.set_ylim(-0.2, 1.0)
plt.savefig('serialcorrexamples.pdf') plt.savefig('serialcorrexamples.pdf')