improved noiseficurves projects

projects/project_noiseficurves/lifspikes.m

@@ -1,9 +1,8 @@
-function spikes = lifspikes( trials, input, tmaxdt, D )
+function spikes = lifspikes(trials, input, tmax, D)
 % Generate spike times of a leaky integrate-and-fire neuron
 % trials: the number of trials to be generated
 % input: the stimulus either as a single value or as a vector
-% tmaxdt: in case of a single value stimulus the duration of a trial
+% tmax: duration of a trial
-%         in case of a vector as a stimulus the time step
 % D: the strength of additive white noise
     
     tau = 0.01;
@@ -14,19 +13,15 @@ function spikes = lifspikes( trials, input, tmaxdt, D )
     vthresh = 10.0;
     dt = 1e-4;
     
-    if length( input ) == 1
+    n = ceil(tmax/dt);
-        input = input * ones( ceil( tmaxdt/dt ), 1 );
+    spikes = cell(trials, 1);
-    else
-        dt = tmaxdt;
-    end
-    spikes = cell( trials, 1 );
     for k=1:trials
         times = [];
         j = 1;
         v = vreset;
-        noise = sqrt(2.0*D)*randn( length( input ), 1 )/sqrt(dt);
+        noise = sqrt(2.0*D)*randn(n, 1)/sqrt(dt);
-        for i=1:length( noise )
+        for i=1:n
-            v = v + ( - v + noise(i) + input(i))*dt/tau;
+            v = v + (- v + noise(i) + input)*dt/tau;
             if v >= vthresh
                 v = vreset;
                 times(j) = i*dt;

projects/project_noiseficurves/noiseficurves.tex

@@ -57,8 +57,10 @@
 
   Measure the tuning curve (also called the intensity-response curve) of the
   neuron.  That is, what is the mean firing rate of the neuron's response
-  as a function of the input $I$. How does this depend on the level of
+  as a function of the input $I$?
-  the intrinsic noise of the neuron?
+
+  How does the intensity-response curve of a neuron depend on the
+  level of the intrinsic noise of the neuron?
 
   The neuron is implemented in the file \texttt{lifspikes.m}.  Call it
   with the following parameters:
@@ -67,39 +69,56 @@ trials = 10;
 tmax = 50.0;
 input = 10.0;  % the input I
 Dnoise = 1.0;  % noise strength
-
+spikes = lifspikes(trials, input, tmax, Dnoise);
-spikes = lifspikes( trials, input, tmax, Dnoise );
     \end{lstlisting}
-    The returned \texttt{spikes} is a cell array with \texttt{trials} elements, each being a vector
+    The returned \texttt{spikes} is a cell array with \texttt{trials}
-    of spike times (in seconds) computed for a duration of \texttt{tmax} seconds.
+    elements, each being a vector of spike times (in seconds) computed
-    The input is set via the \texttt{input} variable, the noise strength via \texttt{Dnoise}.
+    for a duration of \texttt{tmax} seconds.  The input is set via the
+    \texttt{input} variable, the noise strength via \texttt{Dnoise}.
 
-    Think of calling the \texttt{lifspikes()} function as a
+    Think of calling the \texttt{lifspikes()} function as a simple way
-    simple way of doing an electrophysiological experiment. You are
+    of doing an electrophysiological experiment. You are presenting a
-    presenting a stimulus with a constant intensity $I$ that you set. The
+    stimulus with a constant intensity $I$ that you set. The neuron
-    neuron responds to this stimulus, and you record this
+    responds to this stimulus, and you record this response. After
-    response. After detecting the timepoints of the spikes in your
+    detecting the timepoints of the spikes in your recordings you get
-    recordings you get what the \texttt{lifspikes()} function
+    what the \texttt{lifspikes()} function returns. In addition you
-    returns. The advantage over real data is, that you have the
+    can record from different neurons with different noise properties
-    possibility to simply modify the properties of the neuron via the
+    by setting the \texttt{Dnoise} parameter to different values.
-    \texttt{Dnoise} parameter.
 
   \begin{parts}
     \part First set the noise \texttt{Dnoise=0} (no noise). Compute
-    and plot the mean firing rate (number of spikes within the
+    and plot neuron's $f$-$I$ curve, i.e. the mean firing rate (number
-    recording time \texttt{tmax} divided by \texttt{tmax} and averaged
+    of spikes within the recording time \texttt{tmax} divided by
-    over trials) as a function of the input for inputs ranging from 0
+    \texttt{tmax} and averaged over trials) as a function of the input
-    to 20.
+    for inputs ranging from 0 to 20.
 
-    \part Do the same for various noise strength \texttt{Dnoise}. Use $D_{noise} = 1e-3$,
+    How are different stimulus intensities encoded by the firing rate
-    1e-2, and 1e-1. How does the intrinsic noise influence the response curve?
+    of this neuron?
 
