Inserting ion channels

[aburke]

When I try to run my hoc code in Neuron an error appears when it tries inserting the ion mechanisms I wrote using NMODL. This is the error that appears:

Code: Select all

loading membrane mechanisms from /Users/audreyburke/Documents/Grill/Copy (2) of Otsuka STN/i686/.libs/libnrnmech.so
dlopen failed - 
dlopen(/Users/audreyburke/Documents/Grill/Copy (2) of Otsuka STN/i686/.libs/libnrnmech.so, 2): Symbol not found: _ss
  Referenced from: /Users/audreyburke/Documents/Grill/Copy (2) of Otsuka STN/i686/.libs/libnrnmech.so
  Expected in: flat namespace

        1 
/Applications/NEURON-5.9/nrn/i686/bin/nrniv.app/Contents/MacOS/nrniv: parse error
 in start_STN.hoc near line 6
        insert kv3 ek=-90
           ^

Here is my hoc code:

Code: Select all

load_file("nrngui.hoc")

create stn
stn	{nseg=1 diam=5.5 L= 100/(PI*diam)
	Ra= 300 cm=100
	insert kv3 ek=-90
	insert atypeK ek=-90
	insert nacurrent ena=60
	insert LlikeCa
	insert TtypeCa
	insert Calc
	insert caK ek=-90
	insert leak el=-60}
access stn

dt = 0.0125
tstop = 1000

And this is one of my ion mechanisms:

Code: Select all

NEURON {
	SUFFIX kv3
	USEION k READ ek WRITE ik
	RANGE gk_bar, ik, n_inf
}

UNITS {
	(S)  = (siemens)
	(mV) = (millivolts)
	(mA) = (milliamps)
	(mmA) = (microamps)
	(M) = (moles/liter)
}

PARAMETER {
gk_bar      =   57   : mS/cm^2,  K conductance
ek          =  -90   : mV,       K Nernst potential
theta_n     =  -41   : mV,       K half-activation voltage
k_n         =  -14   : mV,       K activation slope
tau0_n      =    0   : ms,       K time constant activation minimum
tau1_n      =   11   : ms,       K time constant activation gain
theta_tau1n =  -40   : mV,       K time constant - half-activation voltage
sigma1_n    =  -40   : mV,       K time constant - activation slope 
theta_tau2n =  -40   : mV,       K time constant + half-activation voltage
sigma2_n    =   50   : mV,       K time constant + activation slope 
}


ASSIGNED {
	v  (mv)
	ik (mA/cm2)
	n_inf
	tau_n
}

STATE { n }

BREAKPOINT {
	SOLVE states METHOD cnexp
	ss()
	ik = gk_bar * n^4 * (v - ek)
}



INITIAL {
	ss()
	n = n_inf
}



DERIVATIVE states {
	n'   = (n_inf - n) / tau_n
}


PROCEDURE ss() {
	n_inf = w_inf(v,theta_n,k_n)
	tau_n = tau_ws(v,tau0_n,tau1_n,theta_tau1n,sigma1_n,theta_tau2n,sigma2_n)
}

FUNCTION w_inf (Vm (mV), theta, kk) {
	UNITSOFF
	w_inf = 1 / (1 + exp((Vm-theta)/kk))
	UNITSON
}

FUNCTION tau_wf (Vm (mv), tau0, tau1, theta_tau, sigma) {
	UNITSOFF
	tau_wf = tau0 + tau1 / (1 + exp((theta_tau-Vm)/sigma))
	UNITSON
}

FUNCTION tau_ws (Vm (mV), tau0, tau1, theta_tau1, sigma1, theta_tau2, sigma2) {
	UNITSOFF
	tau_ws = tau0 + tau1 / (exp((theta_tau1-Vm)/sigma1) + exp((theta_tau2-Vm)/sigma2))
	UNITSON
}

Any help would be greatly appreciated. Thanks!

[ted]

The clue is in the error message

Code: Select all

loading membrane mechanisms from /Users/audreyburke/Documents/Grill/Copy (2) of Otsuka STN/i686/.libs/libnrnmech.so
dlopen failed - 
dlopen(/Users/audreyburke/Documents/Grill/Copy (2) of Otsuka STN/i686/.libs/libnrnmech.so, 2): Symbol not found: _ss
  Referenced from: /Users/audreyburke/Documents/Grill/Copy (2) of Otsuka STN/i686/.libs/libnrnmech.so

