making .mod files CVODE compatable

[mikey]

I wish to get this model to run faster.

http://senselab.med.yale.edu/senselab/m ... odel=17664

I want to do this by using CVODE. But the .mod files are not appropriate for using with CVODE - get a warning message when compiling them that they are not compatable with CVODE. So, I write to ask how to change them such that they will work with CVODE. There are a lot of .mod files in the model - but our task is simplified by the fact that lots of these channels have the same .mod form. I have cut and pasted a representative .mod mechanism below:

Code: Select all

TITLE Fast sodium current
 
COMMENT
  from "An Active Membrane Model of the Cerebellar Purkinje Cell
        1. Simulation of Current Clamp in Slice"
ENDCOMMENT
 
UNITS {
        (mA) = (milliamp)
        (mV) = (millivolt)
}
 
NEURON {
        SUFFIX NaF
	USEION na WRITE ina
        RANGE  gnabar, gna, minf, hinf, mexp, hexp
} 
 
INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}
 
PARAMETER {
        v (mV)
        celsius = 37 (degC)
        dt (ms)
        gnabar	= 7.5 (mho/cm2)
        ena	= 45 (mV)

}
 
STATE {
        m h
}
 
ASSIGNED {
        ina (mA/cm2)
        gna minf hinf mexp hexp 
}
 
BREAKPOINT {
        SOLVE states
        gna = gnabar *m*m* m*h 
	ina = gna* (v-ena)
}
 
UNITSOFF
 
INITIAL {
	rates(v)
	m = minf
	h = hinf
}

PROCEDURE states() {  :Computes state variables m, h
        rates(v)      :             at the current v and dt.
        m = m + mexp*(minf-m)
        h = h + hexp*(hinf-h)
}
 
PROCEDURE rates(v) {  :Computes rate and other constants at current v.
                      :Call once from HOC to initialize inf at resting v.
        LOCAL  q10, tinc, alpha, beta, sum
        TABLE minf, mexp, hinf, hexp DEPEND dt, celsius FROM -100 TO 100 WITH 200
        q10 = 3^((celsius - 37)/10)
        tinc = -dt * q10
                :"m" sodium activation system
        alpha = 35/exp((v+5)/(-10))
        beta =  7/exp((v+65)/20)
        sum = alpha + beta
        minf = alpha/sum
        mexp = 1 - exp(tinc*sum)
                :"h" sodium inactivation system
        alpha = 0.225/(1+exp((v+80)/10))
        beta = 7.5/exp((v-3)/(-18))
        sum = alpha + beta
        hinf = alpha/sum
        hexp = 1 - exp(tinc*sum)
}

 
UNITSON

------------

The channels that do not have the above form are those that are calcium regulated as well as voltage regulated. I have cut and paste a representative of this class below:

Code: Select all

 
TITLE BK calcium-activated potassium current
: Calcium activated K channel.
COMMENT
  from "An Active Membrane Model of the Cerebellar Purkinje Cell
        1. Simulation of Current Clamp in Slice"
ENDCOMMENT

UNITS {
	(molar) = (1/liter)
}

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


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

NEURON {
	SUFFIX KC3
	USEION ca READ cai
	USEION k WRITE ik
	RANGE gkbar,gk,zinf,ik
}


PARAMETER {
	celsius=37	(degC)
	v		(mV)
	gkbar=.08	(mho/cm2)	: Maximum Permeability
	cai = .04e-3	(mM)
	ek  = -85	(mV)
	dt		(ms)
}


ASSIGNED {
	ik		(mA/cm2)
	minf
	mexp
	zinf
	zexp
	gk
}

STATE {	m z }		: fraction of open channels

BREAKPOINT {
	SOLVE state
:	gk = gkbar*m*z*z
	ik = gkbar*m*z*z*(v - ek)
}
:UNITSOFF
:LOCAL fac

:if state_cagk is called from hoc, garbage or segmentation violation will
:result because range variables won't have correct pointer.  This is because
: only BREAKPOINT sets up the correct pointers to range variables.
PROCEDURE state() {	: exact when v held constant; integrates over dt step
	rate(v, cai)
	m = m + mexp*(minf - m)
	z = z + zexp*(zinf - z)
	VERBATIM
	return 0;
	ENDVERBATIM
}

