TITLE SDH Neuron - Hybrid Firing with AdEx 

COMMENT
AdEx + IL, Iadapt, Iadapt2, Isub, IK_A, IK_lt

Implemented by L Medlock 
ENDCOMMENT


NEURON {
  SUFFIX Hybrid_AdEx
  USEION k READ ek WRITE ik
  USEION na READ ena WRITE ina
  NONSPECIFIC_CURRENT il, iklt, ika, iadapt, isub, iadapt2
  RANGE gnabar, gkbar, gleak, gkltbar, gkabar, gadaptbar, gadaptslowbar, gsubbar, el, eklt, eadapt, eadaptslow, eka, esub, ina, ik, il, iklt, ika, iadapt, iadapt2, isub
  RANGE betam, gammam, phiw, betaw, gammaw, phiz, betaz, gammaz, phia, phib, betap, gammap, taup, betaq, gammaq, tauq, betapslow, gammapslow, taupslow
  GLOBAL minf, winf, tauw, zinf, tauz, ainf, taua, binf, taub, pinf, qinf, pslowinf, taub_1
  THREADSAFE
}

UNITS {
	(mV) = (millivolt)
	(mA) = (milliamp)
	(S)  = (siemens)
}

PARAMETER {
  v (mV)
  dt (ms)
  
  :Ileak
  gleak  = 0.00005 (S/cm2)
  el = -70 (mV)

  : IK_lt  (low-threshold non-inactivating (Kv1-type) potassium)
  gkltbar = 0.001 (S/cm2) : Changes per neuron type  
  eklt = -100 (mV)  
  phiz = 0.15
  betaz = -21 (mV)
  gammaz = 15 (mV)

  : IK_A (inactivating (A-type) potassium)
  gkabar = 0.001 (S/cm2)  : Changes per neuron type 
  eka = -100 (mV) 
  phia = 1.0
  phib = 1.0

  : Iadapt (AHP-like current) from Ratte et al. (2015)
  gadaptbar = 0.1 (S/cm2)       
  eadapt = -100 (mV) 
  betap = 0 (mV)
  gammap = 5 (mV)
  taup = 5 (ms)

  : INa
  gnabar = 0 (S/cm2)   
  ena = 50 (mV) 
  betam  = -1.2 (mV)
  gammam = 18 (mV)

  : IK_dr (delayed-rectifier K+ current)
  gkbar  = 0 (S/cm2)
  ek = -100 (mV)
  phiw = 0.15
  betaw = -10 (mV)
  gammaw = 10 (mV)

  : Iadapt2 (slower AHP) from Ratte et al. (2015)
  gadaptslowbar = 0 (S/cm2)  
  eadaptslow = -100 (mV)       
  betapslow = 0 (mV)
  gammapslow = 5 (mV)
  taupslow = 40 (ms)

  : Isub (subthreshold inward current) from Ratte et al. (2015)
  gsubbar = 0 (S/cm2)       
  esub = 50 (mV) 
  betaq = -40 (mV)
  gammaq = 10 (mV)
  tauq = 2 (ms)
}

ASSIGNED {
  ina (mA/cm2)
	ik (mA/cm2)
  il (mA/cm2)
  iklt (mA/cm2)
  ika (mA/cm2)
  iadapt (mA/cm2)
  isub (mA/cm2)
  iadapt2 (mA/cm2)
  gk (S/cm2)  
  gna (S/cm2) 
  gklt (S/cm2)
  gka (S/cm2)
  gadapt (S/cm2)
  gsub (S/cm2)
  gadaptslow (S/cm2)
	minf
	winf
	tauw (ms)
  zinf
  tauz (ms)
  ainf
  taua (ms)
  binf
  taub (ms)
  taub_1 (ms)
  pinf
  qinf
  pslowinf
}

STATE {
	m   FROM 0 TO 1
  w   FROM 0 TO 1
  z   FROM 0 TO 1
  a   FROM 0 TO 1
  b   FROM 0 TO 1
  p   FROM 0 TO 1
  q   FROM 0 TO 1
  pslow FROM 0 TO 1
}

INITIAL {
  rates(v)
  m = minf
  w = winf
  z = zinf
  a = ainf
  b = binf
  p = pinf
  q = qinf
  pslow = pslowinf
}

BREAKPOINT {
	SOLVE states METHOD cnexp
  gna = gnabar * m
  gk  = gkbar * w  
  gklt = gkltbar * z
  gka = gkabar * a^4 * b
  gadapt = gadaptbar * p
  gadaptslow = gadaptslowbar * pslow
  gsub = gsubbar * q

	ik  = gk * (v - ek)
  ina = gna * (v - ena)
  il  = gleak * (v - el)
  iklt = gklt * (v - eklt)
  ika = gka * (v - eka)
  iadapt = gadapt * (v - eadapt)
  iadapt2 = gadaptslow * (v - eadaptslow)
  isub = gsub * (v - esub)
}

DERIVATIVE states {
	rates(v)
  : Sodium Current (INa)
  m = minf
  : Delayed Rectifier K+ Current (IK_dr)
  w' = phiw*(winf-w)/tauw
  : Low-Threshold K+ Current (IK_lt)
  z' = phiz*(zinf-z)/tauz
  : A-type K+ Current (IK_A)
  a' = phia*(ainf-a)/taua
  b' = phib*(binf-b)/taub
  : Fast AHP Current (Iadapt)
  p' = (pinf-p)/taup
  : Slow AHP Current (Iadapt2)
  pslow' = (pslowinf-pslow)/taupslow
  : Inward Subthreshold Current (Isub)
  q' = (qinf-q)/tauq
}


PROCEDURE rates( v (mV) ) {
  TABLE minf, winf, tauw, zinf, tauz, ainf, taua, binf, taub_1, pinf, taub, pslowinf, qinf FROM -100 TO 100 WITH 200
  : Fast Activating Inward Current (INa):
  minf = 0.5*(1 + tanh((v-betam)/gammam))
  : Delayed Rectifier K+ Current (IK_dr)
  winf = 0.5*(1 + tanh((v-betaw)/gammaw))
  tauw = 1/(cosh((v-betaw)/(2*gammaw)))
  : Low-Threshold K+ Current (IK_lt)
  zinf = 0.5*(1 + tanh((v-betaz)/gammaz))
  tauz = 1/(cosh((v-betaz)/(2*gammaz)))
  : A-type K+ Current (IK_A)
  ainf = (1 / (1 + exp(-((v+65)/8.5))))
  taua = 1/(exp((v+35.82)/19.69)+exp(-(v+79.69)/12.7)) + 0.37
  binf = (1 / (1 + exp((v+78)/6)))
  taub_1 = 1/((exp((v+46.05)/5)+exp(-(v+238.4)/37.45)))
  if (v<-63) {
		taub = taub_1
		}
	else {
		taub = 19
		}
  : Fast AHP Current (Iadapt)
  pinf = (1 / (1 + exp((betap-v)/gammap)))
  : Slow AHP Current (Iadapt2)
  pslowinf = (1 / (1 + exp((betapslow-v)/gammapslow)))
  : Subthreshold Current (Isub)
  qinf = (1 / (1 + exp((betaq-v)/gammaq)))
 
}