126899/Rothman_Manis_2003.py at master · ModelDBRepository/126899 · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
"""
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()