"""
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 for Brian 1 from their Neuron implementation by Romain Brette,
updated for Brian 2 with additional functionality by Thomas McCall.
"""
from brian import *
import time

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

'''
Simulation parameters: choose current amplitude and neuron type
(from type1c, type1t, type12, type 21, type2, type2o)
'''

start = time.time()

neuron_type='type1c'
Ipulse=250 * pA
ncell = 1

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(ncell,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')

finish = time.time()

print "Runtime: ", finish-start

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