% This model is the voltage-dependence of gating variables 
% in the Connors-Stevens model, similar to Hodgkin-Huxley, but
% more like neurons in the cortex, being type-I. 
% See Dayan and Abbott pp 166-172 then pp.196-198 and p.224.
% A slowly inactivating calcium T-type conductance is also included
clear
dV = 0.0001;
Vmin = -0.100
Vmax = 0

V=Vmin:dV:Vmax;    % voltage vector

n_inf=zeros(size(V));   % n: potassium activation gating variable
tau_n=zeros(size(V));
m_inf=zeros(size(V));   % m: sodium activation gating variable
tau_m=zeros(size(V));
h_inf=zeros(size(V));   % h: sodim inactivation gating variable
tau_h=zeros(size(V));

a_inf=zeros(size(V));   % A-current activation gating variable
tau_a=zeros(size(V));
b_inf=zeros(size(V));   % A-current inactivation gating variable
tau_b=zeros(size(V));

mca_inf=zeros(size(V));   % A-current activation gating variable
tau_mca=zeros(size(V));
hca_inf=zeros(size(V));   % A-current inactivation gating variable
tau_hca=zeros(size(V));


for i = 1:length(V); % now see how things change as a function of V
    
    Vm = V(i)*1000; % converts voltages to mV as needed in the equations on p.224 of Dayan/Abbott
    
    % Sodium and potassium gating variables are defined by the
    % voltage-dependent transition rates between states, labeled alpha and
    % beta. Written out from Dayan/Abbott, units are 1/ms.
    
    if ( Vm == -29.7 ) 
        alpha_m = 0.38/0.1;
    else 
        alpha_m = 0.38*(Vm+29.7)/(1-exp(-0.1*(Vm+29.7)));
    end
    beta_m = 15.2*exp(-0.0556*(Vm+54.7));

    alpha_h = 0.266*exp(-0.05*(Vm+48));
    beta_h = 3.8/(1+exp(-0.1*(Vm+18)));
    
    if ( Vm == -45.7 ) 
       alpha_n = 0.02/0.1;
    else
        alpha_n = 0.02*(Vm+45.7)/(1-exp(-0.1*(Vm+45.7)));
    end
    beta_n = 0.25*exp(-0.0125*(Vm+55.7));
    
    % From the alpha and beta for each gating variable we find the steady
    % state values (_inf) and the time constants (tau_) for each m,h and n.
    
    tau_m(i) = 1e-3/(alpha_m+beta_m);      % time constant converted from ms to sec
    m_inf(i) = alpha_m/(alpha_m+beta_m);
    
    tau_h(i) = 1e-3/(alpha_h+beta_h);      % time constant converted from ms to sec
    h_inf(i) = alpha_h/(alpha_h+beta_h);
    
    tau_n(i) = 1e-3/(alpha_n+beta_n);      % time constant converted from ms to sec
    n_inf(i) = alpha_n/(alpha_n+beta_n);   
    
    % For the A-type current gating variables, instead of using alpha and
    % beta, we just use the steady-state values a_inf and b_inf along with 
    % the time constants tau_a and tau_b that are found empirically
    % (Dayan-Abbott, p. 224)
    
    a_inf(i) = (0.0761*exp(0.0314*(Vm+94.22))/(1+exp(0.0346*(Vm+1.17))))^(1/3.0);
    tau_a(i) = 0.3632*1e-3 + 1.158e-3/(1+exp(0.0497*(Vm+55.96)));
    
    b_inf(i) = (1/(1+exp(0.0688*(Vm+53.3))))^4;
    tau_b(i) = 1.24e-3 + 2.678e-3/(1+exp(0.0624*(Vm+50)));
    
    % for the Ca_T current gating variables are given by formulae for the 
    % steady states and time constants (Dayan-Abbott p.225):
    
    mca_inf(i) = 1/(1+exp(-(Vm+57)/6.2));
    tau_mca(i) = 0.612e-3 + 1e-3/(exp(-(Vm+132)/16.7)+exp((Vm+16.8)/18.2));
    
    hca_inf(i) = 1/(1+exp((Vm+81)/4));
    if ( Vm < -80 ) 
        tau_hca(i) = 1e-3*exp((Vm+467)/66.6);
    else
        tau_hca(i) = 28e-3+1e-3*exp(-(Vm+22)/10.5);
    end
    

       
end