Two channels don't like each other

NMODL and the Channel Builder.
Post Reply
MBeining
Posts: 35
Joined: Thu Apr 03, 2014 8:18 am

Two channels don't like each other

Post by MBeining »

Hi there,
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)
}
ted
Site Admin
Posts: 5810
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Two channels don't like each other

Post by ted »

These problems can be context sensitive. Do you have a toy model that reproduces the error?
Post Reply