"""
Cochlear neuron model of Rothman & Manis
----------------------------------------
Rothman JS, Manis PB (2003) The roles potassium currents play in
regulating the electrical activity of ventral cochlear nucleus neurons.
J Neurophysiol 89:3097-113.

All model types differ only by the maximal conductances.

Adapted from their Neuron implementation by Romain Brette
"""
from brian import *

#defaultclock.dt=0.025*ms # for better precision

'''
Simulation parameters: choose current amplitude and neuron type
(from type1c, type1t, type12, type 21, type2, type2o)
'''
neuron_type='type1c'
Ipulse=250*pA

C=12*pF
Eh=-43*mV
EK=-70*mV # -77*mV in mod file
El=-65*mV
ENa=50*mV
nf = 0.85 # proportion of n vs p kinetics
zss = 0.5 # steady state inactivation of glt
celsius=22. # temperature
q10 = 3.**((celsius - 22)/10.)
# hcno current (octopus cell)
frac=0.0
qt=4.5**((celsius-33.)/10.)

# Maximal conductances of different cell types in nS
maximal_conductances=dict(
type1c=(1000,150,0,0,0.5,0,2),
type1t=(1000,80,0,65,0.5,0,2),
type12=(1000,150,20,0,2,0,2),
type21=(1000,150,35,0,3.5,0,2),
type2=(1000,150,200,0,20,0,2),
type2o=(1000,150,600,0,0,40,2) # octopus cell
)
gnabar,gkhtbar,gkltbar,gkabar,ghbar,gbarno,gl=[x*nS for x in maximal_conductances[neuron_type]]

# Classical Na channel
eqs_na="""
ina = gnabar*m**3*h*(ENa-v) : amp
dm/dt=q10*(minf-m)/mtau : 1
dh/dt=q10*(hinf-h)/htau : 1
minf = 1./(1+exp(-(vu + 38.) / 7.)) : 1
hinf = 1./(1+exp((vu + 65.) / 6.)) : 1
mtau =  ((10. / (5*exp((vu+60.) / 18.) + 36.*exp(-(vu+60.) / 25.))) + 0.04)*ms : ms
htau =  ((100. / (7*exp((vu+60.) / 11.) + 10.*exp(-(vu+60.) / 25.))) + 0.6)*ms : ms
"""

# KHT channel (delayed-rectifier K+)
eqs_kht="""
ikht = gkhtbar*(nf*n**2 + (1-nf)*p)*(EK-v) : amp
dn/dt=q10*(ninf-n)/ntau : 1
dp/dt=q10*(pinf-p)/ptau : 1
ninf =   (1 + exp(-(vu + 15) / 5.))**-0.5 : 1
pinf =  1. / (1 + exp(-(vu + 23) / 6.)) : 1
ntau =  ((100. / (11*exp((vu+60) / 24.) + 21*exp(-(vu+60) / 23.))) + 0.7)*ms : ms
ptau = ((100. / (4*exp((vu+60) / 32.) + 5*exp(-(vu+60) / 22.))) + 5)*ms : ms
"""

# Ih channel (subthreshold adaptive, non-inactivating)
eqs_ih="""
ih = ghbar*r*(Eh-v) : amp
dr/dt=q10*(rinf-r)/rtau : 1
rinf = 1. / (1+exp((vu + 76.) / 7.)) : 1
rtau = ((100000. / (237.*exp((vu+60.) / 12.) + 17.*exp(-(vu+60.) / 14.))) + 25.)*ms : ms
"""

# KLT channel (low threshold K+)
eqs_klt="""
iklt = gkltbar*w**4*z*(EK-v) : amp
dw/dt=q10*(winf-w)/wtau : 1
dz/dt=q10*(zinf-z)/wtau : 1
winf = (1. / (1 + exp(-(vu + 48.) / 6.)))**0.25 : 1
zinf = zss + ((1.-zss) / (1 + exp((vu + 71.) / 10.))) : 1
wtau = ((100. / (6.*exp((vu+60.) / 6.) + 16.*exp(-(vu+60.) / 45.))) + 1.5)*ms : ms
ztau = ((1000. / (exp((vu+60.) / 20.) + exp(-(vu+60.) / 8.))) + 50)*ms : ms
"""

# Ka channel (transient K+)
eqs_ka="""
ika = gkabar*a**4*b*c*(EK-v): amp
da/dt=q10*(ainf-a)/atau : 1
db/dt=q10*(binf-b)/btau : 1
dc/dt=q10*(cinf-c)/ctau : 1
ainf = (1. / (1 + exp(-(vu + 31) / 6.)))**0.25 : 1
binf = 1. / (1 + exp((vu + 66) / 7.))**0.5 : 1
cinf = 1. / (1 + exp((vu + 66) / 7.))**0.5 : 1
atau =  ((100. / (7*exp((vu+60) / 14.) + 29*exp(-(vu+60) / 24.))) + 0.1)*ms : ms
btau =  ((1000. / (14*exp((vu+60) / 27.) + 29*exp(-(vu+60) / 24.))) + 1)*ms : ms
ctau = ((90. / (1 + exp((-66-vu) / 17.))) + 10)*ms : ms
"""

# Leak
eqs_leak="""
ileak = gl*(El-v) : amp
"""

# h current for octopus cells
eqs_hcno="""
ihcno = gbarno*(h1*frac + h2*(1-frac))*(Eh-v) : amp
dh1/dt=(hinfno-h1)/tau1 : 1
dh2/dt=(hinfno-h2)/tau2 : 1
hinfno = 1./(1+exp((vu+66.)/7.)) : 1
tau1 = bet1/(qt*0.008*(1+alp1))*ms : ms
tau2 = bet2/(qt*0.0029*(1+alp2))*ms : ms
alp1 = exp(1e-3*3*(vu+50)*9.648e4/(8.315*(273.16+celsius))) : 1
bet1 = exp(1e-3*3*0.3*(vu+50)*9.648e4/(8.315*(273.16+celsius))) : 1 
alp2 = exp(1e-3*3*(vu+84)*9.648e4/(8.315*(273.16+celsius))) : 1
bet2 = exp(1e-3*3*0.6*(vu+84)*9.648e4/(8.315*(273.16+celsius))) : 1
"""

eqs="""
dv/dt=(ileak+ina+ikht+iklt+ika+ih+ihcno+I)/C : volt
vu = v/mV : 1 # unitless v
I : amp
"""
eqs+=eqs_leak+eqs_ka+eqs_na+eqs_ih+eqs_klt+eqs_kht+eqs_hcno

neuron=NeuronGroup(1,eqs,implicit=True)
neuron.v=El

run(50*ms) # Go to rest

M=StateMonitor(neuron,'v',record=0)
neuron.I=Ipulse

run(100*ms,report='text')

plot(M.times/ms,M[0]/mV)
show()