# modified Steve Prescott, Aug 29, 2011
# For Coggan et al. (Explaining pathological changes in axonal excitability through dynamical analysis of conductancebased models. J Neural Eng 2011; 8: 065002)
# Incorporates Na influx and simple decay. No pumps yet implemented in this code.
# Nai dynamically varies and is used to update Ena based on Nernst equation.
# DIFFERENTIAL EQUATIONS
dv/dt = (I1(t)+I2(t)+Idcgna*minf(V)*(VVna)gk*w*(VVK)gl*(VVl)gnap*z*(VVna))/cap
dw/dt = phi*(winf(V)w)/tauw(V)
dz/dt = phi_z*(zinf(V)z)/tauz(V)
dnai/dt = (SAvol*(gna*minf(v)*(VVna)+gnap*z*(VVna)))/F(nai17.5)/tau_na
# if you don't want Nai to vary, comment out above line and define Nai as a fixed parameter,
# Either that, or just define Vna as a parameter rather than using the Nernst equation to calculate it (see ### below)
param tau_na=100
# PARAMETERS AND FUNCTIONS
minf(v)=.5*(1+tanh((vv1)/v2))
winf(v)=.5*(1+tanh((vv3)/v4))
zinf(v)=.5*(1+tanh((vv5)/v6))
tauw(v)=1/cosh((vv3)/(2*v4))
tauz(v)=1/cosh((vv5)/(2*v6))
# v5 and v6 correspond to beta_z and gamma_z  checked, these parameters match PNAS paper
Vna=25*ln(nao/nai)
### see note above.
# param Vna=50
param vk=100,vl=70
param gk=20,gl=2,gna=30
param v1=1.2,v2=18,v3=10,v4=10
param gnap=0.8
param v5=45,v6=10
param F=96485
param phi=.15,phi_z=0.05,cap=2
param nao=138
#nao and nai(0) adjusted to set Vna(0)=50mV, differences from 150 & 15 (resp) divided between them
param r=.005
# param h=.01
# !vol=(4/3)*pi*(r^3)
# !SA=4*pi*(r^2)
!SAvol=2/r
# SA:vol ratio for cylinder without ends
# note because of units, set r to 0.005 to get radius of 0.5 microns. The decimal place is not wrong.
# spherical soma r=7.5 um, cylindrical node r=.5 um, h=1 um
# for a cylinder: vol=pi*h*(r^2), SA=2*pi*r*(r+h), SA/vol=(2/r)+(2/h)
aux vna_=vna
aux vk_=vk
aux gna_=gna*minf(V)
# INITIAL CONDITIONS
z(0)=0
nai(0)=18.67
V(0)=0
w(0)=0.000025
#STIMULUS PARAMETERS
# stimulus, set I_stims to 0 if you want to use noisy stim, slope to control ramp
param idc=0
I1(t)=heav(t>=t_start)*I_stim1
# Control stimulus offset or twostage step, slope to control ramp
I2(t)=heav(t>(t_start+t_run))*I_stim1*(1)
param I_stim1=50,t_start=300,t_run=0
# ALWAYS USE EULER when using noise, but no noise processes are implemented here
@ total=150,dt=.05,xlo=100,xhi=60,ylo=.125,yhi=.6,xp=v,yp=w
@ meth=euler
@ MAXSTOR=1000000,bounds=10000
done