which means NEURON can't find the file that contains the compiled mechanisms in the
usual place (a subdirectory called .i686 that sits inside the directory where your hoc file
is located). This bit

Code: Select all

  Expected in: flat namespace

means you didn't tell NEURON to look anywhere else than the usual place.

Put the necessary mod files in the same directory as the hoc file, compile them there
(on a Mac I think just drag and drop that directory onto the mknrndll icon will do it), and
your hoc file should launch without complaint.

[aburke]

Thanks for getting back to me.

I had compiled my files like you described above. Also, other Neuron code functions properly on my computer. Do you think a problem with my code itself could cause this error message?

[ted]

Well, the mod file for kv3 is full of units inconsistencies, some serious--you need to check
it with modlunit and fix them. I'd bet that the mod files for the other mechanisms are similarly
afflicted, but that wouldn't cause the error message.

If you want to try to diagnose the problem, comment out all of the insert statements in the
hoc file, then see if NEURON will load it without complaint. If so, exit NEURON,
uncomment just one insert statement, and try again. Repeat until you run into a
problem, and when the problem occurs, note the error message.

By the way, are the mod files from a model in ModelDB?

[aburke]

I found an error in one of my MOD files by commenting out all the insert statements, and then adding them back in one by one. I had ss() in one of my MOD files that did not have a PROCEDURE statement. That particular MOD file calcuates the internal calcium concentration and calcium potential as the program runs. As I looked at the error message again, I noticed it read:

Code: Select all

2): Symbol not found: _ss 

I've been using modlunit to check the units in my MOD files. This is my first time using modlunit, so I'm still somewhat unfamiliar with its error messages. When I checked one MOD file using modlunit (the code that calculates eca and cai), the following message appeared:

Code: Select all

Checking units of /Users/audreyburke/Desktop/untitled folder/Calc.mod
The previous primary expression with units: 1000 /m3-sec
is missing a conversion factor and should read:
  (0.0001)*()
 at line 50 in file /Users/audreyburke/Desktop/untitled folder/Calc.mod
cai' = -((2*ica)/(diam*F)) - (K_Ca * cai)<<ERROR>>

Here's the code:

Code: Select all

NEURON {
	SUFFIX Calc
	USEION ca READ ica,cai	WRITE cai, eca
	RANGE cai0, K_Ca
}

UNITS{
	(mV) = (millivolt)
	(mA) = (milliamp)
	(molar) = (1/liter)
	(mM) = (millimolar)
	F = (faraday) (coulomb)
	(um) = (micron)
}

PARAMETER {
	diam (um)
	cao = 2 (mM)
	cai0 = 1E-4 (mM)
	K_Ca = 2 (/ms)
	RTF = 25.7 (mV)
}

ASSIGNED{
ica     (mA/cm2)
eca	(mV)
}

STATE {
	cai (mM)
}

BREAKPOINT {
SOLVE state METHOD cnexp
eca = RTF/2 * log(cao/cai)
}

INITIAL {
cai = cai0
}

DERIVATIVE state {
cai' = -((2*ica)/(diam*F)) - (K_Ca * cai)
}

When the message mentions the "previous primary expression", is it talking about the first term or the second term of the equation? I've tried adding the conversion factor to each side of the equation, but when I check the file with modlunit again, it tells me that another conversion factor is needed.

Thanks for all of your help!

[aburke]

I forgot to mention that these MOD files are not from MODELDB, but I wrote them myself. I am trying to model the STN cell from Otsuka (2004), "Conductance-Based Model of the Voltage-Dependent Generation of a Plateau Potential in Subthalamic Neurons", J Neurophysiol 92: 255-264

[ted]

I forgot to mention that these MOD files are not from MODELDB

Good, because if it had been from ModelDB, I would have been very concerned. Units
errors are the rule for models that are in the development stage, but one hopes that
serious ones are resolved long before publication.

