CVODE and *.mod file

[chinhou]

Hello,

I'm a new user of Neuron. I have read some Neuron documents to start writing code in NEURON. So, this is my first step to write a code in NEURON.
I 'm sure my question is very simple for you, but the answer will help me so much.this is my code:

Code: Select all

create soma, apical_dend,distal_dend
 
//define geometry
insert pass
 insert hh2
 // connect section

//d_lambda rule  
xopen("fixnseg.hoc") //  i found this code in your site
geom_nseg()            // so i have different nseg for my sections

objref cvode
cvode = new CVode()
cvode.active(1)
cvode.atol(1e-4)

when i click on my hoc file i see the following error:

nrniv: hh2 cannot be used with CVODE in C:/neuron1.hoc near line 30
cvode.active(1)
^
CVode[1].active(1 )

hh2 is a mod file. that i found in ModelDB

Code: Select all

TITLE Hippocampal HH channels
:
: Fast Na+ and K+ currents responsible for action potentials
: Iterative equations
:
: Equations modified by Traub, for Hippocampal Pyramidal cells, in:
: Traub & Miles, Neuronal Networks of the Hippocampus, Cambridge, 1991
:
: range variable vtraub adjust threshold
:
: Written by Alain Destexhe, Salk Institute, Aug 1992
:

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

NEURON {
	SUFFIX hh2
	USEION na READ ena WRITE ina
	USEION k READ ek WRITE ik
	RANGE gnabar, gkbar, vtraub
	RANGE m_inf, h_inf, n_inf
	RANGE tau_m, tau_h, tau_n
	RANGE m_exp, h_exp, n_exp
}


UNITS {
	(mA) = (milliamp)
	(mV) = (millivolt)
}

PARAMETER {
	gnabar	= .003 	(mho/cm2)
	gkbar	= .005 	(mho/cm2)

	ena	= 50	(mV)
	ek	= -90	(mV)
	celsius = 36    (degC)
	dt              (ms)
	v               (mV)
	vtraub	= -63	(mV)
}

STATE {
	m h n
}

ASSIGNED {
	ina	(mA/cm2)
	ik	(mA/cm2)
	il	(mA/cm2)
	m_inf
	h_inf
	n_inf
	tau_m
	tau_h
	tau_n
	m_exp
	h_exp
	n_exp
	tadj
}


BREAKPOINT {
	SOLVE states
	ina = gnabar * m*m*m*h * (v - ena)
	ik  = gkbar * n*n*n*n * (v - ek)
}


:DERIVATIVE states {   : exact Hodgkin-Huxley equations
:	evaluate_fct(v)
:	m' = (m_inf - m) / tau_m
:	h' = (h_inf - h) / tau_h
:	n' = (n_inf - n) / tau_n
:}

PROCEDURE states() {	: exact when v held constant
	evaluate_fct(v)
	m = m + m_exp * (m_inf - m)
	h = h + h_exp * (h_inf - h)
	n = n + n_exp * (n_inf - n)
	VERBATIM
	return 0;
	ENDVERBATIM
}

UNITSOFF
INITIAL {
	m = 0
	h = 0
	n = 0
:
:  Q10 was assumed to be 3 for both currents
:
: original measurements at roomtemperature?

	tadj = 3.0 ^ ((celsius-36)/ 10 )
}

PROCEDURE evaluate_fct(v(mV)) { LOCAL a,b,v2

	v2 = v - vtraub : convert to traub convention

	a = 0.32 * (13-v2) / ( exp((13-v2)/4) - 1)
	b = 0.28 * (v2-40) / ( exp((v2-40)/5) - 1)
	tau_m = 1 / (a + b) / tadj
	m_inf = a / (a + b)

	a = 0.128 * exp((17-v2)/18)
	b = 4 / ( 1 + exp((40-v2)/5) )
	tau_h = 1 / (a + b) / tadj
	h_inf = a / (a + b)

	a = 0.032 * (15-v2) / ( exp((15-v2)/5) - 1)
	b = 0.5 * exp((10-v2)/40)
	tau_n = 1 / (a + b) / tadj
	n_inf = a / (a + b)

	m_exp = 1 - exp(-dt/tau_m)
	h_exp = 1 - exp(-dt/tau_h)
	n_exp = 1 - exp(-dt/tau_n)
}