INITIAL {
	rate(v, cai)
	m = minf
	z = zinf
}

FUNCTION alp(v (mV), ca (mM)) (1/ms) { :callable from hoc
	alp = 400/(ca*1000)
}

FUNCTION bet(v (mV)) (1/ms) { :callable from hoc
	bet = 0.11/exp((v-35)/14.9)
}

PROCEDURE rate(v (mV), ca (mM)) { :callable from hoc
	LOCAL a,b
	a = alp(v,ca)
	zinf = 1/(1+a)
	zexp = (1 - exp(-dt/10))
	b = bet(v)
	minf = 7.5/(7.5+b)
	mexp = (1 - exp(-dt*(7.5+b)))
}
:UNITSON

Would be so so grateful to anyone that can help out.

[ted]

Tell you what: I'll answer your question now. After you read my answer, you'll edit your
original msg so that code formatting is preserved, as described here
https://www.neuron.yale.edu/phpBB2/viewtopic.php?t=493
To help you get going, I already fixed the first big block of code in your original msg.
Sounds like a deal?

In HH-style models of channel gating, the gating states are described by 1st order DEs.
An old trick for speeding up simulations avoided solving these DEs numerically by
treating the voltage-dependent rate constants (alphas and betas) as if they were
constant during each time step, with values that correspond to membrane potential at
the start of the time step. This allowed one to compute the future value of each gating
variable from its present value and its time constant and steady state value at the
present membrane potential.

That's what these mod files are doing. The way to make them compatible with cvode is
to rewrite them so that the gating variables are once again described by DEs. For
example, here's how to fix the description of the NaF mechanism (omitting most of
the stuff that stays unchanged). Stuff that is relevant to cvode is flagged by upper case CVODE

Code: Select all

NEURON {
 . . .
: CVODE--eliminate unused mexp, hexp, add mtau, htau
:        RANGE  gnabar, gna, minf, hinf, mexp, hexp
        RANGE  gnabar, gna, minf, hinf, mtau, htau
}
 . . .
PARAMETER {
 . . .
: CVODE--eliminate dt
:        dt (ms)
 . . . 
}
 . . .
ASSIGNED {
 . . .
: CVODE--eliminate mexp, hexp
:	mexp hexp
: declare mtau, htau
	mtau (ms)
	htau (ms)
}


BREAKPOINT {
: CVODE--replace PROCEDURE states with a DERIVATIVE block
:        SOLVE states
: cnexp method for DERIVATIVE block with linear DEs
        SOLVE states METHOD cnexp
 . . .
}
 . . .
: CVODE--replace this PROCEDURE with a DERIVATIVE block
: PROCEDURE states() {  :Computes state variables m, h
:        rates(v)      :             at the current v and dt.
:        m = m + mexp*(minf-m)
:        h = h + hexp*(hinf-h)
: }
DERIVATIVE states {
        rates(v)
        m' = (minf-m)/mtau
        h' = (hinf-h)/htau
}

PROCEDURE rates(v) {  :Computes rate and other constants at current v.
                      :Call once from HOC to initialize inf at resting v.
: CVODE--eliminate tinc
:        LOCAL  q10, tinc, alpha, beta, sum
        LOCAL  q10, alpha, beta, sum
: CVODE--eliminate mexp, hexp, dt, insert mtau, htau
:        TABLE minf, mexp, hinf, hexp DEPEND dt, celsius FROM -100 TO 100 WITH 200
        TABLE minf, mtau, hinf, htau DEPEND celsius FROM -100 TO 100 WITH 200
        q10 = 3^((celsius - 37)/10)
: CVODE
:        tinc = -dt * q10
                :"m" sodium activation system
        alpha = 35/exp((v+5)/(-10))
        beta =  7/exp((v+65)/20)
        sum = alpha + beta
        minf = alpha/sum
: CVODE
:        mexp = 1 - exp(tinc*sum)
	mtau = 1/(sum*q10)
                :"h" sodium inactivation system
        alpha = 0.225/(1+exp((v+80)/10))
        beta = 7.5/exp((v-3)/(-18))
        sum = alpha + beta
        hinf = alpha/sum
: CVODE
:        hexp = 1 - exp(tinc*sum)
	htau = 1/(sum*q10)
}