I am trying to model the STN cell from Otsuka (2004), "Conductance-Based Model of the Voltage-Dependent Generation of a Plateau Potential in Subthalamic Neurons", J Neurophysiol 92: 255-264

It can be pretty difficult to reproduce published models without the authors' source code.
That paper isn't in NEURON's bibliography, so I assume that they implemented the model
with some other simulator or programming language. Even so, having the original code can
resolve many uncertainties that arise from omissions or errors in the published paper. It's
not in ModelDB; is it available from the journal as supplementary material, have they
posted it on their WWW site, or do you think they might send it to you if you asked?

[ted]

This is my first time using modlunit, so I'm still somewhat unfamiliar with its error messages.

Familiarity isn't necessarily beneficial; they can be almost impenetrably terse.

When the message mentions the "previous primary expression", is it talking about the first term or the second term of the equation? I've tried adding the conversion factor to each side of the equation, but when I check the file with modlunit again, it tells me that another conversion factor is needed.

Rather than try to decipher these oracular utterances, it can be much easier to simplify
the expression in order to deal with problems one at a time.

First iteration: change the body of the DERIVATIVE block to

Code: Select all

: cai' = -((2*ica)/(diam*F)) - (K_Ca * cai)
cai' = -((2*ica)/(diam*F))

This eliminates the 2nd term on the RHS of the equation.
The error message becomes

Code: Select all

The previous primary expression with units: 1+07 /m3-sec
is missing a conversion factor and should read:
  (10000)*()
 at line 44 in file calc.mod
cai' = -((2*ica)/(diam*F))<<ERROR>>

Changing the simplified equation to

Code: Select all

cai' = -(10000)*((2*ica)/(diam*F))

gets rid of the error msg.

Second iteration: restore the 2nd term to the RHS of the equation, which becomes

Code: Select all

cai' = -(10000)*((2*ica)/(diam*F)) - (K_Ca * cai)

modlunit now has no complaints.

One comment: generally there is no need to calculate eca in the BREAKPOINT block,
or to WRITE it in the NEURON block. Whenever you insert a mechanism that WRITEs
the internal or external concentration of an ionic species, NEURON will automatically
compute the equilibrium potential. I suppose it's OK to do as you have done if you are
trying to capture every last aspect of a model that was implemented with some other
software than NEURON, but you should check to make sure that eca is what you
expect, and not what NEURON computes from its own assumptions about celsius,
R, and F.

[aburke]

I've fixed all the units in my MOD files. All the files compile succesfully and load successfully when I run the HOC file. However, when I try running the program, the following error appears:

Code: Select all

loading membrane mechanisms from /Users/audreyburke/Desktop/FinalSTN/i686/.libs/libnrnmech.so
Additional mechanisms from files
 Calc.mod LlikeCa.mod TtypeCa.mod atypeK.mod caK.mod kv3.mod leak.mod nacurrent.mod
        1 
oc>exp(813.75) out of range, returning exp(700)
A math function was called that returned an out of range value
errno=34 during call to mechanism LlikeCa
/Applications/NEURON-5.9/nrn/i686/bin/nrniv.app/Contents/MacOS/nrniv: errno set during call to INITIAL block
 near line 11
 {run()}
        ^
exp(814.625) out of range, returning exp(700)
A math function was called that returned an out of range value
errno=34 during call to mechanism caK
/Applications/NEURON-5.9/nrn/i686/bin/nrniv.app/Contents/MacOS/nrniv: errno set during call to INITIAL block
 near line 11
 ^
No more errno warnings during this execution

Here's the HOC code:

Code: Select all

// This file loads the appropriate parameters for the STN model

create soma

load_file("nrngui.hoc")

dt = 0.0125
tstop = 1000

soma{
diam = 5.5
L = 100/(PI*diam)
nseg = 1
Ra = 300
cm = 100
Celsius = 25
insert kv3 ek = -90
insert atypeK ek = -90
insert nacurrent ena = 60
insert LlikeCa
insert TtypeCa
insert Calc
insert caK ek = -90
insert leak el = -60
}

and the MOD file that the error references:

Code: Select all

NEURON {
	SUFFIX LlikeCa
	USEION ca READ eca,cai WRITE ica
	RANGE gL_bar, c_inf, d1_inf, d2_inf
}

