Simnon Programs Used in the Text,

INTRODUCTION to DYNAMIC MODELING of NEURO-SENSORY SYSTEMS

Presented in the order in which they appear.  

Note: all programs are Text Only (*.txt) in Word for Windows format.  They must be saved with the *.t extension rather than *.txt in order to run as Simnon programs. Quotes (") tell Simnon to ignore symbols to their right.  I.e., they denote comments. Simnon is case-insensitive, so for example, Jin = jin.


continuous system HODHUX  " THE original Hodgkin-Huxley equations. 6/06/99. 
"(From Sec. 1.4.1.)
" Run w/ Euler integration w/ dt = 0.00001.
STATE v m n h             " v in mV.  v is depolarization if < 0.
DER dv dm dn dh           " v = Vmo - Vm.  Vm is actual transmembrane V.
"                           Vmo = resting potential = - 70 mV.
TIME t                    " t in ms.
"
" HH membrane patch ODE. (1 cm^2)
dv = Jin/Cm - Jk/Cm - Jna/Cm - Jl/Cm  " ODE from Kirchoff's current law.
"
" Ionic current densities.  
Jk = gko*(n^4)*(v - Vk)			   
Jna = gnao*(m^3)*h*(v - Vna)
Jl = glo*(v - Vl)
Jc = Cm*dv
Jnet = Jc + Jk + Jna + Jl  " Microamps/cm^2
"
" Ionic conductances:
gk = gko*(n^4)		   " Potassium conductance
gna = gnao*(m^3)*h	   " Sodium conductance
gl = glo			   " Leakage conductance
gnet = gk + gna + gl	   " Net membrane conductance
negv = -v
"
dn = - n*(an + bn) + an    " K+ activation parameter ODE.
dm = - m*(am + bm) + am    " Na+ activation parameter ODE.
dh = - h*(ah + bh) + ah    " Na+ inactivation parameter ODE.
"
"  VOLTAGE-DEPENDENT PARAMETER FUNCTIONS: 
an = 0.01*(v + 10)/(exp(0.1*v + 1) - 1)
bn = 0.125*exp(v/80)
am = 0.1*(v + 25)/(exp(0.1*v + 2.5) - 1)
bm = 4*exp(v/18)
ah = 0.07*exp(v/20)
bh = 1/(exp(0.1*v + 3) + 1)
"
Vm = -(v + 70)
Jksc = Jk/10   " Scaled current densities
Jnasc = Jna/10
Jcsc = 10*Jc   " Scaled capacitor current density.
Vmsc = Vm/70   " Scaled membrane voltage.
"
" CONSTANTS:
zero:0 
Vk:12         " mV
Vna:-115      " mV
VL:-10.613    " mV
"
glo:0.3       " milliS/cm^2.
gko:36        "   "
gnao:120      "   "
Cm:.001       " milliF/cm^2 = 1 microF/cm^2. 
f:0.10        " Hz
pi:3.14159
"
" INPUTS:
"Jin1 = IF t > to THEN K*(t - to) ELSE 0   " Delayed current ramp input.
"Jin2 = IF t > t2 THEN -K*(t - t2) ELSE 0  " Jin in microamps/cm^2.
t2:10                                      " End of ramp
to:1                                       " Jin < 0 depolarizes membrane. 
K:-0.3                                     " Jin inward is + charges inward
Jin1 = IF t > to THEN Jino ELSE 0          " Delayed current pulse input.
Jin2 = IF t > (to + delt) THEN -Jino ELSE 0 
delt:5					   	 " Duration of current pulse, ms.
Jino:-2                                    " Microamps.
Jinac:1
" Jin1 = IF t > to THEN Jino ELSE 0        " Delayed current step input.
" Jin2 = Jinac*sin(2*pi*f*t)               " Sinusoidal current input.
Jin = Jin1 + Jin2 
"
" INITIAL CONDS. (Also final values for v = 0.)
m:0.052932 
n:0.31768
h:0.59612
v:0
"
END_


continuous system HHVCLAMP  " 6/04/99  In Sec. 1.4.3.
" Run w/ Euler integration w/ dt = 0.0001.
STATE v n m h           " v in mV.  v is depolarization if < 0.
DER dv dn dm dh         " v = Vmo - Vm.  Vm is actual transmembrane V.
"                         Vmo = resting potential = - 70 mV.
TIME t                  " t in ms.
"
" HH membrane patch ODE.
dv = Jin/Cm - Jl/Cm - Jk/Cm - Jna/Cm   
"
" Ionic current densities.  
Jk = gko*(n^4)*(v - Vk)			   
Jna = gnao*(m^3)*h*(v - Vna)
Jl = glo*(v - Vl)
Jc = Cm*dv
Jnet = Jc + Jl + Jk + Jna  " Microamps/cm^2
"
dn = - n*(an + bn) + an    " K+ activation parameter ODE.
dm = - m*(am + bm) + am    " Na+ activation parameter ODE.
dh = - h*(ah + bh) + ah    " Na+ inactivation parameter ODE.
"
"  VOLTAGE-DEPENDENT PARAMETER FUNCTIONS: 
an = 0.01*(v + 10)/(exp(0.1*v + 1) - 1)
bn = 0.125*exp(v/80)
am = 0.1*(v + 25)/(exp(0.1*v + 2.5) - 1)
bm = 4*exp(v/18)
ah = 0.07*exp(v/20)
bh = 1/(exp(0.1*v + 3) + 1)
"
Vm = -(v + 70)
Jksc = Jk/10   " Scaled current densities
Jnasc = Jna/10
Jcsc = 10*Jc   " Scaled capacitor current density.
Vmsc = Vm/70   " Scaled membrane voltage.
"
" VOLTAGE CLAMP CONTROL ALGORITHM:
Jin = Kc*(Vs - v - 70) 
"
" CONSTANTS:
zero:0
Kc:1.E4	  " Controller gain, Kc = (1/2)Kv*Gc
Vk:12         " mV
Vna:-115      " mV
VL:-10.613    " mV
glo:0.3       " mS/cm^2.
gko:36        "   "
gnao:120      "   "
Cm:1          " mF/cm^2. 
"
" INPUTS: Jin < 0 depolarizes membrane. Jin inward is + charges inward
"Vs1 = IF t > to THEN Vso ELSE 0  " Delayed voltage pulse input.
"Vs2 = IF t > (to + delt) THEN -Vso ELSE 0 
delt:14					 " Duration of voltage pulse, ms.
to:1					 	 " ms.                                       
Vso:-20 					 " mV.
"Vs = Vs1 + Vs2 + 70			 " -Vso sets Vm
Vs = 70 + ramp
ramp = IF t > to THEN Kr*(t - to) ELSE 0
Kr:-3 
"
" INITIAL CONDS. (Also final values for v = 0.)
m:0.052932 
n:0.31768
h:0.59612
v:0
"
END


continuous system ADTHRESH  " 7/02/99. System uses NFB to raise phi2 of 
" RPFM sensory neuron to reduce false positives. (In Sec.2.1.4)
STATE v1 v2 r1 r2 p2 q2 Vmn in ein
DER dv1 dv2 dr1 dr2 dp2 dq2 dVmn din dein
TIME t
"  Use Euler Integration w/ dt = tau.
"  
"  IPFM VPC GENERATES RANDOM QUM EVENTS:
dv1 = rin - z1
z1 = y1*phi1/tau
w1 = IF v1 > phi1 THEN 1 ELSE 0
s1 = DELAY(w1, tau)
x1 = w1 - s1
y1 = IF x1 > 0 THEN x1 ELSE 0
e1 = z1*Do1
"
"  THE RPFM SGL
dv2 = -c2*v2 + c2*ein + vmn - z2    " Vmn added to v2 = Vm.
z2 = y2*phi2/tau
w2 = IF v2 > phi2 THEN 1 ELSE 0
s2 = DELAY(w2, tau)
x2 = w2 - s2
y2 = IF x2 > 0 THEN x2 ELSE 0
"
"  NOISE GENERATION:
dr1 = -b1*r1 + SD1*NORM(t) 		" Membrane noise, 2-pole filtered.
dVmn = -b2*Vmn + r1
"
dr2 = -b3*r2 + SD2*NORM(t + Td) 	" Input noise = in drives input events.
din = -b4*in + r2				
rin = IF in > 0 THEN in ELSE 0	" Rectified input noise to IPFM random 
"						  event generator.
dp2 = -a2*p2 + z2*K2	" 2-pole ballistic LPF for threshold feedback	
dq2 = -a3*q2 + p2 
"
dein = -a1*ein + e1	" LPF to condition random event impulses.
"
PHI2 = phi20 + Kf*q2	" Sensory neuron SGL threshold is sum of 
"				  membrane noise + slow NFB term.
plty1 = 0.3*y1 + O1
plty2 = O2 + 0.3*y2
"
"  CONSTANTS:
phi1:2.5 
phi20:0.75
Do1:15.
K2:.2   	" Adjusts gain of BF.
Kf:1	 	" Kf = 0 turns off NFB.
a1:1.
a2:.1		" a,b,c units are radians/s.
a3:.05
b1:1   
b2:2
b3:1.
b4:1. 
c2:.05
Td:100
tau:.001
SD1:50
SD2:8.
O1:-1.5
O2:-1 
zero:0
"
END


Continuous system ARENcns3   " v. 2/27/99 (from Appendix 1)
"  Use EULER integration with delT = tau.
"  Neural model to simulate enhancement of positional orientation 
"  by protocerebral PPI units, as recorded by Willey, 1981.
"  There are 47 states.
STATE r1 r2 r3 r4 
STATE v1 v2 v3 v4 v5 v6 v7 
STATE p1 q1 p2 q2 p3 q3 p4 q4 p5 q5 p6 q6 p7 q7 p8 q8 p9 q9 
STATE p10 q10 p11 q11 p12 q12 p13 q13 p14 q14 p15 q15  
STATE p16 q16 p17 q17 p18 q18
"
DER dr1 dr2 dr3 dr4 
DER dv1 dv2 dv3 dv4 dv5 dv6 dv7 
DER dp1 dq1 dp2 dq2 dp3 dq3 dp4 dq4 dp5 dq5 dp6 dq6 dp7 dq7 dp8 dq8 dp9 dq9 
DER dp10 dq10 dp11 dq11 dp12 dq12 dp13 dq13 dp14 dq14 dp15 dq15
DER dp16 dq16 dp17 dq17 dp18 dq18
"
TIME t   " t is in ms. 
"
"  ANALOG DRIVES FOR 4 VNC PIs:  (Theta can range 0 to 360 deg.)
"
frcpi =(maxF1/2)*(1 + cos(2*(theta - 45)/R)) "Choose theta in deg. rcpi is +
rcpi1 = IF THETA > 135 THEN 1 ELSE 0	    " Supresses rcpi for theta	
rcpi2 = IF THETA < 315 THEN 1 ELSE 0	    " between 135 - 315o.
rcpi3 = rcpi1*rcpi2
rcpi = IF rcpi3 > 0 THEN 0 ELSE frcpi
"	
fripi = (maxF2/2)*(1 + cos(2*(theta - 135)/R))  " cos arg must be in rads.
ripi1 = IF theta > 225 THEN 1 ELSE 0
ripi2 = IF theta < 45 THEN 1 ELSE 0
ripi3 = ripi1 + ripi2
ripi = IF ripi3 > 0 THEN 0 ELSE fripi
"
flipi = (maxF2/2)*(1 + cos(2*(theta - 225)/R))
lipi1 = IF theta > 315 THEN 1 ELSE 0
lipi2 = IF theta < 135 THEN 1 ELSE 0
lipi3 = lipi1 + lipi2
lipi = IF lipi3 > 0 THEN 0 ELSE flipi
"
"
flcpi = (maxF1/2)*(1 + cos(2*(theta - 315)/R))
lcpi1 = IF theta > 45 THEN 1 ELSE 0
lcpi2 = IF theta < 225 THEN 1 ELSE 0
lcpi3 = lcpi1*lcpi2
lcpi = IF lcpi3 > 0 THEN 0 ELSE flcpi 
"
"  IPFM VFCs TO GENERATE INPUT SPIKES: Driven by analog f(theta)s above.
"
dr1 = rcpi - zi1 
wi1 = IF r1 > phi2 THEN 1 ELSE 0
si1 = DELAY(wi1, tau)
xi1 = wi1 - si1
b1 = IF xi1 > 0 THEN xi1 ELSE 0
zi1 = b1*phi2/tau 
uRCPI = zi1*Do	" Pulse train emulating VNC positional interneuron, RCPI.
"
dr2 = ripi - zi2
wi2 = IF r2 > phi2 THEN 1 ELSE 0
si2 = DELAY(wi2, tau)
xi2 = wi2 - si2
b2 = IF xi2 > 0 THEN xi2 ELSE 0
zi2 = b2*phi2/tau
uRIPI = zi2*Do	" Pulse train emulating VNC positional interneuron, RIPI.
"
dr3 = lipi - zi3
wi3 = IF r3 > phi2 THEN 1 ELSE 0
si3 = DELAY(wi3, tau)
xi3 = wi3 - si3
b3 = IF xi3 > 0 THEN xi3 ELSE 0
zi3 = b3*phi2/tau
uLIPI = zi3*Do	" Pulse train emulating VNC positional interneuron, LIPI.
"
dr4 = lcpi - zi4
wi4 = IF r4 > phi2 THEN 1 ELSE 0
si4 = DELAY(wi4, tau)
xi4 = wi4 - si4
b4 = IF xi4 > 0 THEN xi4 ELSE 0
zi4 = b4*phi2/tau
uLCPI = zi4*Do	" Pulse train emulating VNC positional interneuron, LCPI.
"
"  THE 3 OUTPUT RPFM NEURONS
"
dv1 = -co*v1 + co*e1 - z1   " RPFM Output Neuron # 1 @45o. z1 resets v1 to 0.
z1 = y1*phi/tau			" e1 is sum of e- & i- psps to N1.
w1 = IF v1 > phi THEN 1 ELSE 0	
s1 = DELAY(w1, tau)
x1 = w1 - s1
y1 = IF x1 > 0 THEN x1 ELSE 0  " + pulse generator when v1 > phi.
u1 = y1*Doo/tau
e1 = ge1*q1 - gi1*q2	       " Sum of epsps - ipsps for N1.
"
dv2 = -co*v2 + co*e2 - z2      " Output Neuron 2 @ 90o. 
z2 = y2*phi/tau
w2 = IF v2 > phi THEN 1 ELSE 0
s2 = DELAY(w2, tau)
x2 = w2 - s2
y2 = IF x2 > 0 THEN x2 ELSE 0
u2 = y2*Doo/tau
e2 = ge2*(q7 + q8) - gi2*(q6 + q9)
"
dv3 = -co*v3 + co*e3 - z3     " Output Neuron 3 @ 135o. 
z3 = y3*phi/tau
w3 = IF v3 > phi THEN 1 ELSE 0
s3 = DELAY(w3, tau)
x3 = w3 - s3
y3 = IF x3 > 0 THEN x3 ELSE 0
u3 = y3*Doo/tau
e3 = ge3*q14 - gi3*q13
"
"  4 INHIBITORY INTERNEURONS:
"
dv4 = -ci*v4 + ci*e4 -z4
z4 = y4*phi/tau
w4 = IF v4 > phi THEN 1 ELSE 0
s4 = DELAY(w4, tau)
x4 = w4 - s4
y4 = IF x4 > 0 THEN x4 ELSE 0
u4 = y4*Doi/tau
e4 = ge4*(q3 + q17 + q18)
"
dv5 = -ci*v5 + ci*e5 - z5
z5 = y5*phi/tau
w5 = IF v5 > phi THEN 1 ELSE 0
s5 = DELAY(w5, tau)
x5 = w5 - s5
y5 = IF x5 > 0 THEN x5 ELSE 0
u5 = y5*Doi/tau
e5 = ge5*(q4 + q5)
"
dv6 = -ci*v6 + ci*e6 - z6
z6 = y6*phi/tau
w6 = IF v6 > phi THEN 1 ELSE 0
s6 = DELAY(w6, tau)
x6 = w6 - s6
y6 = IF x6 > 0 THEN x6 ELSE 0
u6 = y6*Doi/tau
e6 = ge6*(q10 + q11)
"
dv7 = -ci*v7 + ci*e7 - z7
z7 = y7*phi/tau
w7 = IF v7 > phi THEN 1 ELSE 0
s7 = DELAY(w7, tau)
x7 = w7 - s7
y7 = IF x7 > 0 THEN x7 ELSE 0
u7 = y7*Doi/tau
e7 = ge7*(q12 + q15 + q16)
"
"  SYNAPTIC BALLISIC FILTER ODEs. (There are 14 synapses)
"
c1 = uRCPI
dp1 = -ae*p1 + c1		" c1 is spike input to BF.		
dq1 = -ae*q1 + p1		" q1 is BF analog output, > 0.
"
c2 = u4
dp2 = -ai*p2 + c2		
dq2 = -ai*q2 + p2
"
c3 = u2
dp3 = -ae*p3 + c3
dq3 = -ae*q3 + p3
"
c4 = u1
dp4 = -ae*p4 + c4
dq4 = -ae*q4 + p4
"
c5 = uLCPI
dp5 = -ae*p5 + c5
dq5 = -ae*q5 + p5
"
c6 = u5
dp6 = -ai*p6 + c6
dq6 = -ai*q6 + p6
"
c7 = uRIPI
dp7 = -ae*p7 + c7
dq7 = -ae*q7 + p7
"
c8 = uRCPI
dp8 = -ae*p8 + c8
dq8 = -ae*q8 + p8
"
c9 = u6
dp9 = -ai*p9 + c9
dq9 = -ai*q9 + p9
"
c10 = uLIPI
dp10 = -ae*p10 + c10
dq10 = -ae*q10 + p10
"
c11 = u3
dp11 = -ae*p11 + c11
dq11 = -ae*q11 + p11
"
c12 = u2
dp12 = -ae*p12 + c12
dq12 = -ae*q12 + p12
"
c13 = u7
dp13 = -ai*p13 + c13
dq13 = -ai*q13 + p13
"
c14 = uRIPI
dp14 = -ae*p14 + c14
dq14 = -ae*p14 + p14
"
c15 = uRCPI
dp15 = -ae*p15 + c15
dq15 = -ae*q15 + p15
"
c16 = uLIPI
dp16 = -ae*p16 + c16
dq16 = -ae*q16 + p16
"
c17 = uRIPI
dp17 = -ae*p17 + c17
dq17 = -ae*q17 + p17
"
c18 = uLCPI
dp18 = -ae*p18 + c18
dq18 = -ae*q18 + p18
"
"  OFFSET SCALED SPIKE OUTPUTS FOR PLOTTING:
"
bo1 = b1/5 + .1   " VNC INPUTS set by theta.
bo2 = b2/5 + .4
bo3 = b3/5 + .7
bo4 = b4/5 + 1.0
"
yo1 = y1/4 + 1.3  " Positional unit outputs: #1 @ 45o, #2 @ 90o, #3 @ 135o.
yo2 = y2/4 + 1.6
yo3 = y3/4 + 1.9
"
yo4 = y4/4 + 2.2  " Inhib interneuron outputs
yo5 = y5/4 + 2.5
yo6 = y6/4 + 2.8
yo7 = y7/4 + 3.1
angle = theta/100
"
theta = t  " Plot over 0 to 225 sec.= degrees.
"
"  CONSTANTS:
pi:3.14159
ge1:1      "N1 epsp input weight.
gi1:3      "N1 ipsp input weight.
ge2:1
gi2:3
ge3:1
gi3:3
ge4:1
ge5:1
ge6:.3      "N6 epsp input weight.
ge7:1
ae:.5
ai:.5
co:1
ci:1
phi:2.5 
phi2:3
tau:0.01 "(sec)
Do:1.
Doo:.5
Doi:2
maxF1:1
maxF2:1
R:57.296
" 
END


continuous system LinSynBF  " v. 3/10/99  A Linear ballistic filter
STATE x1 x2				" From Sec. 4.1.1.
DER dx1 dx2
TIME t
"
dx1 = -a*x1 + Vs
dx2 = -b*x2 + x1*a*b		" ab factor causes all h(t) to have same 
"					  area under curve.
"
Vs1 = IF t > dt THEN 0 ELSE Vso	" Makes three input impulses.
Vs2 = IF t > t2 THEN Vso ELSE 0
Vs2n = IF t > (t2 + dt) THEN -Vso ELSE 0
Vs3 = IF t > t3 THEN Vso ELSE 0
Vs3n = IF t > (t3 + dt) THEN -Vso ELSE 0
Vs = Vs1 + Vs2 + Vs2n + Vs3 + Vs3n
"
" CONSTANTS: 
a:1
b:1
Vso:100
dt:.01
t2:3
t3:6
"
END


Continuous System nicochan  " 3/11/99  System to emulate nicotinic ACh
STATE x y                   " channel dynamics. EPSP V proportional to x.
DER dx dy                   " Note 2 ACh molecules required/channel to
TIME t  " t in ms.            open it and cause epsp.
"                             y = ACh conc in synaptic cleft.
"					From Sec. 4.1.2.
"
dy = alphad - k4*y + k2*2*x - k1*(y^2)*(N - x) 
dx = k1*(y^2)*(N -x) - (k2 + k3)*x  " x = Density of open channels on SSM.
"						  vm epsp is proportional to x.
alpha = IF t < to THEN ao ELSE 0    " Pulse input rate of ACh.
alphad = DELAY(alpha/to, D)
"
" CONSTANTS:
D:.3
k1:10
k2:1
k3:1
k4:1
ao:1
N:1	" Max fraction of ion channels. (N - x)=fraction of unbound channels.
to:.1
"
END


continuous system FACILBF   " 3/10/99  Use Euler integration w/ DT = .001.
STATE x1 x2 f g		    " From Sec. 4.1.3.
DER dx1 dx2 df dg
TIME t
"
dx1 = -a*x1 + Vs		" Primary BF
dx2 = -b*x2 + a*b*x1
"
Vsd = DELAY(Vs, DEL)
df = -c*f + Vsd		" Facilitating BF
dg = -d*g + c*d*f 
ff = (1 + g*k)		" Facilitation factor
xo = x2*ff			" Facilitated epsp
"
Vs1 = IF t > dt THEN 0 ELSE Vso	" 3 input pulses.
Vs2 = IF t > t2 THEN Vso ELSE 0
Vs2n = IF t > (t2 + dt) THEN -Vso ELSE 0
Vs3 = IF t > t3 THEN Vso ELSE 0
Vs3n = IF t > (t3 + dt) THEN -Vso ELSE 0
Vs = Vs1 + Vs2 + Vs2n + Vs3 + Vs3n 
"
a:1
b:1
c:.3
d:.3
k:7	" Facilitation gain
DEL:.2
dt:.001
t2:3
t3:6
Vso:1000
"
END


continuous system ANTIFAC   " 3/10/99  Use Euler integration w/ DT = .001.
STATE x1 x2 f g		    " Antifacilitating psp.
DER dx1 dx2 df dg		    " From Sec. 4.1.4.
TIME t
"
dx1 = -a*x1 + Vs		" Primary BF
dx2 = -b*x2 + a*b*x1
"
df = -c*f + Vs		" Anti-Facilitating BF
dg = -d*g + c*d*f 
ff = exp(-k*g)		" Anti-Facilitation factor
xo = x2*ff			" Anti-Facilitated psp
"
Vs1 = IF t > dt THEN 0 ELSE Vso
Vs2 = IF t > t2 THEN Vso ELSE 0
Vs2n = IF t > (t2 + dt) THEN -Vso ELSE 0
Vs3 = IF t > t3 THEN Vso ELSE 0
Vs3n = IF t > (t3 + dt) THEN -Vso ELSE 0
Vs = Vs1 + Vs2 + Vs2n + Vs3 + Vs3n 
"
"  PARAMETERS
a:1  	" radians/ms
b:1
c:.3
d:.3
k:1	" Antifacilitating factor.
dt:.001
t2:3
t3:6
Vso:1000
"
END


continuous system ClInh1   " 3/12/99 Chloride inhibitory synapse model.
STATE x1 x2 x3 x4 v1 f g	" From Sec. 4.1.5.
DER dx1 dx2 dx3 dx4 dv1 df dg
TIME t  " ms.
"
dx1 = -a*x1 + Vs1		   " BF for axodendritic synapse 1.
dx2 = -b*x2 + a*b*x1
Ve1 = DELAY(x2, D1)
"
dx3 = -a*x3 + Vs2		   " BF for axodendritic synapse 2.
dx4 = -b*x4 + a*b*x3
Ve2 = DELAY(x4, D2)
"
df = -c*f + Vi		   " BF for axosomatic inhibitory synapse
dg = -d*g + c*d*f
" 
Hi = 1/(1 + g*Ki)          " Hill function for Cl- inhibitory synapse.
Vex = (Ve1*k1 + Ve2*k2)*Hi " Sum 2 epsps; inhib with Cl- synapse.
"
dv1 = -c1*v1 + c1*Vex - z1	"RPFM SGL
w1 = IF v1 > phi THEN 1 ELSE 0
s1 = DELAY(w1, tau)
q1 = w1 - s1
y1 = IF q1 > 0 THEN q1 ELSE 0
z1 = y1*phi/tau
"
"  INPUTS:
Vs1 = IF t < dt THEN Vso ELSE 0
Vs2 = IF t < dt THEN Vso ELSE 0
vi1 = IF t > ti THEN Vso ELSE 0
vi2 = IF t > (ti + dt) THEN -Vso ELSE 0
Vi = vi1 + vi2

" PARAMETERS:
tau:.001  " ms.
ti:1
dt:.001
Vso:1000
phi:.1
c1:.3
a:1
b:1
c:.5
d:.5
D1:0.3
D2:0.1
k1:.5
k2:.9
Ki:5
"
END


To simulate an IPFM SGL using Simnon(, we can use the following subroutine:
" (From Sec. 4.3.1)
"
dv  =  e - z    " Integrator with Ki = 1. v is a state, dv is its derivative.
w  =  IF v > phi THEN 1 ELSE 0 "(   
s  =  DELAY(w, tau)            "( (  Pulse generator.
x  =  w - s                    "(
y  =  IF x > 0 THEN x ELSE 0   " Pulses are > 0.
z  =  y*phi/tau                " Pulse resets integrator from phi to 0. 
u  =  y*Do/tau	             " Pulse train output, pulse areas = Do.
"


A Simnon model for an RPFM SGL is given below: (Note that only the first line of code differs from the IPFM model.) From Sec. 4.3.2.
"   
dv = -c*v + c*e - z		" Analog LPF with inputs e and z.
w = IF v > phir THEN 1 ELSE 0
s = DELAY(w, tau)
x = w - s
y = IF x > 0 THEN x ELSE 0	" y(t) is triangular: height =1, base = 2tau.
z = y*phir/tau			" The pulse z resets v to 0.
u = y*Do/tau			" RPFM SGL output pulses
"


A SGL model showing neural fatigue: From Sec. 4.3.4.
"
dv = -c*v + c*e - z		" Analog LPF with inputs e and z.
w = IF v > phi THEN 1 ELSE 0
s = DELAY(w, tau)
x = w - s
y = IF x > 0 THEN x ELSE 0	" y(t) is triangular: height =1, base = 2tau.
z = y*phi/tau			" The pulse z resets v to 0.
u = y*Do/tau			" RPFM SGL output pulses
dp = -a*p + k*u			" Feedback LPF
dq = -b*q + p
phi = phi0 + q			" Phi is raised prop.to firing rate.
"


continuous system TNeuron  " 4/22/99. Use EULER integration with Dt = .001.
"					From Sec. 4.4.1.
STATE x1 x2 x3 x4 v	   
DER dx1 dx2 dx3 dx4 dv 
TIME t  " msec. 
"
dx1 = -a*x1 + Vs1      " 2-pole, synaptic LPF
dx2 = -b*x2 + a*b*x1   " x2 is output epsp
"
dx3 = -a*x3 + Vs2	     " 2-pole, synaptic LPF
dx4 = -b*x4 + a*b*x3
"
Vex = x2 + x4	     " Two epsps summed to form generator potential.
"
"  RPFM T-NEURON
dv = -c*v + c*Vex - z 
w = IF v > phir THEN 1 ELSE 0
d = DELAY(w, tau)
q = w - d
y = IF q > 0 THEN q ELSE 0
z = y*phir/tau
"
"  MAKE INPUT PULSES:
Vs1 = IF t < dt THEN Vso/dt ELSE 0
Vs21 = IF t > t2 THEN Vso/dt ELSE 0
Vs22 = IF t > (t2 + dt) THEN -Vso/dt ELSE 0
Vs2 = Vs21 + Vs22
"
RPFMout = y/8 + .5
"
"  PARAMETERS:
tau:.001
dt:.001
Vso:1.
phir:.2
a:1	" Units radians/ms.
b:1	"   "		"
c:1	"   "		"
t2:1
"
END



Continuous system BDsupr1 	" V. 3/04/99. From Sec. 4.4.2
"  Use EULER integration with delT = tau. There are 9 states.		
"	This system is Reiss' BD with harmonic suppression.
"
STATE noise v1 v3 p1 q1 p2 q2 p3 q3 
DER dnoise dv1 dv3 dp1 dq1 dp2 dq2 dp3 dq3 
TIME t   " t is in ms. 
"
dnoise = -wo*noise + SD*NORM(t) 	" BW limiting ODE for noise.
Vin = noise + Kr*t + Vo			" Noise + ramp + const. drive for N1
"							  IPFM VFC.
" THE IPFM SPIKE SOURCE: 
"
dv1 = Vin - z1                     	" IPFM neuron N1: Integrator Ki = 1. 
w1 = IF v1 > phi1 THEN 1 ELSE 0	" N1 behaves as an ideal VFC.      
s1 = DELAY(w1, tau) 			" Pulse generator.
x1 = w1 - s1                    
y1 = IF x1 > 0 THEN x1 ELSE 0     	" Pulses are > 0.
z1 = y1*phi1/tau                  	" Pulse resets integrator. 
u1 = y1*Do1/tau	             	" Pulse train output, pulse areas = Do1.
"
y2 = DELAY(y1, D2) 		" Excit. delay D2 ms. to N3 synapse. 
u2 = Do2*y2/tau			" Thru excit. 1:1 interneuron (xparent).		
"
yi = DELAY(y1, D1)		" Inhibitory delay (thru xparent 1:1 neuron). 
ui = Do3*yi/tau
"
"  SYNAPTIC BALLISIC FILTER ODEs. (There are 3 synapses on N3)
dp1 = -a1*p1 + a1*u1
dq1 = -a1*q1 + a1*p1
"
dp2 = -a2*p2 + a2*u2		
dq2 = -a2*q2 + a2*p2
"
dp3 = -a3*p3 + a3*ui    " BF for ipsp.  Determines effective delay, I.
dq3 = -a3*q3 + a3*p3
"
" OUTPUT RPFM T-NEURON, N3:
dv3 = -c3*v3 + c3*e3 - z3   	
z3 = y3*phi3/tau
w3 = IF v3 > phi3 THEN 1 ELSE 0
s3 = DELAY(w3, tau)
x3 = w3 - s3
y3 = IF x3 > 0 THEN x3 ELSE 0
e3 = g1*q1 + g2*q2 - g3*q3
"
"  OFFSET SCALED SPIKE OUTPUTS FOR PLOTTING:
b0 = y1/4 + .1
b1 = y2/4 + .4
b2 = yi/4 + .7
b3 = y3/4 + 1.0
"
"  CONSTANTS:
a1:1 " r/ms.
a2:1
a3:.55
g1:1
g2:1      
g3:.5
c3:7
phi1:3 
phi3:.22
tau:0.001 "(millisec)
Do1:1
Do2:1
Do3:1.2
wo:1.5
SD:0
Kr:.0125
Vo:0
D1:0.333
D2:2
" 
END



Continuous system RECIPIN3	" V. 3/16/99. From Sec.4.5.1.
"  The Basic 4 Neuron Reciprocal Inhibition model. 
" 
STATE v1 v2 v3 v4 p1 p2 p3 p4 p5 p6 va 
DER dv1 dv2 dv3 dv4 dp1 dp2 dp3 dp4 dp5 dp6 dva 
TIME t	" t is in ms.
"
"  IPFM VFC TO GENERATE INPUT SPIKES:
dva = Ea - za				" Ea is analog input to VFC. 
wa = IF va > phia THEN 1 ELSE 0
sa = DELAY(wa, tau)
xa = wa - sa
ya = IF xa > 0 THEN xa ELSE 0
za = ya*phia/tau
ua = ya*Doa/tau
"
"  THE RPFM RI PAIR:
dv1 = -c1*v1 + c1*E1 - z1            " Output Neuron 1. RPFM model.
w1 = IF v1 > phi1 THEN 1 ELSE 0
s1 = DELAY(w1, tau)
x1 = w1 - s1
y1 = IF x1 > 0 THEN x1 ELSE 0
z1 = y1*phi1/tau
u1 = y1*Do1/tau
"
dv2 = -c2*v2 + c2*E2 - z2            " Output Neuron 2. RPFM model.
w2 = IF v2 > phi2 THEN 1 ELSE 0
s2 = DELAY(w2, tau)
x2 = w2 - s2
y2 = IF x2 > 0 THEN x2 ELSE 0
z2 = y2*phi2/tau
u2 = y2*Do2/tau
"
" 2 RPFM INHIBITORY INTERNEURONS:
dv3 = -c3*v3 + c3*e3 - z3             " Ihibitory interneuron 3.
w3 = IF v3 > phi3 THEN 1 ELSE 0
s3 = DELAY(w3, tau)
x3 = w3 - s3
y3 = IF x3 > 0 THEN x3 ELSE 0
z3 = y3*phi3/tau
u3 = y3*Do3/tau
"
dv4 = -c4*v4 + c4*e4 - z4             " Inhibitory interneuron 4.
w4 = IF v4 > phi4 THEN 1 ELSE 0
s4 = DELAY(w4, tau)
x4 = w4 - s4
y4 = IF x4 > 0 THEN x4 ELSE 0
z4 = y4*phi4/tau
u4 = y4*Do4/tau
"
" 1 TIME-CONSTANT, SYNAPTIC BFs:
" 
dp1 = -a1*p1 + ua		" Synapse 1 ballistic filter (ipsp)
"
dp2 = -a2*p2 + ua		" Synapse 2 ballistic filter. q2 is output.
"
dp3 = -a3*p3 + u1		 
"
dp4 = -a4*p4 + u2
"
dp5 = -a5*p5 + u3
"
dp6 = -a6*p6 + u4
"
" Inputs to RPFM neurons
E1 = p1 - p6      " epsp - ipsp input to output neuron 1.			
E2 = p2 - p5      " epsp - ipsp input to output neuron 2.
E3 = p3 		" epsp input to inhibitory interneuron 3
E4 = p4      	" epsp input to inhibitory interneuron 4.
"
" Offset outputs for plotting.
oa = ya/5 + 1
o1 = y1/5 + 1.3	
o2 = y2/5 + 1.6			
o3 = y3/5 + 1.9
o4 = y4/5 + 2.2
"
"  INPUTS TO IPFM VFCs:
Ea = A	" Nonzero input to IPFM VFCa.
"
" PARAMETERS: 
A:2.
tau:0.001	" ms
a1:1		" All natural frequencies are in r/ms.
a2:1
a3:.5
a4:.5
a5:.2
a6:.2
c1:1.
c2:1
c3:1
c4:1
phia:1
phi1:1
phi2:1
phi3:1
phi4:1
Doa:1.
Do1:0.75
Do2:0.75
Do3:1
Do4:1
"
END


Szentagothai's RI circuit.  See Sec. 4.5.2.

Continuous system SZGRIsys " V. 3/14/99. From Appendix 2.
"  4 Neuron RI model after Szentagothai in B&H, Fig. 5.16. 14 states.
" 
STATE v1 v2 v3 v4 p1 p2 p3 p4 p5 p6 p7 p8 va vb           
DER dv1 dv2 dv3 dv4 dp1 dp2 dp3 dp4 dp5 dp6 dp7 dp8 dva dvb
TIME t
"
"  IPFM VFCa TO GENERATE INPUT SPIKES:
dva = Ea - za				 
wa = IF va > phia THEN 1 ELSE 0
sa = DELAY(wa, tau)
xa = wa - sa
ya = IF xa > 0 THEN xa ELSE 0
za = y*phia/tau
ua = z*Doa
"
" IPFM VFCb TO GENERATE INPUT SPIKES
dvb = Eb - zb 
wb = IF vb > phib THEN 1 ELSE 0
sb = DELAY(wb, tau)
xb = wb - sb
yb = IF xb > 0 THEN xb ELSE 0
zb = y*phib/tau
ub = z*Dob

"  THE RPFM RI SYSTEM:
dv1 = -c1*v1 + c1*E1 - z1            " Output Neuron 1. 
w1 = IF v1 > phi1 THEN 1 ELSE 0
s1 = DELAY(w1, tau)
x1 = w1 - s1
y1 = IF x1 > 0 THEN x1 ELSE 0
z1 = y1*phi1/tau
"
dv2 = -c2*v2 + c2*E2 - z2            " Output Neuron 2. 
w2 = IF v2 > phi2 THEN 1 ELSE 0
s2 = DELAY(w2, tau)
x2 = w2 - s2
y2 = IF x2 > 0 THEN x2 ELSE 0
z2 = y2*phi2/tau
"
dv3 = -c3*v3 + c3*E3 - z3             " Ihibitory interneuron 3.
w3 = IF v3 > phi3 THEN 1 ELSE 0
s3 = DELAY(w3, tau)
x3 = w3 - s3
y3 = IF x3 > 0 THEN x3 ELSE 0
z3 = y3*phi3/tau
u3 = y3*Do3/tau
"
dv4 = -c4*v4 + c4*E4 - z4             " Inhibitory interneuron 4.
w4 = IF v4 > phi4 THEN 1 ELSE 0
s4 = DELAY(w4, tau)
x4 = w4 - s4
y4 = IF x4 > 0 THEN x4 ELSE 0
z4 = y4*phi4/tau
u4 = y4*Do4/tau
"
" EXCITATORY SYNAPSE BFs:
"
dp1 = -a1*p1 + ua		" Synapse 1 ballistic filter (ipsp)
"
dp2 = -a2*p2 + ub		" Synapse 2 ballistic filter. q2 is output.
"
dp3 = -a3*p3 + ua
"
dp4 = -a4*p4 + ub
"
" INHIBITORY SYNAPSE BFs:
"
dp5 = -a5*p5 + u4
"
dp6 = -a6*p6 + u3
"
dp7 = -a7*p7 + u3
"
dp8 = -a8*p8 + u4
"
" Inputs to RPFM neurons
E1 = p1 - p5          " epsp - ipsp input to output neuron 1.			
E2 = p2 - p6          " epsp - ipsp input to output neuron 2.
E3 = p3 - p8          " epsp - ipsp input to inhibitory interneuron 3
E4 = p4 - p7	    " epsp - ipsp input to inhibitory interneuron 4.
"
" Offset outputs for plotting.
oa = ya/5 + 1
ob = yb/5 + 1.3
o1 = y1/5 + 1.6	
o2 = y2/5 + 1.9			
o3 = y3/5 + 2.2
o4 = y4/5 + 2.5
"
" ANALOG INPUTS TO IPFM VFCs:
Ea = A*(1 - cos(wo*t))	" Nonzero input to IPFM VFCa.
Eb = B*(1 + cos(wo*t))  "  Nonzero input to IPFM VFCb.
wo = 6.28*fo
"
" PARAMETERS
A:2.
B:2 
fo:. 05
tau:0.001
a1:1
a2:1
a3:1
a4:1
a5:1
a6:1
a7:1
a8:1
c1:1.
c2:1
c3:1
c4:1
phia:1
phib:1
phi1:1
phi2:1
phi3:1
phi4:1
Doa:1.
Dob:1		
Do3:1
Do4:1
"
END

Parameters used:

[SZGRIsys]  " Data file RI3.t  For burst sharpening by RI. 3/16/99.
v1:0.
v2:0.
v3:0.
v4:0.
p1:0.
p2:0.
p3:0.
p4:0.
p5:0.
p6:0.
p7:0.
p8:0.
va:0.
vb:0.
phia:1.
tau:1.E-3
Doa:1.
phib:1.
Dob:1.
c1:1.
phi1:1.
c2:1.
phi2:1.
c3:1.
phi3:1.
Do3:1.
c4:1.
phi4:1.
Do4:1.
a1:1.
a2:1.
a3:1.
a4:1.
a5:1.
a6:1.
a7:0.3
a8:0.3
A:0.8
B:0.8
fo:0.05



Continuous system BURSTNM1 " v. 3/16/99. From Sec. 4.5.3.
"  4 Neuron RI model by RBN. 
" 
STATE v1 v2 v3 v4 p1 p2 p3 p4 p5 p6 va 
DER dv1 dv2 dv3 dv4 dp1 dp2 dp3 dp4 dp5 dp6 dva 
TIME t
"
"  IPFM VFCa TO GENERATE INPUT SPIKES:
dva = Ea - za 
wa = IF va > phia THEN 1 ELSE 0
sa = DELAY(wa, tau)
xa = wa - sa
ya = IF xa > 0 THEN xa ELSE 0
za = ya*phia/tau
ua = ya*Doa/tau
"
"  THE RPFM RI PAIR:
dv1 = -c1*v1 + c1*E1 - z1            " Output Neuron 1. RPFM model.
w1 = IF v1 > phi1 THEN 1 ELSE 0
s1 = DELAY(w1, tau)
x1 = w1 - s1
y1 = IF x1 > 0 THEN x1 ELSE 0
z1 = y1*phi1/tau
u1 = y1*Do1/tau
"
dv2 = -c2*v2 + c2*E2 - z2            " Output Neuron 2. RPFM model
w2 = IF v2 > phi2 THEN 1 ELSE 0
s2 = DELAY(w2, tau)
x2 = w2 - s2
y2 = IF x2 > 0 THEN x2 ELSE 0
z2 = y2*phi2/tau
u2 = y2*Do2/tau
"
"  THE INHIBITORY INTERNEURONS:
dv3 = -c3*v3 + c3*e3 - z3             " Ihibitory interneuron 3.
w3 = IF v3 > phi3 THEN 1 ELSE 0
s3 = DELAY(w3, tau)
x3 = w3 - s3
y3 = IF x3 > 0 THEN x3 ELSE 0
z3 = y3*phi3/tau
u3 = y3*Do3/tau
u3d = DELAY(u3, D)
"
dv4 = -c4*v4 + c4*e4 - z4             " Inhibitory interneuron 4.
w4 = IF v4 > phi4 THEN 1 ELSE 0
s4 = DELAY(w4, tau)
x4 = w4 - s4
y4 = IF x4 > 0 THEN x4 ELSE 0
z4 = y4*phi4/tau
u4 = y4*Do4/tau
"
" 1 TC SYNAPTIC BFs: 
dp1 = -a1*p1 + ua		" Synapse 1 ballistic filter (ipsp)
"
dp2 = -a2*p2 + ua		" Synapse 2 ballistic filter. q2 is output.
"
dp3 = -a3*p3 + u3d	" Delayed FB around N1
"
dp4 = -a4*p4 + u4
"
dp5 = -a5*p5 + u1
"
dp6 = -a6*p6 + u1
"
" Inputs to RPFM neurons
E1 = p1 - p3          " epsp - ipsp input to output neuron 1.			
E2 = p2 - p4          " epsp - ipsp input to output neuron 2.
E3 = p6 		    " epsp input to inhibitory interneuron 3
E4 = p5      	    " epsp input to inhibitory interneuron 4.
"
" Offset outputs for plotting.
oa = ya/5 + 1
o1 = y1/5 + 1.3	
o2 = y2/5 + 1.6			
o3 = y3/5 + 1.9
o4 = y4/5 + 2.2
"
"
" PARAMETERS (Adjusted to get 180o phase between N1 and N1 bursts)
A:2.
fo:.025
tau:0.001
D:3	" msec.
a1:1
a2:1
a3:.5
a4:5
a5:5
a6:.5
c1:1.
c2:1
c3:1
c4:1
phia:1
phi1:0.5
phi2:0.7
phi3:2.1
phi4:0.2
Doa:1.
Do1:2
Do2:1
Do3:.5
Do4:3
"
END

The parameter file, RINH1.t:

[BURSTNM1]  " 3/16/99. Pars for BURSTNM1.t
v1:0.
v2:0.
v3:0.
v4:0.
p1:0.
p2:0.
p3:0.
p4:0.
p5:0.
p6:0.
va:0.
phia:1.
tau:1.E-3
Doa:1.
c1:1.
phi1:0.5
Do1:2.
c2:1.
phi2:0.7
Do2:1.
c3:1.
phi3:2.1
Do3:0.5
D:3.
c4:1.
phi4:0.2
Do4:3.
a1:1.
a2:1.
a3:0.5
a4:5.
a5:5.
a6:0.5
A:2.
fo:0.025
a7:1.
a8:1.



continuous system RINGOSC1  " 3/16/99 Uses PFB around 4 RPFM neurons
" Synaptic BFs have 2 TCs.	From Sec. 4.5.4.
"
STATE v1 v2 v3 v4 p1 p2 p3 p4 p5 q1 q2 q3 q4 q5 
DER dv1 dv2 dv3 dv4 dp1 dp2 dp3 dp4 dp5 dq1 dq2 dq3 dq4 dq5
TIME t
"
dv1 = -c1*v1 + c1*E1 - z1            " Neuron 1. RPFM model.
w1 = IF v1 > phi1 THEN 1 ELSE 0
s1 = DELAY(w1, tau)
x1 = w1 - s1
y1 = IF x1 > 0 THEN x1 ELSE 0
z1 = y1*phi1/tau
u1 = y1*Do1/tau
"
dv2 = -c2*v2 + c2*E2 - z2            " Neuron 2. RPFM model
w2 = IF v2 > phi2 THEN 1 ELSE 0
s2 = DELAY(w2, tau)
x2 = w2 - s2
y2 = IF x2 > 0 THEN x2 ELSE 0
z2 = y2*phi2/tau
u2 = y2*Do2/tau
"
dv3 = -c3*v3 + c3*E3 - z3             " N3.
w3 = IF v3 > phi3 THEN 1 ELSE 0
s3 = DELAY(w3, tau)
x3 = w3 - s3
y3 = IF x3 > 0 THEN x3 ELSE 0
z3 = y3*phi3/tau
u3 = y3*Do3/tau
"
dv4 = -c4*v4 + c4*E4 - z4             " N 4.
w4 = IF v4 > phi4 THEN 1 ELSE 0
s4 = DELAY(w4, tau)
x4 = w4 - s4
y4 = IF x4 > 0 THEN x4 ELSE 0
z4 = y4*phi4/tau
u4 = y4*Do4/tau
"
" SYNAPSES:
"
dp1 = -a1*p1 + u1
dq1 = -b1*q1 + p1*a1*b1
"
dp2 = -a2*p2 + u2
dq2 = -b2*q2 + p2*a2*b2
"
dp3 = -a3*p3 + u3
dq3 = -b3*q3 + p3*a3*b3
"
dp4 = -a4*p4 + u4
dq4 = -b4*q4 + p4*a4*b4
"
dp5 = -a5*p5 + uo
dq5 = -b5*q5 + p5*a5*b5
"
"  NEURON INPUTS:
"
E1 = q4 + q5
E2 = q1
E3 = q2
E4 = q3
uo = IF t < tau THEN Ao ELSE 0
"
"  PLOTTING VARIABLES:
o1 = y1/5 + .1
o2 = y2/5 + .4
o3 = y3/5 + .7
o4 = y4/5 + 1.0
"
"  CONSTANTS:
Ao:1000
tau:.001
c1:1
c2:1
c3:1
c4:1
Do1:1
Do2:1
Do3:1
Do4:1
phi1:1
phi2:1
phi3:1
phi4:1
a1:1
b1:1
a2:1
b2:1
a3:1
b3:1
a4:1
b4:1
a5:1
b5:1
"
END



Continuous system RING4	" V. 3/21/99. From Sec. 4.5.4
"  A Basic 4 Neuron RING model with NFB. 1 input to N1. 
"  Uses IPFM neurons and single TC synaptic BFs. System unstable for
"  hi gain; does burst oscs, but not 2-phase. 
"
STATE v1 v2 v3 v4 p1 p3 p4 p5 p6 va 
DER dv1 dv2 dv3 dv4 dp1 dp3 dp4 dp5 dp6 dva 
TIME t
"
"  IPFM VFCa TO GENERATE INPUT SPIKES:
dva = Ea - za 
wa = IF va > phia THEN 1 ELSE 0
sa = DELAY(wa, tau)
xa = wa - sa
ya = IF xa > 0 THEN xa ELSE 0
za = ya*phia/tau
ua = ya*Doa/tau
"
"  THE IPFM SGLs:
"
dv1 = E1 - z1           		" Output Neuron 1. 
w1 = IF v1 > phi1 THEN 1 ELSE 0
s1 = DELAY(w1, tau)
x1 = w1 - s1
y1 = IF x1 > 0 THEN x1 ELSE 0
z1 = y1*phi1/tau
u1 = y1*Do1/tau
"
dv2 = E2 - z2           		" Interneuron 2. 
w2 = IF v2 > phi2 THEN 1 ELSE 0
s2 = DELAY(w2, tau)
x2 = w2 - s2
y2 = IF x2 > 0 THEN x2 ELSE 0
z2 = y2*phi2/tau
u2 = y2*Do2/tau
"
dv3 = E3 - z3             		" Interneuron 3.
w3 = IF v3 > phi3 THEN 1 ELSE 0
s3 = DELAY(w3, tau)
x3 = w3 - s3
y3 = IF x3 > 0 THEN x3 ELSE 0
z3 = y3*phi3/tau
u3 = y3*Do3/tau
"		
dv4 = E4 - z4            		 " Inhibitory interneuron 4.
w4 = IF v4 > phi4 THEN 1 ELSE 0
s4 = DELAY(w4, tau)
x4 = w4 - s4
y4 = IF x4 > 0 THEN x4 ELSE 0
z4 = y4*phi4/tau
u4 = y4*Do4/tau
"
" 1 TC SYNAPTIC BFs: 
dp1 = -a1*p1 + a1*ua	" Synapse 1 BF (epsp). dc gain = 1
"
dp3 = -a3*p3 + a3*u1		
"
dp4 = -a4*p4 + A4*u2
"
dp5 = -a5*p5 + a5*u3
"
dp6 = -a6*p6 + a6*u4	" Synapse 6 BF. Generates ipsp.
"
" Inputs to IPFM neurons
E1 = p1 - p6      " epsp - ipsp input to neuron 1.			
E2 = p5      	" epsp input to neuron 2.
E3 = p3 	    	" epsp input to interneuron 3
E4 = p4      	" epsp input to inhibitory interneuron 4.
"
" Offset outputs for plotting.
oa = ya/5 + 1
o1 = y1/5 + 1.3	
o2 = y2/5 + 1.6			
o3 = y3/5 + 1.9
o4 = y4/5 + 2.2
"
"  INPUT TO IPFM VFC:
Ea = A "*(1 - cos(wo*t))	" Nonzero input to IPFM VFCa.
wo = 6.28*fo
"
" PARAMETERS: 
A:0.75
fo:.025
tau:0.001
a1:1.
a3:.2
a4:.2
a5:.2
a6:.2
phia:1
phi1:1
phi2:1
phi3:1
phi4:1
Doa:1.
Do1:1
Do2:2
Do3:2
Do4:1
"
END


Program used in home problem 1.7:

CONTINUOUS SYSTEM epsp	" 02/05/2000 Use Euler integration with delT = tau.  
STATE v psNT bNT Nbar Ns NT	" Linear kinetic model for epsp generation.
DER dv dpsNT dbNT dNbar dNs dNT
TIME t	" t in milliseconds.
"
" IPFM SOURCE OF PERIODIC INPUT PULSES:
dv = Vin - z	
w = IF v > phi THEN 1 ELSE 0
s = DELAY(w, tau)
x = w - s
y = IF x > 0 THEN x ELSE 0
yo = y/10 - .15	" yo for plotting
z = y*phi/tau
Pin = y*Do/tau		" Input pulses.
"				  Ko is fraction of NT released at every pulse.
dpsNT = -Kb*psNT + DELAY(Ko*NT*Pin, dT)  " Conc. NT @ SSM. Kb is mono-
"				molecular rate const for NT binding to receptors.
dbNT = Kb*psNT - Kh*bNT		" Relative density of NT bound to (XS)
"			  	receptors. Kh is hydrolysis rate constant.
dNbar = Kh*bNT - Kd*Nbar	" Conc. inactive NT in cleft.
"				  Kd is diffusion rate constant.
dNs = Kd*Nbar - Ku*Ns		" Relative density of Nbar bound to bouton 
"		        	membrane. Ku is rate of active Nbar uptake by bouton.
dNT = c + Ku*Ns - Ko*NT*Pin	" Conc. stored neurotransmitter in bouton.
"					Note that NT is recycled with rate Ku.
Vm = Kv*bNT				" epsp is proportional to bNT.
"
Vin = Vo + Kr*t			" Causes the freq. of y to ramp up.
"
" CONSTANTS:
phi:1
tau:.001
Vo:.10
Do:8
Ko:.05
Kb:1.
Kh:2.
Kd:0.5
Kv:10
Ku:5
Kr:.00333
c:0
dT:0.2 " ms.
zero:0
"
" IC:
NT:.5
" 
END



Program used in home problem 4.9:

Continuous system HP49      " v. 3/13/99. rev. 01/25/00
" Use EULER integration with delT = tau. There are 7 states.
" This system is Reiss' Band Suppressor.
STATE noise v1 v3 p1 p2 q1 q2  
DER dnoise dv1 dv3 dp1 dp2 dq1 dq2   
TIME t   " t is in ms. 
"
dnoise = -wo*noise + SD*NORM(t) " BW limiting ODE for noise.
Vin = noise + Kr*t + Vo		  " Analog noise drive for N1
"
"  THE IPFM INPUT VFC, N1: 
dv1 = Vin - z1                 	 
w1 = IF v1 > phi1 THEN 1 ELSE 0	      
s1 = DELAY(w1, tau) 		 " (Pulse generator.
x1 = w1 - s1                    
y1 = IF x1 > 0 THEN x1 ELSE 0  " Pulses are > 0.
z1 = y1*phi1/tau       		 " Pulse resets integrator. 
u1 = y1*Do1/tau	             " Pulse train output, pulse areas = Do1.
"
y2 = DELAY(y1, D) 		 " Inhib. delay D ms to N3 synapse. 
u2 = Do2*y2/tau			 " Thru excit. 1:1 interneuron (xparent).
"
"  2-POLE, SYNAPTIC BALLISTIC FILTERs: (There are 2 synapses on N3)
dq1 = -a1*q1 + a1*u1
dp1 = -a1*p1 + a1*q1		 " Direct excitatory synapse BF.
"
dq2 = -a2*q2 + a2*u2
dp2 = -a2*p2 + a2*q2		 " Inhibitory synapse BF.		
"
e3 = p1 - p2			 " Generator potential for N3:
"
" RPFM T-NEURON, N3:
dv3 = -c3*v3 + c3*e3 - z3 	
w3 = IF v3 > phi3 THEN 1 ELSE 0
s3 = DELAY(w3, tau)
x3 = w3 - s3
y3 = IF x3 > 0 THEN x3 ELSE 0
z3 = y3*phi3/tau
"
"  OFFSET SCALED SPIKE OUTPUTS FOR PLOTTING:
yo1 = y1/5 - .75
yo2 = y2/5 - .5
yo3 = y3/5 - .25
"
"  CONSTANTS:
a1:0.9
a2:3
c3:3.33
phi1:1 
phi3:0.30
tau:0.001 	"(millisec)
Do1:1
Do2:1.1
wo:1.5	" Cutoff freq. of noise LPF in r/ms.
SD:0		" Standard deviation of broad-band Gaussian noise.
Kr:6.5E-3
Vo:.20
D:2.5 	" ms.
" 
END







Program used in home prob. 4.11:

CONTINUOUS SYSTEM HP411	"02/01/00. Use Euler integration with delT =.001.
STATE v phi
DER dv dphi
TIME t			" t in ms.
"
" RPFM SGL WITH EXPONENTIAL REFRACTORY PERIOD:
dv = -c*v + c*Vin - z
w = IF v > phir THEN 1 ELSE 0
s = DELAY(w, dT)
x = w - s
y = IF x > 0 THEN x ELSE 0	" RPFM output.
z = y*phir/dT
dphi = -a*phi + a*Do*z	" Refractory increase in firing threshold.
phir = phio + phi		" Firing threshold, mV.
Vin = Vo + Kr*t		" SGL generator potential, mV.
"
yo = y/1.25 - 1		" Offset y for plotting
"
" PARAMETERS:
phio:10 " mV.
c:0.2	  " r/ms.
a:0.5	  " r/ms.
Vo:7.5  " mV.
Kr:.5
dT:.001 " ms.
Do:2
"
END