UNITSON

Although the big trouble is wrong way of using cvode(i guess), but let my ask my question in order:

1- Does my total nseg mean (for example secname.nseg=2 , 3, 4 then total nseg=9 ) number of compartment or as i defined in code i have 3 compartment due to 3 section?

2- From hh2 code, "vtraub" reperesnt the threshold for spike, right?
but which parameter define the maximum value of membrane potential during a action potential ? (i'm not sure even this question is correct!)

3-and finally, why i have this error?



In advance, i grateful for your kindly advice!

[ted]

when i click on my hoc file i see the following error:

Code: Select all

nrniv: hh2 cannot be used with CVODE in C:/neuron1.hoc near line 30
cvode.active(1)
                ^
CVode[1].active(1        )

The error message occurs because the time step dt is an explicit variable in the source
code for the hh2 mechanism (see PROCEDURE evaluate_fct). This is sometimes called
"explicit exponential integration." Mechanisms that make explicit use of dt are not compatible with CVODE. Please see
http://www.neuron.yale.edu/neuron/stati ... cvode.html
especially the DESCRIPTION.

Two questions for you:
1. Why are you using CVODE? You're just starting to learn how to use NEURON. Why
add extra complications at the very beginning?
2. Do you really need to use the hh2 mechanism, or did it just "seem interesting"? If you
really need that particular mechanism, and you absolutly must also use CVODE, then
you're going to have to make hh2 CVODE-compatible, or maybe it would make sense
to see if a CVODE-compatible implementation of that mechanism is already available.

1- Does my total nseg mean (for example secname.nseg=2 , 3, 4 then total nseg=9 ) number of compartment or as i defined in code i have 3 compartment due to 3 section?

Imagine a toy model with 3 sections

Code: Select all

create a,b,c
a.nseg = n1
b.nseg = n2
c.nseg = n3

Total number of compartments in toy model is n1+n2+n3
To learn more about nseg see
http://www.neuron.yale.edu/neuron/stati ... .html#nseg
Also go to NEURON's FAQ page
http://www.neuron.yale.edu/neuron/faq/general-questions
and read the answers to these questions
Why should I use an odd value for nseg?
What's a good strategy for specifying nseg?

2- From hh2 code, "vtraub" reperesnt the threshold for spike, right?

No. In this source code, vtraub shifts the voltage-dependence of the voltage-gated
currents, and that will have an effect on the amount of excitation necessary to trigger a
spike, but vtraub is not spike threshold. HH-style models of the action potential do not
have a spike threshold.

which parameter define the maximum value of membrane potential during a action potential ?

There isn't one. The peak amplitude of an action potential depends on the equilibrium
potentials and conductances or permeabilities of all ions that can pass through the
membrane.

[Raj]

By coincidence I was just working on the cells given by Destexhe:
http://senselab.med.yale.edu/senselab/m ... model=3817

The code already contains commented out code which does not use the exact solution which relies on a fixed time step. Removing these comments and commenting out the fixed time step makes the model suitable for use with CVode. (See code below)

For simulations without CVode I think all the mod-files can switch to the recommended numerical integrator for mechanism whose states satisfy a linear differential equation: cnexp.

So in the other mod files

Code: Select all

SOLVE states METHOD euler 

can be replaced by

Code: Select all

SOLVE states METHOD cnexp 

Running the code from the original model with these changes I get very similar results. Using cnexp there are no visible changes. When using CVode there are considerable changes in spike timing, but the behaviour is still qualitatively the same.

Code: Select all

TITLE Hippocampal HH channels
:
: Fast Na+ and K+ currents responsible for action potentials
: Iterative equations
:
: Equations modified by Traub, for Hippocampal Pyramidal cells, in:
: Traub & Miles, Neuronal Networks of the Hippocampus, Cambridge, 1991
:
: range variable vtraub adjust threshold
:
: Written by Alain Destexhe, Salk Institute, Aug 1992
:

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

NEURON {
	SUFFIX hh2:cvode
	USEION na READ ena WRITE ina
	USEION k READ ek WRITE ik
	RANGE gnabar, gkbar, vtraub
	RANGE m_inf, h_inf, n_inf
	RANGE tau_m, tau_h, tau_n
	RANGE m_exp, h_exp, n_exp
}


UNITS {
	(mA) = (milliamp)
	(mV) = (millivolt)
}

PARAMETER {
	gnabar	= .003 	(mho/cm2)
	gkbar	= .005 	(mho/cm2)

	ena	= 50	(mV)
	ek	= -90	(mV)
	celsius = 36    (degC)
	dt              (ms)
	v               (mV)
	vtraub	= -63	(mV)
}

STATE {
	m h n
}

ASSIGNED {
	ina	(mA/cm2)
	ik	(mA/cm2)
	il	(mA/cm2)
	m_inf
	h_inf
	n_inf
	tau_m
	tau_h
	tau_n
	m_exp
	h_exp
	n_exp
	tadj
}


BREAKPOINT {
	SOLVE states METHOD cnexp
	ina = gnabar * m*m*m*h * (v - ena)
	ik  = gkbar * n*n*n*n * (v - ek)
}


DERIVATIVE states {   : exact Hodgkin-Huxley equations
	evaluate_fct(v)
	m' = (m_inf - m) / tau_m
	h' = (h_inf - h) / tau_h
	n' = (n_inf - n) / tau_n
}

:PROCEDURE states() {	: exact when v held constant
:	evaluate_fct(v)
:	m = m + m_exp * (m_inf - m)
:	h = h + h_exp * (h_inf - h)
:	n = n + n_exp * (n_inf - n)
:	VERBATIM
:	return 0;
:	ENDVERBATIM
:}

UNITSOFF
INITIAL {
	m = 0
	h = 0
	n = 0
:
:  Q10 was assumed to be 3 for both currents
:
: original measurements at roomtemperature?

	tadj = 3.0 ^ ((celsius-36)/ 10 )
}

PROCEDURE evaluate_fct(v(mV)) { LOCAL a,b,v2

	v2 = v - vtraub : convert to traub convention

	a = 0.32 * (13-v2) / ( exp((13-v2)/4) - 1)
	b = 0.28 * (v2-40) / ( exp((v2-40)/5) - 1)
	tau_m = 1 / (a + b) / tadj
	m_inf = a / (a + b)

	a = 0.128 * exp((17-v2)/18)
	b = 4 / ( 1 + exp((40-v2)/5) )
	tau_h = 1 / (a + b) / tadj
	h_inf = a / (a + b)

	a = 0.032 * (15-v2) / ( exp((15-v2)/5) - 1)
	b = 0.5 * exp((10-v2)/40)
	tau_n = 1 / (a + b) / tadj
	n_inf = a / (a + b)

:	m_exp = 1 - exp(-dt/tau_m)
:	h_exp = 1 - exp(-dt/tau_h)
:	n_exp = 1 - exp(-dt/tau_n)
}

UNITSON

[ted]

which parameter define the maximum value of membrane potential during a action potential ?

There isn't one. The peak amplitude of an action potential depends on the equilibrium
potentials and conductances or permeabilities of all ions that can pass through the
membrane.

A perturbation analysis of spike peak amplitude will show that it depends most on the
equilibrium potential of the current that is associated with the maximum conductance--
which is ena for sodium spikes. A plot of the peak depolarization vs. ena will be almost
linear with a slope close to 1 for a wide range of ena.

[chinhou]

Hello,

Thank you so much for useful comments.

when i run my code, in voltage graph during current clamp, i see spikes,but the amplitude of first and last spike is a bit bigger , about 5 mv.

However, thanks again.

I have other question about IV curve topics, that i will send on that post.
Hopefully, i get other exact answer.