Now you know how to deal with the code for KC3.

A footnote: there is no guarantee that using cvode will accelerate simulations. Depending
on the error tolerance that you choose, and the detailed kinetics of the mechanisms that
are in the model, simulations may even execute more slowly, because the integrator is
quite fastidious about controlling error and handling every event at the time it occurs.
However, for a given execution time, cvode simulations with probably be more accurate
than fixed dt simulations.

Further comments about this mod file:

Note that the USEION statement omits READ ena. This has several consequences:
ena, which is declared in the PARAMETER block, is truly a parameter of this
mechanism. That is, its value is independent of the hoc-level ena (which WILL affect
almost every other HH-style na channel mechanism that you will ever encounter).
Its value will also be independent of intra- and extracellular sodium concentration.

v and celsius are declared in the PARAMETER block, but these are not really
parameters of this mechanism. These variables will get their values from
outside this mechanism, and they should really be declared in the ASSIGNED block.
The assignment statement celsius = 37, no matter where it occurs in this mod file,
will have no effect; if you want celsius to be anything other than 6.3, you need to
execute an assignment statement at the hoc level.

The units of gna should be declared (in the ASSIGNED block).

[mikey]

Thanks so much Ted for your kind reply. I have edited my initial post with the

Code: Select all

 tags. Using your worked out example I have been able to edit practially all the model .mod files to make them CVODE compatable. But a few remain (3) that I have tried and failed to adjust. The reason for all of these is that they are just so different in the PROCEDURE rates(v) block from your worked example that I cannot figure out how to manipulate them in the correct fashion. I have put the original code for these below (using [code] tags to keep it all clean). Although I have cut and pase the full code for all of these - it is only the PROCEDURE rates(v) block in each that I have the problem with. I cannot see how I can effectively replace the mexp entity with mtau (or hexp for htau, or zexp for ztau) in these cases. For instance, in one I am confused because no alpha and beta equations are displayed, whereas the worked example did make use of such equations. 

EDIT: I am having real problem with getting the code to format properly in this post for some reason. I have cut and paste the code to a later post.

[mikey]

Dear Ted, I tried to get the code posted in the proper fashion using the

Code: Select all

 and 

tags - in fact you can see the tags in my post, as they haven't worked in formatting the code in the nice white background. I dont understand why not as I have BBcode ON in my settings. And I did not tick it off or anything like that when I poseted.

EDIT - I have solved this by just putting all the code in a later post - it has formatted good in this later post. Please see below.

[mikey]

------------
1
------------

Code: Select all

TITLE Anomalous rectifier
 
COMMENT
  from "An Active Membrane Model of the Cerebellar Purkinje Cell
        1. Simulation of Current Clamp in Slice"
  
  "Anomalous Rectification in Neurons from Cat Sensorimotor Cortex
    In Vitro" W.J.SPAIN, P.C.SCHWINDT, and W.E.CRILL
           JOURNAL OF NEUROPHYSIOLOGY Vol 57, No.5
ENDCOMMENT
 
UNITS {
        (mA) = (milliamp)
        (mV) = (millivolt)
}
 
NEURON {
        SUFFIX Kh
	USEION k WRITE ik
        RANGE  gkbar, gk, minf, mexp, nexp, ik
} 
 
INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}
 
PARAMETER {
        v (mV)
        celsius = 37 (degC)
        dt (ms)
        gkbar	= .0003 (mho/cm2)
        ek	= -30 (mV)

}
 
STATE {
        m 
}
 
ASSIGNED {
        ik (mA/cm2)
        gk minf  mexp nexp
}
 
BREAKPOINT {
        SOLVE states
        gk = gkbar *m
	ik = gk* (v-ek)
}
 
UNITSOFF
 
INITIAL {
	rates(v)
	m = minf
}

