Recording range variable in Vector

[regger]

Hi,

I'm having a problem recording a range variable into a Vector. Specifically, I'm using the Mainen et al. 1995 spike initiation model (http://senselab.med.yale.edu/modeldb/Sh ... model=8210) and want to look at the different state variables, for example of the sodium channels.

Just to demonstrate what I mean I put together this little example script:

Code: Select all

import neuron
import numpy as np
import matplotlib.pyplot as plt
h = neuron.h

Rm = 15000.0
cm = 1.0
Ra = 150.0
ePas = -65.0

soma = h.Section()
soma.L = 15.0
soma.diam = 15.0
soma.nseg = 1

soma.cm = cm
soma.Ra = Ra
soma.insert('pas')
for seg in soma:
    seg.pas.g = 1/Rm
    seg.pas.e = ePas
soma.insert('na3')
soma.insert('kd3')

rec = {}
for lbl in 't','vSoma':
   rec[lbl] = h.Vector()
rec['t'].record(h._ref_t)
rec['vSoma'].record(soma(0.5)._ref_v)
stateVars = {}
for var in 'm', 'h', 'minf', 'hinf', 'mtau', 'htau':
    stateVars[var] = h.Vector()
stateVars['m'].record(soma(0.5).na3._ref_m, sec=soma)
stateVars['h'].record(soma(0.5).na3._ref_h, sec=soma)
stateVars['mtau'].record(soma(0.5).na3._ref_mtau, sec=soma)
stateVars['htau'].record(soma(0.5).na3._ref_htau, sec=soma)
stateVars['minf'].record(soma(0.5).na3._ref_minf, sec=soma)
stateVars['hinf'].record(soma(0.5).na3._ref_hinf, sec=soma)

somaIC = h.IClamp(0.5, sec=soma)
somaIC.delay = 10.0
somaIC.dur = 90.0
somaIC.amp = 0.01

h.load_file('stdrun.hoc')
h.dt = 0.025
h.celsius = 37.0
h.finitialize(-65.0)
h.tstop = 100.0
h.run()

tArray = np.array(rec['t'])
vArray = np.array(rec['vSoma'])
plt.figure(1)
plt.subplot(221)
plt.plot(tArray, vArray, 'b')
plt.ylabel('Voltage [mV]')
plt.xlabel('Time [ms]')
plt.subplot(222)
plt.plot(tArray, stateVars['m'], 'b', label='m')
plt.plot(tArray, stateVars['h'], 'r', label='h')
plt.ylabel('Na state variables m/h')
plt.xlabel('Time [ms]')
plt.legend()
plt.subplot(223)
plt.plot(tArray, stateVars['mtau'], 'b', label='$\\tau_{m}$')
plt.plot(tArray, stateVars['htau'], 'r', label='$\\tau_{h}$')
plt.ylabel('Na state variables $\\tau_m$/$\\tau_h$')
plt.xlabel('Time [ms]')
plt.legend()
plt.subplot(224)
plt.plot(tArray, stateVars['minf'], 'b', label='$m_{\infty}$')
plt.plot(tArray, stateVars['hinf'], 'r', label='$h_{\infty}$')
plt.ylabel('Na state variables $m_{\infty}$/$h_{\infty}$')
plt.xlabel('Time [ms]')
plt.legend()

plt.show()

Recording m, h, minf and hinf works. However, recording htau and mtau does not work. The Vector is just filled with zeros after the simulation run.

Here's the mod file:

Code: Select all

COMMENT

na3h5.mod

Sodium channel, Hodgkin-Huxley style kinetics.  

Kinetics were fit to data from Huguenard et al. (1988) and Hamill et
al. (1991)

qi is not well constrained by the data, since there are no points
between -80 and -55.  So this was fixed at 5 while the thi1,thi2,Rg,Rd
were optimized using a simplex least square proc

voltage dependencies are shifted approximately +5mV from the best
fit to give higher threshold

use with kd3h5.mod

Author: Zach Mainen, Salk Institute, 1994, zach@salk.edu

ENDCOMMENT

INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}

NEURON {
	SUFFIX na3
	USEION na READ ena WRITE ina
	RANGE m, h, gna, gbar, vshift
	GLOBAL tha, thi1, thi2, qa, qi, qinf, thinf
	RANGE minf, hinf, mtau, htau
	GLOBAL Ra, Rb, Rd, Rg
	GLOBAL q10, temp, tadj, vmin, vmax
}

