-
Notifications
You must be signed in to change notification settings - Fork 1
/
TypeIISomaCablePairFig2.ode
87 lines (63 loc) · 2.16 KB
/
TypeIISomaCablePairFig2.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
# Soma is Hodgkin Huxley
p Ek=-77,Ena=50, El=-54.4
p gl=.3,gkdr=36,gna=120
p i=11
p C=1.
p amphi=.1,amhalf=-40,amwidth=10
p bmphi=4,bmhalf=-65,bmwidth=18
p ahphi=0.07,ahhalf=-65,ahwidth=20
p bhphi=1,bhhalf=-35,bhwidth=10
p anphi=.01,anhalf=-55,anwidth=10
p bnphi=.125,bnhalf=-65,bnwidth=80
# functions
am(v)=amphi*(v-amhalf)/(1-exp(-(v-amhalf)/amwidth))
bm(v)=bmphi*exp(-(v-bmhalf)/bmwidth)
ah(v)=ahphi*exp(-(v-ahhalf)/ahwidth)
bh(v)=bhphi/(1+exp(-(v-bhhalf)/bhwidth))
an(v)=anphi*(v-anhalf)/(1-exp(-(v-anhalf)/anwidth))
bn(v)=bnphi*exp(-(v-bnhalf)/bnwidth)
#currents
ina(v,m,h)=gna*m^3*h*(v-Ena)
ikdr(v,n)=gkdr*n^4*(v-Ek)
il(v)=gl*(v-El)
Isyn(v,y)=gsyn*y*(v-Esyn)
#diff. equ.
v1'=(i-(ina(v1,m1,ha)+ikdr(v1,n1)+il(v1)+p0*Isyn(v1,y2))+eps*(ua1-v1)/dx)/C
v2'=(i-(ina(v2,m2,hb)+ikdr(v2,n2)+il(v2)+p0*Isyn(v2,y1))+eps*(ub1-v2)/dx)/C
m1'=am(v1)*(1-m1)-bm(v1)*m1
m2'=am(v2)*(1-m2)-bm(v2)*m2
n1'=an(v1)*(1-n1)-bn(v1)*n1
n2'=an(v2)*(1-n2)-bn(v2)*n2
ha'=ah(v1)*(1-ha)-bh(v1)*ha
hb'=ah(v2)*(1-hb)-bh(v2)*hb
## synapse
parameter taur=1,taud=3,thresh=-30
x1'=(-x1+.5*(1+tanh((v1-thresh)/3.0)))/taur
x2'=(-x2+.5*(1+tanh((v2-thresh)/3.0)))/taur
y1'=(-y1+x1)/taud
y2'=(-y2+x2)/taud
init x1=.001, y1=.001,x2=.001, y2=.001
p gsyn=.1, Esyn=0
p Vp=-50, Vsp=9, gnad=0.02, gld=.1, taupna=10
# !!!! the cable is passive if gnad=0 !!!!
pinfd(V)=1/(1+exp(-(V-Vp)/Vsp))
Ih(V,y)=gnad*y*(V-Ena)/gld
Ild(V)=V-El
ha[1..50]'=(pinfd(ua[j])-ha[j])/taupna
hb[1..50]'=(pinfd(ub[j])-hb[j])/taupna
# NOT TO CONFUSE WITH h GATE IN SOMA!!
# cable equation
ua1'=((lambda/dx)^2*(ua2-2*ua1+v1)-Ild(ua1)-Ih(ua1,ha1)-p1*Isyn(ua1,y2)/gld)/tau
ua[2..50]'= ((lambda/dx)^2*(ua[j+1]-2*ua[j]+ua[j-1])-Ild(ua[j])-Ih(ua[j],ha[j]) -p[j]*Isyn(ua[j],y2)/gld )/tau
ua51=(c1+b1*ua50/dx)/(a1+b1/dx)
ub1'=((lambda/dx)^2*(ub2-2*ub1+v2)-Ild(ub1)-Ih(ub1,hb1)-p1*Isyn(ub1,y1)/gld)/tau
ub[2..50]'= ((lambda/dx)^2*(ub[j+1]-2*ub[j]+ub[j-1])-Ild(ub[j])-Ih(ub[j],hb[j])-p[j]*Isyn(ub[j],y1)/gld )/tau
ub51=(c1+b1*ub50/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=.05
par t0=14.45
aux prc=t0-t
p p[0..50]=0
@ total=300,xlo=0,xhi=300,ylo=-100,yhi=60,dt=0.05,bounds=10000000
d