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