From 1362b39fd640a83d16805d98e67ce000a155edfb Mon Sep 17 00:00:00 2001
From: Jan Benda <jan.benda@uni-tuebingen.de>
Date: Tue, 21 Jan 2020 19:07:36 +0100
Subject: [PATCH] [projects] solution for activation curve

---
 projects/README                               |  1 +
 .../solution/ivcurve.m                        | 16 ++++++
 .../project_activation_curve/solution/main.m  | 50 +++++++++++++++++++
 .../solution/plotcurrents.m                   | 13 +++++
 .../solution/reversalpotential.m              |  4 ++
 5 files changed, 84 insertions(+)
 create mode 100644 projects/project_activation_curve/solution/ivcurve.m
 create mode 100644 projects/project_activation_curve/solution/main.m
 create mode 100644 projects/project_activation_curve/solution/plotcurrents.m
 create mode 100644 projects/project_activation_curve/solution/reversalpotential.m

diff --git a/projects/README b/projects/README
index 665350c..0f688a4 100644
--- a/projects/README
+++ b/projects/README
@@ -23,6 +23,7 @@ Projects
 
 1) project_activation_curve
 medium
+also normalize activation curve to maximum.
 
 2) project_adaptation_fit
 OK, medium
diff --git a/projects/project_activation_curve/solution/ivcurve.m b/projects/project_activation_curve/solution/ivcurve.m
new file mode 100644
index 0000000..0431d4f
--- /dev/null
+++ b/projects/project_activation_curve/solution/ivcurve.m
@@ -0,0 +1,16 @@
+function [vsteps, peakcurrents] = ivcurve(vsteps, time, currents, tmax)
+
+peakcurrents = zeros(1, length(vsteps));
+for k = 1:length(peakcurrents)
+    c = currents((time>0.0)&(time<tmax), k);
+    minc = min(c);
+    maxc = max(c);
+    if abs(minc) > maxc
+        peakcurrents(k) = minc;
+    else
+        peakcurrents(k) = maxc;
+    end
+end
+
+end
+
diff --git a/projects/project_activation_curve/solution/main.m b/projects/project_activation_curve/solution/main.m
new file mode 100644
index 0000000..f2bd03f
--- /dev/null
+++ b/projects/project_activation_curve/solution/main.m
@@ -0,0 +1,50 @@
+%% plot data:
+x = load('../data/WT_01.mat');
+wtdata = x.data;
+plotcurrents(wtdata.t, wtdata.I);
+
+x = load('../data/A1622D_01.mat');
+addata = x.data;
+plotcurrents(addata.t, addata.I);
+
+%% I-V curve:
+[wtsteps, wtpeaks] = ivcurve(wtdata.steps, wtdata.t, wtdata.I, 100.0);
+
+[adsteps, adpeaks] = ivcurve(addata.steps, addata.t, addata.I, 100.0);
+
+figure();
+plot(wtsteps, wtpeaks, '-b');
+hold on;
+plot(adsteps, adpeaks, '-r');
+hold off;
+
+%% reversal potential:
+wtE = reversalpotential(wtsteps, wtpeaks);
+adE = reversalpotential(adsteps, adpeaks);
+
+%% activation curve:
+wtg = wtpeaks./(wtsteps - wtE);
+adg = adpeaks./(adsteps - adE);
+
+wtinfty = wtg(wtsteps<40.0)/mean(wtg((wtsteps>=20.0)&(wtsteps<=40.0)));
+adinfty = adg(adsteps<40.0)/mean(adg((adsteps>=20.0)&(adsteps<=40.0)));
+wtsteps = wtsteps(wtsteps<40.0);
+adsteps = adsteps(adsteps<40.0);
+
+figure();
+plot(wtsteps, wtinfty, '-b');
+hold on;
+plot(adsteps, adinfty, '-r');
+
+%% boltzmann fit:
+bf = @(p, v) 1.0./(1.0+exp(-p(1)*(v - p(2))));
+p = lsqcurvefit(bf, [1.0, -40.0], wtsteps, wtinfty);
+wtfit = bf(p, wtsteps);
+p = lsqcurvefit(bf, [1.0, -40.0], adsteps, adinfty);
+adfit = bf(p, adsteps);
+
+plot(wtsteps, wtfit, '-b');
+plot(wtsteps, adfit, '-r');
+hold off;
+
+
diff --git a/projects/project_activation_curve/solution/plotcurrents.m b/projects/project_activation_curve/solution/plotcurrents.m
new file mode 100644
index 0000000..e2c510c
--- /dev/null
+++ b/projects/project_activation_curve/solution/plotcurrents.m
@@ -0,0 +1,13 @@
+function plotcurrents(time, currents)
+
+figure();
+hold on;
+for k = 1:size(currents, 2)
+    plot(time, currents(:, k))
+end
+hold off;
+xlabel('Time [ms]')
+ylabel('Current')
+
+end
+
diff --git a/projects/project_activation_curve/solution/reversalpotential.m b/projects/project_activation_curve/solution/reversalpotential.m
new file mode 100644
index 0000000..686b667
--- /dev/null
+++ b/projects/project_activation_curve/solution/reversalpotential.m
@@ -0,0 +1,4 @@
+function E = reversalpotential(vsteps, currents)
+    p = polyfit(vsteps((vsteps>=20.0)&(vsteps<50.0)), currents((vsteps>=20.0)&(vsteps<50.0)), 1);
+    E = -p(2)/p(1);
+end
\ No newline at end of file