-    \part Show some interspike interval histograms for some
+    \part Compute the $f$-$I$ curves of neurons with various noise
-    interesting values of the input and the noise strength.
+    strengths \texttt{Dnoise}. Use $D_{noise} = 1e-3$, $1e-2$, and
+    $1e-1$. 
+
+    How does the intrinsic noise influence the response curve?
+
+    How is the encoding of stimuli influenced by increasing intrinsic
+    noise?
+
+    What are possible sources of this intrinsic noise?
+
+    \part Show spike raster plots and interspike interval histograms
+    of the responses for some interesting values of the input and the
+    noise strength. For example, you might want to compare the
+    responses of the four different neurons to the same input, or by
+    the same resulting mean firing rate.
 
     \part How does the coefficient of variation $CV_{isi}$ (standard
     deviation divided by mean) of the interspike intervalls depend on
     the input and the noise level?
+
+    \part Based o your results, discuss how intrinsic noise might
+    improve and how it might deteriote the encoding of different
+    stimulus intensities.
     
 
  \end{parts}

projects/project_noiseficurves/solution/ficurve.m

@@ -0,0 +1,8 @@
+function rates = ficurve(trials, inputs, tmax, D)
+% compute f-I curve.
+rates = zeros(length(inputs), 1);
+for k=1:length(inputs)
+	spikes = lifspikes(trials, inputs(k), tmax, D);
+	rates(k) = firingrate(spikes, 0.0, tmax);
+end
+end

projects/project_noiseficurves/solution/firingrate.m

@@ -0,0 +1,9 @@
+function rate = firingrate(spikes, tmin, tmax)
+% mean firing rate between tmin and tmax.
+rates = zeros(length(spikes), 1);
+for i = 1:length(spikes)
+	times= spikes{i};
+	rates(i) = length(times((times>=tmin)&(times<=tmax)))/(tmax-tmin);
+end
+rate = mean(rates);
+end

projects/project_noiseficurves/solution/isih.m

@@ -0,0 +1,11 @@
+function isih(spikes, bins)
+isis = [];
+for i = 1:length(spikes)
+	times= spikes{i};
+	isis = [isis; diff(times(:))];
+end
+[h, b] = hist(isis, bins);
+h = h / sum(h) / (bins(2)-bins(1));
+bar(1000.0*b, h);
+xlim([0 1000.0*b(end)])
+end

projects/project_noiseficurves/solution/lifspikes.m

@@ -0,0 +1,33 @@
+function spikes = lifspikes(trials, input, tmax, D)
+% Generate spike times of a leaky integrate-and-fire neuron
+% trials: the number of trials to be generated
+% input: the stimulus either as a single value or as a vector
+% tmax: duration of a trial
+% D: the strength of additive white noise
+    
+    tau = 0.01;
+    if nargin < 4
+        D = 1e0;
+    end
+    vreset = 0.0;
+    vthresh = 10.0;
+    dt = 1e-4;
+    
+    n = ceil(tmax/dt);
+    spikes = cell(trials, 1);
+    for k=1:trials
+        times = [];
+        j = 1;
+        v = vreset;
+        noise = sqrt(2.0*D)*randn(n, 1)/sqrt(dt);
+        for i=1:n
+            v = v + (- v + noise(i) + input)*dt/tau;
+            if v >= vthresh
+                v = vreset;
+                times(j) = i*dt;
+                j = j + 1;
+            end
+        end
+        spikes{k} = times;
+    end
+end

projects/project_noiseficurves/solution/noiseficurves.m

@@ -0,0 +1,26 @@
+%% general settings for the model neuron:
+trials = 10;
+tmax = 50.0;
+
+%% f-I curves:
+figure()
+Ds = [0, 0.001, 0.01, 0.1];
+for j = 1:length(Ds)
+	D = Ds(j);
+	inputs = 0.0:0.5:20.0;
+	rates = ficurve(trials, inputs, tmax, D);
+	plot(inputs, rates);
+	hold on;
+end
+hold off;
+
+%% spike raster and CVs
+input = 12.0;
+for j = 1:length(Ds)
+	D = Ds(j);
+	spikes = lifspikes(trials, input, tmax, D);
+	subplot(4, 2, 2*j-1);
+	spikeraster(spikes, 0.0, 1.0);
+	subplot(4, 2, 2*j);
+	isih(spikes, [0:0.001:0.04]);
+end

projects/project_noiseficurves/solution/spikeraster.m

@@ -0,0 +1,30 @@
+function spikeraster(spikes, tmin, tmax)
+% Display a spike raster of the spike times given in spikes.
+%
+% spikeraster(spikes, tmax)
+%   spikes: a cell array of vectors of spike times in seconds
+%   tmin: plot spike raster starting at tmin seconds
+%   tmax: plot spike raster upto tmax seconds
+
+ntrials = length(spikes);
+for k = 1:ntrials
+    times = spikes{k};
+    times = times((times>=tmin) & (times<=tmax));
+    if tmax < 1.5
+        times = 1000.0*times; % conversion to ms
+    end
+    for i = 1:length( times )
+      line([times(i) times(i)],[k-0.4 k+0.4], 'Color', 'k'); 
+    end
+end
+if (tmax-tmin) < 1.5
+    xlabel('Time [ms]');
+    xlim([1000.0*tmin 1000.0*tmax]);
+else
+    xlabel('Time [s]');
+    xlim([tmin tmax]);
+end
+ylabel('Trials');
+ylim([0.3 ntrials+0.7 ]);
+end
+