pump kinetic scheme

jamie · Post by **jamie** » Thu Jan 04, 2007 12:48 pm

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)
	FARADAY = (faraday)	 (10000 coulomb)
	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
voli = PI*diam*diam/4       : not sure about this decleration of voli here
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
: ~ nai << (-ina/(FARADAY)*PI*diam*(1e4))                          
: 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)
ina_pmp = 2*FARADAY*(f_flux - b_flux)/parea
}

Raj · Post by **Raj** » Thu Jan 04, 2007 2:56 pm

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 · Post by **jamie** » Thu Jan 04, 2007 5:06 pm

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 · Post by **Raj** » Fri Jan 05, 2007 3:35 am

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 · Post by **jamie** » Fri Jan 05, 2007 7:53 am

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 · Post by **jamie** » Sat Jan 06, 2007 11:09 am

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 · Post by **jamie** » Sat Jan 06, 2007 2:23 pm

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 · Post by **Raj** » Sun Jan 07, 2007 4:24 am

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 · Post by **jamie** » Tue Jan 09, 2007 10:16 am

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:
http://senselab.med.yale.edu/senselab/m ... nadifl.mod

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 {
	SUFFIX nadifl
	USEION na READ ina WRITE nai
	RANGE D
}

UNITS {
	(mM) = (milli/liter)
	(um) = (micron)
	FARADAY = (faraday) (coulomb)
	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}
	~ nai << (-ina/(FARADAY)*PI*diam*(1e4))
}

Raj · Post by **Raj** » Tue Jan 09, 2007 11:25 am

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 · Post by **menica** » Fri Oct 20, 2017 10:45 am

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 
	tadj
}
 
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
	tadj = 3.0^((celsius-23.5)/10)
	:"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)
   FARADAY = 9.6485e4 (coulombs)
   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 · Post by **Raj** » Fri Oct 20, 2017 12:07 pm

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 · Post by **menica** » Mon Oct 23, 2017 4:56 am

Dear Raj,
thanks for your answer.
I hope someone will reply to my post.
Best
Menica

Post by **ted** » Mon Oct 23, 2017 10:59 am

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 · Post by **menica** » Tue Oct 24, 2017 1:58 pm

Dear Ted,
thanks for your input.
I read about control theory (new to me) and I understood the steady-state problem of the proportional controllers.
One solution to this problem would be adding an integral term (over the change of nai) to the controller, but I don't know how to do this in nmodl.
An another way will be deactivating the Na channel at rest (as a voltage-gated ion channel is supposed to work) but it is not clear to me how to achieve this. when I wrote ina=inapump and inapump is never 0.

www.neuron.yale.edu

pump kinetic scheme

pump kinetic scheme

Re: pump kinetic scheme

Re: pump kinetic scheme

Re: pump kinetic scheme

Re: pump kinetic scheme

Re: pump kinetic scheme