TITLE  potassium delayed rectifier, and potassium A channels
: Anna's version adding some shift values
: Kv4 from Tarfa et al 2017
: sodium from Yu et al 2014
 
UNITS {
        (molar) = (1/liter)
        (S) = (siemens)
        (mA) = (milliamp)
        (mV) = (millivolt)
         F = (faraday) (coulomb)
         R = (mole k) (mV-coulomb/degC)
        (mM) =  (millimolar)
}
 
NEURON {
        SUFFIX hhb
        USEION k WRITE ik
        RANGE  nspeed,gkhhbar,gkabar,ikhh,ika,ik,gk,gka,nshift,bkcor, ashift, nslope, taukv4,asshift,fptog, qfix,qinf,pinf,ninf,ntau, qtau,ptau,kchip,fchip,q1,q2
        :GLOBAL minf,hinf,ninf
}
 
INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}
 
PARAMETER {
        v   (mV)
        dt  (ms)
		celsius = 35.0 (degC)
        gkhhbar = 665.0e-6 (S/cm2)
        gkabar  = 266.0e-6  (S/cm2) :kv4 conductance should be ~80% of gkhh
        ek  = -90.0  (mV)
        qv = 66.0 (mV)
        qs = 10.0  (1)
        ashift = 0.0 (mV)
        nshift = 0 (mV)
        nslope = 12
        bkcor = 0.0 (ms)
        taukv4 = 1.0
        asshift = 0
        fptog = 0
        qfix = 0
        nspeed = 2
        kchip=1
        fchip=0
 	
}
 
STATE {
         n <1e-4> p <1e-4> qone <1e-4> qtwo <1e-4> p2 <1e-4>
}
 
ASSIGNED {
        ik (mA/cm2)
        ika (mA/cm2)
        ikhh (mA/cm2)
		ninf qinf pinf hsinf p2inf q2inf
		qtau ptau ntau q2tau p2tau
		gk gka
		q1 q2
}
 
BREAKPOINT {
        SOLVE states METHOD cnexp
  	gk = n*n*n
  	q1 = qone*kchip
  	q2 = qtwo*fchip
	if(fptog){
	  gka = p*p*p*qfix
	} else {
	  gka = p*p*p*q1+p2*p2*p2*q2 : single activation state
	}
    ikhh = gkhhbar*gk*(v - ek)      
    ika = gkabar*gka*(v - ek)      
    ik = ika + ikhh

}
 
UNITSOFF
 


: a-type modified to model from Latorre 2012
: this model features better inactivation kinetics
PROCEDURE rates(vf){
LOCAL ahf, bhf, ahs, bhs, am, bm

        
        ninf = boltz(vf,-20,nslope)

        
        
        p2inf = boltz(vf,-40.8,7.73) : fit to data, there will be some discrepency around window current
        :p2inf = boltz(vf,-25,8.0) : 25 activation, with n=4, this makes p4inf at -37 mV (28 originally)
        qinf = boltz(vf,-68.0+ashift,-6.44+asshift) : inactivation
        q2inf = boltz(vf,-68.0+ashift,-6.44+asshift)
        p2inf = pow(p2inf,1.0/3.0)
        
        :pinf = boltz(vf,-25.0,8.0)
        :pinf = pow(pinf,1.0/3.0)
        pinf = p2inf


        

        if (v > -40-nshift){
          ntau = 1 + 5.0*exp(-(log(1.0+0.05*(vf+40+nshift))/0.05*log(1.0+0.05*(vf+40+nshift))/0.05)/300) : peak tau at -40 mV, tau ~5ms at 0
        } else {
          :ntau = 2.0 + 8.0*exp(-((vf+40+nshift))*((vf+40+nshift))/18.0) : peak tau at -40 mV, tau ~5ms at 0 :5 5 9
          ntau = nspeed +(6-nspeed)*exp(-((vf+40+nshift))*((vf+40+nshift))/8.0) : peak tau at -40 mV, tau ~5ms at 0 :5 5 9
          :ntau = 10.0
        }
        ntau =1.0*ntau
        p2tau = (1.029 + 4.83*boltz(vf,-56.7,-6.22))
        ptau= p2tau : slower activation
        :ptau = p2tau
        qtau = (20.0+(taukv4*230)*boltz(vf,-60,5.0)) : recovers at comparable rate to fast components
        q2tau = (20.0+25.0*boltz(vf,-60,5.0)) : faster may be needed for pacing 


}
INITIAL {
        rates(v)

        n = ninf
        p = pinf
        qone =qinf
        qtwo = q2inf
        p2 = p2inf

}

DERIVATIVE states {  :Computes state variables m, h, and n 
        rates(v)

        n' = 1.0*(ninf-n)/ntau
        p' = (pinf-p)/ptau
        qone' = (qinf-qone)/qtau
        qtwo' = (q2inf-qtwo)/q2tau
        p2' = (p2inf-p2)/p2tau
}
 
 
 
FUNCTION boltz(x,y,z) {
                boltz = 1/(1 + exp(-(x - y)/z))
}
 
UNITSON