calcium dependent inactivation function

[GTR]

Dear Ted,
I am trying to make a CaL.mod file of a high voltage activated calcium channel CaL with the current equation of the form:
ilca=glcabar*q(v)*h(cai)*ghk with q and h the activation and inactivation functions respectively.
In order to calculate the intracellular calcium concentration (cai) i have a ccanl.mod (like this in Migliore et al. or Santhakumar et al. of ModelDB) that calculates cai concentration:

Code: Select all

NEURON {
	SUFFIX ccanl
USEION nca READ ncai, inca, enca WRITE enca, ncai VALENCE 2
USEION lca READ lcai, ilca, elca WRITE elca, lcai VALENCE 2
USEION tca READ tcai, itca, etca WRITE etca, tcai VALENCE 2
RANGE caiinf, catau, cai, ncai, lcai,tcai, eca, elca, enca, etca
}
...................................................

so my first problem is how can i access the calculated cai concentration to my CaL mechanism?

I read the topic about communication between different NMODL mechanisms so i tried to implement the useion method but unfortunately I didn't make it.So what to add in the two .mod files?And how this cai concentration I need to
access should be stated in the CaL mechanism ?in parameter section?as a RANGE variable?

and one more question:
so in my CaL.mod the activation function is voltage-dependent q(v) and the
inactivation function is calcium-dependent h(cai) as well as the steady-state hinf(cai) and the time constant htau(cai) so my problem is how and where to declare them?
I imagine sth like this in this part of my CaL code:

Code: Select all

.........................................................................
BREAKPOINT {
    SOLVE states METHOD cnexp
    glca=glcabar*q*h(cai)
    ilca = glca*ghk(v,cai,cao)

DERIVATIVE state {	
	rate(cai(mM))
h'=(hinf-h)/htau
}

DERIVATIVE state {	: exact when v held constant
	rate(v(mV))
q'=(qinf-q)/qtau
}
............................................................
UNITSOFF
PROCEDURE rate(v(mV)) 
        qinf=1/(1+exp((-24.6-v)/11.3))
        qtau=1.25/cosh(-0.03*(v+37.1))	
}
PROCEDURE rate(cai(mM)) 
        hinf=0.47/(1+exp(cai)-0.7)
        htau=1220	
}
UNITSON

can I treat this h(cai) as a function?
Are you aware of any model that might help me?

Excuse my boring description but I want to make my problems clear to you.
regards,
George.

[ted]

I am trying to make a CaL.mod file of a high voltage activated calcium channel CaL with the current equation of the form:
ilca=glcabar*q(v)*h(cai)*ghk with q and h the activation and inactivation functions respectively.
In order to calculate the intracellular calcium concentration (cai) i have a ccanl.mod (like this in Migliore et al. or Santhakumar et al. of ModelDB) that calculates cai concentration:

Code: Select all

NEURON {
	SUFFIX ccanl
USEION nca READ ncai, inca, enca WRITE enca, ncai VALENCE 2
USEION lca READ lcai, ilca, elca WRITE elca, lcai VALENCE 2
USEION tca READ tcai, itca, etca WRITE etca, tcai VALENCE 2
RANGE caiinf, catau, cai, ncai, lcai,tcai, eca, elca, enca, etca
}

so my first problem is how can i access the calculated cai concentration to my CaL mechanism?

I presume that the model by Migliore et al., to which you refer, is their parallel network
remake of the model by Santhakumar et al. (2005). So let's focus on the latter.

Several of the mod files supplied by Santhakumar et al. contain complexities whose
utility escapes me. For example, the Ca accumulation mechanism, which is specified
in ccanl.mod, fractionates calcium current into three different components inca, ilca,
and itca, each with its own corresponding concentration and equilibrium potential
(ncai, enca, lcai, elca, tcai, etca). It integrates each of the three current components
separately, then pools the three resulting concentrations into cai, computes eca from
cai, and finally assigns the value of eca to enca, elca, and etca. However, there is
no mechanism that WRITEs cai, so there is no cai to be accessed from hoc.

As far as I can tell, all these manipulations are pointless. Neither the NMODL nor the
hoc code makes any use of this decomposition of Ca currents and concentrations.
Maybe they're just reusing code from another model in which this complexity played
some useful role, but in this particular model it is baroque.

I read the topic about communication between different NMODL mechanisms

Don't even bother. Your best bet is to eliminate the problem at its source, i.e. write your
own model descriptions so that you don't have to wrestle with this balkanization of
calcium currents and concentrations.

