calling a function from NMODL

[Nooshin]

Hi Ted,
I have a question about functions built into NMODL.
I am trying to incorporate channel noise in the Hodgkin-Huxley 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.

[ramcdougal]

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.

[hines]

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.

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
}

[hines]

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

[Nooshin]

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.

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^((celsius-22 (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' = (ninf-n)/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)
}

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.

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*alphan-betan, 0, 0, 0),
        (0, 3*alphan, -2*alphan-2*betan, 3*betan, 0),
        (0, 0, 2*alphan, -alphan-3*betan, 4*betan),
        (0, 0, 0, alphan, -4*betan) ))
    
    SK = lambda alphan, betan ,X: mysqrtm(DK(alphan, betan ,X))
    
    KBar = array(((1-n)**4, 4*n*(1-n)**3, 6*n**2*(1-n)**2, 4*n**3*(1-n), 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)
    

Thank you for your help.

[hines]

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.

[Nooshin]

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.

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, -alphan-3*betan, 4*betan),
            (0, 0, 0, alphan, -4*betan) ))
     
        X = array(((1-n)**4, 4*n*(1-n)**3, 6*n**2*(1-n)**2, 4*n**3*(1-n), 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()

[hines]

I see you have a

Code: Select all

def wrapper(alphan, betan, n, noise)

but in your

Code: Select all

recording_callback = (wrapper, (seg.alphan_Kv1 , seg.betan_Kv1, seg.n_Kv1 , 0 ))

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 .

[Nooshin]

Thanks a lot. Appreciate it.

[Nooshin]

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

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

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.

[hines]

I hope you are doing a sanity check when you have problems by printing some details. e.g

Code: Select all

  print(seg.sec.name(), seg.x, seg.v, seg.Kv1.alphan)

Anyway, one problem I see is that

Code: Select all

for seg in sec.allseg():

includes in the iteration sec(0) and sec(1) which are zero area nodes and do not contain Kv1. The line should read

Code: Select all

for seg in sec:

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

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.