Modifying standard Hodgkin-Huxley dynamics

[pascal]

I am new to NEURON, and I decided to learn how to modify mechanisms by trying to replicate the Wang-Buzsaki neuron (from the 1996 paper "Gamma Oscillation by Synaptic Inhibition in a Hippocampal Interneuronal Network Model"). The model is very simple, using a single compartment with Hodgkin-Huxley-like dynamics, but for some reason when I try to replicate it I get very unstable dynamics.

My approach was to modify the file hh.mod found in nrn/src/nrnoc. Here is the revised mod file, hhbuz.mod:

Code: Select all

UNITS {
        (mA) = (milliamp)
        (mV) = (millivolt)
	(S) = (siemens)
}
 
? interface
NEURON {
        SUFFIX hhbuz
        USEION na READ ena WRITE ina
        USEION k READ ek WRITE ik
        NONSPECIFIC_CURRENT il
        RANGE gnabar, gkbar, gl, el, gna, gk
	   : do not need mtau
        : GLOBAL minf, hinf, ninf, mtau, htau, ntau 
	   GLOBAL minf, hinf, ninf, htau, ntau
	THREADSAFE : assigned GLOBALs will be per thread
}
 
PARAMETER {
	   : modify these values
        :gnabar = .12 (S/cm2)	<0,1e9>
        :gkbar = .036 (S/cm2)	<0,1e9>
        :gl = .0003 (S/cm2)	<0,1e9>
        :el = -54.3 (mV)

        gnabar = .035 (S/cm2)	<0,1e9>
        gkbar = .009 (S/cm2)	<0,1e9>
        gl = .0001 (S/cm2)	<0,1e9>
        el = -65 (mV)
}
 
STATE {
	: do not need m, since we are just using minf
      :  m h n
	h n
}
 
ASSIGNED {
        v (mV)
        celsius (degC)
        ena (mV)
        ek (mV)

	gna (S/cm2)
	gk (S/cm2)
        ina (mA/cm2)
        ik (mA/cm2)
        il (mA/cm2)
        minf hinf ninf
	: do not need mtau
	: mtau (ms) htau (ms) ntau (ms)
	mtau (ms) htau (ms) ntau (ms)
}
 
? currents
BREAKPOINT {
        SOLVE states METHOD cnexp
        : gna = gnabar*m*m*m*h
	gna = gnabar*minf*minf*minf*h 
	ina = gna*(v - ena)
        gk = gkbar*n*n*n*n
	ik = gk*(v - ek)      
        il = gl*(v - el)
}
 
 
INITIAL {
	rates(v)
:	m = minf  we do not have an 'm' gating variable, since using m_inf
	h = hinf
	n = ninf
}

? states
DERIVATIVE states {  
        rates(v)
	:   m' =  (minf-m)/mtau
        h' = 5*(hinf-h)/htau
        n' = 5*(ninf-n)/ntau
}
 
:LOCAL q10


? rates
PROCEDURE rates(v(mV)) {  :Computes rate and other constants at current v.
                      :Call once from HOC to initialize inf at resting v.
        LOCAL  alpha, beta, sum, q10
        TABLE minf, mtau, hinf, htau, ninf, ntau DEPEND celsius FROM -100 TO 100 WITH 200

:modified value according to Buzsaki model
UNITSOFF
     :   q10 = 3^((celsius - 6.3)/10)
                :"m" sodium activation system
        alpha = .1 * vtrap(-(v+35),10)
        beta =  4 * exp(-(v+60)/18)
        sum = alpha + beta
	:mtau = 1/(q10*sum)
        minf = alpha/sum
                :"h" sodium inactivation system
        alpha = .07 * exp(-(v+58)/20)
        beta = 1 / (exp(-(v+28)/10) + 1)
        sum = alpha + beta
	htau = 1/(q10*sum)
        hinf = alpha/sum
                :"n" potassium activation system
        alpha = .01*vtrap(-(v+34),10) 
        beta = .125*exp(-(v+44)/80)
	sum = alpha + beta
        ntau = 1/(q10*sum)
        ninf = alpha/sum
}
 
FUNCTION vtrap(x,y) {  :Traps for 0 in denominator of rate eqns.
        if (fabs(x/y) < 1e-6) {
                vtrap = y*(1 - x/y/2)
        }else{
                vtrap = x/(exp(x/y) - 1)
        }
}
 
UNITSON

And here is my hoc file, which simply plots the voltage trace of the neuron for 125 ms, with the neuron receiving input current from a current clamp (0.001 nA) for the last 100 ms:

Code: Select all

create soma 
access soma 
// I copied and pasted this from Buzsaki's hoc file on modeldb
  soma {
    Ra = 100		// geometry 
    nseg = 1
    cm=1
    diam = 10.
    L = 10./PI
  }

insert hhbuz
soma ek=-90
soma ena=55

//insert hh
//soma ek=-77
//soma ena=50

psection()

//set up current to  be injected into model neuron
objref clamp
soma clamp = new IClamp(0.5)
clamp.del = 25
clamp.dur = 100
clamp.amp = 0.001

//graphical display
objref g
g = new Graph()
g.size(0,125,-80,40)
g.addvar("soma.v(0.5)")

//simulation control
dt=0.025
tstop=125
v_init=-70

proc initialize() {
  finitialize(v_init)
  fcurrent()
}

proc integrate() {
  g.begin()

  while(t<tstop) {
    fadvance()
    g.plot(t)
  }
  g.flush()
}

proc go() {
  initialize()
  integrate()
}

Here is the plot of the resultant voltage trace (using NEURON 7.2 on a Windows 7 machine):




Uploaded with ImageShack.us


In the above hoc file, if I comment the lines insert hhbuz, soma ek=-90, and soma ena=55, and then uncomment the lines insert hh, soma ek=-77, and soma ena=50, I get the following plot:



Uploaded with ImageShack.us

So the code works just fine for the normal HH model, but when I revise it to try to replicate the Wang-Buzsaki neuron, it fails miserably. Any ideas about what I did wrong would be much appreciated. Thanks.

[ted]

A reasonable strategy for debugging voltage-gated channel mechanisms is to start by examining the voltage-dependence of the rate constants or time constants, and steady-state values. The steady state values seem to have reasonable voltage dependence, but the time constants htau and ntau are both infinite. The cause is that q10 is used in the calculation of htau and ntau, but q10 itself is never assigned a value.

[pascal]

Thanks for the help, Ted. It definitely works now that I uncommented the q10 assignment. How do I investigate the voltage dependence of htau and ntau? I tried to record their values throughout the simulation, but that did not work, and I also I couldn't simply print them out using the 'print' command, like you can do with the voltage...

[ted]

Look in the Forum's "Hot tips" area for the topic
Using the Grapher to plot voltage-dependent gating variables

While you're there, see what else might interest you.