implementing recovery from inactivation

NMODL and the Channel Builder.
Post Reply
fabien tell
Posts: 46
Joined: Mon Mar 25, 2013 1:36 pm
Location: france
Contact:

implementing recovery from inactivation

Post by fabien tell »

Hello,

I'm a new user of Neuron. I tried to make a model for a transient potassium current (IKV4).
My model seems to give decent data when I generate activation and inactivation curves or when i inject it in current clamp.
However, this current is characterized by two constants of inactivation and by a constant for recovery from inactivation , I called tauinac1=50 (ms), tauinac2=200 (ms) and taurecov=20 (ms).

I tried to "force" the model to set the constant for recovery from inactivation when the membrane potential is more negative than -80 mV. The mod file seems OK in terms of units (checked by modlunit) and is compiled without any problem.
However when I run a protocol for measuring the recovery from inactivation, I found values that are always dependent on the constant for inactivation ( tauinac1 and tauinac2) but not on recovery from inactivation.

Thanks a lot for you help

Here is my mod file: (I put comments only in this mail)

Code: Select all

NEURON {
  SUFFIX kaNTS6
  USEION k READ ek WRITE ik
  RANGE gbar, g, i
}

UNITS {
  (S) = (siemens)
  (mV) = (millivolt)
  (mA) = (milliamp)

}

PARAMETER {
  gbar = 0.006 (S/cm2)
    eK = -90 (mV)
Vmid_ac = -45 (mV)
k_ac = 13 (mV) 
Vmid_ina = -73 (mV) 
k_ina = -4.66 (mV)
tauinac1=50 (ms)
tauinac2=200 (ms)
taurecov=20 (ms)
A1=0.5     : (A1 A2 : respective ratio of the fast and slow inactivation)
A2=0.5
m=3
h=1
: parametres des constantes 

}

ASSIGNED {
  v	(mV)
  ek	(mV)
  ik 	(mA/cm2)
  i 	(mA/cm2)
  g	(S/cm2)
atau (ms)
btau (ms)
b2tau (ms)
}

STATE {a b b2} : creation variables


BREAKPOINT {
  SOLVE states METHOD cnexp
  g = (A1*gbar*(a^m)*(b^h)) + (A2*gbar*(a^m)*(b2^h))  : I sum two currents to simulate two inactivation constants
  i = g*(v-eK)
  ik = i
}
  
INITIAL {
  b = binf(v)
b2= binf(v)
  a = ainf(v)
atau = a_tau(v)
btau=b_tau (v)  
b2tau=b2_tau(v)

}

DERIVATIVE states {
 b'= (binf(v)-b)/btau
 a' = (ainf(v)-a)/atau
b2'= (binf(v)-b2)/b2tau

}

FUNCTION ainf (V (mV)) () {

  UNITSOFF
    ainf = 1/(1+exp(-(V-Vmid_ac)/k_ac)) : values obtained from fit
  UNITSON
}


FUNCTION binf (V (mV)) () {

  UNITSOFF
    binf = 1/(1+exp(-(V-Vmid_ina)/k_ina))   : values obtained from fit
  UNITSON
}

FUNCTION b2inf (V (mV)) () {

  UNITSOFF
    b2inf = 1/(1+exp(-(V-Vmid_ina)/k_ina))  : values obtained from fit
  UNITSON
}


FUNCTION a_tau (V (mV)) (ms) {

 UNITSOFF
a_tau= 1+ ((5.5 *exp((-0.022)^2) * ((v+65)^2))) : function for activation rate I wonder if I should use V instead of v
 UNITSON

}

FUNCTION b_tau (V (mV)) (ms) {

 UNITSOFF
if (V>-80) {b_tau=tauinac1    : function for inactivation rate 
}else{
b_tau=taurecov}
 UNITSON

}

FUNCTION b2_tau (V (mV)) (ms) {
UNITSOFF
if (V>-80) {b2_tau=tauinac2     : function for inactivation rate 
}else{
b2_tau=taurecov}
UNITSON
}

PROCEDURE rates (ms) {

atau=a_tau(v)
btau=b_tau(v)        : I donit know if this part it 's really 
b2tau=b2_tau(v)
}
ted
Site Admin
Posts: 5795
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: implementing recovery from inactivation

Post by ted »

The nmodl file looks like it should work. Can you describe the "protocol" you are using to examine recover from inactivation?
fabien tell
Posts: 46
Joined: Mon Mar 25, 2013 1:36 pm
Location: france
Contact:

Re: implementing recovery from inactivation

Post by fabien tell »

Hello Ted,

