blockhh3.mod

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, tfast
        RANGE taukv4,asshift,fptog, qfix,qinf,pinf,ninf,ntau, qtau,ptau,kchip,fchip, qfast, qslow,gka2,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
        qfast = 0
        qslow = 0
        tfast = 25
 	
}
 
STATE {
         n <1e-3> p <1e-3> qone <1e-3> qtwo <1e-3> p2 <1e-3>
}
 
ASSIGNED {
        ik (mA/cm2)
        ika (mA/cm2)
        ikhh (mA/cm2)
		ninf qinf pinf hsinf p2inf q2inf
		qtau ptau ntau q2tau p2tau
		gk gka gka2
		q1 q2
}
 
BREAKPOINT {
        SOLVE states METHOD cnexp
  	gk = n*n*n

	if(fabs(fptog-1.0)<1e-6){ : neuron can be wonky with boolean-esque values
	  q2 = qfast*fchip
	  q1 = qslow*kchip
	} else {
	  q1 = qone*kchip
	  q2 = qtwo*fchip
	}
	gka = p*p*p*q1+p2*p2*p2*q2 : single activation state
	gka2 = gka*gkabar
    ikhh = gkhhbar*gk*(v - ek)      
    ika = gka2*(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, va

        
        ninf = boltz(vf,-20,nslope)

        
        if(vf>-40.8+2*ashift){
			p2inf = boltz(vf,-40.8+2*ashift,7.73) : fit to data, there will be some discrepency around window current
		} else {
			p2inf = boltz(vf,-40.8+2*ashift,5.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 + 2.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
          ntau = 1.0 +9.0*exp(-((vf+40+nshift))*((vf+40+nshift))/15.0) : peak tau at -40 mV, tau ~5ms at 0 :5 5 9

        } 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 +(10.0-nspeed)*exp(-((vf+40+nshift))*((vf+40+nshift))/15.0) : peak tau at -40 mV, tau ~5ms at 0 :5 5 9
          :ntau = 3.0
        }
        ntau =1.0*ntau
        p2tau = (1.029 + 4.83*boltz(vf,-56.7+ashift,-6.22))
        ptau= p2tau : slower activation
        :ptau = p2tau
        qtau = (50.0+(taukv4*250-25)*boltz(vf,-65+ashift,4.0)) : saturates after 100 ms per Tarfa et al
        q2tau = tfast:(10.0+65.0*boltz(vf,-65+ashift,4.0)) : faster may be needed for pacing : fast component inactivation


}
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