RANGE variables not outputting data

[fred]

I'm trying to get data out of my model by declaring RANGE variables, which should make available key inner workings of my mechanism. When I try to graph these RANGE variables, there is no plot. Where have they gone?

The RANGE variables are defined at the end of the BREAKPOINT block, and they depend only on PARAMETERs and STATEs. From reading "Expanding NEURON's Repertoire of Mechanisms with NMODL" (Hines and Carnevale), my code should work.

All of the RANGE variables in this mod file disappear in the same way. I also tried moving the equations defining these variables into the KINETIC block with no success. I've had this problem in both NEURON 5.9 and the latest alpha version of NEURON. This is my code:

Code: Select all

:Accum_Ca

NEURON {
    SUFFIX Accum_Ca
    USEION ca READ cai, ica WRITE cai
    USEION cass READ cassi, icass WRITE cassi VALENCE 2
    RANGE Jrel, Jup, Jleak, Jtr, Jxfer
}

UNITS { 
    (S) = (siemens)
    (mV) = (millivolt)
    (mA) = (milliamp)
    (mM) = (milli/liter)
    (pA) = (picoamp)
    (uF) = (microfarad)
    (pF) = (picofarad)
}

PARAMETER {
    :Table 2 - Cell geometry params
        Acap = 1.534e-4 (cm2)
        Vmyo = 25.84e-6 (microliter)
        VJSR = 0.12e-6 (microliter)
        VNSR = 2.098e-6 (microliter)
        Vss = 1.485e-9 (microliter)

    F = 96.5 (kilocoulombs)     :C/mmol  NOT C/mM !!

    iCaL_max = 7.0e-3 (mA/cm2)

    :table 6- Buffering parameters
:        LTRPN_tot = 0.0700 (mM)
:        HTRPN_tot = 0.1400 (mM)
        k_htrpn_pos = 2.37 (/mM-ms)
        k_htrpn_neg = 3.2e-5 (/ms)
        k_ltrpn_pos = 32.7 (/mM-ms)
        k_ltrpn_neg = 0.0196 (/ms)
        CMDN_tot = 0.0500 (mM)
        CSQN_tot = 15.0000 (mM)
        K_CMDN_m = 0.000238 (mM)
        K_CSQN_m = 0.8000 (mM)

    :table 4- SR parameters
        new1 = 4.5 (/ms)
        new2 = 1.74e-5 (/ms)
        new3 = 0.00045 (mM/ms)
        Km_up = 0.0005 (mM)
        tau_tr = 20.0 (ms)
        tau_xfer = 8.0 (ms)
        ka_pos = 0.006075e12 (/mM4-ms)
        ka_neg = 0.07125 (/ms)
        kb_pos = 0.00405e9 (/mM3-ms)
        kb_neg = 0.965 (/ms)
        kc_pos = 0.009 (/ms)
        kc_neg = 0.0008 (/ms)
        n = 4   :these 2 (n,m) are hard coded-  there's a conflict with the unit-checking for 'm' (meters)
        m = 3
}

ASSIGNED { 
    v (mV)
    ica (mA/cm2)
    icass (mA/cm2)
    Bi
    Bss
    Bjsr
    Jrel (mM/ms)
    Jxfer (mM/ms)
    Jtr (mM/ms)
    Jup (mM/ms)
    Jleak (mM/ms)
}
    
STATE { 
    cai (mM) <1e-10>
    cassi (mM) <1e-10>
    cajsr (mM) <1e-10>
    cansr (mM) <1e-10>
    HTRPNCa (mM) <1e-10>
    LTRPNCa (mM) <1e-10>
    HTRPN (mM) <1e-10>
    LTRPN (mM) <1e-10>

    :Ryanodine
        PO1 PC1 PO2 PC2 P_RyR      
}

BREAKPOINT { 
    SOLVE nain METHOD sparse 
    :Calc'd here for validation

        Jxfer = (cassi - cai)/tau_xfer
        Jrel = new1*(PO1+PO2)*(cajsr-cassi)*P_RyR
        Jtr = (cansr - cajsr)/tau_tr
        Jleak = new2*(cansr-cai)
        Jup = new3*cai^2/(Km_up^2+cai^2)
}

INITIAL {
    LTRPNCa = 0.0112684 (mM)
    HTRPNCa = 0.125290 (mM)
    LTRPN = .0587316 :0.0700 - 0.0112684 (mM)     :LTRPN_tot - LTRPNCa
    HTRPN = 0.01471  :0.1400 - 0.125290 (mM)
    PC1 = 0.999817
    PC2 = 0.167740e-3
    PO1 = 1.149102e-4
    PO2 = 9.951726e-10
    P_RyR = 0.0
    cajsr = 1.29950 (mM)
    cansr = 1.29950 (mM)
}

