Skip to content

Commit db31eea

Browse files
tony-azevedotony-azevedo
committed
HH and plotting and SongSoundBox
1 parent 6bc8cac commit db31eea

File tree

7 files changed

+474
-1
lines changed

7 files changed

+474
-1
lines changed

‎plotting/savePDFandFIG.m

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,13 @@
1+
function s = savePDFandFIG(fig,savedir,auxdir,name)
2+
3+
if ~isdir(fullfile(savedir,auxdir))
4+
mkdir(fullfile(savedir,auxdir))
5+
end
6+
7+
fn = fullfile(savedir,auxdir,[name '.pdf']);
8+
figure(fig)
9+
eval(['export_fig ' regexprep(fn,'\sAzevedo',''' Azevedo''') ' -pdf -transparent']);
10+
11+
set(fig,'paperpositionMode','auto');
12+
saveas(fig,fullfile(savedir,auxdir,name),'fig');
13+
s = 1;

‎sandbox/HH0.m

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@
3232
end
3333

3434

35-
gNa = 12; ENa=115; % gNa = 120; ENa=115;
35+
gNa = 120; ENa=115; % gNa = 120; ENa=115;
3636
gK = gNa*36/120; EK=-12; % gK = 36; EK=-12;
3737
gL= gNa*.3/120; ERest=100; % gL=0.3; ERest=10.6;
3838

‎sandbox/HHGSandbox.m

Lines changed: 138 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,138 @@
1+
% Playing with HH model
2+
3+
%% Standard time vector
4+
5+
t = -100:.01:200;
6+
V = -65*ones(size(t));
7+
V(t>=0&t<100) = 30;
8+
9+
plot(t,V)
10+
11+
%% Play with conductance ratios, where does the neuron come to rest?
12+
13+
% Cool! this is very sensitive to these parameters
14+
close all
15+
16+
gNa0 = 120; % 120;
17+
gKratio = .10; % 0.3;
18+
gLratio0 = .008; % 0.0025;
19+
20+
m0 = 0.0529;
21+
h0 = 0.5961;
22+
n0 = 0.3177;
23+
24+
[I_na,I_k,I_l,m,h,n,t,V] = HH_IClamp(0,t,'V0',-70,'m0',m0,'h0',h0,'n0',n0,'gNa',gNa0,'gKratio',gKratio,'gLratio',gLratio0);
25+
HHPlot(t,V,I_na,I_k,I_l,[],m,h,n);
26+
27+
28+
29+
%% allow neuron to settle
30+
V = -65*ones(size(t));
31+
32+
[I_na,I_k,I_l,m,h,n,t,V] = HH_VClamp(V,t,'m0',0.05,'h0',0.54,'n0',0.34,'gNa',120,'gKratio',0.1,'gLratio',0.01);
33+
HHPlot(t,V,I_na,I_k,I_l,[],m,h,n);
34+
35+
m0 = m(end);
36+
h0 = h(end);
37+
n0 = n(end);
38+
39+
%% settled
40+
V = -65*ones(size(t));
41+
42+
[I_na,I_k,I_l,m,h,n,t] = HH_VClamp(V,t,'m0',m0,'h0',h0,'n0',n0,'gNa',120,'gKratio',0.3,'gLratio',0.0025);
43+
HHPlot(t,V,I_na,I_k,I_l,[],m,h,n);
44+
45+
m0 = m(end);
46+
h0 = h(end);
47+
n0 = n(end);
48+
49+
%% Vclamp
50+
V = -65*ones(size(t));
51+
V(t>=0&t<100) = -30;
52+
[I_na,I_k,I_l,m,h,n,t] = HH_VClamp(V,t,'m0',m0,'h0',h0,'n0',n0,'gNa',120,'gKratio',0.3,'gLratio',0.0025);
53+
HHPlot(t,V,-I_na,-I_k,-I_l,[],m,h,n);
54+
55+
56+
%% allow neuron to settle
57+
V = -65*ones(size(t));
58+
59+
[I_na,I_k,I_l,m,h,n,t] = HH_VClamp(V,t,'m0',0.05,'h0',0.54,'n0',0.34,'gNa',120,'gKratio',0.3,'gLratio',0.0025);
60+
HHPlot(t,V,I_na,I_k,I_l,[],m,h,n);
61+
62+
63+
%% Start with settled
64+
[V,~,~,~,~] = HH0(0,t,'V0',V0,'m0',m0(end),'h0',h0(end),'n0',n0(end));
65+
figure(1);
66+
plot(t,V);
67+
68+
%% Inject steady
69+
I = zeros(size(t));
70+
I(t>100) = 1;
71+
[V,~,~,~,t] = HH0(I,t,'V0',V0,'m0',m0(end),'h0',h0(end),'n0',n0(end));
72+
figure(1);
73+
plot(t,V);
74+
75+
%% Inject steady
76+
I = zeros(size(t));
77+
I(t>100) = 30;
78+
[V,~,~,~,t] = HH0(I,t,'V0',V0,'m0',m0(end),'h0',h0(end),'n0',n0(end));
79+
figure(1);
80+
plot(t,V);
81+
xlabel('ms')
82+
ylabel('mV')
83+
84+
85+
%% Inject steady
86+
I = zeros(size(t));
87+
I(t>100) = 100;
88+
[V,~,~,~,t] = HH0(I,t,'V0',V0,'m0',m0(end),'h0',h0(end),'n0',n0(end));
89+
figure(1);
90+
plot(t,V);
91+
92+
%% Inject steady
93+
I = zeros(size(t));
94+
I(t>100) = 20;
95+
[V,m,h,n,t] = HH0(I,t,'V0',V0,'m0',m0(end),'h0',h0(end),'n0',n0(end));
96+
figure(1);
97+
plot(t,V);
98+
xlabel('ms')
99+
ylabel('mV')
100+
101+
VD = V(end);
102+
nD = n(end);
103+
mD = m(end);
104+
hD = h(end);
105+
106+
%% Inject oscillating
107+
I_off = 20;
108+
ramp = ones(size(t));
109+
dt = t(2)-t(1);
110+
ramp(t<100) = t(t<100)/100;
111+
ramp(t>t(end)-100) = flipud(t(t<=100)/100);
112+
113+
f = [25, 50, 100, 200, 400]; %Hz
114+
a = [1 2 4 8];
115+
116+
figure(2)
117+
ylims = [Inf,-Inf];
118+
for ii = 1:length(a)
119+
for jj = 1:length(f)
120+
f_ = f(jj)/1000; % cyc/ms
121+
I = a(ii)*sin(2*pi*f_ * t).*ramp + I_off;
122+
%I(:) = 170;
123+
124+
[V,n,m,h,t] = HH0(I,t,'V0',VD,'m0',mD,'h0',hD,'n0',nD);
125+
126+
ax = subplot(length(a),length(f),(ii-1)*length(f)+jj);
127+
plot(t,V);
128+
yl = get(ax,'ylim');
129+
ylims = [min(ylims(1),yl(1)),max(ylims(2),yl(2))];
130+
if ii==1
131+
title(ax,[num2str(f(jj)) ' Hz']);
132+
end
133+
end
134+
end
135+
set(get(2,'children'),'ylim',ylims)
136+
xlabel('ms')
137+
138+

