Tekaetal2016_SingleGates.m

% 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. 

%There are 6 sections in this file that compute the simulations and
%analyze the results for Figures 1 and 2 of our paper.
%Section 1-3 simulate voltage clamp on the individual gates having 
%power-law dynamcis. They correspond to figure 1 in the paper. 
% Section 4-6 analyze the data produced by section 1-3.

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulation section 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Section 1
% Simulationg power-law dynamics in the N gate of the Hodgking-Huxley model
% corresponding to the potassium conductance.

clear
homedir='.'; %Your home direcoty
cd(homedir)

f2save='fractoinalNgateSweep'; %Name of file

tstop   = 110; %ms
tstep   = 1e-3; %ms
mVinit  = -65; %mV
Iinjamp = -45; %nA


Ncells  = 1; %Model only 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

%Setup the modeling structure NetProp
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; %noise added to the input current

alphavalue=[1 0.8 0.4 0.2 ];

%This is the sweep variable. 
Vamp=[-110:10:120];

V=zeros(length(t),length(Vamp));
for a=1:length(Vamp)
    V(:,a)=Vamp(a)*ones(length(t),Ncells);
    V(1:30/dt,a)=mVinit;
end

save(f2save); %Save all the structures used to run the simulations

%if you have parallel matlab then use this, otherwise comment out next line
matlabpool('open',12) 

parfor b=1:length(Vamp)
    out02(b)=fractionalNgate(NetProp,V(:,b),t,0.2); 
end
save(f2save,'out02','-append')
clear out02;
fprintf('N donee 02\n');

parfor b=1:length(Vamp)
    out04(b)=fractionalNgate(NetProp,V(:,b),t,0.4);
end
save(f2save,'out04','-append')
clear out04;
fprintf('N done 04\n');

parfor b=1:length(Vamp)
    out06(b)=fractionalNgate(NetProp,V(:,b),t,0.6);
end
save(f2save,'out06','-append')
clear out06;
fprintf('N done 06\n');

parfor b=1:length(Vamp)
    out08(b)=fractionalNgate(NetProp,V(:,b),t,0.8);
end
save(f2save,'out08','-append')
clear out08;
fprintf('N done 08\n');

parfor b=1:length(Vamp)
    out10(b)=fractionalNgate(NetProp,V(:,b),t,1); % main output
end
save(f2save,'out10','-append')
clear out10;
fprintf('N done 10\n');

matlabpool('close') %only if you have parallel Matlab
save fractoinalNgateSweep

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Section 2 - Power-law dynamics of the H gate (Sodium conductance)

clear
homedir='/home/fidel/Wondimu';
cd(homedir)
tstop  = 130;
tstep  = 1e-3;
mVinit = -65;
Iinjamp = -45;
alphav = 1;

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

Iinjamplitude=Iinjamp; % nA (or microA/cm2) injected current amplitude,
V=Iinjamplitude*ones(length(t),Ncells);
%V(1:60/dt,Ncells)=mVinit;

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;

Vamp=[-110:5:-90 -88:2:-40 -35:5:60];
f2save='fractoinalHgateSweep';

V=zeros(length(t),length(Vamp));
for a=1:length(Vamp)
    V(:,a)=Vamp(a)*ones(length(t),Ncells);
    V(1:30/dt,a)=mVinit;
end


save(f2save);
matlabpool('open',12)

tic
parfor b=1:length(Vamp)
    out02(b)=fractionalHgate(NetProp,V(:,b),t,0.2); % main output
end
save(f2save,'out02','-append')
clear out02;
fprintf('H donee 02\n');
toc

tic
parfor b=1:length(Vamp)
    out04(b)=fractionalHgate(NetProp,V(:,b),t,0.4);
end
save(f2save,'out04','-append')
clear out04;
fprintf('H donee 04\n');
toc

tic
parfor b=1:length(Vamp)
    out06(b)=fractionalHgate(NetProp,V(:,b),t,0.6);
