COMMENT ---------------------------------------------------------------- This a stochastic version of the na3h5.mod of Z. Mainen in Mainen & Sejnowski 95 (Zach Mainen, Salk Institute, 1994, zach@salk.edu). It is a sodium channel, Hodgkin-Huxley like kinetics, based on the gates model, assuming stochastic opening and closing. The rates function are adapted directly from the na3h5.mod by Zach Mainen, available through: http://www.cnl.salk.edu/~zach/initdemo.html The stochastic model is similar to that used by Schneidman, Freedman, Segev, 1998 and is as follows (ascii-art): = 3alpha_m => = 2alpha_m => = alpha_m => [m0h1] [m1h1] [m2h1] <[m3h1]> <= beta_m = <= 2beta_m = <= 3beta_m = ^ ^ ^ ^ beta_h| | beta_h| | beta_h| | beta_h| | | |alpha_h | |alpha_h | | alpha_h | | alphah | | | | | | | | V = 3alpha_m => V = 2alpha_m => V = alpha_m => V [m0h0] [m1h0] [m2h0] [m3h0] <= beta_m = <= 2beta_m = <= 3beta_m = The model keeps track of the number of channels in each state and uses a binomial rng to update this number. The model requires that the RNG mechanism be inserted somewhere in the simulation in order to provide the BnlDev_RNG function. Jan 1999 Mickey London, Hebrew University, 1998, mikilon@lobster.ls.huji.ac.il Peter N. Steinmetz, Caltech, peter@klab.caltech.edu 14 Sep 99 PNS. Added deterministic flag. 01 Sep 02 K. Diba changed deterministic to RANGE variable ---------------------------------------------------------------- ENDCOMMENT INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)} NEURON { SUFFIX sna USEION na READ ena WRITE ina GLOBAL minf, hinf, mtau, htau,am,bm,ah,bh RANGE N, reff, eta, gamma, deterministic, gna GLOBAL P_am, P_bm, P_ah, P_bh, hinf_curve GLOBAL tha, thi1, thi2, qa, qi, qinf, thinf, vshift GLOBAL Ra, Rb, Rd, Rg GLOBAL vmin, vmax, q10, orig_temp, wflag, tadj GLOBAL DONT_VECTORIZE : prevent vectorization to agree with RNG.mod } UNITS { (mA) = (milliamp) (mV) = (millivolt) (pS) = (picosiemens) (um) = (micron) } PARAMETER { v (mV) dt (ms) area vshift = -10 (mV) : voltage shift (affects all) tha = -35 (mV) : v 1/2 for act (-42) qa = 9 (mV) : act slope Ra = 0.182 (/ms) : open (v) Rb = 0.124 (/ms) : close (v) thi1 = -50 (mV) : v 1/2 for inact thi2 = -75 (mV) : v 1/2 for inact qi = 5 (mV) : inact tau slope hinf_curve = 1 : 0 if use normal relation for hinf : otherwise use curve specified by following : two parameters thinf = -65 (mV) : inact inf slope qinf = 6.2 (mV) : inact inf slope Rg = 0.0091 (/ms) : inact (v) Rd = 0.024 (/ms) : inact recov (v) gamma = 20 (pS) eta = 50 (1/um2) : channel density celsius (degC) orig_temp = 23 (degC) : original temperature q10 = 2.3 : temperature sensitivity deterministic = 0 : if non-zero, use deterministic variables vmin = -120 (mV) : range to construct tables for vmax = 100 (mV) DONT_VECTORIZE : required declaration } ASSIGNED { ina (mA/cm2) gna (pS/um2) ena (mV) am (/ms) bm (/ms) ah (/ms) bh (/ms) minf hinf mtau (ms) htau (ms) tadj N reff (pS/um2) scale_dens (pS/um2) P_am : probability of one channel making am transition P_bm P_ah P_bh wflag } STATE { m h : deterministic variables m0h0 m0h1 m1h0 m1h1 m2h0 m2h1 m3h0 m3h1 : states m0h0_m1h0 m1h0_m2h0 m2h0_m3h0 m0h1_m1h1 m1h1_m2h1 m2h1_m3h1 m3h0_m2h0 m2h0_m1h0 m1h0_m0h0 m3h1_m2h1 m2h1_m1h1 m1h1_m0h1 m0h0_m0h1 m0h1_m0h0 m1h0_m1h1 m1h1_m1h0 m2h0_m2h1 m2h1_m2h0 m3h0_m3h1 m3h1_m3h0 } : ---------------------------------------------------------------- : Initialization. INITIAL { trates(v+vshift) wflag = 1 : only give a warning once! m = minf h = hinf scale_dens = gamma/area reff = gamma*eta N = floor(eta*area + 0.5) : round to nearest number of channels m1h0 = floor(3*m*(1-m)*(1-m)*(1-h)*N + 0.5) m2h0 = floor(3*m*m*(1-m)*(1-h)*N + 0.5) m3h0 = floor(m*m*m*(1-h)*N + 0.5) m0h1 = floor((1-m)*(1-m)*(1-m)*h*N + 0.5) m1h1 = floor(3*m*(1-m)*(1-m)*h*N + 0.5) m2h1 = floor(3*m*m*(1-m)*h*N + 0.5) m3h1 = floor(m*m*m*h*N + 0.5) : put the rest of the channels in the non-conducting & inactivated state m0h0 = N - (m1h0 + m2h0 + m3h0 + m0h1 + m1h1 + m2h1 + m3h1) m0h0_m1h0=0 m1h0_m2h0=0 m2h0_m3h0=0 m0h1_m1h1=0 m1h1_m2h1=0 m2h1_m3h1=0 m3h0_m2h0=0 m2h0_m1h0=0 m1h0_m0h0=0 m3h1_m2h1=0 m2h1_m1h1=0 m1h1_m0h1=0 m0h0_m0h1=0 m0h1_m0h0=0 m1h0_m1h1=0 m1h1_m1h0=0 m2h0_m2h1=0 m2h1_m2h0=0 m3h0_m3h1=0 m3h1_m3h0=0 } : ---------------------------------------------------------------- : Breakpoint for each integration step BREAKPOINT { SOLVE states if (deterministic) { if (deterministic-1){ gna = m*m*m*h*reff*tadj } else { gna = floor(m*m*m*h* N + 0.5) * scale_dens *tadj} } else{ gna = strap(m3h1) * scale_dens * tadj } ina = (1e-4) * gna * (v - ena) } : ---------------------------------------------------------------- : states - compute state variables PROCEDURE states() { VERBATIM extern double BnlDev_RNG(); ENDVERBATIM trates(v+vshift) : deterministic versions of state variables : integrated by relaxing toward the steady state value m = m + (1 - exp(-dt/mtau)) * (minf-m) h = h + (1 - exp(-dt/htau)) * (hinf-h) P_am = strap(am*dt) P_bm = strap(bm*dt) : check that will represent probabilities when used ChkProb( 3.0 * P_am) ChkProb( 3.0 * P_bm) ChkProb( P_bm/(1.0-2.0*P_am) ) ChkProb( 2.0 * P_bm/(1.0-P_am) ) : m gate transitions : if (deterministic) { : m0h0_m1h0 = 3.0*P_am*m0h0 : m1h0_m2h0 = 2.0*P_am*m1h0 : m1h0_m0h0 = P_bm/(1.0-2.0*P_am)*(m1h0 - m1h0_m2h0) : m2h0_m3h0 = P_am*m2h0 : m2h0_m1h0 = 2.0*P_bm/(1.0-P_am)*(m2h0 - m2h0_m3h0) : m3h0_m2h0 = 3.0*P_bm*m3h0 : m0h1_m1h1 = 3.0*P_am*m0h1 : m1h1_m2h1 = 2.0*P_am*m1h1 : m1h1_m0h1 = P_bm/(1.0-2.0*P_am)*(m1h1 - m1h1_m2h1) : m2h1_m3h1 = P_am*m2h1 : m2h1_m1h1 = 2.0*P_bm/(1.0-P_am)*(m2h1 - m2h1_m3h1) : m3h1_m2h1 = 3.0*P_bm*m3h1 : } : else { m0h0_m1h0 = BnlDev_RNG(3.0*P_am,m0h0) m1h0_m2h0 = BnlDev_RNG(2.0*P_am,m1h0) m1h0_m0h0 = BnlDev_RNG(P_bm/(1.0-2.0*P_am), m1h0 - m1h0_m2h0) m2h0_m3h0 = BnlDev_RNG(P_am,m2h0) m2h0_m1h0 = BnlDev_RNG(2.0*P_bm/(1.0-P_am), m2h0 - m2h0_m3h0) m3h0_m2h0 = BnlDev_RNG(3.0*P_bm,m3h0) m0h1_m1h1 = BnlDev_RNG(3.0*P_am, m0h1) m1h1_m2h1 = BnlDev_RNG(2.0*P_am, m1h1) m1h1_m0h1 = BnlDev_RNG(P_bm/(1.0-2.0*P_am), m1h1 - m1h1_m2h1) m2h1_m3h1 = BnlDev_RNG(P_am,m2h1) m2h1_m1h1 = BnlDev_RNG(2.0*P_bm/(1.0-P_am), m2h1 - m2h1_m3h1) m3h1_m2h1 = BnlDev_RNG(3.0*P_bm,m3h1) : } : new numbers in each state after the m gate transitions m0h0 = m0h0 - m0h0_m1h0 + m1h0_m0h0 m1h0 = m1h0 - m1h0_m2h0 - m1h0_m0h0 + m2h0_m1h0 + m0h0_m1h0 m2h0 = m2h0 - m2h0_m3h0 - m2h0_m1h0 + m3h0_m2h0 + m1h0_m2h0 m3h0 = m3h0 - m3h0_m2h0 + m2h0_m3h0 m0h1 = m0h1 - m0h1_m1h1 + m1h1_m0h1 m1h1 = m1h1 - m1h1_m2h1 - m1h1_m0h1 + m2h1_m1h1 + m0h1_m1h1 m2h1 = m2h1 - m2h1_m3h1 - m2h1_m1h1 + m3h1_m2h1 + m1h1_m2h1 m3h1 = m3h1 - m3h1_m2h1 + m2h1_m3h1 : probabilities of making h gate transitions P_ah = strap(ah*dt) P_bh = strap(bh*dt) ChkProb(P_ah) ChkProb(P_bh) : number making h gate transitions : if (deterministic) { : m0h0_m0h1 = P_ah*m0h0 : m0h1_m0h0 = P_bh*m0h1 : m1h0_m1h1 = P_ah*m1h0 : m1h1_m1h0 = P_bh*m1h1 : m2h0_m2h1 = P_ah*m2h0 : m2h1_m2h0 = P_bh*m2h1 : m3h0_m3h1 = P_ah*m3h0 : m3h1_m3h0 = P_bh*m3h1 : } : else { m0h0_m0h1 = BnlDev_RNG(P_ah,m0h0) m0h1_m0h0 = BnlDev_RNG(P_bh,m0h1) m1h0_m1h1 = BnlDev_RNG(P_ah,m1h0) m1h1_m1h0 = BnlDev_RNG(P_bh,m1h1) m2h0_m2h1 = BnlDev_RNG(P_ah,m2h0) m2h1_m2h0 = BnlDev_RNG(P_bh,m2h1) m3h0_m3h1 = BnlDev_RNG(P_ah,m3h0) m3h1_m3h0 = BnlDev_RNG(P_bh,m3h1) : } m0h0 = m0h0 - m0h0_m0h1 + m0h1_m0h0 m1h0 = m1h0 - m1h0_m1h1 + m1h1_m1h0 m2h0 = m2h0 - m2h0_m2h1 + m2h1_m2h0 m3h0 = m3h0 - m3h0_m3h1 + m3h1_m3h0 m0h1 = m0h1 - m0h1_m0h0 + m0h0_m0h1 m1h1 = m1h1 - m1h1_m1h0 + m1h0_m1h1 m2h1 = m2h1 - m2h1_m2h0 + m2h0_m2h1 m3h1 = m3h1 - m3h1_m3h0 + m3h0_m3h1 } : ---------------------------------------------------------------- : trates - compute rates, using table if possible PROCEDURE trates(vm) { TABLE minf, mtau, hinf, htau, am, bm, ah, bh, tadj DEPEND dt,Ra,Rb,Rd,Rg,tha,thi1,thi2,qa,qi,qinf,q10,orig_temp,celsius, hinf_curve FROM vmin TO vmax WITH 199 tadj = q10^((celsius-orig_temp)/10) : m activation variable am = SigmoidRate(vm,tha,Ra,qa) am = am * tadj bm = SigmoidRate(-vm,-tha,Rb,qa) bm = bm * tadj mtau = 1/(am+bm) minf = am*mtau : h inactivation variable ah = SigmoidRate(vm,thi1,Rd,qi) ah = ah * tadj bh = SigmoidRate(-vm,-thi2,Rg,qi) bh = bh * tadj htau = 1/(ah+bh) : hinf_curve is non-zero if using explicit fit, zero otherwise if (hinf_curve == 0) { hinf = ah*htau } else { hinf = 1/(1+exp((vm-thinf)/qinf)) : recompute rates based on hinf ah = hinf/htau bh = 1/htau - ah } } : ---------------------------------------------------------------- : SigmoidRate - Compute a sigmoid rate function given the : 50% point th, the slope q, and the amplitude a. FUNCTION SigmoidRate(v,th,a,q) { if (fabs(v-th) > 1e-6) { SigmoidRate = a * (v - th) / (1 - exp(-(v - th)/q)) } else { SigmoidRate = a * q } } : ---------------------------------------------------------------- : sign trap - trap negative numbers and replace with zero FUNCTION strap(x) { if (x < 0) { strap = 0 if (wflag){ VERBATIM fprintf (stderr,"sna.mod:strap: negative value for state"); ENDVERBATIM wflag = 0} } else { strap = x } } : ---------------------------------------------------------------- : ChkProb - Check that number represents a probability PROCEDURE ChkProb(p) { if (p < 0.0 || p > 1.0) { if (wflag){ VERBATIM fprintf(stderr, "sna.mod:ChkProb: argument not a probability.\n"); ENDVERBATIM wflag =0} } }