KINETIC nain {
    :COMPARTMENT changes LHS to mass/time in explicit flux statement
    :   I choose not to use this method, using intensive values in the calc.
    :   This is clearer, I think.
    :
    :From Ted, "I always hate wrestling with that stuff.  Your approach is fine, as
    :long as you never have to change nseg or the dimensions (diameter or
    :length) of your model cell."
    :
    :   I may want to increase resolution- in this case, compartments will need to
    :       be changed to COMPARTMENT type.
        
        Bi = (1+CMDN_tot*K_CMDN_m/(K_CMDN_m+cai)^2)^-1
        Bss = (1+CMDN_tot*K_CMDN_m/(K_CMDN_m+cassi)^2)^-1
        Bjsr = (1+CSQN_tot*K_CSQN_m/(K_CSQN_m+cajsr)^2)^-1

    :Explicit Fluxes
        ~ cai << (Bi*(-(ica)*Acap/Vmyo/F/2 + (cassi - cai)/tau_xfer + new2*(cansr-cai) - new3*cai^2/(Km_up^2+cai^2)))     :these units work
        ~ cassi << (Bss*(new1*(PO1+PO2)*(cajsr-cassi)*P_RyR*VJSR/Vss-(cassi - cai)/tau_xfer*Vmyo/Vss-icass*Acap/2/Vss/F))
        ~ cajsr << (Bjsr*((cansr - cajsr)/tau_tr-new1*(PO1+PO2)*(cajsr-cassi)*P_RyR))
        ~ cansr << ((new3*cai^2/(Km_up^2+cai^2)-new2*(cansr-cai))*Vmyo/VNSR-(cansr - cajsr)/tau_tr*VJSR/VNSR)

    :Buffering in Vmyo- this replaces eqtns A14,A16-17 in Bond2004
        ~ cai + LTRPN <-> LTRPNCa (k_ltrpn_pos, k_ltrpn_neg)
        ~ cai + HTRPN <-> HTRPNCa (k_htrpn_pos, k_htrpn_neg)

    :Ryanodine Kinetics
        ~ PO1 <-> PC1 (ka_neg, ka_pos*cassi^4)
        ~ PO1 <-> PO2 (kb_pos*cassi^3, kb_neg)
        ~ PO1 <-> PC2 (kc_pos, kc_neg)
        CONSERVE PO1 + PO2 + PC1 + PC2 = 1
        ~ P_RyR << (-0.04(/ms)*P_RyR - 0.1(/ms)*icass/iCaL_max*exp(-(v-5.0)^2/648.0(mV2)))
}

[ted]

Trying to compile this with nrnivmodl (under Linux) generates this error msg:
Method sparse can't be used with Block nain
I suspect that the cause of the problem lies in the use of "explicit flux statements" that
should be recast in the form of coupled reactions. It's fine to use a statement of the form
~ cai << (-ica*k)
to represent the effect of transmembrane ca current on cai, but exchange between
intracellular compartments must be represented with NMODL's reaction syntax. Thus,
instead of

Code: Select all

~ cai << (Bi*(-(ica)*Acap/Vmyo/F/2 + (cassi - cai)/tau_xfer + new2*(cansr-cai) - new3*cai^2/(Km_up^2+cai^2)))

one should write something like

Code: Select all

~ cai << (Bi*(-(ica)*Acap/Vmyo/F/2))
~ cai <-> cassi  (Bi/tau_xfer, Bi/tau_xfer)
~ cai <-> cansr  (Bi*new2, Bi*new2)

plus some kind of reaction for the fourth term (what does it represent?). The three other
explicit flux expressions probably also need to be broken down into families of reactions.

With all these reactions, I don't see how it will be possible to avoid using
COMPARTMENT statements to ensure consistency of units.

Bi looks like an "instantaneous buffering with saturation" factor; since it involves cai and
appears in the rate constants for transitions between states, it may degrade stability.
In that case, it may be necessary to represent the buffering reaction--just make up a
new pair of STATEs for the bound and free buffer, and give the rate constants of the
buffering reaction appropriate values to achieve the correct affinity whlie also making
them 3 or 4 orders of magnitude faster than everything else.

I made an initial stab at refactoring this NMODL code by making a duplicate copy of
the KINETIC block and surrounding it with COMMENT . . . ENDCOMMENT. Then in
the rest of the NMODL code, I commented out all STATEs, except for cai, and all of the
expressions that are invalidated by doing that (modlunit helped me find those). This
preserved your carefully entered notation in a way that allows easy reuse when the
time comes.

Next I made a copy of the commented-out explicit flux statement for cai, uncommented it,
and changed it to