end
save(f2save,'out06','-append')
clear out06;
fprintf('H donee 06\n');
toc

tic
parfor b=1:length(Vamp)
    out08(b)=fractionalHgate(NetProp,V(:,b),t,0.8);
end
save(f2save,'out08','-append')
clear out08;
fprintf('H donee 08\n');
toc

tic
parfor b=1:length(Vamp)
    out10(b)=fractionalHgate(NetProp,V(:,b),t,1); % main output
end
save(f2save,'out10','-append')
clear out10;
fprintf('H donee 10\n');
toc

matlabpool('close')

save fractoinalHgateSweep

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Section 3 - Power-law dynamics of the M gate (Sodium conductance)

clear
homedir='/home/fidel/Wondimu';
cd(homedir)
tstop  = 1130;
tstep  = 0.5e-3;
mVinit = -65;
Iinjamp = -45;
alphav = 1;

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

Iinjamplitude=Iinjamp; % nA (or microA/cm2) injected current amplitude,
V=Iinjamplitude*ones(length(t),Ncells);
V(1:20/dt,Ncells)=mVinit;

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;

Vamp=[-110:5:60];
f2save='fractoinalMgateSweep';

V=zeros(length(t),length(Vamp));
for a=1:length(Vamp)
    V(:,a)=Vamp(a)*ones(length(t),Ncells);
    V(1:30/dt,a)=mVinit;
end


save(f2save);

matlabpool('open',12)

% These simulation were unstable and were not performed for the project.
% See text for details.
% parfor b=1:length(Vamp)
%     out02(b)=fractionalMgate(NetProp,V(:,b),t,0.2); % main output
% end
% save(f2save,'out02','-append')
% clear out02;
% fprintf('M donee 02\n');
% toc

tic
parfor b=1:length(Vamp)
    out04(b)=fractionalMgate(NetProp,V(:,b),t,0.4);
end
save(f2save,'out04','-append')
clear out04;
fprintf('M donee 04\n');
toc

tic
parfor b=1:length(Vamp)
    out06(b)=fractionalMgate(NetProp,V(:,b),t,0.6);
end
save(f2save,'out06','-append')
clear out06;
fprintf('M donee 06\n');
toc

tic
parfor b=1:length(Vamp)
    out08(b)=fractionalMgate(NetProp,V(:,b),t,0.8);
end
save(f2save,'out08','-append')
clear out08;
fprintf('M donee 08\n');
toc

tic
parfor b=1:length(Vamp)
    out10(b)=fractionalMgate(NetProp,V(:,b),t,1); % main output
end
save(f2save,'out10','-append')
clear out10;
fprintf('M donee 10\n');
toc

matlabpool('close')
save fractoinalMgateSweep


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% After running the simulations of voltage clamp for each gate we then 
% analyze the reults.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Section 4 analyze fractional N kinetics
clear
homedir='.';%type your directory
cd(homedir)

load fractoinalNgateSweep

%dual exponential function to fit the gate response.
f=fittype('ninf*(1-exp(-x*tau))+n2*(1-exp(-x*tau2))'); 
f0=fitoptions(f);
g=[];
cc1=1;
%the suffix for the simulations with different values of fractional
%exponent
fs={'02','04','06','08','10'};
V2u=[1:length(Vamp)]; %Vamp was defined in the particular simulations


