% Simulation of power-law dynamic gate variables in the Hodgkin-Huxley
% model using fractional order derivatives.
% Teka W, Stockton D, Santamaria F. "Power-law dynamics of membrane
% conductances increase spiking diversity in a Hdgkin-Huxley model" PLoS
% Computational Biology, in press, 2016.
% If you use this software please reference our paper.
% Similarly to the Tekaetaal_SignleGate, sections 7-9 generate a few action potentials.
%Sections 10-12 analyze these data which corresponds to Figure 2.
% The other figures were generate by running sections 7-9 for very long
% periods of time. They require super-computer resources and multiple days
%of simulation time. Data files can be available by request.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The full Hodgkin-Huxley model with individual gates having power-law
% dynamics. We model only a few miliseconds to classify the spike shapes.
% Sections 7, 8, and 9 run the models for power-law N, H, and M. Section
% 10, 11, 12 analyze the respective data.
% Data files can be available by request.
% In sections 7-10 we vary the current (0-24 nA) and the values of the fractional
% order derivative (0.2-1.0).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Section 7 - Modeling Spiking N
clear
homedir='.'; %this is your directory
cd(homedir)
tstop = 100; %ms
tstep = 1e-3; %ms
mVinit = -65; %mV
Ncells = 1; %only modeling one cell
Cm=1; %microF/cm^2
dt =tstep;
t=0:dt:tstop; %ms
v0 = mVinit*ones([1 Ncells]); % mV initial value
vrest=v0(1,1);% mV the resting potential.
gK=36; gNa=120; gL=0.3; % channel conductances: mS/cm2 -- (micro A/mV)/cm^2
EK=-12 + v0(1,1); ENa=115 + v0(1,1); EL=10.6 + v0(1,1); % channel reversal potentials: mV
m=0.0529;h=0.596; n=0.3177; % Initial vlaues: steady state for 0 input n=0.3177, m=0.0529, h=0.596, v=-64.9997
NetProp.Ncells=Ncells;
NetProp.dt=dt;
NetProp.Cm=Cm;
NetProp.v0=v0;
NetProp.vrest=vrest;
NetProp.gK=gK;
NetProp.gNa=gNa;
NetProp.gL=gL;
NetProp.EK=EK;
NetProp.ENa=ENa;
NetProp.EL=EL;
NetProp.m=m;
NetProp.h=h;
NetProp.n=n;
NetProp.Noise=0;
Iamp_1=[1:12];
Iamp_2=[13:24];
Iamp=[Iamp_1 Iamp_2];
I=ones(length(t),Ncells);
I(1:10/dt,Ncells)=0;
f2safe='HH_fracN_firingrateVI';
save(f2safe)
matlabpool('open',12) %use if you have parallel matlab
tic
parfor b=1:length(Iamp)
out02(b)=runNetworkderivativeHHFractionalPotassium(NetProp,Iamp(b)*I,t,0.2);
end
toc
fprintf('Done with 02 \n')
save(f2safe,'out02','-append');
clear out02
tic
parfor b=1:length(Iamp)
out04(b)=runNetworkderivativeHHFractionalPotassium(NetProp,Iamp(b)*I,t,0.4);
end
fprintf('Done with 04 \n')
save(f2safe,'out04','-append');
clear out04
toc
tic
parfor b=1:length(Iamp)
out06(b)=runNetworkderivativeHHFractionalPotassium(NetProp,Iamp(b)*I,t,0.6);
end
fprintf('Done with 06 \n')
save(f2safe,'out06','-append');
clear out06
toc
tic
parfor b=1:length(Iamp)
out08(b)=runNetworkderivativeHHFractionalPotassium(NetProp,Iamp(b)*I,t,0.8);
end
fprintf('Done with 08 \n')
save(f2safe,'out08','-append');
clear out08
toc
% we don't run this because we used the one from the H gate. This is the
% classic Hodgkin-Huxley model
% parfor b=1:length(Iamp)
% out10(b)=runNetworkderivativeHHFractionalPotassium(NetProp,Iamp(b)*I,t,1.0);
%
% end
% fprintf('done\n');
% save(f2safe,'out10','-append');
% clear out10
matlabpool('close')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Section 8- Modeling spiking fractional M
clear
homedir='/home/fidel/Wondimu';
cd(homedir)
tstop = 100; %ms
tstep = 1e-3; %ms
mVinit = -65; %mv
Ncells = 1;
Cm=1; %microF/cm^2
dt =tstep;
t=0:dt:tstop; %ms
v0 = mVinit*ones([1 Ncells]); % mV initial value
vrest=v0(1,1);% mV the resting potential.
gK=36; gNa=120; gL=0.3; % channel conductances: mS/cm2 -- (micro A/mV)/cm^2
EK=-12 + v0(1,1); ENa=115 + v0(1,1); EL=10.6 + v0(1,1); % channel reversal potentials: mV
m=0.0529;h=0.596; n=0.3177; % Initial vlaues: steady state for 0 input n=0.3177, m=0.0529, h=0.596, v=-64.9997
NetProp.Ncells=Ncells;
NetProp.dt=dt;
NetProp.Cm=Cm;
NetProp.v0=v0;
NetProp.vrest=vrest;
NetProp.gK=gK;
NetProp.gNa=gNa;
NetProp.gL=gL;
NetProp.EK=EK;
NetProp.ENa=ENa;
NetProp.EL=EL;
NetProp.m=m;
NetProp.h=h;
NetProp.n=n;
NetProp.Noise=0;
Iamp_1=[1:12];
Iamp_2=[13:24];
Iamp=[Iamp_1 Iamp_2];
I=ones(length(t),Ncells);
I(1:10/dt,Ncells)=0;
f2safe='HH_fracM_firingrateVI';
save(f2safe)
matlabpool('open',12)
tic
parfor b=1:length(Iamp)
out02(b)=runNetworkderivativeHHFractionalNa_m(NetProp,Iamp(b)*I,t,0.2);
end
toc
save(f2safe,'out02','-append');
clear out02
fprintf('Done with 02 \n')
parfor b=1:length(Iamp)
out04(b)=runNetworkderivativeHHFractionalNa_m(NetProp,Iamp(b)*I,t,0.4);
end
fprintf('Done with 04 \n')
save(f2safe,'out04','-append');
clear out04
parfor b=1:length(Iamp)
out06(b)=runNetworkderivativeHHFractionalNa_m(NetProp,Iamp(b)*I,t,0.6);
end
fprintf('Done with 06 \n')
save(f2safe,'out06','-append');
clear out06
parfor b=1:length(Iamp)
out08(b)=runNetworkderivativeHHFractionalNa_m(NetProp,Iamp(b)*I,t,0.8);
end
fprintf('Done with 08 \n')
save(f2safe,'out08','-append');
clear out08
% parfor b=1:length(Iamp)
% out10(b)=runNetworkderivativeHHFractionalNa_m(NetProp,Iamp(b)*I,t,1.0);
%
% end
%
% fprintf('done\n');
% save(f2safe,'out10','-append');
% clear out10
matlabpool('close')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Section 9 - Modeling Spiking fractional H
clear
homedir='.';
cd(homedir)
tstop = 100;
tstep = 1e-3;
mVinit = -65;
Ncells = 1;
Cm=1; %microF/cm^2
dt =tstep;
t=0:dt:tstop; %ms
v0 = mVinit*ones([1 Ncells]); % mV initial value
vrest=v0(1,1);% mV the resting potential.
gK=36; gNa=120; gL=0.3; % channel conductances: mS/cm2 -- (micro A/mV)/cm^2
EK=-12 + v0(1,1); ENa=115 + v0(1,1); EL=10.6 + v0(1,1); % channel reversal potentials: mV
m=0.0529;h=0.596; n=0.3177; % Initial vlaues: steady state for 0 input n=0.3177, m=0.0529, h=0.596, v=-64.9997
NetProp.Ncells=Ncells;
NetProp.dt=dt;
NetProp.Cm=Cm;
NetProp.v0=v0;
NetProp.vrest=vrest;
NetProp.gK=gK;
NetProp.gNa=gNa;
NetProp.gL=gL;
NetProp.EK=EK;
NetProp.ENa=ENa;
NetProp.EL=EL;
NetProp.m=m;
NetProp.h=h;
NetProp.n=n;
NetProp.Noise=0;
alphavalue=[1 0.8 0.4 0.2 ];
Iamp_1=[1:12];
Iamp_2=[13:24];
Iamp=[Iamp_1 Iamp_2];
b=1;
I=ones(length(t),Ncells);
I(1:10/dt,Ncells)=0;
f2safe='HH_fracH_firingrateVI';
save(f2safe)
matlabpool('open',12)
tic
parfor b=1:length(Iamp)
out02(b)=runNetworkderivativeHHFractionalNa_h(NetProp,Iamp(b)*I,t,0.2);
end
toc
save(f2safe,'out02','-append');
clear out02
fprintf('Done with 02 \n')
parfor b=1:length(Iamp)
out04(b)=runNetworkderivativeHHFractionalNa_h(NetProp,Iamp(b)*I,t,0.4);
end
fprintf('Done with 04 \n')
save(f2safe,'out04','-append');
clear out04
parfor b=1:length(Iamp)
out06(b)=runNetworkderivativeHHFractionalNa_h(NetProp,Iamp(b)*I,t,0.6);
end
fprintf('Done with 06 \n')
save(f2safe,'out06','-append');
clear out06
parfor b=1:length(Iamp)
out08(b)=runNetworkderivativeHHFractionalNa_h(NetProp,Iamp(b)*I,t,0.8);
end
fprintf('Done with 08 \n')
save(f2safe,'out08','-append');
clear out08
parfor b=1:length(Iamp)
out10(b)=runNetworkderivativeHHFractionalNa_h(NetProp,Iamp(b)*I,t,1.0);
end
fprintf('done\n');
save(f2safe,'out10','-append');
clear out10
matlabpool('close')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Section 10 - Analyze full HH model with power-law N
clear
homedir='.';
cd(homedir)
load HH_fracN_firingrateVI
%the alpha=1 is the same for all the simulations, it is the classic HH
%model
load HH_fracH_firingrateVI out10
fs={'02','04','06','08','10'};
V2u=[1:length(Iamp)];
for cc1=1:length(fs);
eval(['o2a=out' fs{cc1} ';']);
spSh{cc1}=[];cc2=1;
for b=V2u %go through the currents
t=o2a(b).t;
mTh=(o2a(b).v>-20); %detect if there's a spike
sp=find(diff(mTh)==-1);
spN(b,cc1)=numel(sp);
fr(b,cc1)=(numel(sp))/90e-3;
if numel(sp)>1 %record the first evoked spike for each input current
if (sp(2)+7000)<=length(o2a(b).v)
[~,mpr]=max(o2a(b).v(sp(2)+[-2500:7000]));
maxPeak=sp(1)+mpr-2000;
spSh{cc1}(:,cc2)=o2a(b).v(maxPeak+[-2500:7000])';
cc2=cc2+1;
end
end
end
end
clf
% Plot the spike shape of the second generated action potential for all
% values of the fractional exponent at the input current values just above
% rheobase.
cv=['bgrcm'];
etav=[0.2:0.2:1]';
ppv=[1:5];
tsp=[-2500:7000]*dt;
for a=1:5
subplot(4,4,1)
%h=plot(tsp,mean(spSh{a},2),cv(a));
h=plot(tsp,spSh{a}(:,1),cv(a));
formatFig(h,1)
hold on
axis([-2 7 -80 40])
end
%Plot the current at which there was at least 2 action potentials
%generated.
subplot(4,4,5)
spTh=(spN.*~(spN==1));
for a=1:5
dummy=find(spTh(:,a));
Ithv(a)=Iamp(dummy(1));
plot(etav(a),Ithv(a),[cv(a) 'o'])
hold on
end
h=plot([0.2:0.2:1],Ithv,'k--');
formatFig(h,1)
axis([0 1 0 15]);
nst.etavI=etav;
nst.ITh=Ithv;
% now the action potential threshold using phase plane
pcV=[9:16];
for cc1=1:length(fs);
eval(['o2a=out' fs{cc1} ';']);
cc2=1;
for b=V2u
t=o2a(b).t;
mTh=(o2a(b).v>-20);
sp=find(diff(mTh)==-1);
spN(b,cc1)=numel(sp);
fr(b,cc1)=(numel(sp))/90e-3;
if numel(sp)>1
data=o2a(b).v(sp(2)-4000:end);
v=data(1:end-1);
dvdt=diff(data)/1e-3;
w20=((dvdt<21) .* (dvdt>19));
vthw=(v<-30);
th1=(w20.*vthw);
be=find(diff(th1));
be2=diff(reshape(be,2,length(be)/2));
be3=be(1:2:end)+be2';
vthV=v(be3);
dvthV=dvdt(be3);
figure(1)
subplot(4,4,9);%pcV(cc1))
% h=plot(v(1:10:end),dvdt(1:10:end),'k');
h=plot(spSh{cc1}(1:50:(end-50),1),diff(spSh{cc1}(1:50:end,1),1,1)/dt/50,cv(cc1));
formatFig(h,1)
hold on
plot(vthV,dvthV,'rs')
mvthV(cc2)=mean(vthV);
cc2=cc2+1;
end
end
meanspTh(cc1)=mean(mvthV);
end
axis([-80 50 -200 700])
%Plot the voltage threshold as a function of the fractional order
%derivative
subplot(4,4,13)
for a=1:5
h=plot(etav(a),meanspTh(a),[cv(a) 'o']);
hold on
end
h=plot(etav,meanspTh,'k');
formatFig(h,1)
axis([0 1.0 -55 -40])
nst.etav=etav;
nst.meanspTh=meanspTh;
%save VThCompare nst -append %just when you want to compare to the other
%analyses.
print -depsc2 FracN_APshape.eps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Section 11 - Analyse full Hodgkin-Huxley model with power-law M gate
clear
homedir='.';
cd(homedir)
load HH_fracM_firingrateVI
load HH_fracH_firingrateVI out10
%remember that the simulation at fractional order 0.2 were unstable
fs={'04','06','08','10'};
V2u=[1:length(Iamp)];
for cc1=1:length(fs);
eval(['o2a=out' fs{cc1} ';']);
spSh{cc1}=[];cc2=1;
for b=V2u
t=o2a(b).t;
mTh=(o2a(b).v>0);
sp=find(diff(mTh)==-1);
spN(b,cc1)=numel(sp);
fr(b,cc1)=(numel(sp))/90e-3;
if numel(sp)>1
[~,mpr]=max(o2a(b).v(sp(2)+[-2000:7000]));
maxPeak=sp(2)+mpr-2000;
spSh{cc1}(:,cc2)=o2a(b).v(maxPeak+[-2000:7000])';
cc2=cc2+1;
end
end
end
clf
cv=['bgrcm'];
etav=[0.4:0.2:1]';
ppv=[1:5];
tsp=[-2000:7000]*dt;
for a=1:4
subplot(4,4,1)
h=plot(tsp,mean(spSh{a},2),cv(a));
formatFig(h,1)
hold on
axis([-2 7 -80 50])
end
subplot(4,4,5)
spTh=(spN.*~(spN==1));
for a=1:4
dummy=find(spTh(:,a));
Ithv(a)=Iamp(dummy(1));
plot(etav(a),Ithv(a),[cv(a) 'o'])
hold on
end
h=plot([0.4:0.2:1],Ithv,'k-');
formatFig(h,1)
axis([0 1 0 15]);
mst.etavI=etav;
mst.ITh=Ithv;
% now the action potential threshold using phase plane
pcV=[9:16];
for cc1=1:length(fs);
eval(['o2a=out' fs{cc1} ';']);
cc2=1;
for b=V2u
t=o2a(b).t;
mTh=(o2a(b).v>-20);
sp=find(diff(mTh)==-1);
spN(b,cc1)=numel(sp);
fr(b,cc1)=(numel(sp))/90e-3;
if numel(sp)>1
data=o2a(b).v(sp(2)-4000:end);
v=data(1:end-1);
dvdt=diff(data)/1e-3;
w20=((dvdt<21) .* (dvdt>19));
vthw=(v<-30);
th1=(w20.*vthw);
be=find(diff(th1));
be2=diff(reshape(be,2,length(be)/2));
be3=be(1:2:end)+be2';
vthV=v(be3);
dvthV=dvdt(be3);
figure(1)
subplot(4,4,9);%pcV(cc1))
% h=plot(v(1:10:end),dvdt(1:10:end),'k');
h=plot(spSh{cc1}(1:10:(end-10),1),diff(spSh{cc1}(1:10:end,1),1,1)/dt/10,cv(cc1));
formatFig(h,1)
hold on
plot(vthV,dvthV,'rs')
mvthV(cc2)=mean(vthV);
cc2=cc2+1;
end
end
meanspTh(cc1)=mean(mvthV);
end
axis([-80 50 -200 700])
subplot(4,4,13)
for a=1:4
h=plot(etav(a),meanspTh(a),[cv(a+1) 'o']);
hold on
end
h=plot(etav,meanspTh,'k');
formatFig(h,1)
axis([0 1.0 -55 -40])
mst.etav=etav;
mst.meanspTh=meanspTh;
%save VThCompare mst -append
print -depsc2 FracM_APshape.eps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Section 12 - Analize full Hodgkin-Huxley model with power-law H gate
clear
homedir='.';
cd(homedir)
load HH_fracH_firingrateVI
fs={'02','04','06','08','10'};
V2u=[1:length(Iamp)];
for cc1=1:length(fs);
eval(['o2a=out' fs{cc1} ';']);
spSh{cc1}=[];cc2=1;
for b=V2u
t=o2a(b).t;
mTh=(o2a(b).v>0);
sp=find(diff(mTh)==-1);
spN(b,cc1)=numel(sp);
fr(b,cc1)=(numel(sp))/90e-3;
if numel(sp)>0
[~,mpr]=max(o2a(b).v(sp(1)+[-2000:7000]));
maxPeak=sp(1)+mpr-2000;
spSh{cc1}(:,cc2)=o2a(b).v(maxPeak+[-2000:7000])';
cc2=cc2+1;
end
end
end
clf
cv=['bgrcm'];
etav=[0.2:0.2:1]';
ppv=[1:5];
tsp=[-2000:7000]*dt;
for a=1:5
subplot(4,4,1)
h=plot(tsp,mean(spSh{a},2),cv(a));
formatFig(h,1)
hold on
axis([-2 7 -80 40])
end
subplot(4,4,5)
spTh=(spN.*~(spN==1));
for a=2:5
dummy=find(spTh(:,a));
Ithv(a)=Iamp(dummy(1));
plot(etav(a),Ithv(a),[cv(a) 'o'])
hold on
end
h=plot(etav,Ithv,'k--');
formatFig(h,1)
axis([0 1 0 15]);
hst.etavI=etav;
hst.ITh=Ithv;
% now the action potential threshold using phase plane
pcV=[9:16];
for cc1=2:length(fs);
eval(['o2a=out' fs{cc1} ';']);
cc2=1;
for b=V2u
t=o2a(b).t;
mTh=(o2a(b).v>0);
sp=find(diff(mTh)==-1);
spN(b,cc1)=numel(sp);
fr(b,cc1)=(numel(sp))/90e-3;
if numel(sp)>1
data=o2a(b).v(sp(2)-4000:end);
v=data(1:end-1);
dvdt=diff(data)/1e-3;
w20=((dvdt<21) .* (dvdt>19));
vthw=(v<-30);
th1=(w20.*vthw);
be=find(diff(th1));
be=be(1:2*floor(length(be)/2));
be2=diff(reshape(be,2,length(be)/2));
be3=be(1:2:end)+be2';
vthV=v(be3);
dvthV=dvdt(be3);
figure(1)
subplot(4,4,9);%pcV(cc1))
% h=plot(v(1:10:end),dvdt(1:10:end),'k');
h=plot(spSh{cc1}(1:10:(end-10),1),diff(spSh{cc1}(1:10:end,1),1,1)/dt/10,cv(cc1));
formatFig(h,1)
hold on
plot(vthV,dvthV,'rs')
mvthV(cc2)=mean(vthV);
cc2=cc2+1;
end
end
meanspTh(cc1)=mean(mvthV);
end
axis([-80 50 -200 700])
subplot(4,4,13)
for a=2:5
h=plot(etav(a),meanspTh(a),[cv(a) 'o']);
hold on
end
h=plot(etav,meanspTh,'k');
formatFig(h,1)
axis([0 1.0 -55 -40])
hst.etav=etav;
hst.meanspTh=meanspTh;
save VThCompare hst -append
print -depsc2 FracH_APshape.eps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Comparing thresholds
%This makes a graph that compares the current and voltage thresholds from
%sections 8-10
load VThCompare
clf
subplot(4,4,13)
for a=1:5
h=plot(hst.etav(a),hst.ITh(a),[cv(a) 'o']);
hold on
h=plot(nst.etav(a),nst.ITh(a),[cv(a) 'o']);
end
for a=1:4
h=plot(mst.etav(a),mst.ITh(a),[cv(a+1) 'o']);
hold on
end
h=plot(hst.etav,hst.ITh,'k');
h=plot(mst.etav,mst.ITh,'k');
h=plot(nst.etav,nst.ITh,'k');
formatFig(h,1)
subplot(4,4,14)
for a=1:5
h=plot(nst.etav(a),nst.meanspTh(a),[cv(a) 'o']);
hold on
end
for a=2:5
h=plot(hst.etav(a),hst.meanspTh(a),[cv(a) 'o']);
end
for a=1:4
h=plot(mst.etav(a),mst.meanspTh(a),[cv(a+1) 'o']);
hold on
end
h=plot(hst.etav(2:end),hst.meanspTh(2:end),'k');
h=plot(mst.etav,mst.meanspTh,'k');
h=plot(nst.etav,nst.meanspTh,'k');
formatFig(h,1)
h2=get(h,'parent');
set(h2(1),'xtick',[0.2 0.4 0.8 1])
print -depsc2 ThrehsoldCompare.eps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%