%   hvcmodel_ode.m

%   This function contains the right hand sides for the odes, and is called
%   by the various single neuron codes.

function dydt= hvcmodel_ode(t, y, parameters)

    Iap= parameters(1);
	
    % Nerst potentials in mV
    VL= parameters(2);
	VK = parameters(3);
	VNa= parameters(4);
    VH=parameters(5);

    
    % Conductances in nS for all the neurons
    gL   = parameters(6);
    gNap = parameters(7);
    gCaL = parameters(8);
    
    % Conductances in nS for specific the neurons
    % HVCx neuron
    gNa  = parameters(9);
    gK   = parameters(10);
    gA   = parameters(11);
    gM   = parameters(12);
    gCaT = parameters(13);
    gSK  = parameters(14);
    gH   = parameters(15);
    kr   = parameters(16);
       C = parameters(17);
    
    % Time constants in mS
    tauNbar  = parameters(18);
    
    tauHpbar = parameters(19);
    tauE     = parameters(20);
    tauZ     = parameters(21);
    tauR0    = parameters(22);
    tauR1    = parameters(23);
    
    
	thetaM  = parameters(24);
	thetaN  = parameters(25);
	thetaMp = parameters(26);
    thetaHp = parameters(27);
    thetaA  = parameters(28);
    thetaE  = parameters(29);
    thetaZ  = parameters(30);
    thetaS  = parameters(31);
    thetaAT  = parameters(32);
    thetaBT  = parameters(33);
    thetaRT  = parameters(34);
    thetaRRT  = parameters(35);
    thetaRf  = parameters(36);
    thetaRs  = parameters(37);
    
    sigmaM   = parameters(38);
    sigmaN   = parameters(39);
    sigmaMp  = parameters(40);
    sigmaHp  = parameters(41);
    sigmaA   = parameters(42);
    sigmaE   = parameters(43);
    sigmaZ   = parameters(44);
    sigmaS   = parameters(45);
    sigmaAT   = parameters(46);
    sigmaBT   = parameters(47);
    sigmaRT   = parameters(48);
    sigmaRRT   = parameters(49);
    sigmaRf   = parameters(50);
    sigmaRs   = parameters(51);
    
    % Other parameters
    Ca_ex = parameters(52);
    RTF   = parameters(53); 
    ks    = parameters(54); 
    f     = parameters(55); 
    eps   = parameters(56); 
    kCa   = parameters(57); 
    bCa   = parameters(58); 
    prf   = parameters(59);
    ton   = parameters(60);
    toff  = parameters(61);
    WtauH  = parameters(62);
    
	V  = y(1);
	n  = y(2);
    h  = y(3);
    hp = y(4);
    e  = y(5);
    z  = y(6);
    rT = y(7);
    Cai=y(8); 
    rf=y(9);
    rs=y(10);
    
    % Leak current
    iL = gL*(V-VL);

    % K+ Current (iK)   
    ninf = 1./(1+exp((V-thetaN)./sigmaN));
    tauN = tauNbar./cosh((V-thetaN)./(2*sigmaN));
    iK = gK*(n^4)*(V-VK);
    
    % Na+ Current (iNa)  
    minf =  1./(1+exp((V-thetaM)./sigmaM));
    alphah =  0.128*exp((V+50)./-18);
    betah  =  4./(1+exp((V+27)./-5));
    hinf   =  alphah ./(alphah+betah);
    tauH   = (1.2)./(alphah+betah);
    iNa    =  gNa*(minf^3)*h*(V-VNa);

    % Nap+ current (iNap)
    mpinf = 1./(1+exp((V-thetaMp)./sigmaMp));
    hpinf = 1./(1+exp((V-thetaHp)./sigmaHp));
    tauHp = tauHpbar./cosh((V-thetaHp)./(2*sigmaHp));
    iNap  = gNap*mpinf*hp*(V-VNa);
    
    % A+ Current (iA)   
    ainf = 1./(1+exp((V-thetaA)./sigmaA));
    einf = 1./(1+exp((V-thetaE)./sigmaE));
    iA   = gA*ainf*e*(V-VK);
    
    % M+ Current (iM)   
    zinf = 1./(1+exp(-(V-thetaZ)./sigmaZ));
    iM   = gM*z*(V-VK);
    
    % High-threshold L-type Ca 2+ Current (iCaL)
    sinf = 1./(1+exp((V-thetaS)./sigmaS));
    iCaL = gCaL*(sinf^2)*V*(Ca_ex./(1-exp((2*V)./RTF)));
    
    % Low-threshold T-type Ca 2+ current (iCaT)
    aTinf = 1./(1+exp((V-thetaAT)./sigmaAT));
    bTinf = 1./(1+exp((rT-thetaBT)./sigmaBT))-1./(1+exp(-thetaBT./sigmaBT));
    rTinf = 1./(1+exp((V-thetaRT)./sigmaRT));
    tauRT = (tauR0 + tauR1./(1+exp((V-thetaRRT)./sigmaRRT)));
    iCaT  = gCaT*(aTinf^3)*(bTinf^3)*V*(Ca_ex./(1-exp((2*V)./RTF)));
    
    % Ca 2+ dependent K+ current (iSK)
    kinf  = (Cai^2)./(Cai^2 + ks^2);
    iSK   = gSK*kinf*(V-VK);
    
    % Hyperpolarization-activated inward Current (iH)
    rfinf = 1./(1+exp((V-thetaRf)./sigmaRf));
    rsinf = 1./(1+exp((V-thetaRs)./sigmaRs));
    tauRf = prf./(-7.4*(V+70)./(exp(-(V+70)./0.8)-1)+65*exp(-(V+56)./23));
    tauRs = exp((V+289.71)./33.3);
    iH    = gH*(kr*rf+(1-kr)*rs)*(V-VH);
      
    if t>ton && t<toff
        Iapp=Iap;
    else
        Iapp=0;
    end
    
    dydt(size(y,1),1)= 0;
	dydt(1)= +((-1./C)*(iL+iK+iNa+iNap+iA+iM+iCaL+iCaT+iSK+iH-Iapp));
	dydt(2)= +((ninf-n)./(tauN)); 
	dydt(3)= +((hinf-h)./(tauH));
    dydt(4)= +((hpinf-hp)./(tauHp));
    dydt(5)= +((einf-e)./tauE);
    dydt(6)= +((zinf-z)./tauZ);
    dydt(7)= +((rTinf-rT)./tauRT);
    dydt(8)= +(-f*(eps*(iCaL + iCaT)+ kCa*(Cai-bCa)));
    dydt(9)= +((rfinf-rf)./tauRf);
    dydt(10)= +((rsinf-rs)./tauRs);

	