%All this section does is to automate the fitting process by selecting
%appropriate bounds. If the gate decreased value from the initial condition
%the the min-max bound would be different if the gate value increased for a
%different voltage. 
for cc1=1:length(fs);
    eval(['o2a=out' fs{cc1} ';']);
    cc2=1;
    for b=V2u
        g(cc2).tr=(o2a(b).t(logical(o2a(b).t>=30))-30)';
        d2f=o2a(b).nV(logical(o2a(b).t>=30));
        g(cc2).nfinf=d2f(end);
        g(cc2).d2f=d2f-d2f(1);
        g(cc2).d2fraw=d2f;
        if cc1==5
            f0.Start=[(g(cc2).d2f(end))*[1 1] 1 1 ];
        else
            f0.Start=[(g(cc2).d2f(end))*[0.1 0.5 ] 0.002 1 ];
        end
        
        if (max(g(cc2).d2f))>0 && min(g(cc2).d2f)==0
            f0.Upper=[abs(g(cc2).d2f(end))*[1.5 1.5 ] 40 40 ];
            f0.Lower=[0*[1.5 1.5 ] 0 0];
        elseif (max(g(cc2).d2f))>0 && min(g(cc2).d2f)<0
            f0.Upper=[0*[1.5 1.5 ] 40 40 ];
            f0.Lower=[min(g(cc2).d2f)*[1.5 1.5 ] 0 0];
        elseif (max(g(cc2).d2f))==0
            f0.Upper=[0*[1.5 1.5 ] 40 40 ];
            f0.Lower=[-abs(g(cc2).d2f(end))*[1.5 1.5 ] 0 0];
        else 
            disp('WHAaaaat!')
        end

        
        g(cc2).fL=fit(g(cc2).tr,...
            g(cc2).d2f,f,f0);
        tauVdummy=1./[g(cc2).fL.tau g(cc2).fL.tau2];
        ref2E=find(g(cc2).tr>min(tauVdummy));
        if isempty(ref2E)
            ref2E=length(g(cc2).tr);
        end
         [b Vamp(b)]
        clf
        plot(g(cc2).tr,(g(cc2).d2f),'k',...
            g(cc2).tr(ref2E(1)),g(cc2).d2f(ref2E(1)),'r*',...
            g(cc2).tr,g(cc2).fL(g(cc2).tr),'r');
        drawnow
        
        cc2=cc2+1;
    end
    gAll(cc1).g=g;%the selected fit.
end

%Analyze and plot data
clf

%Plot the value of the N gate vs time for V=30 and values of fractional 
%order 0.2, 0.4, 0.6, 0.8, 1.0
subplot(4,4,1)
cla
vamp2p=15;
Vamp(vamp2p)
for a=1:length(gAll)
    trEx4(:,a)=gAll(a).g(vamp2p).tr;
    d2a4r(:,a)=gAll(a).g(vamp2p).d2fraw;
    d2a4(:,a)=gAll(a).g(vamp2p).d2f;
    fL(a).f=gAll(a).g(vamp2p).fL;
    trEx12(:,a)=gAll(a).g(vamp2p).tr;
    d2a12(:,a)=gAll(a).g(vamp2p).d2fraw;
end
h=plot(trEx4(:,1),(d2a4r),'-');
formatFig(h,1)
axis([0 30 0.3 1]);

%Log-log plot of the same data and fitting - Not used in the paper
%Same data as in subplot(4,4,1)
subplot(4,4,5)
cla
h=loglog(trEx4(:,1),(d2a4));
formatFig(h,1)
cla
cv='bgrmc';
for a=1:size(d2a4,2)
    [logfit(a).f, logfit(a).g]=fit(log10(trEx4(2:500,1)),log10(d2a4(2:500,a)),'poly1');
    
    h=plot(log10(trEx4(:,1)),log10(d2a4(:,a)),cv(a));
    hold on
    fracexpfig(a)=logfit(a).f.p1;
end
formatFig(h,1)
axis([-3.2 2 -1 0])


%do an adaptive fit for all the log traces!
for a=1:length(gAll)
    for b=1:length(gAll(a).g)
        g2f=gAll(a).g(b);
        out=adaptiveLogFit(g2f,0.95);
        ftM(a,b)=out.time;% get the average from this.
        t2int(a,b)=10.^((log10(gAll(5).g(b).d2fraw(end))-out.f.p2)./out.f.p1); %he last point with alpha =1
        
    end
end

%plot the exponent of the log fit to the voltage response. This corresponds
% to the fractional order of the derivative - Not used in the paper.
subplot(4,4,9)
h=plot(fracexpfig,'o--k');
formatFig(h,1)

