From 1175dc3ab98aa738c2469c14620e44e6f8f8601c Mon Sep 17 00:00:00 2001
From: Jan Benda <jan.benda@uni-tuebingen.de>
Date: Mon, 23 Jan 2017 13:30:51 +0100
Subject: [PATCH] added lif project
---
projects/project_lif/Makefile | 10 ++
projects/project_lif/lif.tex | 160 ++++++++++++++++++
projects/project_lif/solution/lif.m | 116 +++++++++++++
projects/project_lif/solution/lifspikes.m | 14 ++
.../project_lif/solution/passivemembrane.m | 8 +
5 files changed, 308 insertions(+)
create mode 100644 projects/project_lif/Makefile
create mode 100644 projects/project_lif/lif.tex
create mode 100644 projects/project_lif/solution/lif.m
create mode 100644 projects/project_lif/solution/lifspikes.m
create mode 100644 projects/project_lif/solution/passivemembrane.m
diff --git a/projects/project_lif/Makefile b/projects/project_lif/Makefile
new file mode 100644
index 0000000..6422eb4
--- /dev/null
+++ b/projects/project_lif/Makefile
@@ -0,0 +1,10 @@
+latex:
+ pdflatex *.tex > /dev/null
+ pdflatex *.tex > /dev/null
+
+clean:
+ rm -rf *.log *.aux *.zip *.out auto
+ rm -f `basename *.tex .tex`.pdf
+
+zip: latex
+ zip `basename *.tex .tex`.zip *.pdf *.dat *.mat *.m
diff --git a/projects/project_lif/lif.tex b/projects/project_lif/lif.tex
new file mode 100644
index 0000000..e42d352
--- /dev/null
+++ b/projects/project_lif/lif.tex
@@ -0,0 +1,160 @@
+\documentclass[addpoints,11pt]{exam}
+\usepackage{url}
+\usepackage{color}
+\usepackage{hyperref}
+
+\pagestyle{headandfoot}
+\runningheadrule
+\firstpageheadrule
+\firstpageheader{Scientific Computing}{Project Assignment}{11/05/2014
+ -- 11/06/2014}
+%\runningheader{Homework 01}{Page \thepage\ of \numpages}{23. October 2014}
+\firstpagefooter{}{}{}
+\runningfooter{}{}{}
+\pointsinmargin
+\bracketedpoints
+
+%%%%% listings %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\usepackage{listings}
+\lstset{
+ basicstyle=\ttfamily,
+ numbers=left,
+ showstringspaces=false,
+ language=Matlab,
+ breaklines=true,
+ breakautoindent=true,
+ columns=flexible,
+ frame=single,
+% captionpos=t,
+ xleftmargin=2em,
+ xrightmargin=1em,
+% aboveskip=11pt,
+ %title=\lstname,
+% title={\protect\filename@parse{\lstname}\protect\filename@base.\protect\filename@ext}
+ }
+
+
+%\printanswers
+%\shadedsolutions
+
+
+\begin{document}
+%%%%%%%%%%%%%%%%%%%%% Submission instructions %%%%%%%%%%%%%%%%%%%%%%%%%
+\sffamily
+% \begin{flushright}
+% \gradetable[h][questions]
+% \end{flushright}
+
+\begin{center}
+ \input{../disclaimer.tex}
+\end{center}
+
+%%%%%%%%%%%%%% Questions %%%%%%%%%%%%%%%%%%%%%%%%%
+
+\begin{questions}
+ \question The temporal evolution of the membrane voltage $V(t)$ of a
+ passive neuron is described by the membrane equation
+ \begin{equation}
+ \label{passivemembrane}
+ \tau \frac{dV}{dt} = -V + E
+ \end{equation}
+ where $\tau=10$\,ms is the membrane time constant and $E(t)$ is the
+ reversal potential that also depends on time $t$.
+
+ Such a differential equation can be numerically solved with the Euler method.
+ For this the time is discretized by a time step $\Delta t=0.1$\,ms.
+ The $i$-th time point is then at time $t_i = i \cdot \Delta t$.
+ In matlab we get the time points $t_i$ simply by
+ \begin{lstlisting}
+dt = 0.1;
+tmax = 100.0;
+time = [0.0:dt:tmax]; % t_i
+ \end{lstlisting}
+ When the membrane potential at time $t_0 = 0$ is $V_0$, the so
+ called ``initial condition'', then we can iteratively compute the
+ membrane potentials $V_i$ for successive time points $t_i$ according to
+ \begin{equation}
+ \label{euler}
+ V_{i+1} = V_i + (-V_i + E_i) \frac{\Delta t}{\tau}
+ \end{equation}
+
+ \begin{parts}
+ \part Write a function that computes the time course of the
+ membrane potential of the passive membrane. The function gets as
+ input arguments the initial condition $V_0$, the vector with the
+ time course of $E(t)$, the value of the membrane time-constant
+ $\tau$, and the time step $\Delta t$.
+
+ \part In order to test your function set $V_0=1$\,mV and $E(t)=0$
+ and compute $V(t)$ for $t_{max}=50$\,ms. Plot $V(t)$ and compare it to
+ the expected result of $V(t) = \exp(-t/\tau)$.
+
+ Why is $V=0$ the resting potential of this neuron?
+
+ \part Response of the passive membrane to a step input.
+
+ Set $V_0=0$. Construct a vector for the input $E(t)$ such that
+ $E(t)=0$ for $t<20$\,ms and $t>70$\,ms and $E(t)=10$\,mV for
+ $20$\,ms $<t<70$\,ms. Plot $E(t)$ and the resulting $V(t)$ for
+ $t_{max}=120$\,ms.
+
+ \part Response to sine waves.
+
+ As an input we now use $E(t)=\sin(2\pi f t)$. Compute the time
+ course of the membrane potential in response to this input
+ ($t_{max}=1$\,s). Vary the frequency $f$ between 1 and 100\,Hz. Be
+ careful with the units within the sine function --- $ft$ must be
+ unitless.
+
+ What do you observe?
+
+ \part Transfer function of the passive neuron.
+
+ Measure the amplitude of the voltage responses evoked by the
+ sinusoidal inputs as the maximum of the last 900\,ms of the
+ responses. Plot the amplitude of the response as a function of
+ input frequency. This is the transfer function of the passive neuron.
+
+ How does the transfer function depend on the membrane time
+ constant?
+
+ \part Leaky integrate-and-fire neuron.
+
+ The passive neuron can be turned into a spiking neuron by
+ introducing a fixed voltage threshold. Whenever the computed
+ membrane potential of the passive neuron crosses the voltage
+ threshold a spike is generated and the membrane voltage is set to
+ the reset potential $V_R$ that we here set to zero. ``Generating a
+ spike'' only means that we note down the time of the threshold
+ crossing as a time where an action potential occurred. The
+ waveform of the action potential is not modeled. Here we use a
+ voltage threshold of one.
+
+ Write a function that implements this leaky integrate-and-fire
+ neuron by expanding the function for the passive neuron
+ appropriate. The function returns a vector of spike times.
+
+ Illustrate how this model works by appropriate plot(s) and
+ input(s) $E(t)$, e.g. constant inputs lower and higher than the
+ voltage threshold.
+
+ \part Show the response of the leaky integrate-and-fire neuron to
+ a sine wave $E(t)=A\sin(2\pi ft)$ with $A=2$\,mV and frequency
+ $f=10$, 20, and 30\,Hz.
+
+ \part Compute the firing rate as a function of the frequency of
+ the stimulating sine wave ($A=2$\,mV and frequencies between 5 and
+ 30\,Hz). For a spike train with $n$ spikes at times $t_k$ ($k=1,
+ 2, \ldots n$) the firing rate is
+ \begin{equation}
+ \label{firingrate}
+ r = \frac{n-1}{t_n - t_1}
+ \end{equation}
+
+ What do you observe? Does the firing rate encode the frequency of
+ the stimulus?
+ \end{parts}
+
+\end{questions}
+
+\end{document}
diff --git a/projects/project_lif/solution/lif.m b/projects/project_lif/solution/lif.m
new file mode 100644
index 0000000..f287491
--- /dev/null
+++ b/projects/project_lif/solution/lif.m
@@ -0,0 +1,116 @@
+%% general settings:
+tau = 10.0;
+dt = 0.1;
+
+%% test passive membrane:
+tmax = 50.0;
+time = [0:dt:tmax];
+E = zeros(length(time), 1);
+V = passivemembrane(1.0, E, tau, dt);
+expfun = exp(-time/tau);
+figure()
+plot(time, V, 'b', 'linewidth', 2);
+hold on;
+plot(time, expfun, 'r');
+hold off;
+
+%% step input:
+tmax = 120.0;
+time = [0:dt:tmax];
+E = zeros(length(time), 1);
+E(20/dt:70/dt) = 10;
+V = passivemembrane(0.0, E, tau, dt);
+figure()
+plot(time, E, 'r');
+hold on;
+plot(time, V, 'b', 'linewidth', 2);
+hold off;
+
+%% sine waves:
+tmax = 1000.0;
+time = [0:dt:tmax];
+freqs = [1.0 10.0 30.0 100.0];
+figure();
+for k = 1:length(freqs)
+ f = freqs(k);
+ E = sin(2*pi*0.001*f*time);
+ V = passivemembrane(0.0, E, tau, dt);
+ subplot(4, 1, k);
+ plot(time, E, 'r');
+ hold on;
+ plot(time, V, 'b', 'linewidth', 2);
+ hold off;
+end
+
+
+%% transfer function:
+tmax = 1000.0;
+time = [0:dt:tmax];
+taus = [3.0 10.0 30.0];
+figure();
+for tau = taus
+ freqs = [1.0:1.0:100.0];
+ rates = zeros(length(freqs), 1);
+ for k = 1:length(freqs)
+ f = freqs(k);
+ E = sin(2*pi*0.001*f*time);
+ V = passivemembrane(0.0, E, tau, dt);
+ rates(k) = max(V(100/dt:end));
+ end
+ plot(freqs, rates);
+ hold on;
+end
+hold off;
+
+
+%% leaky IaF:
+tau = 10.0;
+tmax = 50.0;
+time = [0:dt:tmax];
+E = zeros(length(time), 1) + 1.5;
+[spikes, V] = lifspikes(0.0, E, tau, dt);
+figure()
+plot(time, V, 'b', 'linewidth', 2);
+hold on;
+plot(spikes, ones(length(spikes), 1), 'or');
+hold off;
+
+
+%% leaky IaF and sine input:
+tau = 10.0;
+tmax = 500.0;
+time = [0:dt:tmax];
+f = 10.0;
+E = 2.0*sin(2*pi*0.001*f*time);
+[spikes, V] = lifspikes(0.0, E, tau, dt);
+figure()
+subplot(2, 1, 1);
+plot(time, V, 'b', 'linewidth', 2);
+hold on;
+plot(spikes, ones(length(spikes), 1), 'or');
+hold off;
+
+
+%% transfer function of LIF spikes:
+tmax = 1000.0;
+time = [0:dt:tmax];
+taus = [3.0 10.0];
+figure();
+for tau = taus
+ freqs = [5.0:0.1:30.0];
+ rates = zeros(length(freqs), 1);
+ for k = 1:length(freqs)
+ f = freqs(k);
+ E = 2.0*sin(2*pi*0.001*f*time);
+ [spikes, V] = lifspikes(0.0, E, tau, dt);
+ spikes = spikes(spikes>100.0);
+ if length(spikes) > 2
+ rates(k) = 1000.0*(length(spikes)-1)/(spikes(end)-spikes(1));
+ else
+ rates(k) = 0.0;
+ end
+ end
+ plot(freqs, rates);
+ hold on;
+end
+hold off;
\ No newline at end of file
diff --git a/projects/project_lif/solution/lifspikes.m b/projects/project_lif/solution/lifspikes.m
new file mode 100644
index 0000000..328b927
--- /dev/null
+++ b/projects/project_lif/solution/lifspikes.m
@@ -0,0 +1,14 @@
+function [spikes, voltage] = lifspikes(V0, E, tau, dt)
+voltage = zeros(length(E), 1);
+V = V0;
+thresh = 1.0;
+spikes = [];
+for k = 1:length(E)
+ voltage(k) = V;
+ if V > thresh
+ spikes = [spikes; k*dt];
+ V = 0.0;
+ end
+ V = V + (-V+E(k))*dt/tau;
+end
+end
diff --git a/projects/project_lif/solution/passivemembrane.m b/projects/project_lif/solution/passivemembrane.m
new file mode 100644
index 0000000..b6c2a2d
--- /dev/null
+++ b/projects/project_lif/solution/passivemembrane.m
@@ -0,0 +1,8 @@
+function voltage = passivemembrane(V0, E, tau, dt)
+voltage = zeros(length(E), 1);
+V = V0;
+for k = 1:length(E)
+ voltage(k) = V;
+ V = V + (-V+E(k))*dt/tau;
+end
+end