The place to start is with the calcium accumulation mechanism. This should simply
USEION ca READ ica, cai WRITE cai

All calcium currents should
USEION ca READ eca, cai WRITE ica
You already know how to use
RANGE i
in the NEURON block, and
i = g*(v-e)
ica = i
in the BREAKPOINT block if you want to be able to access the individual current
components from hoc.

Code: Select all

BREAKPOINT {
    SOLVE states METHOD cnexp
    glca=glcabar*q*h(cai)
    ilca = glca*ghk(v,cai,cao)
}

Presumably your L current mechanism will have a suffix like lca, so there is no reason
to call its density or net conductance glcabar and glca--why have hiccupping variable
names glcabar_lca and glca_lca at the hoc level, when gbar and g in your mod file will
let you use the perfectly mnemonic gbar_lca and g_lca in hoc?

Put all of the ODEs into a single DERIVATIVE block, and do all of the rate constant
evaluations in a single PROCEDURE, e.g.

Code: Select all

DERIVATIVE state {	
	rate(cai(mM), v(mV))
h'=(hinf-h)/htau
q'=(qinf-q)/qtau
}

UNITSOFF
PROCEDURE rate(v(mV), cai(mM)) {
        qinf=1/(1+exp((-24.6-v)/11.3))
        qtau=1.25/cosh(-0.03*(v+37.1))	
        hinf=0.47/(1+exp(cai)-0.7)
        htau=1220	
}
UNITSON

You might find it useful to look at this model by Lazarewicz et al. (2002)--
http://senselab.med.yale.edu/SenseLab/M ... odel=20007

[GTR]

Dear Ted,
Unfortunately I have to wrestle with all these currents and concentrations (CaL,CaN,CaT) because the subthalamic neuron model specification I want to implement says that :
Intracellular calcium levels were modelled in a sub-membrane shell with buffering and diffusion modelled as an exponential decay.
With the equation (that includes the three currents) of the form:

[cai] ' =(Ilca+Inca+Itca)*c*(cai0-[cai])/tca

c :is the conversion constant
tca:the time constant of the decay
cai0:the basal intracellular calcium level

this is the reason I want to use ccanl.mod file from which I want to calculate cai concentration for my calcium .mod files so to start with the calcium accumulation mechanism I added to Neuron block:

Code: Select all

NEURON {
	SUFFIX ccanl
USEION nca READ ncai, inca, enca WRITE enca, ncai VALENCE 2
USEION lca READ lcai, ilca, elca WRITE elca, lcai VALENCE 2
USEION tca READ tcai, itca, etca WRITE etca, tcai VALENCE 2
RANGE caiinf, catau, cai, ncai, lcai,tcai, eca, elca, enca, etca}

USEION ca READ cai,ica WRITE cai VALENCE 2

and in the DERATIVE block:

Code: Select all