Thanks for your prompt reply.
The protocol is as follows :
a cell body with only pas, HH and IKA
V clamp mode : (HH off)
time duration : 2000
first epoch duration : 1000 at -20 mv to fully inactivate IKA
second epoch: second duration: from 10 to 500 ms (20 ms increase) at -100 mV
third epoch : third duration: 500 ms at -20 mv.

The recovery from inactivation is then measured as the time (which depends on the duration of epoch 2) for the IKA current (in epoch3) to reach 63 % of its maximum value.
Using this protocol, this recovery time constant varies when I increase or decrease the inactivation rate (1 or 2 , to be easier you can put A2 t0 O to have only one time constant of inactivation) BUT not when I change taurecovery. It means that if i set the inactivation rate at 100 ms , I will find with this protocol a recovery around 100 ms and if I change taurecovery at 100 , 200 or else it do not affect it.

To be complete, I do not have in the hoc file a statement to initialize the value (initial block in the mod file) because I did not understand if it was necessary and how to do it.

Thanks again
ted
Site Admin
Posts: 5795
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: implementing recovery from inactivation

Post by ted »

OK, now that I look closer I see many problems, all fixable. But first, a couple of questions for you.

The mod file has FUNCTIONs binf() and b2inf(), but the bodies of these functions are identical, and you're not using b2inf() in the DERIVATIVE or INITIAL blocks. Do you really need both of these FUNCTIONs?

Are you quite sure about FUNCTION a_tau()? Eliminating superfluous parentheses reduces it to this
a_tau= 1+ 5.5*exp(0.022^2) * (v+65)^2
i.e. the (v+65)^2 lies outside of the exp(). This seems quite unlike the mathematical formulation for any voltage-dependent gating variable I have seen before.
fabien tell
Posts: 46
Joined: Mon Mar 25, 2013 1:36 pm
Location: france
Contact:

Re: implementing recovery from inactivation

Post by fabien tell »

Thanks a lot Ted,

Sorry for the mistakes. My answers :


"The mod file has FUNCTIONs binf() and b2inf(), but the bodies of these functions are identical, and you're not using b2inf() in the DERIVATIVE or INITIAL blocks".

A: That's was a mistake in the derivate block : Now one should read b2'= (b2inf(v)-b2)/b2tau.


"Do you really need both of these FUNCTIONs?" "

A: My goal was to model a current with two inactivation rates, a fast and a slow. It was easier two sum two currents even if they had the same parameters apart from tauinac. It may be not the right or the more elegant way to do it but for me it was easier to understand. Do you think it may be a problem?


"Are you quite sure about FUNCTION a_tau()? Eliminating superfluous parentheses reduces it to this
a_tau= 1+ 5.5*exp(0.022^2) * (v+65)^2
i.e. the (v+65)^2 lies outside of the exp(). This seems quite unlike the mathematical formulation for any voltage-dependent gating variable I have seen before."

A: you're right. Now the equation is a_tau= 1+ 5.5 *exp(-0.000484 * ((v+65)^2)) .I've taken it form a previous published article on the same model.

I corrected all mistakes but the recovery is still dictated by the inactivation rates ( tauinac1 or tauniac2).


Thanks again for your time

Regards.



: modified from SCHILD et al. J physiol (1993)


NEURON {
SUFFIX kaNTS7
USEION k READ ek WRITE ik
RANGE gbar, g, i
}

UNITS {
(S) = (siemens)
(mV) = (millivolt)
(mA) = (milliamp)

}

PARAMETER {
gbar = 0.006 (S/cm2)
eK = -90 (mV)
Vmid_ac = -55 (mV)
k_ac = 15 (mV)
Vmid_ina = -73 (mV)
k_ina = -7 (mV)
tauinac1=25 (ms)
tauinac2=200 (ms)
taurecov=20 (ms)
A1=0.75 : (A1 A2 : respective ratio of the fast and slow inactivation)
A2=0.25
m=3
h=1
: parametres des constantes

}

ASSIGNED {
v (mV)
ek (mV)
ik (mA/cm2)
i (mA/cm2)
g (S/cm2)
atau (ms)
btau (ms)
b2tau (ms)
}

STATE {a b b2} : creation variables


BREAKPOINT {
SOLVE states METHOD cnexp
g = (A1*gbar*(a^m)*(b^h)) + (A2*gbar*(a^m)*(b2^h)) : I sum two currents to simulate two inactivation constants
i = g*(v-eK)
ik = i
}

INITIAL {
b = binf(v)
b2= binf(v)
a = ainf(v)
atau = a_tau(v)
btau=b_tau (v)
b2tau=b2_tau(v)

}

DERIVATIVE states {
b'= (binf(v)-b)/btau
a' = (ainf(v)-a)/atau
b2'= (b2inf(v)-b2)/b2tau :fixed

}