PARAMETER {
	gbar = 1200   	(pS/um2)	: 0.12 mho/cm2
	vshift = 0	(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
	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) 

	temp = 23	(degC)		: original temp 
	q10  = 2.3			: temperature sensitivity

	v 		(mV)
	dt		(ms)
	celsius		(degC)
	vmin = -120	(mV)
	vmax = 100	(mV)
}


UNITS {
	(mA) = (milliamp)
	(mV) = (millivolt)
	(pS) = (picosiemens)
	(um) = (micron)
} 

ASSIGNED {
	ina 		(mA/cm2)
	gna		(pS/um2)
	ena		(mV)
	minf 		hinf
	mtau (ms)	htau (ms)
	tadj
}
 

STATE { m h }

INITIAL { 
	trates(v+vshift)
	m = minf
	h = hinf
}

BREAKPOINT {
        SOLVE states
        gna = gbar*m*m*m*h
	ina = (1e-4) * gna * (v - ena)
} 

LOCAL mexp, hexp 

PROCEDURE states() {   :Computes state variables m, h, and n 
        trates(v+vshift)      :             at the current v and dt.
        m = m + mexp*(minf-m)
        h = h + hexp*(hinf-h)
        VERBATIM
        return 0;
        ENDVERBATIM
}

PROCEDURE trates(v) {  
                      
        LOCAL tinc
        TABLE minf, mexp, hinf, hexp
	DEPEND dt, celsius, temp, Ra, Rb, Rd, Rg, tha, thi1, thi2, qa, qi, qinf
	
	FROM vmin TO vmax WITH 199

	rates(v): not consistently executed from here if usetable == 1

        tadj = q10^((celsius - temp)/10)
        tinc = -dt * tadj

        mexp = 1 - exp(tinc/mtau)
        hexp = 1 - exp(tinc/htau)
}


PROCEDURE rates(vm) {  
        LOCAL  a, b

	a = trap0(vm,tha,Ra,qa)
	b = trap0(-vm,-tha,Rb,qa)
	mtau = 1/(a+b)
	minf = a*mtau

		:"h" inactivation 

	a = trap0(vm,thi1,Rd,qi)
	b = trap0(-vm,-thi2,Rg,qi)
	htau = 1/(a+b)
	hinf = 1/(1+exp((vm-thinf)/qinf))
}


FUNCTION trap0(v,th,a,q) {
	if (fabs(v/th) > 1e-6) {
	        trap0 = a * (v - th) / (1 - exp(-(v - th)/q))
	} else {
	        trap0 = a * q
 	}
}

I can get the recording to work when I also create tables of htau and mtau; i.e. I change the procedure trates to the following:

Code: Select all

PROCEDURE trates(v) {  
                      
        LOCAL tinc
        TABLE minf, mexp, hinf, hexp, mtau, htau
	DEPEND dt, celsius, temp, Ra, Rb, Rd, Rg, tha, thi1, thi2, qa, qi, qinf
	
	FROM vmin TO vmax WITH 199

	rates(v): not consistently executed from here if usetable == 1

        tadj = q10^((celsius - temp)/10)
        tinc = -dt * tadj

        mexp = 1 - exp(tinc/mtau)
        hexp = 1 - exp(tinc/htau)
}

So, I'm trying to understand why it works when I create the tables of mtau and htau and otherwise it doesn't. What am I missing? Or am I doing something wrong?

Thanks,

Robert

[ted]

Using the original mod file, I just ran my own test but with model specification and user interface done in hoc, and got a different result: the mtau and htau values, whether displayed via NEURON's InterViews-based GUI or captured to Vectors, weren't zero, but they didn't evolve over the course of a simulation--instead they were 0.188... and 17.4..., respectively.

