-
Notifications
You must be signed in to change notification settings - Fork 1
/
TypeISomaCablewHCNFig3.ode
109 lines (72 loc) · 2.31 KB
/
TypeISomaCablewHCNFig3.ode
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
###########################################################################
# This model has a delayed rectifier, a leak and
# a sodium current based on the famous Traub and Miles model but with changes put in
# by Bard Ermentrout in linearization of F-I curves by adaptation, Neural computation,
# 10:1721-9, 1998
#
# Input Currents are in nA, maximal conductances in mS/cm2, membrane potential
# in mV time in ms, and capacitance in uF/cm2
############################################################################
#adjust defaults for xpp internal variables
@ MAXSTOR=40000
@ YP=v
@ TOTAL=1000
@ DT=0.05
@ BOUND=100000
@ XHI=1000
@ YLO=-100
@ YHI=50
#@ METH=gear
@ BACKGROUND=white
# parameters
par gbarNa=80
par gbarK=80
par gleak=0.2
par C=1., i=1
v' = (i-(gleak*(v-Vl)+gK*(v-vk)+gNa*(v-vna))+eps*(u1-v)/dx)/C+p0*pulse(t-taup)
aux longcur=eps*(u1-v)/dx
#****************************************************
gNa=gbarNa*(m**3)*h
alpham=if (v+54) then (0.32*(V+54)/(1-exp(-(V+54)/4))) else (0.32)
betam=if (V+27) then (0.28*(V+27)/(exp((V+27)/5)-1)) else (0.28)
alphah=0.128*exp(-(V+50)/18)
betah=4/(1+exp(-(V+27)/5))
dm/dt=alpham*(1-m)-betam*m
dh/dt=alphah*(1-h)-betah*h
#**********************************************************
gK=gbarK*(n**4)
alphan=if (V+52) then (0.032*(V+52)/(1-exp(-(V+52)/5))) else (0.032)
betan=0.5*exp(-(V+57)/40)
dn/dt=alphan*(1-n)-betan*n
# auxiliary variables
# sod is sodium current
aux sod=gbarNa*(m**3)*h*(v-vna)
# pot is potassium current
aux pot=gbarK*(n**4)*(v-vk)
#initial conditions
init v=-70
init h=1-.000001,m=.000001,n=.000001
par Vna=50,Vk=-100,Vl=-67,Vh=-65,Vrh=-40,Vsh=6,gld=0.1,gh=0.05,tauh=400
hinf(V)=1/(1+exp((V-Vh)/Vsh))
Ild(V)=V-Vl
Ih(V,y)=gh*y*(V-Vrh)/gld
h[1..50]'=(hinf(u[j])-h[j])/tauh
# NOT TO CONFUSE WITH h GATE IN SOMA!!
# cable equation has HCN current
p p[0..50]=0
u1'=((lambda/dx)^2*(u2-2*u1+v)-Ild(u1)-Ih(u1,h1))/tau +p1*pulse(t-taup)
u[2..50]'= ((lambda/dx)^2*(u[j+1]-2*u[j]+u[j-1])-Ild(u[j])-Ih(u[j],h[j]))/tau +p[j]*pulse(t-taup)
u51=(c1+b1*u50/dx)/(a1+b1/dx)
par lambda=1,tau=10,dx=.1,c1=0,a1=0,b1=1,c0=0,a0=0,b0=1,eps=.025
pulse(t)=heav(t)*heav(sigma-t)
par sigma=.2
par t0=241.3
aux prc=t0-t
taup'=0
i u[1..50]=-60
i h[1..50]=.1
par taur=1,taud=3,thresh=-30, gsyn=.1 Esyn=0
x1'=(-x1+.5*(1+tanh((v-thresh)/3.0)))/taur
y1'=(-y1+x1)/taud
init x1=0,y1=0
d