Code: Select all

        ~ cai << (Bi*(-(ica)*Acap/Vmyo/F/2 ))

At this point the mechanism does nothing but accumulate ca into a single pool (with what
looks like instantaneous saturating buffering). It did compile with nrnivmodl, and I was
able to put it in a single compartment model along with hh and cachan (from the
cachan.mod that is in c:\nrnxx\examples\nrniv\nmodl (also in
nrn-x.x/share/examples/nrniv/nmodl/cachan.mod of the expanded source code for
NEURON). Triggering a spike did produce an increase of cai, as expected, so it was
ready for further revision to include other states/reactions/compartments.

I am concerned about this mechanism, however, because inserting the statement
SOLVE nain STEADYSTATE sparse
into its INITIAL block caused problems--the first init & run succeeded, the second one
produced an incorrect result with a nonsense value for cai (5980.5603), and the third
failed with this error msg

Code: Select all

oc>at line 114 in file cabal.mod:
        SOLVE nain STEADYSTATE sparse
Error at section location soma(0.5)
Convergence not achieved in maximum number of iterations
/usr/local/nrn/i686/bin/nrniv: scopmath library error
 near line 1
 {run()}
        ^
        finitialize(-65)
      init()
    stdinit()
  run()

which may be due to the fact that there is nothing that pumps ca out of the cell. At
some point it will probably be necessary to be able to do a STEADYSTATE initialization
of this mechanism because of the complexity of the various reactions involved.

[fred]

I moved the J defining statements into the KINETIC block, and now their values are correctly calculated. I don't understand why they were not calculated correctly when they were located in the BREAKPOINT block. They depend on states and parameters only, just like other range variables I've used in the past (and as seen in examples in "Expanding NEURON’s Repertoire of Mechanisms with NMODL" by Hines and Carnevale).

[hines]

When I try to graph these RANGE variables, there is no plot. Where have they gone?

Do you mean the values are identically zero? I think it would be best take this offline to email (michael.hines@yale.edu) and send me in a zip file the hoc/ses/mod files necessary to exhibit the problem.
The calculation of flux values within the BREAKPOINT block should indeed work. I see several issues of concern in the KINETIC block but I don't want to prejudge without seeing actual simulation results.

[fred]

The values were identically zero, but in preparing the files for you, I seem to have solved the problem. Perhaps some kind of error had built up in my testing folder, where I had compiled many versions of the NMODL code. This is the second time I've had a problem disappear when I compile into a new folder.

RE: Comments on the KINETIC block:

instead of

Code: Select all

~ cai << (Bi*(-(ica)*Acap/Vmyo/F/2 + (cassi - cai)/tau_xfer + new2*(cansr-cai) - new3*cai^2/(Km_up^2+cai^2)))

one should write something like

Code: Select all

~ cai << (Bi*(-(ica)*Acap/Vmyo/F/2))
~ cai <-> cassi  (Bi/tau_xfer, Bi/tau_xfer)
~ cai <-> cansr  (Bi*new2, Bi*new2)

plus some kind of reaction for the fourth term (what does it represent?). The three other
explicit flux expressions probably also need to be broken down into families of reactions.

The fourth term represents active Ca uptake by the SR. I don't see a way to write this flux without using the explicit form, " A << x*y*C". In my model, I can rewrite all but 1 flux using the preferred reaction syntax, but this doesn't appear to be necessary because my model, as written, appears to run without error (though it's not bug-free yet). Do you expect that, as written, my code might have some hidden error?

[hines]

This is the second time I've had a problem disappear when I compile into a new folder.

MSWIN, LINUX, or MAC? In the old folder with the problem, if LINUX or MAC does the problem go away when you remove the i686 or powerpc subdirectory? If MSWIN does the problem go away when you

Code: Select all

rm  *.c *.o *.dll

I am assuming the problem is a failure to recompile when one of the dependencies has changed. Sometimes that can be traced to a bad date (far in the future) of the dependent file.

The fourth term represents active Ca uptake by the SR. I don't see a way to write this flux without using the explicit form, " A << x*y*C".

Active uptake often means that one introduces an enzyme along with two reactions

Code: Select all

~ s + e <-> se
~ se <-> p + e

but then, when converting to differential equations, people often do some analysis to eliminate the se state and hence the complicated denominator.
Anyway, the only danger of allowing states into the rate constants of a kinetic scheme is the potential for numerical instability due to missing jacobian elements. If it is a relatively slow process then it won't be a problem but if you make one of the coefficients in that flux very very large then I believe a numerical error (such as a negative concentration) will occur. If you need numerical stability over all parameter values and it is too tedious to reverse engineer the equations to form a pure
KINETIC scheme, then you can go the opposite direction and use a DERIVATIVE block with the
SOLVE block METHOD derivimplicit