‎sandbox/HHPlot.m

Lines changed: 58 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
function fig = HHPlot(T,V,I_na,I_k,I_l,I_inj,m,h,n)
2+
3+
fig = figure;
4+
set(fig,'color',[1 1 1],'position',[680 195 1083 783],'name','HH Plot');
5+
6+
pnl = panel(fig);
7+
pnl.margin = [20 20 20 20];
8+
pnl.pack('v',{1/9 4/9 4/9});
9+
pnl.de.margin = [10 10 10 10];
10+
11+
line(T,V,'parent',pnl(1).select(),...
12+
'color',[0 0 0],'displayname','V');
13+
pnl(1).ylabel('mV');
14+
legend(pnl(1).select(),'toggle');
15+
legend(pnl(1).select(),'boxoff');
16+
set(pnl(1).select(),'xtick',[],'xcolor',[1 1 1])
17+
18+
if ~isempty(I_na)
19+
line(T,I_na,'parent',pnl(2).select(),...
20+
'color',[1 0 0],'displayname','I_na');
21+
end
22+
if ~isempty(I_k)
23+
line(T,I_k,'parent',pnl(2).select(),...
24+
'color',[0 0 1],'displayname','I_k');
25+
end
26+
if ~isempty(I_l)
27+
line(T,I_l,'parent',pnl(2).select(),...
28+
'color',[0 1 0],'displayname','I_l');
29+
end
30+
pnl(2).ylabel('pA');
31+
legend(pnl(2).select(),'toggle');
32+
legend(pnl(2).select(),'boxoff');
33+
34+
if ~isempty(I_inj)
35+
line(T,I_inj,'parent',pnl(2).select(),...
36+
'color',[0 0 0],'displayname','I_inj');
37+
else
38+
line(T,I_na+I_k+I_l,'parent',pnl(2).select(),...
39+
'color',[1 1 1]*.8,'displayname','I_t');
40+
end
41+
42+
line(T,m,'parent',pnl(3).select(),...
43+
'color',[.7 0 0],'displayname','m');
44+
line(T,h,'parent',pnl(3).select(),...
45+
'color',[1 .7 .7],'displayname','h');
46+
line(T,m.^3.*h,'parent',pnl(3).select(),...
47+
'color',[1 0 0],'linewidth',2,'displayname','gNa');
48+
line(T,n,'parent',pnl(3).select(),...
49+
'color',[.7 .7 1],'displayname','n');
50+
line(T,n.^4,'parent',pnl(3).select(),...
51+
'color',[0 0 1],'linewidth',2,'displayname','gK');
52+
53+
pnl(3).ylabel('p_0');
54+
pnl(3).xlabel('ms');
55+
legend(pnl(3).select(),'toggle');
56+
legend(pnl(3).select(),'boxoff');
57+
58+
linkaxes([pnl(1).select() pnl(2).select() pnl(3).select()],'x')

