pump kinetic scheme

NMODL and the Channel Builder.
jamie

pump kinetic scheme

Below is a .mod mechanism for a sodium-potassium pump. It is a modification of Carmen Canavier's mechanism which can be found here:

http://senselab.med.yale.edu/senselab/m ... model=6763

(it is to be used with another of Canavier's mechanisms, for intracellular sodium accumulation).

The modification made is that a phenomenological lag term, tau_pump, has been added. Unfortunately this modification makes the mechanism unstable, I think because it promotes/allows negative values of nai. As I understand it, the mechanism would be more robust and secure when using a Kinetic scheme. So, I am now trying to express the below in a kinetic scheme based mechanism.

Code: Select all

``````
TITLE PUMP

UNITS {
(molar) = (1/liter)
(pA) = (picoamp)
(mV) =  (millivolt)
(uS) = (micromho)
(mA) =  (milliamp)
(mM) =  (millimolar)
}

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

NEURON {
SUFFIX pump
USEION na READ nai  WRITE ina
USEION k  WRITE ik
RANGE  inapump,ipumpmax,n,km,ipump_bar
}

PARAMETER {
dt                (ms)
nai               (mM)
ipumpmax  = 0.04 (mA/cm2)
km = 10.0        (mM)
n=1.5

nainit = 4               (mM)
celsius = 35             (degC)
tau_pump = 10000         (ms)
}

ASSIGNED {
ina             (mA/cm2)
ik              (mA/cm2)
ipump_bar       (mA/cm2)
}

BREAKPOINT {
ipump_bar = ipumpmax*(1/(1 + pow(km/nai,n)))
SOLVE state METHOD cnexp
ina = 3.0*inapump
ik = -2.0*inapump
}

STATE {inapump}        : activation variable to be solved in the DEs

DERIVATIVE state {
inapump' = (ipump_bar - inapump) / tau_pump
}

INITIAL {
inapump = ipump_bar
}

``````

So, I have cut and paste my best effort below at reproducing the above in a kinetic scheme based mechanism. It actually compiles OK. But when put into a cell model, causes errors (errnos).

As an aside, I am really not sure what to set the rates to. I am not sure how to use the Hill equation parameter values in the Canavier mechanism to set the rates in this kinetic scheme.

I write here because maybe someone has trodden this road before. I would be very grateful for any advice. Wishing you a happy new year.

Code: Select all

``````

NEURON {
SUFFIX pump
USEION na READ nai, nao, ina  WRITE ina
USEION k READ ko, ki, ik WRITE ik
RANGE ina_pmp, TotalPump
}

UNITS {
(mol) = (1)
(molar) = (1/liter)
(mM)	= (millimolar)
(um)	= (micron)
(mA)	= (milliamp)
PI	= (pi) (1)
}

PARAMETER {
k1 = 1 (/mM-ms)
k2 = 0.005 (/ms)
k3 = 1 (/ms)
k4 = 0.005 (/mM-ms)
TotalPump = 1e-11 (mol/cm2)
}

ASSIGNED {
nai		(mM)
nao (mM)
ki (mM)
ko (mM)
ina_pmp (mA/cm2)
ina_pmp_last (mA/cm2)
parea (um)
diam		(um)
ina		(mA/cm2)
ik  (mA/cm2)
}

: volume extracellular space.value irrellevant as nao is treated constant.
CONSTANT { volo = 1e10 (um2)
}

STATE {
pump (mol/cm2)
nkpump (mol/cm2)
}

BREAKPOINT {
SOLVE state METHOD sparse
ina_pmp_last = ina_pmp
ina = ina_pmp
}

INITIAL{
parea = PI*diam
pump = TotalPump/(1 + (nai*k1/k2))
nkpump = TotalPump - pump
}

KINETIC state { LOCAL voli
COMPARTMENT (1e10)*parea {pump nkpump}
COMPARTMENT voli { nai ki }
: COMPARTMENT voli PI*diam*diam/4 {nai ki}
COMPARTMENT volo { nao ko }
: ~ nai << (-(ina - last_ina_pmp)*vol2surfaceratio/FARADAY
: these lines are for declaring nai accumulation. I think they are taken
: care of in a seperate nai accumulation mechanism so not needed
~ 3 nai + 4 ko + pump <-> nkpump	(k1*parea*(1e10), k2*parea*(1e10))
~ nkpump <-> 3 nao + 4 ki		(k3*parea*(1e10), k4*parea*(1e10))
CONSERVE pump + nkpump = TotalPump * parea * (1e10)
}

``````

Raj
Posts: 219
Joined: Thu Jun 09, 2005 1:09 pm
Location: Hanze University of Applied Sciences
Contact:
I remember a discussion in this forum with David Sterratt about how to add non-linearities to your differential equations. Michael Hines in response to our attempts to do an analysis discussed in some detail what is the problem and what is the way out. You might want to have a look at his post in the following topic:
https://www.neuron.yale.edu/phpBB2/view ... light=#774
and work on the basis of the original file.

jamie
hi raj. thanks so, so much for replying. Although im a bit confused as Im not sure what I can take from that thread as it doesnt refer to KINETIC schemes. Im really sorry but i'm a bit unclear. Would be so great if you could spell it out for me a bit.

Raj
Posts: 219
Joined: Thu Jun 09, 2005 1:09 pm
Location: Hanze University of Applied Sciences
Contact:
My comment regards your non-kinetic solution to the problem and indicates the probable cause of its failing. The thread I directed you to can help you with correctly implementing the change you had in mind without reimplementing the whole mechansim as a kinetic scheme.

It is up to you to decide whether you want to fix the original solution you had in mind or move to kinetic schemes.

jamie
ah. thanks raj. sorry for my misunderstanding. Yes. perhaps it would be better if I can fix the original problem without having to go down the KINETIC route. So, I will read this Sterrat thread with this now in focus and let you know back here how I get on.

jamie
Hi raj. once again thanks so much for your help. So, I have read the Sterrat thread and taken a few things away from it to modify the pump mechanism. Alas, I still dont have a stable mechanism. Compiles ok. but crashes a cell model because nai is still able to go negative. I have cut and paste the latest version below. There is a lot commented out for lines and avenues I went down, that didn't work. Perhaps the biggest advance with the Sterrat thread input is that I am now using derivimplicit.

Code: Select all

``````
TITLE PUMP

UNITS {
(molar) = (1/liter)
(pA) = (picoamp)
(mV) =  (millivolt)
(uS) = (micromho)
(mA) =  (milliamp)
(mM) =  (millimolar)
}

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

NEURON {
SUFFIX pump
USEION na READ nai  WRITE ina
USEION k  WRITE ik
RANGE  inapump,ipumpmax,n,km,ipump_bar
}

PARAMETER {
dt                (ms)
nai               (mM)
ipumpmax  = 0.04 (mA/cm2)
km = 10.0        (mM)
n=1.5

nainit = 4               (mM)
celsius = 35             (degC)
tau_pump = 10000         (ms)
}

ASSIGNED {
ina             (mA/cm2)
ik              (mA/cm2)
ipump_bar       (mA/cm2)
}

BREAKPOINT {
:        ipump_bar = ipumpmax*(1/(1 + pow(km/nai,n)))
:         SOLVE state METHOD cnexp
SOLVE state METHOD derivimplicit
:        ina = 3.0*inapump
:        ik = -2.0*inapump
}

STATE {inapump}        : activation variable to be solved in the DEs

DERIVATIVE state {
ipump_bar = ipumpmax*(1/(1 + pow(km/nai,n)))
ina = 3.0*inapump
ik = -2.0*inapump
inapump' = (ipump_bar - inapump) / tau_pump
: if (inapump <= 0. ) { inapump = 0. } : cannot pump inward
: if (ipump_bar <= 0. ) { ipump_bar = 0. } : cannot pump inward
: if (ina <= 0. ) { ina = 0. }
}

INITIAL {
inapump = ipump_bar
}

``````

jamie
Error message recieved when I use the mechanism in a cell model. when it crashes.

"a math function was called with argument out of domain
errno=33 during call to mechanism pump"

Raj
Posts: 219
Joined: Thu Jun 09, 2005 1:09 pm
Location: Hanze University of Applied Sciences
Contact:
What is your initial value nai? If it is very small or (less likely) large compared to km it might explain your error.

The assignments to the currents do belong in the breakpoint block.

You might try to avoid evaluating pow for small and large km/nai values:

Define what you find large and what small, i.e. define epsilon (eps) and assign a small value to it e.g. 1e-6 and use this to see whether you need to evaluate the pow function by replacing your ipump_bar calculation with the following code:

Code: Select all

``````if(km/nai < eps ){
ipump_bar = ipumpmax
}else if (km/nai > 1/eps ){
ipump_bar = 0
}else{
ipump_bar = ipumpmax*(1/(1 + pow(km/nai,n)))
}``````
If nai can become zero you should test nai/km provided km is non-zero.

jamie
Thanks so much raj. And I have tried your solution and it works in that the simulation no longer crashes. Great! That is the good news. But the bad news is that the voltage trajectory jumps around a bit when we get to the point in the simulation where nai approaches 0. But then it doesnt crash. And that is great news.

I came up with an alternative solution and I would really appreciate knowing what you think about it. As it is a bit of a shot in the dark by me. It does have the benefit though that the voltage trajectory doesnt jump around as nai approaches 0 and there is no crashing.

Instead of adjusting the pump mechanism I decided to make a modification to the sodium accumulation mechanism. I added this line:

if(nai <= 0. ){ nai = 0.0001}

The accumulation mechanism in full, with this modification:

NB. this mechanism is, aside from the modification, basically equivalent to this modeldb entry:

Code: Select all

``````
COMMENT
Longitudinal diffusion of sodium (no buffering)
(equivalent modified euler with standard method and
equivalent to diagonalized linear solver with CVODE )
ENDCOMMENT

NEURON {
USEION na READ ina WRITE nai
RANGE D
}

UNITS {
(mM) = (milli/liter)
(um) = (micron)
PI = (pi) (1)
}

PARAMETER {
D = .6 (um2/ms)
}

ASSIGNED {
ina (milliamp/cm2)
diam (um)
}

STATE {
nai (mM)
}

INITIAL {

nai = 0

}

BREAKPOINT {
SOLVE conc METHOD sparse
if(nai <= 0. ){ nai = 0.0001}
}

KINETIC conc {
COMPARTMENT PI*diam*diam/4 {nai}
LONGITUDINAL_DIFFUSION D {nai}
}

``````

Raj
Posts: 219
Joined: Thu Jun 09, 2005 1:09 pm
Location: Hanze University of Applied Sciences
Contact:
The modification made is that a phenomenological lag term, tau_pump, has been added. Unfortunately this modification makes the mechanism unstable, I think because it promotes/allows negative values of nai.
Which phenomenom are you trying to model with tau_pump? Normally in reaction or reaction diffusion the depletion of one of the substances in the equations will stop the reaction, as you noted by adding your decay this might now happen to late with zero crossings of nai as a result.

Simply fixing your problem by defining lowest value for nai, might be a good thing to do, but that requires that details of the dynamics of nai below a certain value are not relevant to your model and that the dynamics is faithfully represented otherwise and furthermore the fixing of the lowest level should not significantly influence the build up of relevant levels of nai. If one of these conditions does not hold then your solution is probably wrong.

menica
Posts: 71
Joined: Wed Sep 02, 2015 11:02 am

Re: pump kinetic scheme

Dear Raj,
I am dealing with the same problem discussed in this thread.
I need to insert a tau term in order to have control on the pump dynamics.
I am considering only one soma in which I inserted a voltage-activate Na channel, Na pump and Na accumulation mechanism.
Here are the mod files:
Na voltage gated channel:

Code: Select all

``````TITLE dgkin-Huxley like sodium  :mrate>1 fast activation, <1 slowactivation

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

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

NEURON {
SUFFIX nahh
USEION na READ ena WRITE ina
RANGE gnabar,m_inf,h_inf,tau_m,tau_h,ina,mrate,hrate
}

PARAMETER {
v		(mV)
celsius		(degC)
dt		(ms)
gnabar= 0.1	(mho/cm2)
ena		(mV)
mrate =1 (1)
hrate =1 (1)
}

STATE {
m h
}

ASSIGNED {
ina	(mA/cm2)
m_inf h_inf tau_m tau_h
}

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

UNITSOFF
INITIAL {
rates(v)
m = m_inf
h = h_inf
}

PROCEDURE states() {
rates(v)
m = m + (1-exp(-dt/tau_m)) * (m_inf-m)
h = h + (1-exp(-dt/tau_h)) * (h_inf-h)
}

PROCEDURE rates(v) { LOCAL alpha, beta, q10, tinc
:TABLE m_inf, tau_m, h_inf, tau_h DEPEND dt,
: celsius FROM -100 TO 100 WITH 200
:"m" sodium activation system
alpha = .091 * vtrap(v+38,5)
beta =  .062 * vtrap(-(v+38),5)
tau_m = 1 / (alpha+beta) / tadj / mrate
m_inf = alpha/(alpha+beta)
:"h" sodium inactivation system
alpha = .016 * exp(-(v+55)/15)
beta = 2.07 / (1+exp((17-v)/21))
tau_h = 1 / (alpha+beta) / tadj / hrate
h_inf = alpha/(alpha+beta)

}

FUNCTION vtrap( x, b) {
: Traps for 0 in denominator of rate equations
if (fabs(x/b) < 1e-6) {
vtrap = b+x/2 }
else {
vtrap = x / (1-exp(-x/b)) }
}
UNITSON``````
Na accumulation

Code: Select all

``````NEURON
{
SUFFIX naaccum

USEION na READ ina WRITE nai, nao
RANGE nai0, nao0
RANGE rho
}

UNITS
{
(molar) =          (1/liter)
(mM)   =          (millimolar)
(um)   =          (micron)
(mV)    =          (millivolt)
(mA)   =          (milliamp)
R        = 8.3134   (joule/degC)
PI      = (pi)     (1)
}

PARAMETER
{
nai0 = 10   (mM)
nao0 = 140   (mM)

rho  = 0.2    (1) : extra vol : intra vol
diam             (um)
celsius          (degC)

}

ASSIGNED
{
ina (mA/cm2)
nai (mM)
nax (mM)
}

STATE
{
nao  (mM)
}

INITIAL
{
nai = nai0
nao = nao0
nax = nai0 + nao0 * rho
}

BREAKPOINT
{
SOLVE state METHOD cnexp
}

DERIVATIVE state
{
nao' = 4 * (ina)   / (rho*diam*FARADAY)  * (1e4)
nai  = nax - nao * rho
}
``````
Na pump:

Code: Select all

``````UNITS {
(molar) = (1/liter)
(pA) = (picoamp)
(mV) =	(millivolt)
(uS) = (micromho)
(mA) =	(milliamp)
(mM) =	(millimolar)
}

NEURON {
SUFFIX pumpnatau
USEION na READ nai  WRITE ina
RANGE  inapump,nai0,ipumpmax,n,km, ipump_bar, tau_pump , ina
}

PARAMETER {
dt (ms)
nai   (mM)
ipumpmax  = 0.013   (mA/cm2)
km = 10       (mM)
n=1.5
tau_pump = 100         (ms)
nai0 = 10   (mM)
}

ASSIGNED {
ina		(mA/cm2)
ipump_bar  (mA/cm2)
}

BREAKPOINT {
SOLVE state METHOD derivimplicit
}

STATE {inapump}

INITIAL {
nai=nai0
inapump = ipump_bar
}

DERIVATIVE state {

if(nai <= nai0 ){
ipump_bar = 0
}else{
ipump_bar = ipumpmax*(1/(1 + pow(km/nai,n)))
}
inapump' = (ipump_bar - inapump) / tau_pump
ina = inapump
}

``````
First, when I insert a Iclamp current and acvtivate the sodium channel, this leads to a change in the ion concentrations that should be restored by the pump. but the pump does not restore the nai concentration to its inital value.
Second, when I initialize the code I can see that there is a lot of instability.
How can I fix this problems?
Best

Raj
Posts: 219
Joined: Thu Jun 09, 2005 1:09 pm
Location: Hanze University of Applied Sciences
Contact:

Re: pump kinetic scheme

Dear Menica,

I am not an active Neuron user anymore. So unfortunately I cannot help.

I hope somebody else can help out.

Kind regards,
Raj

menica
Posts: 71
Joined: Wed Sep 02, 2015 11:02 am

Re: pump kinetic scheme

Dear Raj,
I hope someone will reply to my post.
Best
Menica

ted
Posts: 5513
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: pump kinetic scheme

The sodium accumulation mechanism should use derivimplicit, not cnexp. That may take care of stability issues. See
Integration methods for SOLVE statements
in the Forum's Hot tops area.
when I insert a Iclamp current and acvtivate the sodium channel, this leads to a change in the ion concentrations that should be restored by the pump. but the pump does not restore the nai concentration to its inital value.
Your pump's ODE describes a proportional controller. As long as the perturbing input (presence of an Na channel that is open at rest) exists, there will be a difference between the steady state sodium concentration and the pump's "set point." That's an inherent property of proportional controllers.

menica
Posts: 71
Joined: Wed Sep 02, 2015 11:02 am

Dear Ted,