PROCEDURE states() {  :Computes state variables m, 
        rates(v)      :             at the current v and dt.
        m = 0.8*(m + mexp*(minf-m))
		+0.2*(m + nexp*(minf-m))
}
 
PROCEDURE rates(v) {  :Computes rate and other constants at current v.
                      :Call once from HOC to initialize inf at resting v.
        LOCAL  q10, tinc
        TABLE minf, mexp,nexp DEPEND dt, celsius FROM -100 TO 100 WITH 200
        q10 = 3^((celsius - 37)/10)
        tinc = -dt * q10
                :"m" potassium activation system
        minf = 1/(1+exp((v+78)/7))
        mexp = 1 - exp(tinc/38)
        nexp = 1 - exp(tinc/319)
}
UNITSON

------------
2
------------

Code: Select all

TITLE Delayed rectifire
 
COMMENT
  from "An Active Membrane Model of the Cerebellar Purkinje Cell
        1. Simulation of Current Clamp in Slice"
ENDCOMMENT
 
UNITS {
        (mA) = (milliamp)
        (mV) = (millivolt)
}
 
NEURON {
        SUFFIX Kdr
	USEION k WRITE ik
        RANGE  gkbar, gk, minf, hinf, mexp, hexp, ik, alpha, beta
} 
 
INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}
 
PARAMETER {
        v (mV)
        celsius = 37 (degC)
        dt (ms)
        gkbar	= .6 (mho/cm2)
        ek	= -85 (mV)

}
 
STATE {
        m h
}
 
ASSIGNED {
        ik (mA/cm2)
        gk minf hinf mexp hexp
}
 
BREAKPOINT {
        SOLVE states
        gk = gkbar *m*m*h
	ik = gk* (v-ek)
}
 
UNITSOFF
 
INITIAL {
	rates(v)
	m = minf
	h = hinf
}

PROCEDURE states() {  :Computes state variables m,h
        rates(v)      :             at the current v and dt.
        m = m + mexp*(minf-m)
	h = h + hexp*(hinf-h)
}

PROCEDURE rates(v) {  :Computes rate and other constants at current v.
                      :Call once from HOC to initialize inf at resting v.
        LOCAL  q10, tinc, tauh, alpha, beta, gamma, zeta, taum
:        TABLE minf, mexp, hinf, hexp DEPEND dt, celsius FROM -100 TO 100 WITH 2:00
        q10 = 3^((celsius - 37)/10)
        tinc = -dt * q10
                :"m" potassium activation system
        alpha = -0.0047*(v-8)/(exp((v-8)/(-12))-0.9999)
:		if(v == 8) {v = 8.0001}
	beta = exp((v+127)/(-30))
	minf = alpha/(alpha+beta)
        gamma = -0.0047*(v+12)/(exp((v+12)/(-12))-0.9999)
	zeta = exp((v+147)/(-30))
	taum = 1/(gamma + zeta)
	mexp = 1 - exp(tinc/taum)
                :"h" potassium activation system
        hinf = 1.0 / (1+exp((v+25)/4))
        if(v<-25) {
	tauh = 1200
	}else{
	tauh = 10
	}
	hexp = 1 - exp(tinc/tauh)
       
}

 
UNITSON

------------
3
------------

Code: Select all

TITLE K2 calcium-activated potassium current
: Calcium activated K channel.
COMMENT
  from "An Active Membrane Model of the Cerebellar Purkinje Cell
        1. Simulation of Current Clamp in Slice"
ENDCOMMENT

UNITS {
	(molar) = (1/liter)
}

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


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

NEURON {
	SUFFIX K23
	USEION ca READ cai
	USEION k WRITE ik
	RANGE gkbar,gk,zinf,ik
}


PARAMETER {
	celsius=37	(degC)
	v		(mV)
	gkbar=.00039	(mho/cm2)	: Maximum Permeability
	cai = .04e-3	(mM)
	ek  = -85	(mV)
	dt		(ms)
}


ASSIGNED {
	ik		(mA/cm2)
	minf
	mexp
	zinf
	zexp
	gk
}

STATE {	m z }		: fraction of open channels