FUNCTION ainf (v (mV)) () {

UNITSOFF
ainf = 1/(1+exp(-(v-Vmid_ac)/k_ac)) : values obtained from fit
UNITSON
}


FUNCTION binf (V (mV)) () {

UNITSOFF
binf = 1/(1+exp(-(v-Vmid_ina)/k_ina)) : values obtained from fit
UNITSON
}

FUNCTION b2inf (v (mV)) () {

UNITSOFF
b2inf = 1/(1+exp(-(v-Vmid_ina)/k_ina)) : values obtained from fit
UNITSON
}


FUNCTION a_tau (v (mV)) (ms) {

UNITSOFF
a_tau= 1+ 5.5 *exp(-0.000484 * ((v+65)^2)) : fixed I used v everywhere
UNITSON

}

FUNCTION b_tau (v (mV)) (ms) {

UNITSOFF
if (v>-80) {b_tau=tauinac1 : function for inactivation rate
}else{
b_tau=taurecov}
UNITSON

}

FUNCTION b2_tau (v (mV)) (ms) {
UNITSOFF
if (v>-80) {b2_tau=tauinac2 : function for inactivation rate
}else{
b2_tau=taurecov}
UNITSON
}

PROCEDURE rates (ms) {

atau=a_tau(v)
btau=b_tau(v) : I dont know if this part it 's really needed
b2tau=b2_tau(v)
}
ted
Site Admin
Posts: 5795
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: implementing recovery from inactivation

Post by ted »

Thanks, that's useful information.

There were several problems with the original implementation of this mechanism. The most important was that the time constants were not being reevaluated on every time step. To fix this, I changed the INITIAL and DERIVATIVE blocks to

Code: Select all

INITIAL {
  rates(v)
  b = binf
  b2= b2inf
  a = ainf
}
DERIVATIVE states {
  rates(v)
  a' = (ainf-a)/atau
  b' = (binf-b)/btau
  b2'= (b2inf-b2)/b2tau
}
so that rates() would be called at each initialization and time step.

I also revised PROCECURE rates so that it takes care of updating both the time constants and the steady state values.

Code: Select all

PROCEDURE rates(V (mV)) {
  atau=a_tau(V)
  ainf=a_inf(V)
  btau=b_tau(V)
  binf=b_inf(V)
  b2tau=b2_tau(V)
  b2inf=b2_inf(V)
}
Of course, this meant that I had to declare the new steady state "intermediate variables" in the ASSIGNED block

Code: Select all