Then I disabled use of tables by executing
usetable_na3 = 0
at the oc> prompt (see
http://www.neuron.yale.edu/neuron/stati ... tml#Table)
(equivalent from Python would probably be h.usetable_na3 = 0) and ran another simulation. This time mtau and htau varied appropriately with time. In both cases, m, minf, h, and hinf all had appropriate time courses.

So now there are two sets of observations that I don't understand--what you saw with Python, and what I saw with hoc.

[regger]

Ok I just ran it again and I also got non-zero, but constant values for htau and mtau this time (0.27... and 0.04... respectively). The values only become zero when I run it repeatedly from ipython in the same session.
Setting h.usetable_na3 = 0 like you suggested also fixes it. So it looks like it works if tables are created for either all or none of the parameters...
I don't know much about the inner structure of NEURON, but this looks like the reference/pointer given to the recording vector might somehow be invalidated/overwritten...

Also, might this have something to do with Zach's cryptic comment in the mod file:

Code: Select all

PROCEDURE trates(v) {  
        ...
	rates(v): not consistently executed from here if usetable == 1
        ...
}

I guess it's ok for now, because I can get it to work. But it's beginning to look a little like a bug to me...

[ted]

Well, there is a bug, but it's a bug in na3h5.mod, and it's there because the model's authors didn't understand the implications of using a TABLE statement. And I was slow to see the bug because it caused such a small error in my test simulation that I ignored what I knew about the use of TABLEs (about which more later). However, the error exists: a slight difference in the time courses of v and other variables generated by simulations executed with usetable_na3 set to 1 or 0. I was expecting the difference to be so large as to be grossly obvious, but it wasn't, despite the fact that, during a spike simulation, mtau_na3 and htau_na3 vary by factors of ~4 and ~100 when usetable_na3 is 0 and remain constant when usetable_na3 is 1.

What lessons may be drawn from this?

As a general observation, this is a nice illustration of the fact that the voltage dependence of certain rate constants may have little or no qualitative effect on simulation results.

With regard to NMODL specifications of voltage-gated channels, the key points are:

As per Michael Hines
"When trates is called during simulation, the table is used and rates is NOT called. Since the table contains only mexp and hexp instead of htau and mtau, the values of the latter remain as they were when the last table entry was constructed."

In this particular case, it was purely accidental that, with usetables_na3==1, mtau_na3 and htau_na3 just happened to end up with values that gave a very reasonable spike threshold and waveform. There is no guarantee that this will always happen. In computational modeling, nothing should be accidental (i.e. not under of the modeler's control), not even stochasticity. It is easy to imagine scenarios in which accidental values of htau_na3 and mtau_na3 would have caused fatal errors. For example, set usetables_na3 to 0 and run a long simulation in which the model is driven by a prolonged depolarizing current. htau and mtau would then have very different values, and a new simulation after usetables_na3=1 would produce very different results.

So this mechanism has a bug and the fix for it is to change the TABLE statement to
TABLE minf, mexp, mtau, hinf, hexp, htau

The kd3h5.mod mechanism has a similar bug, and the fix is to change its TABLE statement to
TABLE ninf, nexp, ntau

Finally, as a general rule, if a TABLE statement is used, then all time constants, steady state values, and "exp" values (if such are present--the use of those is highly deprecated) should be included in the TABLE statement.

[regger]

Thank you for following up on this so quickly.
Your explanation makes sense. Crucial point that's definitely worth remembering for other mechanisms.

All the best,

Robert

[ted]

Thank you for following up on this so quickly.

Actually I now realize that most of what I wrote in that post is a lot of rubbish. I must apologize. The one part that wasn't rubbish was when I was quoting Hines.

mtau and htau are merely intermediate variables that are used in the calculation of minf, mexp, hinf, and hexp. The latter four are the variables that are used by this "analytical solution" channel specification to advance simulations from one step to the next, and they are the only variables that have to appear in the TABLE statement. Furthermore, the slight difference that I saw between results with usetable on or off amounts to an offset of about 1 time step, just enough to be detected yet truly inconsequential. The h and m values used during simulations are properly controlled by the tables of the voltage dependence of minf, mexp, hinf, and hexp. Finally, throw out my statement that the simulation would be unaffected if mtau and htau remained constant with their "accidental" values--I tested that by commenting out the TABLE statement in PROCEDURE trates (i.e. the three lines that begin with TABLE, DEPEND, and FROM), and inserting
mtau = 0.0407
htau = 0.277778
just before the closing bracket of PROCEDURE rates, recompiling the mechanisms, and running a simulation. Result: the spike waveform was markedly different (much faster rise and fall, and much briefer duration).

So there really isn't anything wrong with na3h5.mod OR kd3h5.mod. They are correct as written, and they work properly for the purpose required by the original model authors.

Broader implications? Only the variables that are actually used during an fadvance need to be included in a TABLE statement. If there are intermediate variables whose time course is of interest, but which remain constant when usetable is 1, it is only necessary to set the mechanism's usetable to 0 before running a simulation.