calling a function from NMODL
calling a function from NMODL
Hi Ted,
I have a question about functions built into NMODL.
I am trying to incorporate channel noise in the HodgkinHuxley equations. At each time step, I need to use the singular value decomposition (svd) function. It seems that svd function is not built into NMODL. I am wondering if and how I can use python interfaced with .mod file, so that at each time step svd will be performed in python and the result will be sent to .mod files.
Thank you in advance for your help.
I have a question about functions built into NMODL.
I am trying to incorporate channel noise in the HodgkinHuxley equations. At each time step, I need to use the singular value decomposition (svd) function. It seems that svd function is not built into NMODL. I am wondering if and how I can use python interfaced with .mod file, so that at each time step svd will be performed in python and the result will be sent to .mod files.
Thank you in advance for your help.

 Posts: 229
 Joined: Fri Nov 28, 2008 3:38 pm
 Location: Yale School of Public Health
Re: calling a function from NMODL
It would be cleaner to have the Python function run at every timestep using a CVode.extra_scatter_gather and then just transfer values to NMODL (using pointers or directly; both should work.)
You could probably get the call from NMODL to work, but that would require a VERBATIM block.
You could probably get the call from NMODL to work, but that would require a VERBATIM block.
Re: calling a function from NMODL
I think that CVode.extra_scatter_gather is fine (in the after scatter direction) for your purpose because, since you are dealing with channel noise, you are most likely using only the fixed step method. In general we need to introduce something like Cvode.before_step(python_callback)
which, for variable step as well as fixed step methods, would guarantee only one call to the callback per time step.
The latter would be equivalent to the following mod file that I use to do a per time step callback for LFP calculations.
which, for variable step as well as fixed step methods, would guarantee only one call to the callback per time step.
The latter would be equivalent to the following mod file that I use to do a per time step callback for LFP calculations.
Code: Select all
$ cat beforestep_py.mod
: Python callback from BEFORE STEP
NEURON {
POINT_PROCESS beforestep_callback
POINTER ptr
}
ASSIGNED {
ptr
}
INITIAL {
}
VERBATIM
extern int (*nrnpy_hoccommand_exec)(Object*);
extern Object** hoc_objgetarg(int);
extern int ifarg(int);
extern void hoc_obj_ref(Object*);
extern void hoc_obj_unref(Object*);
ENDVERBATIM
BEFORE STEP {
:printf("beforestep_callback t=%g\n", t)
VERBATIM
{
Object* cb = (Object*)(_p_ptr);
if (cb) {
(*nrnpy_hoccommand_exec)(cb);
}
}
ENDVERBATIM
}
PROCEDURE set_callback() {
VERBATIM
Object** pcb = (Object**)(&(_p_ptr));
if (*pcb) {
hoc_obj_unref(*pcb);
*pcb = (Object*)0;
}
if (ifarg(1)) {
*pcb = *(hoc_objgetarg(1));
hoc_obj_ref(*pcb);
}
ENDVERBATIM
}
Re: calling a function from NMODL
Forgot to show the usage of the mod file... Note that the callback happens at the end of finitialize and at the beginning of a step.
Code: Select all
nrnivmodl beforestep_py.mod
Code: Select all
$ cat temp.py
from neuron import h
#enough model for the point process to have a home
s = h.Section()
# and for something interesting to happen (an hh action potential)
s.L = 10
s.diam = 10
s.insert("hh")
ic = h.IClamp(s(.5))
ic.delay = 0.1
ic.dur = 1.0
ic.amp = 0.5
def callback():
print("callback t=%g"%h.t)
pycb = h.beforestep_callback(s(.5))
pycb.set_callback(callback)
h.finitialize()
h.fadvance()
h.fadvance()
Code: Select all
$ python temp.py
callback t=0
callback t=0
callback t=0.025
Re: calling a function from NMODL
Hi all,
Thank you very much for your responses. I am using the fixed step method.
The code below is my potassium channel mod file. z1 is the noise and it is POINTER variable. At each time step, I want to have Python calculate z1 using alphan, betan, and n.
The following is my python code. The “findnoiseterm” function calculates the noise. I used h.CVode().extra_scatter_gather so that at every time step the “findnoiseterm” function is run. But I have a question about setting the POINTER variable (z1) to the output of “findnoiseterm” function. I used:
h.setpointer(_ref_hocvar, 'POINTER_name', nrn.Mechanism_object)
But I do not know what I need to write for _ref_hocvar. Since the reference is the output of the “findnoiseterm” function not a hoc variable.
Thank you for your help.
Thank you very much for your responses. I am using the fixed step method.
The code below is my potassium channel mod file. z1 is the noise and it is POINTER variable. At each time step, I want to have Python calculate z1 using alphan, betan, and n.
Code: Select all
NEURON {
SUFFIX Kv1
USEION k READ ek WRITE ik
RANGE gk, gbar, ik, ninf, taun
RANGE alphan, betan
POINTER z1
}
UNITS {
(mV) = (millivolt)
(mA) = (milliamp)
(nA) = (nanoamp)
(pA) = (picoamp)
(S) = (siemens)
(nS) = (nanosiemens)
(pS) = (picosiemens)
(um) = (micron)
(molar) = (1/liter)
(mM) = (millimolar)
}
CONSTANT {
q10 = 3
ca = 0.10 (1/ms) : 0.12889 as original
cva = 50 (mV)
cka = 33.90877 (mV)
cb = 0.12889 (1/ms)
cvb = 50 (mV)
ckb = 7.42101 (mV)
}
PARAMETER {
v (mV)
celsius (degC)
gbar = 0.011 (mho/cm2) <0,1e9>
}
ASSIGNED {
ik (mA/cm2)
ek (mV)
gk (mho/cm2)
ninf
taun (ms)
alphan (1/ms)
betan (1/ms)
qt
z1
}
STATE { n }
INITIAL {
qt = q10^((celsius22 (degC))/10 (degC))
rates(v)
n = ninf
z1=0
}
BREAKPOINT {
SOLVE states METHOD cnexp
gk = gbar * (n^4 + z1)
ik = gk * (v  ek)
}
DERIVATIVE states {
rates(v)
n' = (ninfn)/taun
}
PROCEDURE rates(v (mV)) {
alphan = alphanfkt(v)
betan = betanfkt(v)
ninf = alphan/(alphan+betan)
taun = 1/(qt*(alphan + betan))
}
FUNCTION alphanfkt(v (mV)) (1/ms) {
alphanfkt = ca * exp((v+cva)/cka)
}
FUNCTION betanfkt(v (mV)) (1/ms) {
betanfkt = cb * exp((v+cvb)/ckb)
}
h.setpointer(_ref_hocvar, 'POINTER_name', nrn.Mechanism_object)
But I do not know what I need to write for _ref_hocvar. Since the reference is the output of the “findnoiseterm” function not a hoc variable.
Code: Select all
from __future__ import division
from numpy import *
from pylab import *
from scipy import linalg
from neuron import h, gui
h.load_file("mosinit.hoc")
KHat = zeros((5)) # Initial values set to 0
# Computing matrix square roots with SVD
def mysqrtm(D):
u,s,v = linalg.svd(D)
S = u*sqrt(s)*v
return S
def DK(alphan, betan, X):
return 1/20 * array(( (4*alphan*X[0] + betan*X[1],(4*alphan*X[0] + betan*X[1]) , 0, 0, 0),
((4*alphan*X[0]+betan*X[1]),(4*alphan*X[1]+(3*alphan+ betan)*X[1] + 2*betan*X[2]), (2*betan*X[2] + 3*alphan*X[1]), 0,0),
(0, (2*betan*X[2]+3*alphan*X[1]), (3*alphan*X[1]+(2*alphan+2*betan)*X[2]+3*betan*X[3]),(3*betan*X[3]+2*alphan*X[2]), 0),
(0, 0, (3*betan*X[3]+2*alphan*X[2]), (2*alphan*X[2]+(alphan+3*betan)*X[3]+4*betan*X[4]),(4*betan*X[4]+alphan*X[3])),
(0, 0, 0, (4*betan*X[4]+alphan*X[3]), (alphan*X[3]+4*betan*X[4]) ) ))
def findnoiseterm(alphan, betan, n):
global KHat
KNoise = randn(5,1)
AK = array(( (4*alphan, betan, 0, 0, 0),
(4*alphan, 3*alphanbetan, 0, 0, 0),
(0, 3*alphan, 2*alphan2*betan, 3*betan, 0),
(0, 0, 2*alphan, alphan3*betan, 4*betan),
(0, 0, 0, alphan, 4*betan) ))
SK = lambda alphan, betan ,X: mysqrtm(DK(alphan, betan ,X))
KBar = array(((1n)**4, 4*n*(1n)**3, 6*n**2*(1n)**2, 4*n**3*(1n), n**4))
KHat += dt*dot(AK,KHat) + sqrt(dt)*dot(SK(alphan, betan ,KBar),KNoise[:,1])
noiseterm = KHat[1]
return noiseterm
for sec in h.allsec():
if h.ismembrane('Kv1',sec = sec):
for seg in sec:
recording_callback = (findnoiseterm , (seg.alphan_Kv1 , seg.betan_Kv1, seg.n_Kv1 ))
h.CVode().extra_scatter_gather(1, recording_callback)
h.setpointer(???, 'z1', sec().Kv1)
Re: calling a function from NMODL
My opinion is that findnoiseterm should not be the callback but instead the callback should be a wrapper python function that for every instance of
Kv1, the wrapper would call findnoiseterm and assign the return value into z1_Kv1 (which should be a regular RANGE variable instead of a POINTER.)
If the iteration over all z1_Kv1 instances has performance issues,
https://nrn.readthedocs.io/en/latest/py ... ector.html can help.
Kv1, the wrapper would call findnoiseterm and assign the return value into z1_Kv1 (which should be a regular RANGE variable instead of a POINTER.)
If the iteration over all z1_Kv1 instances has performance issues,
https://nrn.readthedocs.io/en/latest/py ... ector.html can help.
Re: calling a function from NMODL
Thank you for your help. I appreciate it. I tried to implement what you suggested, and it works. In the following code the callback is a wrapper python function and at each time step, z1 is calculated and transferred to the .mod file.
The only problem is that at each time step I want my python code to use the updated alphan, betan, and n for z1 calculation. I thought that
recording_callback = (wrapper, (seg.alphan_Kv1 , seg.betan_Kv1, seg.n_Kv1 , 0 ))
h.CVode().extra_scatter_gather(0, recording_callback)
would update alphan, betan and n per time step. But it seems like it does not. I am wondering how the updated alphan, betan and n from the mod file can be transferred to the python code per time step.
Thanks again.
The only problem is that at each time step I want my python code to use the updated alphan, betan, and n for z1 calculation. I thought that
recording_callback = (wrapper, (seg.alphan_Kv1 , seg.betan_Kv1, seg.n_Kv1 , 0 ))
h.CVode().extra_scatter_gather(0, recording_callback)
would update alphan, betan and n per time step. But it seems like it does not. I am wondering how the updated alphan, betan and n from the mod file can be transferred to the python code per time step.
Thanks again.
Code: Select all
from neuron import h, gui
h.load_file("mosinit.hoc")
KHat = zeros((5,5)) # Initial values set to 0
def findnoiseterm(func):
def inner(alphan,betan,n,noise):
global KHat
KNoise = randn(5,5)
AK = array(( (4*alphan, betan, 0, 0, 0),
(4*alphan, 3*alphan, 0, 0, 0),
(0, 3*alphan, 2*alphan, 3*betan, 0),
(0, 0, 2*alphan, alphan3*betan, 4*betan),
(0, 0, 0, alphan, 4*betan) ))
X = array(((1n)**4, 4*n*(1n)**3, 6*n**2*(1n)**2, 4*n**3*(1n), n**4))
DK= 1/20 * array(( (4*alphan*X[0] + betan*X[1],(4*alphan*X[0]) , 0, 0, 0),
((4*alphan*X[0]+betan*X[1]),(4*alphan*X[1]+(3*alphan+ betan)*X[1]), (2*betan*X[2] + 3*alphan*X[1]), 0,0),
(0, (2*betan*X[2]), (3*alphan*X[1]+(2*alphan)*X[2]+3*betan*X[3]),(3*betan*X[3]), 0),
(0, 0, (3*betan*X[3]), (2*alphan*X[2]+(3*betan)*X[3]+4*betan*X[4]),(4*betan*X[4])),
(0, 0, 0, (4*betan*X[4]), (alphan*X[3]) ) ))
u,s,v = linalg.svd(DK)
SK = u*sqrt(s)*v
KHat += h.dt*dot(AK,KHat) + sqrt(h.dt)*dot(SK,KNoise[:,:])
noise = KHat[4,4]
return noise, func(alphan,betan,n, noise)
return inner
@findnoiseterm
def wrapper(alphan,betan,n, noise):
seg.z1_Kv1=noise
for sec in h.allsec():
if h.ismembrane('Kv1',sec = sec):
for seg in sec:
recording_callback = (wrapper, (seg.alphan_Kv1 , seg.betan_Kv1, seg.n_Kv1 , 0 ))
h.CVode().extra_scatter_gather(0, recording_callback)
h.finitialize(65)
h.fadvance()
Re: calling a function from NMODL
I see you have a
but in your
the arguments to wrapper are just constant values computed at the time the recording_callback assignment statement is executed.
In your mod file, alpahn and betan are RANGE varialbles that are updated when rates(v) is called and those rates depend on the value of v and all the mentioned parameters in alphanfkt and betanfkt functions. Those are updated every time step, so just access their values directly when wrapper is called and calculate the noise .
Code: Select all
def wrapper(alphan, betan, n, noise)
Code: Select all
recording_callback = (wrapper, (seg.alphan_Kv1 , seg.betan_Kv1, seg.n_Kv1 , 0 ))
In your mod file, alpahn and betan are RANGE varialbles that are updated when rates(v) is called and those rates depend on the value of v and all the mentioned parameters in alphanfkt and betanfkt functions. Those are updated every time step, so just access their values directly when wrapper is called and calculate the noise .
Re: calling a function from NMODL
Thanks a lot. Appreciate it.
Re: calling a function from NMODL
My code works well for single compartment model. But for a multicompartment model in which more than a section includes Kv1 mechanism, it only calculates noise (i.e. z1) for the last segment of the last compartment. I have
but it does not calculate the noise for each segment of each section that includes the Kv1 mechanism.
I appreciate any help in this regard.
Code: Select all
for sec in h.allsec():
if h.ismembrane('Kv1',sec = sec):
for seg in sec.allseg():
recording_callback = (wrapper, ( 0 ))
h.CVode().extra_scatter_gather(0, recording_callback)
h.fadvance
I appreciate any help in this regard.
Re: calling a function from NMODL
I hope you are doing a sanity check when you have problems by printing some details. e.g
Anyway, one problem I see is that
includes in the iteration sec(0) and sec(1) which are zero area nodes and do not contain Kv1. The line should read
But the fundamental problem (assuming I understand the context which you did not explicitly provide) is that you registered the callback many times
without specifying the segment the callback is supposed to handle. (Did that actually get called many times or just once per time step?). Anyway
I would have expected you to register the callback once, and in the def wrapper, do the iteration over the segments that contain Kv1.
If you have a lot of sections,only some of which contain Kv1, you can make a list of those sections, or segments, or Kv1 instances and pass it as an arg
to wrapper (or just make it a global).
e.g:
Code: Select all
print(seg.sec.name(), seg.x, seg.v, seg.Kv1.alphan)
Code: Select all
for seg in sec.allseg():
Code: Select all
for seg in sec:
without specifying the segment the callback is supposed to handle. (Did that actually get called many times or just once per time step?). Anyway
I would have expected you to register the callback once, and in the def wrapper, do the iteration over the segments that contain Kv1.
If you have a lot of sections,only some of which contain Kv1, you can make a list of those sections, or segments, or Kv1 instances and pass it as an arg
to wrapper (or just make it a global).
e.g:
Code: Select all
list_of_Kv1_instances = [ seg.Kv1 for sec in h.allsec() if h.ismembrane("Kv1", sec=sec) for seg in sec]
recording_callback = (wrapper, list_of_kv1_instances) # I think one arg does not need to be a tuple but I did not try it.