diff --git a/projects/project_input_resistance/input_resistance.tex b/projects/project_input_resistance/input_resistance.tex new file mode 100644 index 0000000..189aefa --- /dev/null +++ b/projects/project_input_resistance/input_resistance.tex @@ -0,0 +1,77 @@ +\documentclass[addpoints,11pt]{exam} +\usepackage{url} +\usepackage{color} +\usepackage{hyperref} + +\pagestyle{headandfoot} +\runningheadrule +\firstpageheadrule +\firstpageheader{Scientific Computing}{Project Assignment}{WS 2016/17} +%\runningheader{Homework 01}{Page \thepage\ of \numpages}{23. October 2014} +\firstpagefooter{}{}{{\bf Supervisor:} Jan Grewe} +\runningfooter{}{}{} +\pointsinmargin +\bracketedpoints + +%\printanswers +%\shadedsolutions + + +\begin{document} +%%%%%%%%%%%%%%%%%%%%% Submission instructions %%%%%%%%%%%%%%%%%%%%%%%%% +\sffamily +% \begin{flushright} +% \gradetable[h][questions] +% \end{flushright} + +\begin{center} + \input{../disclaimer.tex} +\end{center} + +%%%%%%%%%%%%%% Questions %%%%%%%%%%%%%%%%%%%%%%%%% +\section*{Estimating cellular properties of different cell types.} +You will analyze data from intracellular \textit{in vitro} recordings +of pyramidal neurons from two different maps of the electrosensory +lateral line lobe (ELL) of the weakly electric fish +\textit{Apteronotus leptorhynchus}. The membrane resistance and the +membrane capacitance are fundamental properties of a neuron that have +a great influence on the coding properties of the cell. They are +typically estimated by injecting pulses of hyperpolarizing current +into the cell. From the respective responses we can calculate the +membrane resistance by applying Ohm's law ($U = R \cdot I$). To +estimate the membrane capacitance we need to fit an exponential +function of the form $y = a \cdot e^{(-x/\tau)}$ to the response to get the +membrane time-constant $\tau$. With the knowledge of $R$ and $\tau$ we +can estimate the capacitance $C$ from the simple relation $\tau = R +\cdot C$. + +\begin{questions} + \question{} The accompanying dataset (input\_resistance.zip) + contains datasets from cells originating from two different parts of + the ELL, the medial segment (MS) and the centro-medial segment + (CMS). Each mat-file contains four variables. (i) \textit{V} the + average membrane potential of 20 repeated current injections, (ii) + \textit{V\_std} the across-trial standard deviation of the + responses, (iii) \textit{t} a vector representing the recording + time (in ms), and (iv) \textit{I} a vector containing the time-course of the + injected current. + + \begin{parts} + \part{} Create plots of the raw data. Plot the average response as + a function of time. This plot should also show the across-trial + variability. Also plot the time-course of the injected + current. \\[0.5ex] + \part{} Estimate the input resistances of each cell.\\[0.5ex] + \part{} Fit an exponential to the initial few milliseconds of the + current-on response. Use a gradient-descent approach to do + this.\\ It is very important to understand the exponential decay + function. If you are unsure, play with the function and understand + how the parameters influence the decay. (Hint: It might be + necessary to transform the data a bit.)\\[0.5ex] + \part{} Estimate the membrane capacitance of each cell. Compare + $R$, $I$ and $\tau$ between cells of the two segments.\\[0.5ex] + \part{} Optional: use a double exponential and see, if the fit improves. + \end{parts} +\end{questions} + +\end{document} diff --git a/projects/project_input_resistance/resistance.zip b/projects/project_input_resistance/resistance.zip new file mode 100644 index 0000000..7d37fc6 Binary files /dev/null and b/projects/project_input_resistance/resistance.zip differ diff --git a/projects/project_input_resistance/solution/exponentialDecay.m b/projects/project_input_resistance/solution/exponentialDecay.m new file mode 100644 index 0000000..2ee27cf --- /dev/null +++ b/projects/project_input_resistance/solution/exponentialDecay.m @@ -0,0 +1,5 @@ +function y = exponentialDecay(p, x) + + y = p(1) .* exp(-x./p(2)); + + diff --git a/projects/project_input_resistance/solution/gradientDescent.m b/projects/project_input_resistance/solution/gradientDescent.m new file mode 100644 index 0000000..867e863 --- /dev/null +++ b/projects/project_input_resistance/solution/gradientDescent.m @@ -0,0 +1,14 @@ +function [position, errors] = gradientDescent(x, y, position) + +gradient = []; +errors = []; +count = 1; +eps = 0.02; +close all +while isempty(gradient) || norm(gradient) > 0.025 + gradient = lsqGradient(position, x,y); + % disp(gradient) + errors(count) = lsqError(position, x, y); + position = position - eps .* gradient'; + count = count + 1; +end \ No newline at end of file diff --git a/projects/project_input_resistance/solution/lsqError.m b/projects/project_input_resistance/solution/lsqError.m new file mode 100644 index 0000000..1a48a52 --- /dev/null +++ b/projects/project_input_resistance/solution/lsqError.m @@ -0,0 +1,12 @@ +function error = lsqError(parameter, x, y) +% Objective function for fitting a exponential equation to data. +% +% Arguments: parameter, vector containing the parameters for the exp. decay +% x, vector of the input values +% y, vector of the corresponding measured output values +% +% Returns: the estimation error in terms of the mean sqaure error + + y_est = exponentialDecay(parameter, x); + error = meanSquareError(y, y_est); +end diff --git a/projects/project_input_resistance/solution/lsqGradient.m b/projects/project_input_resistance/solution/lsqGradient.m new file mode 100644 index 0000000..44b3881 --- /dev/null +++ b/projects/project_input_resistance/solution/lsqGradient.m @@ -0,0 +1,9 @@ +function gradient = lsqGradient(parameter, x, y) + + h = 1e-1; % stepsize for derivatives + gradient = zeros(length(parameter),1); + for i = 1:length(parameter) + new_p = parameter; + new_p(i) = parameter(i) + h; + gradient(i) = (lsqError(new_p, x, y) - lsqError(parameter, x, y))/ h; + end diff --git a/projects/project_input_resistance/solution/meanSquareError.m b/projects/project_input_resistance/solution/meanSquareError.m new file mode 100644 index 0000000..253d31a --- /dev/null +++ b/projects/project_input_resistance/solution/meanSquareError.m @@ -0,0 +1,10 @@ +function error = meanSquareError(y, y_est) +% Mean squared error between observed and predicted values. +% +% Arguments: y, vector of observed values. +% y_est, vector of predicted values. +% +% Returns: the error in the mean-squared-deviation sense. + + error = mean((y - y_est).^2); +end diff --git a/projects/project_input_resistance/solution/project_input_resistance.m b/projects/project_input_resistance/solution/project_input_resistance.m new file mode 100644 index 0000000..5e52441 --- /dev/null +++ b/projects/project_input_resistance/solution/project_input_resistance.m @@ -0,0 +1,43 @@ +clear +maps = {'CMS', 'MS'}; + +taus = zeros(length(dir('data/resistance/resistance*.mat')), 1); +cell_types = cell(size(taus)); +count = 1; +size(taus) +for j = 1:length(maps) + files = dir(strcat('data/resistance/resistance_', maps{j}, '*.mat')); + for i = 1:length(files) + load(strcat('data/resistance/', files(i).name)); + disp(files(i).name) + x = t(t>0.5 & t<15); + x = x - x(1); + y = V(t>0.5 & t<15); + y = y - y(end); + initial_params = [1, 2.5]; + [params, errors] = gradientDescent(x, y, initial_params); + disp(count) + taus(count) = params(2); + cell_types{count} = maps{j}; + count = count + 1; + plot_fit(x, y, errors, params); + end +end + + +function plot_fit(x, y, errors, parameter) + figure() + subplot(2,1,1) + hold on + scatter(x,y, 'displayname', 'data') + f_x = exponentialDecay(parameter, x); + plot(x, f_x, 'displayname', 'fit') + xlabel('Input') + ylabel('Output') + grid on + legend show + subplot(2,1,2) + plot(errors) + xlabel('optimization steps') + ylabel('error') +end