DERIVATIVE integrate {
ncai' = -(inca)/depth/FARADAY * (1e7) + (caiinf/3 - ncai)/catau
lcai' = -(ilca)/depth/FARADAY * (1e7) + (caiinf/3 - lcai)/catau
tcai' = -(itca)/depth/FARADAY * (1e7) + (caiinf/3 - tcai)/catau

cai' = ncai'+lcai'+tcai'

My problem is that in the CaL.mod I'm still getting compilation error:
nocmodl CaL
Translating CaL.mod into CaL.c
cai used as both variable and function at line 94 in file CaL.mod
^
make: *** [CaL.o] Error 1

There was an error in the process of creating nrnmech.dll

Code: Select all

TITLE CaL High-voltage activated calcium channel 
: L-Type Calcium Channel with Goldman- Hodgkin-Katz permeability GHK

UNITS {
    (mV) = (millivolt)
    (mA) = (milliamp)
}
NEURON {
    SUFFIX CaL
    USEION ca READ cai,cao WRITE ica VALENCE 2
    RANGE gbar,g, ica 
    GLOBAL qinf,qtau
}
UNITS {
	:FARADAY = 96520 (coul)
	FARADAY = (faraday) (coulomb)
        :R = 8.3134 (joule/degC)
	R = (k-mole) (joule/degC)
        (molar) = (1/liter)
	(mM) = (millimolar)
        (mS)=(millisiemens)

}

PARAMETER {
    v(mV)
    gbar  = 9.50e-4 (mS/cm2)
    cai	(mM)
    cao	(mM)
celsius (degC)

 
}

ASSIGNED {
    
    ica(mA/cm2)
    qinf
    qtau(ms)
    g(mS)

    
}

STATE {
    q h
}
INITIAL {
        rate(v)
	q=qinf
         
        
}

BREAKPOINT {
    SOLVE state METHOD cnexp
    g=gbar*q(v)*q(v)*h(cai)
    ica = g*ghk(v,cai,cao)
}

DERIVATIVE state {    
   rate(cai(mM), v(mV)) 
h'=(hinf-h)/htau 
q'=(qinf-q)/qtau
} 

FUNCTION ghk(v(mV), ci(mM), co(mM)) (.001 coul/cm3)
  {
	LOCAL z, eci, eco
	z = (1e-3)*2*FARADAY*v/(R*(celsius+273.15))
	eco = co*efun(z)
	eci = ci*efun(-z)
	:high cao charge moves inward
	:negative potential charge moves inward
	ghk = (.001)*2*FARADAY*(eci - eco)
}
FUNCTION efun(z) {
	if (fabs(z) < 1e-4) {
		efun = 1 - z/2
	}else{
		efun = z/(exp(z) - 1)
	}
}

UNITSOFF
 
PROCEDURE rate (cai(mM),v(mV)) { 
        qinf=1/(1+exp((-24.6-v)/11.3)) 
        qtau=1.25/cosh(-0.03*(v+37.1))    
        hinf=0.53+0.47/(1+exp((cai-0.7)/0.15))
        htau=1220 
           
} 
UNITSON

Why does it indicate the problem to be in the last line of my code(line94)?

thank you
George.

[ted]

Intracellular calcium levels were modelled in a sub-membrane shell with buffering and diffusion modelled as an exponential decay.
With the equation (that includes the three currents) of the form:

[cai] ' =(Ilca+Inca+Itca)*c*(cai0-[cai])/tca

c :is the conversion constant
tca:the time constant of the decay
cai0:the basal intracellular calcium level

Are you quite sure about this equation? It bears no resemblence to a
mass balance equation, or to the equations in the NMODL code you posted.
Does it contain a typographical error, or is such an error in the source that
you cite?

[GTR]

I am sending you via e-mail the pdf file that contains the model I want to implement in NEURON to see how it specifies intracellular calcium concentration.

[ted]

No need to send the PDF. The citation will be sufficient. I probably already have the
paper, and if not, I will almost certainly be able to retrieve it for myself. I really prefer not
to have my mailbox fill up with emails that have giant attachments.

[GTR]

I was going to send you only 4 pages about the model specification.Anyway the paper is from Andrew Gillies and David Willshaw Membrane channel interactions underlying rat subthalamic projection neuron rhythmic and bursting activity J.Neurophysiology(2005).The equation we talk about is on page 31 .

I appreciate your help,
George.

[ted]

Andrew Gillies and David Willshaw Membrane channel interactions underlying rat subthalamic projection neuron rhythmic and bursting activity J.Neurophysiology(2005)

You mean this one?
Andrew Gillies and David Willshaw
Membrane Channel Interactions Underlying Rat Subthalamic Projection Neuron
Rhythmic and Bursting Activity
J Neurophysiol, Apr 2006; 95: 2352 - 2365

In which case, you are referring to Eq. A11 at the bottom of page 2364. This looks like
just another typographical error--there should be a + between the ) and the c. Think
about it: the equation, as it stands, dictates that

• 1. if cai is at its resting level cai0, then calcium current, no matter how large, will have
  NO effect on calcium concentration
  2. if the total calcium current is 0, cai will NEVER converge toward cai0, no matter how
  big the difference between cai and cai0

Both of these are absurd.

[GTR]

yes this one.me I have a copy from A.P.S., you are right when cai0 equals to cai we have a total zero,cai'=0.It doesn't stand so I'll assume that there is a '+' between.I thought that all these papers are thorougly examined before being published.But I think that this ccanl.mod might fit my needs if properly modified unless you can think for another one.

So what about my previous question:why do I get this compilation error in the CaL.mod?
nocmodl CaL
Translating CaL.mod into CaL.c
cai used as both variable and function at line 94 in file CaL.mod
^
make: *** [CaL.o] Error 1

[ted]

Unfortunately I have to wrestle with all these currents and concentrations (CaL,CaN,CaT)

Follow the recommendations stated above--see my msg that contains