UNITS {
	(S)  = (siemens)
	(mS) = (millisiemens)
	(mV) = (millivolts)
	(mA) = (milliamps)
	(molar) = (1/liter)
	(mM) = (millimolar)
}

PARAMETER {
eca	(mV)
cai     (mM)
gL_bar       =    57   (mS/cm2) : L-like conductance
theta_c      = -30.6   (mV)     : L-like half-activation voltage 
k_c          =    -5   (mV)     : L-like activation slope
tau0_c       =    45   (ms)     : L-like time constant activation minimum
tau1_c       =    10   (ms)     : L-like time constant activation gain
theta_tau1c  =   -27   (mV)     : L-like time constant - half-activation voltage
sigma1_c     =   -20   (mV)     : L-like time constant - activation slope 
theta_tau2c  =   -50   (mV)     : L-like time constant + half-activation voltage
sigma2_c     =    15   (mV)     : L-like time constant + activation slope
theta_d1     =   -60   (mV)     : L-like half-inactivation voltage
k_d1         =   7.5   (mV)     : L-like inactivation slope
tau0_d1      =   400   (ms)     : L-like time constant inactivation minimum
tau1_d1      =   500   (ms)     : L-like time constant inactivation gain
theta_tau1d1 =   -40   (mV)     : L-like time constant - half-inactivation voltage
sigma1_d1    =   -15   (mV)     : L-like time constant - inactivation slope
theta_tau2d1 =   -20   (mV)     : L-like time constant + half-inactivation voltage
sigma2_d1    =    20   (mV)     : L-like time constant + inactivation slope
theta_d2     =   0.1   (mV)     : L-like Ca-dependent half-inactivation voltage
k_d2         = -0.08   (mV)     : L-like Ca-dependent inactivation slope
tau_d2       =   130   (ms)     : L-like Ca-dependent inactivation time constant
}


ASSIGNED {
	v  (mV)
	ica (mA/cm2)
	c_inf
	d1_inf
	d2_inf
	tau_c (ms)
	tau_d1 (ms)
}

STATE { c d1 d2 }

BREAKPOINT {
	SOLVE states METHOD cnexp
	ss()
	ica = gL_bar * c^2 * d1 * d2 * (v - eca) * (0.001)
}



INITIAL {
	ss()
	c = c_inf
	d1 = d1_inf
	d2 = d2_inf
}



DERIVATIVE states {
	c'   = (c_inf - c) / tau_c
	d1'  = (d1_inf - d1)/ tau_d1
	d2'  = (d2_inf - d2)/ tau_d2
}



PROCEDURE ss() {
	
	c_inf = 1 / (1 + exp((v-theta_c)/k_c))
	d1_inf = 1 / (1 + exp((v-theta_d1)/k_d1))
	d2_inf = 1 / (1 + exp((v-theta_d2)/k_d2))
	tau_c = tau0_c + tau1_c / (exp((theta_tau1c-v)/sigma1_c) + exp((theta_tau2c-v)/sigma2_c))
	tau_d1 = tau0_d1 + tau1_d1 / (exp((theta_tau1d1-v)/sigma1_d1) + exp((theta_tau2d1-v)/sigma2_d1))
}

Any idea why this error might be occuring? I tried reducng DT, but the error still appeared.

Thanks!

[ted]

First, something that isn't the cause of the problem you're seeing, but which should be
rectified nonetheless:
ss() should be called at the start of the DERIVATIVE block, not after the SOLVE statement
in the BREAKPOINT block--unless there is some reason for deferring calculation of the
steady state values and rate constants that are to be used to advance to the next time step,
until _after_ the advance has been taken.

A good way to diagnose numerical overflow problems is to embed printf statements that
report calculated values. The most attractive target in this case is right at the end of
PROCEDURE ss(), e.g. like so

Code: Select all

printf("v %g, c_inf %g, d1_inf %g, d2_inf %g\n", v, c_inf, d1_inf, d2_inf)
printf("      tau_c %g, tau_d1 %g\n", tau_c, tau_d1)

Then use the GUI to create a single compartment model that has pas and just this
particular mechanism, and click on the Init button, which reveals . . .