I have some issue with two of my channel models.
When I introdruce one of the two channels below, my model works fine, however when I introdruce both channels to my model, cvode gives this error 7 (the corrector
convergence failed repeatedly or with |h| = hmin.) or without cvode, the voltage is a NaN value.
Anyone sees why they do not like each other?
Thank you in advance,
Marcel Beining
Code: Select all
Ca channels (T,N,L-type)
NEURON {
SUFFIX Caold
USEION ca WRITE ica, cai
RANGE ica, a, b, c, d, e, gtcabar, gncabar, glcabar, gtca, gnca, glca, e_ca
GLOBAL ca0, cao, tau, celsius
}
UNITS {
(molar) = (1/liter)
(mM) = (millimolar)
(mV) = (millivolt)
(mA) = (milliamp)
(S) = (siemens)
B = .26 (mM-cm2/mA-ms)
F = (faraday) (coulomb)
R = (k-mole) (joule/degC)
}
PARAMETER {
ca0 = .00007 (mM) : initial calcium concentration inside
cao = 2 (mM) : calcium concentration outside
tau = 9 (ms)
gtcabar = .01 (S/cm2) : maximum permeability
gncabar = .01 (S/cm2)
glcabar = .01 (S/cm2)
}
ASSIGNED {
v (mV)
e_ca (mV)
ica (mA/cm2)
gtca (S/cm2)
gnca (S/cm2)
glca (S/cm2)
celsius (degC)
}
STATE { cai (mM) a b c d e}
BREAKPOINT {
SOLVE state METHOD cnexp
e_ca = (1000)*(celsius+273.15)*R/(2*F)*log(cao/cai)
gtca = gtcabar*a*a*b
gnca = gncabar*c*c*d
glca = glcabar*e*e
ica = (gtca+gnca+glca)*(v - e_ca)
}
DERIVATIVE state { : exact when v held constant; integrates over dt step
cai' = -B*ica-(cai-ca0)/tau
a' = alphaa(v)*(1-a)-betaa(v)*a
b' = alphab(v)*(1-b)-betab(v)*b
c' = alphac(v)*(1-c)-betac(v)*c
d' = alphad(v)*(1-d)-betad(v)*d
e' = alphae(v)*(1-e)-betae(v)*e
}
INITIAL {
cai = ca0
a = alphaa(v)/(alphaa(v)+betaa(v))
b = alphab(v)/(alphab(v)+betab(v))
c = alphac(v)/(alphac(v)+betac(v))
d = alphad(v)/(alphad(v)+betad(v))
e = alphae(v)/(alphae(v)+betae(v))
}
FUNCTION alphaa(v (mV)) (/ms) {
alphaa = f(2,0.1,v,29.26) :old 19.26
}
FUNCTION betaa(v (mV)) (/ms) {
betaa = exponential(0.009,-0.045393,v,10) :old 0
}
FUNCTION alphab(v (mV)) (/ms) {
alphab = exponential(1e-6,-0.061501,v,10) :old 0
}
FUNCTION betab(v (mV)) (/ms) {
betab = logistic(1,-0.1,v,39.79) :old 29.79
}
FUNCTION alphac(v (mV)) (/ms) {
alphac = f(1.9,0.1,v,29.88) :old 19.88
}
FUNCTION betac(v (mV)) (/ms) {
betac = exponential(0.046,-0.048239,v,10) :old 0
}
FUNCTION alphad(v (mV)) (/ms) {
alphad = exponential(1.6e-4,-0.020661,v,10) :old 0
}
FUNCTION betad(v (mV)) (/ms) {
betad = logistic(1,-0.1,v,49) :old 39
}
FUNCTION alphae(v (mV)) (/ms) {
alphae = f(156.9,0.1,v,91.5) :old 81.5
}
FUNCTION betae(v (mV)) (/ms) {
betae = exponential(0.29,-0.092081,v,10) :old 0
}
FUNCTION f(A, k, v (mV), D) (/ms) {
LOCAL x
UNITSOFF
x = k*(v-D)
if (fabs(x) > 1e-6) {
f = A*x/(1-exp(-x))
}else{
f = A/(1-0.5*x)
}
UNITSON
}
FUNCTION logistic(A, k, v (mV), D) (/ms) {
UNITSOFF
logistic = A/(1+exp(k*(v-D)))
UNITSON
}
FUNCTION exponential(A, k, v (mV), D) (/ms) {
UNITSOFF
exponential = A*exp(k*(v-D))
UNITSON
}
Code: Select all
: Eight state kinetic sodium channel gating scheme
: Modified from k3st.mod, chapter 9.9 (example 9.7)
: of the NEURON book
: 12 August 2008, Christoph Schmidt-Hieber
:
: accompanies the publication:
: Schmidt-Hieber C, Bischofberger J. (2010)
: Fast sodium channel gating supports localized and efficient
: axonal action potential initiation.
: J Neurosci 30:10233-42
NEURON {
SUFFIX na8st
USEION na READ ena WRITE ina
GLOBAL vShift, vShift_inact, maxrate
RANGE vShift_inact_local
RANGE g, gbar
RANGE a1_0, a1_1, b1_0, b1_1, a2_0, a2_1
RANGE b2_0, b2_1, a3_0, a3_1, b3_0, b3_1
RANGE bh_0, bh_1, bh_2, ah_0, ah_1, ah_2
}
UNITS { (mV) = (millivolt) }
: initialize parameters
PARAMETER {
gbar = 33 (millimho/cm2)
a1_0 = 5.142954478051616e+01 (/ms)
a1_1 = 7.674641248142576e-03 (/mV)
b1_0 = 9.132202467321037e-03 (/ms)
b1_1 = 9.342823457307300e-02 (/mV)
a2_0 = 7.488753944786941e+01 (/ms)
a2_1 = 2.014613733367395e-02 (/mV)
b2_0 = 6.387047323688771e-03 (/ms)
b2_1 = 1.501806374396736e-01 (/mV)
a3_0 = 3.838866325780059e+01 (/ms)
a3_1 = 1.253027842782742e-02 (/mV)
b3_0 = 3.989222258297797e-01 (/ms)
b3_1 = 9.001475021228642e-02 (/mV)
bh_0 = 1.687524670388565e+00 (/ms)
bh_1 = 1.210600094822588e-01
bh_2 = 6.827857751079400e-02 (/mV)
ah_0 = 3.800097357917129e+00 (/ms)
ah_1 = 4.445911330118979e+03
ah_2 = 4.059075804728014e-02 (/mV)
vShift = 12 (mV) : shift to the right to account for Donnan potentials
: 12 mV for cclamp, 0 for oo-patch vclamp simulations
vShift_inact = 10 (mV) : global additional shift to the right for inactivation
: 10 mV for cclamp, 0 for oo-patch vclamp simulations
vShift_inact_local = 0 (mV) : additional shift to the right for inactivation, used as local range variable
maxrate = 8.00e+03 (/ms) : limiting value for reaction rates
: See Patlak, 1991
}
ASSIGNED {
v (mV)
ena (mV)
g (millimho/cm2)
ina (milliamp/cm2)
a1 (/ms)
b1 (/ms)
a2 (/ms)
b2 (/ms)
a3 (/ms)
b3 (/ms)
ah (/ms)
bh (/ms)
}
STATE { c1 c2 c3 i1 i2 i3 i4 o }
BREAKPOINT {
SOLVE kin METHOD sparse
g = gbar*o
ina = g*(v - ena)*(1e-3)
}
INITIAL { SOLVE kin STEADYSTATE sparse }
KINETIC kin {
rates(v)
~ c1 <-> c2 (a1, b1)
~ c2 <-> c3 (a2, b2)
~ c3 <-> o (a3, b3)
~ i1 <-> i2 (a1, b1)
~ i2 <-> i3 (a2, b2)
~ i3 <-> i4 (a3, b3)
~ i1 <-> c1 (ah, bh)
~ i2 <-> c2 (ah, bh)
~ i3 <-> c3 (ah, bh)
~ i4 <-> o (ah, bh)
CONSERVE c1 + c2 + c3 + i1 + i2 + i3 + i4 + o = 1
}
: FUNCTION_TABLE tau1(v(mV)) (ms)
: FUNCTION_TABLE tau2(v(mV)) (ms)
PROCEDURE rates(v(millivolt)) {
LOCAL vS
vS = v-vShift
a1 = a1_0*exp( a1_1*vS)
a1 = a1*maxrate / (a1+maxrate)
b1 = b1_0*exp(-b1_1*vS)
b1 = b1*maxrate / (b1+maxrate)
a2 = a2_0*exp( a2_1*vS)
a2 = a2*maxrate / (a2+maxrate)
b2 = b2_0*exp(-b2_1*vS)
b2 = b2*maxrate / (b2+maxrate)
a3 = a3_0*exp( a3_1*vS)
a3 = a3*maxrate / (a3+maxrate)
b3 = b3_0*exp(-b3_1*vS)
b3 = b3*maxrate / (b3+maxrate)
bh = bh_0/(1+bh_1*exp(-bh_2*(vS-vShift_inact-vShift_inact_local)))
bh = bh*maxrate / (bh+maxrate)
ah = ah_0/(1+ah_1*exp( ah_2*(vS-vShift_inact-vShift_inact_local)))
ah = ah*maxrate / (ah+maxrate)
}