%Plot the value of N_inf. In the paper we call this the long term response
%fuction
subplot(4,4,2)
cla
for cc1=1:length(gAll);
    tg=gAll(cc1).g;
    for b=1:length(tg)
        mginf(cc1,b)=tg(b).nfinf;
    end
end
h=plot(Vamp(V2u),mginf','-');
formatFig(h,1);
axis([-110 20 0 1])

%now calculate the slope at the half max value - Not used in the paper
for a=1:size(mginf,1)
    
    V05=(find(diff((mginf(a,:)<0.5))));
    dX=mginf(a,V05+[-5:5]);
    dV=Vamp(V05+[-5:5]);
    flXinf(a).f=fit(dV',dX','poly1');
    slopeXinf(a)=flXinf(a).f.p1;
   % hold on
   % plot(Vamp,flXinf(a).f(Vamp),'k')
end
subplot(4,4,6)
h=plot([0.2:0.2:1],slopeXinf,'o--k');
formatFig(h,1)


%Now calculate the value of Tau_N. Here we fit a slow and a fast time
%constant. When the fractional order is 1 then both time constant are the
%same and equal to the classic Hodgkin-Huxley model
for cc1=1:length(gAll);
    tg=gAll(cc1).g;
    for b=1:length(tg)
        tauV=sort(1./[tg(b).fL.tau tg(b).fL.tau2 ]);
        tauL(cc1,b)=tauV(1);
        tauE(cc1,b)=tauV(2);
        mginf(cc1,b)=tg(b).nfinf;
    end
end

subplot(4,4,10)
cla
v2p=Vamp([1:end] );
tau2p=tauL(1:end,[1:end])';
tau2E=tauE(:,[1:end])';
h=plot(v2p,tau2p);
formatFig(h,1);
 hold on
 h=plot(v2p,tau2E,'--');
 formatFig(h,1);
axis([-110 120 0 35]);

print -depsc Frac_N_Analysis.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Section 5 -analysis of power-law dynamics H gate
clear
homedir='.';%type your directory
cd(homedir)

load fractoinalHgateSweep

f=fittype('ninf*(1-exp(-x*tau))+n2*(1-exp(-x*tau2))');
f0=fitoptions(f);
g=[];
cc1=1;
fs={'02','04','06','08','10'};
V2u=[1:length(Vamp)];

for cc1=1:length(fs);
    eval(['o2a=out' fs{cc1} ';']);
    cc2=1;
    for b=V2u
        g(cc2).tr=(o2a(b).t(logical(o2a(b).t>=30))-30)';
        d2f=o2a(b).nV(logical(o2a(b).t>=30));
        g(cc2).nfinf=d2f(end);
        g(cc2).d2f=d2f-d2f(1);
        g(cc2).d2fraw=d2f;
        if cc1==5
            f0.Start=[(g(cc2).d2f(end))*[1 1] 1 1 ];
        else
            f0.Start=[(g(cc2).d2f(end))*[0.1 0.5 ] 0.002 1 ];
        end
        
        if (max(g(cc2).d2f))>0 && min(g(cc2).d2f)==0
            f0.Upper=[abs(g(cc2).d2f(end))*[1.5 1.5 ] 40 40 ];
            f0.Lower=[0*[1.5 1.5 ] 0 0];
        elseif (max(g(cc2).d2f))>0 && min(g(cc2).d2f)<0
            f0.Upper=[0*[1.5 1.5 ] 40 40 ];
            f0.Lower=[min(g(cc2).d2f)*[1.5 1.5 ] 0 0];
        elseif (max(g(cc2).d2f))==0
            f0.Upper=[0*[1.5 1.5 ] 40 40 ];
            f0.Lower=[-abs(g(cc2).d2f(end))*[1.5 1.5 ] 0 0];
        else 
            disp('WHAaaaat')
        end

        
        g(cc2).fL=fit(g(cc2).tr,...
            g(cc2).d2f,f,f0);
        tauVdummy=1./[g(cc2).fL.tau g(cc2).fL.tau2];
        ref2E=find(g(cc2).tr>min(tauVdummy));
        if isempty(ref2E)
            ref2E=length(g(cc2).tr);
        end
        clf
        plot(g(cc2).tr,(g(cc2).d2f),'k',...
            g(cc2).tr(ref2E(1)),g(cc2).d2f(ref2E(1)),'r*',...
            g(cc2).tr,g(cc2).fL(g(cc2).tr),'r');
        drawnow
%        input('r')
        
        cc2=cc2+1;
    end
    gAll(cc1).g=g;
end

clf
%Plot the value of the H gate vs time for V=-70  and values of fractional 
%order 0.2, 0.4, 0.6, 0.8, 1.0
Vamp(15)
subplot(4,4,1)
cla
for a=1:length(gAll)
    trEx4(:,a)=gAll(a).g(15).tr;
    d2a4r(:,a)=gAll(a).g(15).d2fraw;
    d2a4(:,a)=gAll(a).g(15).d2f;
    fL(a).f=gAll(a).g(15).fL;
    trEx12(:,a)=gAll(a).g(15).tr;
    d2a12(:,a)=gAll(a).g(15).d2fraw;
end
h=plot(trEx4(:,1),(d2a4r),'-');
formatFig(h,1)
axis([0 100 0.6 0.8]);

subplot(4,4,5)
cla
%Log-log plot of the same data and fitting - Not used in the paper
%Same data as in subplot(4,4,1)
cv='bgrmc';
for a=1:size(d2a4,2)
    logfit(a).f=fit(log10(trEx4(2:500,1)),log10(d2a4(2:500,a)),'poly1');
    
    h=plot(log10(trEx4(:,1)),log10(d2a4(:,a)),cv(a));
    hold on
    fracexpfig(a)=logfit(a).f.p1;
end
axis([-2 2 -2 -0.5])

%do an adaptive log fit for all the traces! Not used in the paper
for a=1:length(gAll)
    for b=1:length(gAll(a).g)
        g2f=gAll(a).g(b);
        out=adaptiveLogFit(g2f,0.95);
        ftM(a,b)=out.time;% get the average from this.
    end
end

%plot the exponent of the log fit to the voltage response - Not used in the
%paper
subplot(4,4,9)
h=plot(fracexpfig,'o--k');
formatFig(h,1)

%Plot the value of H_inf. In the paper we call this the long term response
%fuction
subplot(4,4,2)
cla

for cc1=1:length(gAll);
    tg=gAll(cc1).g;
    for b=1:length(tg)
        mginf(cc1,b)=tg(b).nfinf;
    end
end

h=plot(Vamp(V2u),mginf','-');
formatFig(h,1);
axis([-110 20 0 1])

%now calculate the slope at the half max value - Not used in the paper
for a=1:size(mginf,1)   
    V05=(find(diff((mginf(a,:)<0.5))));
    dX=mginf(a,V05+[-5:5]);
    dV=Vamp(V05+[-5:5]);
    flXinf(a).f=fit(dV',dX','poly1');
    slopeXinf(a)=flXinf(a).f.p1;
end
subplot(4,4,6)
h=plot([0.2:0.2:1],slopeXinf,'o--k');
formatFig(h,1)

%Now calculate the value of Tau_H. Here we fit a slow and a fast time
%constant. When the fractional order is 1 then both time constant are the
%same and equal to the classic Hodgkin-Huxley model
%h tau
for cc1=1:length(gAll);
    tg=gAll(cc1).g;
    for b=1:length(tg)
        tauV=sort(1./[tg(b).fL.tau tg(b).fL.tau2 ]);
        tauL(cc1,b)=tauV(1);
        tauE(cc1,b)=tauV(2);
        mginf(cc1,b)=tg(b).nfinf;
    end
end

subplot(4,4,10)
cla
v2p=Vamp([1:end] );
tau2p=tauL(1:end,[1:end])';
tau2E=tauE(:,[1:end])';
h=plot(v2p,tau2p);
formatFig(h,1);
 hold on
 h=plot(v2p,tau2E,'--');
 formatFig(h,1);
axis([-110 20 0 50]);

print -depsc Frac_H_Analysis.eps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Section 6 - Analysis of power-law dynamics M gate

clear
load fractoinalMgateSweep

f=fittype('ninf*(1-exp(-x*tau))'); % here we used a single time constant

f0=fitoptions(f);
g=[];
cc1=1;
fs={'02','04','06','08','10'};
V2u=[1:length(Vamp)];

for cc1=1:length(fs);
    eval(['o2a=out' fs{cc1} ';']);
    cc2=1;
    for b=V2u
        if sum(isnan(o2a(b).nV))==0
            g(cc2).tr=(o2a(b).t(logical(o2a(b).t>=20))-20)';
            d2f=o2a(b).nV(logical(o2a(b).t>=20));
            g(cc2).nfinf=d2f(end);
            g(cc2).d2f=d2f-d2f(1);
            g(cc2).d2fraw=d2f;
            if cc1==5
                f0.Start=[(g(cc2).d2f(end))*[1] 1  ];
            else
                f0.Start=[(g(cc2).d2f(end))*[0.1 ] 0.2  ];
            end
            
            if sum(g(cc2).d2f)==0
                f0.Upper = [1  40];
                f0.Lower = [-1  0 ];
            elseif (max(g(cc2).d2f))>0 && min(g(cc2).d2f)==0
                f0.Upper=[abs(g(cc2).d2f(end))*[ 1.5 ]  40 ];
                f0.Lower=[0*[1.5 ]  0];
            elseif (max(g(cc2).d2f))>0 && min(g(cc2).d2f)<0
                f0.Upper=[0*[ 1.5 ]  40 ];
                f0.Lower=[min(g(cc2).d2f)*[ 1.5 ]  0];
            elseif (max(g(cc2).d2f))==0
                f0.Upper=[0*[ 1.5 ]  40 ];
                f0.Lower=[-abs(g(cc2).d2f(end))*[ 1.5 ]  0];
            else
                disp('WHAaaaat')
            end
            
            if sum(f0.Upper(1:2)==f0.Lower(1:2))==2
                f0.Upper(1:2)=0.2;
                f0.Lower(1:2)=-0.2;
               
            end
            
            g(cc2).fL=fit(g(cc2).tr,...
                g(cc2).d2f,f,f0);
            tauVdummy=1./[g(cc2).fL.tau ];%g(cc2).fL.tau2];
            ref2E=find(g(cc2).tr>min(tauVdummy));
            if isempty(ref2E)
                ref2E=length(g(cc2).tr);
            end
            [b Vamp(b)]
            clf
            plot(g(cc2).tr,(g(cc2).d2f),'k',...
                g(cc2).tr(ref2E(1)),g(cc2).d2f(ref2E(1)),'r*',...
                g(cc2).tr,g(cc2).fL(g(cc2).tr),'r');
            %        input('r') %just in case you can to see the results
        else 
            g(cc2).tr=[];
            g(cc2).nfinf=[];
            g(cc2).d2f=[];
            g(cc2).d2fraw=[];
            g(cc2).fL=[];
        end 
        cc2=cc2+1;
    end
    gAll(cc1).g=g;
end

clf

%Plot the value of the M gate vs time for V=-55  and values of fractional 
%order  0.4, 0.6, 0.8, 1.0

subplot(4,4,1)
cla
for a=1:length(gAll)
   if ~isempty(gAll(a).g(12).tr)
    trEx4(:,a)=gAll(a).g(12).tr;
    d2a4r(:,a)=gAll(a).g(12).d2fraw;
    d2a4(:,a)=gAll(a).g(12).d2f;
    fL(a).f=gAll(a).g(12).fL;
    trEx12(:,a)=gAll(a).g(12).tr;
    d2a12(:,a)=gAll(a).g(12).d2fraw;
   end
end
h=plot(trEx4(:,1),(d2a4r),'-');
formatFig(h,1)
axis([0 20 0 0.2]);

subplot(4,4,5)
cla
h=loglog(trEx4(:,1),(d2a4));
formatFig(h,1)
cla

%Log-log plot of the same data and fitting - Not used in the paper
%Same data as in subplot(4,4,1)
cv='bgrmc';
for a=1:size(d2a4,2)
    logfit(a).f=fit(log10(trEx4(2:500,1)),log10(d2a4(2:500,a)),'poly1');
    
    h=plot(log10(trEx4(:,1)),log10(d2a4(:,a)),cv(a));
    hold on
   % h=plot(log10(trEx4(2:end,1)),(logfit(a).f(log10(trEx4(2:end,1)))),cv(a));
    fracexpfig(a)=logfit(a).f.p1;
end
formatFig(h,1)
axis([-2 1 -1.5 -0.8])

%plot the exponent of the log fit to the voltage response - Not used in the
%paper.
for a=1:length(gAll)
    for b=1:length(gAll(a).g)
        g2f=gAll(a).g(b);
        out=adaptiveLogFit(g2f,0.95);
        ftM(a,b)=out.time;% get the average from this.
    end
end

%plot the exponent of the log fit to the voltage response
subplot(4,4,9)
h=plot(fracexpfig,'o--k');
formatFig(h,1)

%Plot the value of M_inf. In the paper we call this the long term response
%fuction
subplot(4,4,2)
cla
for cc1=1:length(gAll);
    tg=gAll(cc1).g;
    for b=1:length(tg)
        if ~isempty(tg(b).nfinf)
          mginf(cc1,b)=tg(b).nfinf;
        else
            mginf(cc1,b)=0;
        end
    end
end
%some fits got confused because the values were aroung zero, so we deleted them
h=plot(Vamp([1:14 16:end]),mginf(:,[1:14 16:end])','-');
formatFig(h,1);
axis([-110 60 0 1])

%now calculate the slope at the half max value - Not used in the paper

for a=1:size(mginf,1)
    
    V05=(find(diff((mginf(a,:)<0.5))));
    dX=mginf(a,V05(1)+[-5:5]);
    dV=Vamp(V05(1)+[-5:5]);
    flXinf(a).f=fit(dV',dX','poly1');
    slopeXinf(a)=flXinf(a).f.p1;
   % hold on
   % plot(Vamp,flXinf(a).f(Vamp),'k')
end
subplot(4,4,6)
h=plot([0.2:0.2:1],slopeXinf,'o--k');
formatFig(h,1)
axis([0 1 0 0.04])

%Now calculate the value of Tau_M. Here we fit a slow and a fast time
%constant. When the fractional order is 1 then both time constant are the
%same and equal to the classic Hodgkin-Huxley model.
for cc1=1:length(gAll);
    tg=gAll(cc1).g;
    for b=1:length(tg)
        if ~isempty(tg(b).nfinf)
        tauV=sort(1./[tg(b).fL.tau]);% tg(b).fL.tau2 ]);
        tauL(cc1,b)=tauV(1);
        tauE(cc1,b)=tauV(1);
        else
            tauL(cc1,b)=0;
            tauE(cc1,b)=0;
        end
        
    end
end

subplot(4,4,10)
cla
v2p=Vamp([1:9 12:14 16:end] );
tau2p=tauL(1:end,[1:9 12:14 16:end])';
tau2E=tauE(:,[1:9 12:14 16:end])';
h=plot(v2p,tau2p);
formatFig(h,1);
hold on
h=plot(v2p,tau2E,'-');
formatFig(h,1);
axis([-110 120 0 0.6]);


print -depsc Frac_M_Analysis.eps