BREAKPOINT {
	SOLVE state
:	gk = gkbar*m*z*z
	ik = gkbar*m*z*z*(v - ek)
}
:UNITSOFF
:LOCAL fac

:if state_cagk is called from hoc, garbage or segmentation violation will
:result because range variables won't have correct pointer.  This is because
: only BREAKPOINT sets up the correct pointers to range variables.
PROCEDURE state() {	: exact when v held constant; integrates over dt step
	rate(v, cai)
	m = m + mexp*(minf - m)
	z = z + zexp*(zinf - z)
	VERBATIM
	return 0;
	ENDVERBATIM
}

INITIAL {
	rate(v, cai)
	m = minf
	z = zinf
}

FUNCTION alp(v (mV), ca (mM)) (1/ms) { :callable from hoc
	alp = 20/(ca*1000)
}

FUNCTION bet(v (mV)) (1/ms) { :callable from hoc
	bet = 0.075/exp((v+5)/10)
}

PROCEDURE rate(v (mV), ca (mM)) { :callable from hoc
	LOCAL a,b
	a = alp(v,ca)
	zinf = 1/(1+a)
	zexp = (1 - exp(-dt/10))
	b = bet(v)
	minf = 25/(25+b)
	mexp = (1 - exp(-dt*(25+b)))
}
:UNITSON

[ted]

Thanks for using the code formatting tags, which really do improve readability. They also
allow users to select and copy code, with formatting intact, into their own text editors.

The Kh mechanism is too strange to start with. Let's deal with Kdr first.

Code: Select all

PROCEDURE states() {  :Computes state variables m,h
        rates(v)      :             at the current v and dt.
        m = m + mexp*(minf-m)
   h = h + hexp*(hinf-h)
}

The derivative forms for m and h are just
m' = (minf - m) / taum
h' = (hinf - h) / tauh
PROCEDURE rates calculates minf, taum, hinf, and tauh, so no
problem there. Just comment out anything that has to do with dt, tinc,
mexp and hexp. Incorporate temperature effects by dividing the time
constants by q10.

As for K23,

Code: Select all

PROCEDURE state() {   : exact when v held constant; integrates over dt step
   rate(v, cai)
   m = m + mexp*(minf - m)
   z = z + zexp*(zinf - z)
   VERBATIM
   return 0;
   ENDVERBATIM
}

forget about that VERBATIM . . . ENDVERBATIM stuff. That's irrelevant
to the DERIVATIVE block you're going to write. By now you can write the
derivative forms for m and z "by inspection."

"But where are the time constants for z and m?" you may ask. Look at PROCEDURE rate.

Code: Select all

zexp = (1 - exp(-dt/10))
                    ^^ that's tauz
mexp = (1 - exp(-dt*(25+b)))
                     ^^^^ that's 1/taum

Finally we turn to the Kh mechanism.

I had to go back to ModelDB just to make sure PROCEDURE states in
your post was the same as the authors' original, because it looks strange.

Code: Select all

PROCEDURE states() {  :Computes state variables m,
        rates(v)      :             at the current v and dt.
        m = 0.8*(m + mexp*(minf-m)) + 0.2*(m + nexp*(minf-m))
}
where mexp = 1 - exp(tinc/38)
and nexp = 1 - exp(tinc/319)

But the expression for m is just the analytic solution of this differential
equation:
m' = (minf - m)*(0.8/tau1 + 0.2/tau2)
where tau1 is 38/q10 and tau2 is 319/q10
so the BREAKPOINT block is

Code: Select all

BREAKPOINT {
  SOLVE states METHOD cnexp
  gk = gkbar *m
  ik = gk* (v-ek)
}

and the DERIVATIVE block is

Code: Select all

DERIVATIVE states {  :Computes state variables m, 
  rates(v)      :             at the current v and dt.
  m' = (minf-m)*(0.8/tau1 + 0.2/tau2)
}

Any doubts about this are eradicated by comparing simulations of the
cvode-compatible and fixed dt implementations side by side, i.e. insert
them both into a model cell, subject it to voltage clamp, and plot the
values of m and gk. I'll send you the mod files plus a hoc and ses file that
demonstrate this point.

[mikey]