The place to start is with the calcium accumulation mechanism.

This

All calcium currents should
USEION ca READ eca, cai WRITE ica

ensures that when your ca accumulation mechanism READs ica, the value it gets is
the sum of the contributions from all mechanisms that WRITE ica.

Why bother with this

Code: Select all

DERIVATIVE integrate {
ncai' = -(inca)/depth/FARADAY * (1e7) + (caiinf/3 - ncai)/catau
lcai' = -(ilca)/depth/FARADAY * (1e7) + (caiinf/3 - lcai)/catau
tcai' = -(itca)/depth/FARADAY * (1e7) + (caiinf/3 - tcai)/catau

when

Code: Select all

DERIVATIVE integrate {
cai' = -ica/depth/FARADAY * (1e7) + (caiinf - cai)/tau

does the same thing at 1/3 the computational effort? By the way, why use caiinf instead
of cai0?

in the CaL.mod I'm still getting compilation error:
nocmodl CaL
Translating CaL.mod into CaL.c
cai used as both variable and function at line 94 in file CaL.mod

If simply looking at the code doesn't tell you what has caused a conflict in name space,
it's time to do some diagnostic experiments. Look for occurrences of the offending string
(cai in this case), pick one that seems likely, and subsitute something else for it. My
guess is the problem lies in PROCEDURE rate. You might try rewriting it as
PROCEDURE rate (Cai(mM),v(mV)) etc. (i.e. substitute Cai for cai within the procedure
definition). If that doesn't work, maybe the problem is with the way you are calling this
procedure. In the DERIVATIVE block, you might insert a space between the cai and
the (.

I thought that all these papers are thorougly examined before being published.

They are, but it's very easy to miss typos, especially in equations and numeric values.
The authors know what they meant the equation to say, so that's what they perceive
every time they look at it. The reviewers have so much other work to do, that they
can't be preoccupied with every little thing in a paper. The typesetters may have no
idea what the equation is supposed to be, and if they make a mistake, it's up to someone
else to catch it.

That's one of the reasons why it is basically impossible to reproduce most published
models from the information presented in the paper itself. A careful reader who thinks
about everything can discover that this particular error is an error, because of its
absurd consequences. But there are many opportunities for errors that are much
more subtle, that nobody would be able to look at and say, "Aha, this is clearly wrong."
And that's why model code, like sequence data, should be deposited in a searchable
public archive when a paper has been accepted for publication. Otherwise,
computational modeling is not reproducible, and without reproducibility, how can
it be science?

[GTR]

Dear Ted,

I think I found partially the solution to my problem.As far as my KCa.mod is concerned (nrnmech.dll was built successfully).The idea of setting a local function z in the PROCEDURE definition:

Code: Select all

PROCEDURE rates(cai(mM)) {
   LOCAL z
  z=exp(-log(cai)+0.3)        
qt = q10^((celsius - temp)/10)
         wtau =40
	winf =0.81/((1+z)/0.46)
}     

instead of:

Code: Select all

PROCEDURE rates(cai(mM)) {         
qt = q10^((celsius - temp)/10)
        
        wtau =40
	winf =0.81/(1+exp(-log(cai)+0.3))/0.46)

 }   

solved my problem concerning the compilation error I got if you remember :

nocmodl CaL
Translating CaL.mod into CaL.c
cai used as both variable and function at line 94 in file CaL.mod.

So in general my question is:

Why is it difficult for the compiler to read the mathematical expression for winf as it is without setting a local z?I noticed that this technique of 'setting a local function' is followed in the majority of mod files in modelDB. Why has the compiler problem to read as a whole
winf =0.81/(1+exp(-log(cai)+0.3))/0.46)?Because of my false arithmetical operation priority?

And secondly I made the same modification in my CaL.mod file setting LOCALs w,a,b


but now a new compilation error occured:
nocmodl CaL
Translating CaL.mod into CaL.c
v used as both variable and function at line 103 in file CaL.mod

I want to believe that now my problem doesn't lie again in substituting V for v within the procedure definition or insering spaces between ( and v as with cai if you remember?Do you have any idea to suggest?

I appreciate your help,
George.

[ted]

Sorry I didn't examine the code you posted on April 17 in any detail prior to today.
There were far more fundamental issues at hand then, e.g. the form of the mass
balance equation for calcium, whether it made sense to split the intracellular pool of
calcium into different subsets according to which channels it happened to flow
through, etc..

The NMODL code from April 17 contained many errors, most of which you have fixed
in your most recent message. There are also some instances of poor text formatting
that IMO should have generated error messages (for example, v(mV) which needs
whitespace between the v and the opening parenthesis). Unfortunately, the NMODL
compiler does not provide truly helpful clues to most of these, and it has a lazy habit
of complaining about only one error and then abandoning compilation, even when
many errors are present. So you have to be a bit of a detective.

By the way, v, cai, cao, and celsius aren't really parameters. Their values are assigned
by NEURON from elements of the cell model that lie outside the CaL mechanism, so it
makes more sense to declare them in the ASSIGNED block.
Also, it's a good idea to consolidate all units declarations into a single UNITS block.

It is usually helpful to feed NMODL code to modlunit, which at least is faster than
doing a full compile.

Now we're starting to get into serious errors.
The first invocation of modlunit (modlunit cal.mod if you're using Linux) yields
Warning: Cai undefined. (declared within VERBATIM?)
Warning: mM undefined. (declared within VERBATIM?)
Warning: mV undefined. (declared within VERBATIM?)
v used as both variable and function at line 100 in file cal.mod

The first msg brings our attention to
rates(v(mV),Cai(mM))
in the BREAKPOINT block. Three errors here.
1. Cai is local to PROCEDURE rates.
2. v is not a FUNCTION
3. neither is cai
The fix is
rates(v,cai)

It is also worthwhile to check all occurrences of Cai. This should tip you off to the error
in PROCEDURE rates, where
w=exp((cai-0.7)/0.15)
should really be
w=exp((Cai-0.7)/0.15)

Check with modlunit and see that all 4 error msgs are gone, but a new one appears:
too few arguments at line 51 in file cal.mod
rates(v)<<ERROR>>
The fix is obvious.

Less obvious: the INITIAL block fails to initialize h. Another obvious fix.

Running modlunit again reports
units: 10 sec-coul2/m4-kg

units: 0.001 sec-coul2/m2-kg
The units of the previous two expressions are not conformable
at line 60 in file cal.mod
g=gbar*q*q*h<<ERROR>> : maximum channel permeability

The cause?
The units returned by FUNCTION ghk are incompatible with a conductance-based
BREAKPOINT block forumulation
g=gbar*q*q*h
ica = g*ghk(v,cai,cao)
where gbar and g are in (mS/cm2). You need to switch to a permeability-based
formulation, i.e.
open = q*q*h
ica = p*open*ghk(v,cai,cao)
where p is a PARAMETER with units (cm/s), and open is an ASSIGNED that
represents the fraction of channels that are in the open state. If you insist on a
conductance-based formulation then you might want to emulate the style adopted by
Lazarewicz et al. in their CAlM95 mechanism
http://senselab.med.yale.edu/SenseLab/M ... CAlM95.mod
(and pray that an orthodox biophysicist never sees your code)

Other errors may also be present that I have not discovered.

[GTR]

So I attempted to switch to a permeability-based formulation that you had indicated me so as to achieve units consistency in current equation but when I a tried to run a voltage simulation of my neuron in resting state , graph went out of range.

A math function was called that returned an out of range value
errno=34 during call to mechanism CaL
nrniv: errno set during call to INITIAL block
near line 1710
{run()}
^
exp(Inf) out of range, returning exp(700)
exp(Inf) out of range, returning exp(700)
exp(Inf) out of range, returning exp(700)
No more errno warnings during this execution

What's the error in the INITIAL block of CaL.mod?

I understand that I probably have additional mistakes that play important role.

[ted]

Here's how to diagnose the problem yourself: embed a printf statement just after each
line that involves exp. The syntax is identical to C's printf. For example, right after the
line that calculates w

Code: Select all

w=exp((Cai-0.7)/0.15)

insert the following line

Code: Select all

printf("w is %g\n", w)

This will help you locate the cause of the problem. Once you have found the offending line,
insert anoter printf just before it that reports the values of the relevant variables. In this
example, the diagnostic statement would be

Code: Select all

printf("Cai is %g\n", Cai)

Other things to think about:
Have you checked this file with modlunit?
What is the value of cao?
Why calculate qt if you're not going to use it?
If you are going to use qt, why calculate it from a LOCAL q10, which will know nothing
about the q10 that you defined in the PARAMETER block?
Are you absolutely sure that the expression for w should be w=exp((Cai-0.7)/0.15) ? How big do you expect Cai to be?