‎sandbox/HH_IClamp.m

Lines changed: 105 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,105 @@
1+
%Implement the Hodgkin-Huxley equations.
2+
%See Gerstner and Kistler, Spiking Neuron Models, 2002, Section 2.2.
3+
%You'll see I've scaled the voltage by 65 in the equation that updates V
4+
%and the auxillary functions. Hodgkin and Huxley set the resting voltage
5+
%of the neuron to 0 mV and we've set it here to -65 mV (the value accepted
6+
%today).
7+
8+
%INPUTS
9+
% I0 = input current.
10+
% T0 = total time to simulate (in [ms]).
11+
%
12+
%OUTPUTS
13+
% V = the voltage of neuron.
14+
% m = activation variable for Na-current.
15+
% h = inactivation variable for Na-current.
16+
% n = activation variable for K-current.
17+
% t = the time axis of the simulation (useful for plotting).
18+
19+
function [I_na,I_k,I_l,m,h,n,T,V] = HH_IClamp(I,T,varargin)
20+
21+
dt = 0.01; % ms (dt = 0.01ms)
22+
if length(T) == 1
23+
T = ceil(T/dt); %T = ceil(T0/dt);
24+
T = (1:T-1)';
25+
T = T*dt;
26+
else
27+
T = T(:);
28+
dt = T(2)-T(1);
29+
end
30+
if length(I) == 1
31+
I = ones(size(T)) * I; %T = ceil(T0/dt);
32+
end
33+
34+
V = zeros(size(T));
35+
I_na = zeros(size(T));
36+
I_k = zeros(size(T));
37+
I_l = zeros(size(T));
38+
m = zeros(size(T));
39+
h = zeros(size(T));
40+
n = zeros(size(T));
41+
42+
p = inputParser;
43+
p.addParamValue('V0',-70,@isnumeric);
44+
p.addParamValue('m0',0.05,@isnumeric);
45+
p.addParamValue('h0',0.54,@isnumeric);
46+
p.addParamValue('n0',0.34,@isnumeric);
47+
p.addParamValue('gNa',120,@isnumeric);
48+
p.addParamValue('gKratio',0.3,@isnumeric);
49+
p.addParamValue('gLratio',0.0025,@isnumeric);
50+
51+
parse(p,varargin{:});
52+
53+
V(1)= p.Results.V0; % V(1)=-70.0;
54+
m(1)= p.Results.m0;
55+
h(1)= p.Results.h0;
56+
n(1)= p.Results.n0;
57+
58+
gNa = p.Results.gNa; ENa=115; % gNa = 120; ENa=115;
59+
gK = gNa*p.Results.gKratio; EK=-12; % gK = 36; EK=-12;
60+
gL = gNa*p.Results.gLratio; ERest=100; % gL=0.3; ERest=10.6;
61+
62+
63+
for i=1:length(T)-1
64+
V(i+1) = V(i) + dt*(gNa*m(i)^3*h(i)*(ENa-(V(i)+65)) + gK*n(i)^4*(EK-(V(i)+65)) + gL*(ERest-(V(i)+65)) + I(i));
65+
I_na(i) = gNa*m(i)^3*h(i)*(ENa-(V(i)+65));
66+
I_k(i) = gK*n(i)^4*(EK-(V(i)+65));
67+
I_l(i) = gL*(ERest-(V(i)+65));
68+
69+
m(i+1) = m(i) + dt*(alphaM(V(i))*(1-m(i)) - betaM(V(i))*m(i));
70+
h(i+1) = h(i) + dt*(alphaH(V(i))*(1-h(i)) - betaH(V(i))*h(i));
71+
n(i+1) = n(i) + dt*(alphaN(V(i))*(1-n(i)) - betaN(V(i))*n(i));
72+
end
73+
74+
I_na(i+1) = gNa*m(i+1)^3*h(i+1)*(ENa-(V(i+1)+65));
75+
I_k(i+1) = gK*n(i+1)^4*(EK-(V(i+1)+65));
76+
I_l(i+1) = gL*(ERest-(V(i+1)+65));
77+
78+
79+
end
80+
81+
%Below, define the AUXILIARY FUNCTIONS alpha & beta for each gating variable.
82+
83+
function aM = alphaM(V)
84+
aM = (2.5-0.1*(V+65)) ./ (exp(2.5-0.1*(V+65)) -1);
85+
end
86+
87+
function bM = betaM(V)
88+
bM = 4*exp(-(V+65)/18);
89+
end
90+
91+
function aH = alphaH(V)
92+
aH = 0.07*exp(-(V+65)/20);
93+
end
94+
95+
function bH = betaH(V)
96+
bH = 1./(exp(3.0-0.1*(V+65))+1);
97+
end
98+
99+
function aN = alphaN(V)
100+
aN = (0.1-0.01*(V+65)) ./ (exp(1-0.1*(V+65)) -1);
101+
end
102+
103+
function bN = betaN(V)
104+
bN = 0.125*exp(-(V+65)/80);
105+
end

0 commit comments

Comments
 (0)