RE: "I'll send you the mod files plus a hoc and ses file that
demonstrate this point."

Thanks so much Ted. I'm guessing you'll send them to the email account in my profile? I will look out for them. Once again, thanks so, so much. I just hope after all this, the CVODE modification will give us a bit of improved performance rate wise.

[mikey]

When compiling K23.mod I get this warning message:

"Warning: rates undefined. (declared within VERBATIM?)"

Is this any reason to worry? I have cut and paste my modified K23.mod below:

Code: Select all

TITLE K2 calcium-activated potassium current
: Calcium activated K channel.
COMMENT
  from "An Active Membrane Model of the Cerebellar Purkinje Cell
        1. Simulation of Current Clamp in Slice"
ENDCOMMENT

UNITS {
	(molar) = (1/liter)
}

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


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

NEURON {
	SUFFIX K23
	USEION ca READ cai
	USEION k WRITE ik
	RANGE gkbar,gk,zinf,ik
}


PARAMETER {
	celsius=37	(degC)
	v		(mV)
	gkbar=.00039	(mho/cm2)	: Maximum Permeability
	cai = .04e-3	(mM)
	ek  = -85	(mV)
}


ASSIGNED {
	ik		(mA/cm2)
	minf
	zinf
	gk
        mtau (ms)
        ztau (ms)
}

STATE {	m z }		: fraction of open channels

BREAKPOINT {
	SOLVE states METHOD cnexp 
:	gk = gkbar*m*z*z
	ik = gkbar*m*z*z*(v - ek)
}
:UNITSOFF
:LOCAL fac


DERIVATIVE states { 
        rates(v, cai) 
        m' = (minf-m)/mtau 
        z' = (zinf-z)/ztau 
} 


INITIAL {
	rate(v, cai)
	m = minf
	z = zinf
}

FUNCTION alp(v (mV), ca (mM)) (1/ms) { :callable from hoc
	alp = 20/(ca*1000)
}

FUNCTION bet(v (mV)) (1/ms) { :callable from hoc
	bet = 0.075/exp((v+5)/10)
}

PROCEDURE rate(v (mV), ca (mM)) { :callable from hoc
	LOCAL a,b
	a = alp(v,ca)
	zinf = 1/(1+a)
	ztau = 10
	b = bet(v)
	minf = 25/(25+b)
	mtau = 1/(25+b)
}
:UNITSON

[mikey]

The problem I have now is that all my .mod files compile ok (apart from K23.mod above - has that warning message). But when I actually run them in a simulation it doesn't work. The thing breaks down at the first inserted mechanim (which is NAF.mod) in the .hoc file. Get this error message:

syntax error in purkinje.hoc near line 290
insert NaF gnabar_NaF = 10
^

I have cut and paste the code for NaF.mod below. Would be so grateful if you could take a look and see where I have perhaps gone wrong.

Code: Select all

TITLE Fast sodium current
 
COMMENT
  from "An Active Membrane Model of the Cerebellar Purkinje Cell
        1. Simulation of Current Clamp in Slice"
ENDCOMMENT
 
UNITS {
        (mA) = (milliamp)
        (mV) = (millivolt)
}
 
NEURON {
        SUFFIX NaF
	USEION na WRITE ina
        RANGE  gnabar, gna, minf, hinf, mtau, htau
} 
 
INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}
 
PARAMETER {
        v (mV)
        celsius = 37 (degC)
        gnabar	= 7.5 (mho/cm2)
        ena	= 45 (mV)

}
 
STATE {
        m h
}
 
ASSIGNED {
        ina (mA/cm2)
        gna minf hinf 
        mtau (ms)
        htau (ms)
}
 
BREAKPOINT {
        SOLVE states METHOD cnexp
        gna = gnabar *m*m* m*h 
	ina = gna* (v-ena)
}
 
UNITSOFF
 
INITIAL {
	rates(v)
	m = minf
	h = hinf
}

DERIVATIVE states { 
        rates(v) 
        m' = (minf-m)/mtau 
        h' = (hinf-h)/htau 
} 
 