ASSIGNED {
  ainf (1)
  binf (1)
  b2inf (1)
 . . .
Other changes:

The declaration line of FUNCTION a_tau uses V but the body of the function uses v; the v should really be changed to V.

For testing purposes, it is most convenient to expose parameters and intermediate variables to hoc as range variables. I did this by adding these statements

Code: Select all

  RANGE atau, btau, b2tau
  RANGE ainf, binf, b2inf
  RANGE tauinac1, tauinac2, taurecov
to the NEURON block. This allows the inactivation and recovery time constants to be changed by simple hoc assignments, instead of having to edit the mod file and recompile.

I'm emailing you the fixed mod file, plus a "virtual experimental rig" that you can use to explore the effects of parameter changes on the dynamics of the gating state variables a, b, and b2. Compile kants6.mod, then execute

nrngui test.hoc

test.hoc is really a ses file for a GUI that I set up with the NEURON Main Menu toolbar, but there's no harm in giving it a different extension--just don't expect its contents to look like code that a human would write. test.hoc implements a single compartment model that has nothing but kaNTS6 and membrane capacitance (think of this as the computational equivalent to a Xenopus oocyte with kaNTS6 channels). The cell is voltage clamped, and the clamp uses a 100 ms duration depolarizing step to induce inactivation. The "Parameters" panel allows you to change the inactivation and recovery time constants, so you can see how these affect the gating variables.

The resulting mod file will work with fixed time step integration, but it will not work properly with adaptive integration because b_tau and b2_tau change abruptly at t == -80 mv. This is a problem because an adaptive integrator is likely to choose a large integration step that overshoots the point of the discontinuity. It's very easy to change that--instead of implementing an abrupt change between two time constants, just use a sigmoid function like
y = y0 + (y1-y0)/(1+e(-x))
where
y0 = time constant when v == -infinity
y1 = time constant when v == + infinity
x = (v - vhalf)*k
vhalf = voltage at which y == (y1+y0)/2
k is a constant > 0 that specifies the steepness of the relationship between x and y at v == vhalf
I leave that up to you.

Here is the simple (abrupt change of time constants) completed mod file for anyone who is dying to see it:

Code: Select all

NEURON {
  SUFFIX kaNTS6
  USEION k READ ek WRITE ik
  RANGE gbar, g, i
  RANGE atau, btau, b2tau
  RANGE ainf, binf, b2inf
  RANGE tauinac1, tauinac2, taurecov
}

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

PARAMETER {
  gbar = 0.006 (S/cm2)
  eK = -90 (mV)
  Vmid_ac = -45 (mV)
  k_ac = 13 (mV) 
  Vmid_ina = -73 (mV) 
  k_ina = -4.66 (mV)
  tauinac1=50 (ms)
  tauinac2=200 (ms)
  taurecov=20 (ms)
  A1=0.5     : (A1 A2 : respective ratio of the fast and slow inactivation)
  A2=0.5
  m=3
  h=1
}

ASSIGNED {
  v	(mV)
  ek	(mV)
  ik 	(mA/cm2)
  i 	(mA/cm2)
  g	(S/cm2)
  atau (ms)
  btau (ms)
  b2tau (ms)
  : more new stuff
  ainf (1)
  binf (1)
  b2inf (1)
}

STATE {a b b2}

BREAKPOINT {
  SOLVE states METHOD cnexp
  : I sum two currents to simulate two inactivation constants
  g = (A1*gbar*(a^m)*(b^h)) + (A2*gbar*(a^m)*(b2^h))
  i = g*(v-eK)
  ik = i
}

COMMENT
INITIAL {
  b = binf(v)
  b2= binf(v)
  a = ainf(v)
  atau = a_tau(v)
  btau=b_tau (v)  
  b2tau=b2_tau(v)
}
ENDCOMMENT

INITIAL {
  rates(v)
  b = binf
  b2= b2inf
  a = ainf
}

DERIVATIVE states {
:  b'= (binf(v)-b)/btau
:  a' = (ainf(v)-a)/atau
:  b2'= (binf(v)-b2)/b2tau
  rates(v)
  a' = (ainf-a)/atau
  b' = (binf-b)/btau
  b2'= (b2inf-b2)/b2tau
}

COMMENT
FUNCTION ainf (V (mV)) () {
UNITSOFF
  ainf = 1/(1+exp(-(V-Vmid_ac)/k_ac))
UNITSON
}

FUNCTION binf (V (mV)) () {
UNITSOFF
  binf = 1/(1+exp(-(V-Vmid_ina)/k_ina))
UNITSON
}

FUNCTION b2inf (V (mV)) () {
UNITSOFF
  b2inf = 1/(1+exp(-(V-Vmid_ina)/k_ina))
UNITSON
}
ENDCOMMENT

FUNCTION a_inf (V (mV)) () {
  a_inf = 1/(1+exp(-(V-Vmid_ac)/k_ac))
}

FUNCTION b_inf (V (mV)) () {
  b_inf = 1/(1+exp(-(V-Vmid_ina)/k_ina))
}

FUNCTION b2_inf (V (mV)) () {
  b2_inf = 1/(1+exp(-(V-Vmid_ina)/k_ina))
}

FUNCTION a_tau (V (mV)) (ms) {
UNITSOFF
:  a_tau= 1+ ((5.5 *exp((-0.022)^2) * ((v+65)^2)))
: after eliminating superfluous parentheses--
:  a_tau= 1+ 5.5*exp(0.022^2) * (v+65)^2 : really?
  a_tau= 1+ 5.5*exp((0.022*(v+65))^2)
UNITSON
}

FUNCTION b_tau (V (mV)) (ms) {
  if (V>-80) {
    b_tau=tauinac1
  }else{
    b_tau=taurecov
  }
}

FUNCTION b2_tau (V (mV)) (ms) {
  if (V>-80) {
    b2_tau=tauinac2
  }else{
    b2_tau=taurecov
  }
}

COMMENT
PROCEDURE rates (ms) {
  atau=a_tau(v)
  btau=b_tau(v)
  b2tau=b2_tau(v)
}
ENDCOMMENT

PROCEDURE rates(V (mV)) {
  atau=a_tau(V)
  ainf=a_inf(V)
  btau=b_tau(V)
  binf=b_inf(V)
  b2tau=b2_tau(V)
  b2inf=b2_inf(V)
}
fabien tell
Posts: 46
Joined: Mon Mar 25, 2013 1:36 pm
Location: france
Contact:

Re: implementing recovery from inactivation

Post by fabien tell »

Thanks a lot Ted,

Now it works fine and I implemented the inactivation rates as you suggested to avoid problems with integration procedure.

Thanks again for the time spent on this and Bravo for the neuron team !!
ted
Site Admin
Posts: 5795
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: implementing recovery from inactivation

Post by ted »

Thanks for using NEURON in your research. Please be sure to tell us when you publish anything that reports work that was done with NEURON, so we can add it to the Bibliography
http://www.neuron.yale.edu/neuron/stati ... ednrn.html
Post Reply