PROCEDURE rates(v) {  :Computes rate and other constants at current v.
                      :Call once from HOC to initialize inf at resting v.
        LOCAL  q10, alpha, beta, sum
        TABLE minf, mtau, hinf, htau DEPEND celsius FROM -100 TO 100 WITH 200 
        q10 = 3^((celsius - 37)/10)
                :"m" sodium activation system
        alpha = 35/exp((v+5)/(-10))
        beta =  7/exp((v+65)/20)
        sum = alpha + beta
        minf = alpha/sum
       mtau = 1/(sum*q10) 
                :"h" sodium inactivation system
        alpha = 0.225/(1+exp((v+80)/10))
        beta = 7.5/exp((v-3)/(-18))
        sum = alpha + beta
        hinf = alpha/sum
       htau = 1/(sum*q10)
}

 
UNITSON

[ted]

When compiling K23.mod I get this warning message:

"Warning: rates undefined. (declared within VERBATIM?)"

Is this any reason to worry?

Yes. It means that you share the common tendency to making typograpical errors and then
not being able to find them because your brain knows what it thinks is on the printed page,
so that's what you perceive to be there. As a rule, when you see an error message that
says "x is undefined" it means you should use your text editor to look for occurrences of
"x" in your code.

[ted]

[quote="mikey"]The problem I have now is that all my .mod files compile ok (apart from K23.mod above - has that warning message). But when I actually run them in a simulation it doesn't work. The thing breaks down at the first inserted mechanim (which is NAF.mod) in the .hoc file. Get this error message:

syntax error in purkinje.hoc near line 290
insert NaF gnabar_NaF = 10
^{/quote]
The oriiginal hoc syntax error message probably looked like this in NEURON's xterm:

Code: Select all

syntax error in purkinje.hoc near line 290
 insert NaF gnabar_NaF = 10
                     ^

hoc's syntax error messages tend to be quite informative because they usually point to
the offending statement. The problem here is that you have put two statements on one
line. The rule is that you can only do this if the line lies within a compound statement, i.e.
within a pair of curly brackets. That is, the parser would not complain about this
{ insert NaF gnabar_NaF = 10 }

BUT this is a bad practice because code written like that is hard to read, therefore
hard to debug. Only experts who can write hoc with half their brains tied behind their
backs should risk writing multiple statements on a single line.

[mikey]

THE OUTCOME: I have been debugging. The simulation runs all ok with all my new .mod files, except using the original versions for k23.mod and kC3.mod. Thus, the problems must be in these 2 .mod files. Now - these 2 .mod files can be treated as one and the same for our purposes as they both have the same form. Why are they causing problems? Well, we get a warning when we compile them (although they do compile ok) - this warning may provide a clue:

"Warning: rates undefined. (declared within VERBATIM?)"

It was just a typo with the "rates" statement. Some of them were "rate" and others were "rates". So, solved it now. And the whole simulation works with CVODE!! That is the good news. The bad news, is that it doesn't bring much increase in speed. But that was maybe expected given that the simulation does fast spiking behaviour. I will try playing with the error value - see if can increase the speed.

[ted]

the whole simulation works with CVODE!! That is the good news. The bad news, is that it doesn't bring much increase in speed.

Yep. See this

A footnote: there is no guarantee that using cvode will accelerate simulations.

in my reply to your first message in this thread.

You are quite right to suspect that fast spiking could be the principal cause of the lack of
simulation speedup. The presence of a calcium accumulation mechanism will also slow
simulations.

I will try playing with the error value - see if can increase the speed.

Your best bet is to use the GUI's tool for automatically adjusting the error tolerances of
individual variables.
1. NEURON Main Menu / Tools / VariableStepControl
2. Make sure that the VariableTimeStep panel's "Use variable dt" box is checked, then
click on the VariableTimeStep panel's "Atol Scale Tool" button.
3. In the Absolute Tolerance Scale Factors panel, click on the "Analysis Run" button.
Wait for the simulation to finish, then click on that panel's "Rescale" button.
4. Finally, try a new simulation with the RunControl panel's Init & Run button (or by
executing run() ).
5. You will want to click on the VariableTimeStep panel's "Hints" button and read the hints.