www.neuron.yale.edu

Posted: **Tue Jun 10, 2014 1:14 pm**

I would like to efficiently access membrane currents for each segment in the cell at each time step to multiply them by some matrix.

Currently I populate a numpy array on each time step before performing matrix multiplication:

Code: Select all

             jseg = 0    # iseg for the entire cell
             for sec in cell.hoc.all:  # 
                 for seg in sec:    
                     iMemSeg[jseg] = seg.i_membrane*h.area(seg.x)*1E-5 # get membrane current for the segment in (uA)
                     jseg += 1   # increment cell's segment index

This works, but is inefficient and unneccessary. I think It should be possible to generate an array (or list) which keeps objects of seg.i_membrane, which would be updated at each time step within this python array/list, so I would not have to repopulate iMemSeg at each time step. I tried creating function setMemCurArray() in my Cell class:

Code: Select all

class Cell(object):

    def __init__(self,hocTemplate):
        '''
        Constructor
        '''
    
        self.hoc = hocTemplate
        self.setMemCurArray()

        
    def setMemCurArray(self):
    
        self.iMemSeg = []
        for sec in self.hoc.all:  # for now use only a single cell for CSD calculations
            for seg in sec:    
                self.iMemSeg.append(seg.i_membrane) # get membrane current for the segment in (uA)

but this produces me zero currents at all time, so there is something wrong with my understanding of interface between NEURON and Python.

I also tried saving reference to the object in Cell class' method

Code: Select all

        def setMemCurArray(self):
       
            self.iMemSeg = []
            for sec in self.hoc.all:  # for now use only a single cell for CSD calculations
                for seg in sec:   
                    self.iMemSeg.append(seg._ref_i_membrane) # get membrane current for the segment in (uA)

but then I am not sure how to acess values in cell.iMemSeg as now each element contains a memory address of seg.i_membrane not its value.

Thanks in advance for your suggestions.

Posted: **Wed Jun 11, 2014 6:11 am**

I believe that a PtrVector will work nicely for this problem.
http://www.neuron.yale.edu/neuron/stati ... ector.html
In conjuction with Vector.as_numpy(), all the values for the i_membrane will appear in the numpy array, ready for multiplication by a numpy matrix.
http://www.neuron.yale.edu/neuron/stati ... r.as_numpy

Of course, you still need to execute setMemCurArray every time step so that that can execute PtrVec.gather(destvec)

Posted: **Wed Jun 11, 2014 8:21 pm**

Thanks Michael, however I am still not able to make it work. Below is the code for my Cell class which basically only has functions for setting and getting i_membrane: getMemCurArrSlow() is a slow implementation where I loop through the tree at each time step. Following your suggestion I used the PtrVector and created two functions: setMemCurArrEff() to set pointers to seg.i_membrane and another one getMemCurArrEff() to gather vector of seg.i_membrane at each time step. As you can see the setMemCurArrEff() is called only in constructor.

Code: Select all

class Cell(object):
    '''
    Class for cell initialized from hoc Template
    '''

    def __init__(self,hocTemplate):
        '''
        Constructor
        '''
        self.hoc = hocTemplate

        self.getNseg()
        self.iMemSeg = zeros((self.nseg,1))

        self.pVec = h.PtrVector(self.nseg)  # pointer vector
        self.setMemCurArrEff() # this should be called once only during the constructor

    
    def getNseg(self):
        
        self.nseg = 0
        for sec in self.hoc.all:
            self.nseg += sec.nseg # get the total # of segments in the cell

    def getMemCurArrSlow(self):   # slow because it loops over dendrites on each time step

        jseg = 0
        for sec in self.hoc.all:  
            for seg in sec:    
                self.iMemSeg[jseg] = seg.i_membrane 
                jseg += 1   # increment cell's segment index
        return self.iMemSeg

    def setMemCurArrEff(self): 

        jseg = 0
        for sec in self.hoc.all:  
            for seg in sec:    
                self.pVec.pset(jseg,seg._ref_i_membrane)
                jseg += 1  

    def getMemCurArrEff(self): # get membrane currents at each time step from PtrVector (no loop)

        iMemSegfromPtr = h.Vector(self.nseg)    # place here values from pointer to i_membrane
        self.pVec.gather(iMemSegfromPtr)     

        return iMemSegfromPtr.as_numpy()

I call get funcions from calcCSD() belloning to the Population class. calcCSD() is called on each time step from proc advance() after proc fadvance(). I stripped my calcCSD() from much detail just to see whether using PtrVecor works:

Code: Select all

    def calcCSD(self):  # calculate CSD:
  
        iMemSeg = cell.getMemCurArrSlow()   # get membrane currents
            
        iMemArr = cell.getMemCurArrEff()    # must be efficient but does not work 
              
        print ' iMemSeg[67]',iMemSeg[67],' iMemArr[67]',iMemArr[67]

the print statement results in iMemArr being all constant (it is a small number ~1E-316), i.e. the vector does not update, whereas iMemSeg is caluclated correctly, but very inefficiently. Do you have any idea for what is wrong with my implementation of the code?

At this point I wanted to ask you about your last comment:

Of course, you still need to execute setMemCurArray every time step so that that can execute PtrVec.gather(destvec)

This is confusing, if I still need to execute setMemCurArray (which is now actually called getMemCurArrSlow()), then what would be the point of using PtrVector? I realize that I do need to perform a PtrVec.gather(destvec). But can I avoid looping through dendritic tree at each time step?

Posted: **Thu Jun 12, 2014 7:05 am**

Forget my mention of setMemCurArray. I don't see any problem with the way you are settiing up the PtrVector once and using it to gather every time step. (although I would also setup a permanent
self.pVec and a shadow numpy array for it that you can return instead of freeing and allocating a pVec every time step

Anyway, in regard to the fact that it is not working, I would need to diagnose using the code itself. Perhaps, after cell creation, the i_membrane memory is being reallocated if you are requesting
cache efficiency.

Posted: **Thu Jun 12, 2014 3:45 pm**

I do have self.pVec initialized in the costructor and I ruturn numpy array as you recommended.

Posted: **Thu Jun 12, 2014 9:27 pm**

Thanks for sending me the code.

My guess about h.cvode.cache_efficient(1) in test.py was lucky. Comment that out and your code
works. cvode.cache_efficient reallocates all model memory and any pointers previously specified
have to be updated to the new memory locations. So it is a necessary to call Cell.setMemCurArrEff()
after the call to h.cvode.cache_efficient(1).

Posted: **Mon Jun 16, 2014 2:34 am**

Great, it works now. I would have never figured that out on my own.

Could you please tell whether h.cvode.cache_efficient(1) works when using a constant time step?

I am still puzzled with the conceptual issue of why simply populating a numpy array or (a python list) with segment currents as shown below is not sufficient to get an access to segment currents at successive time steps.

Code: Select all

def getMemCur(self):   # called in constructor

        jseg = 0
        for sec in self.hoc.all: 
            for seg in sec:   
                self.iMemSeg[jseg] = seg.i_membrane
                jseg += 1   # increment cell's segment index
        return self.iMemSeg

If I call this function only once during the consturctor, then this assignment should survive on each successive time step and elements of self.iMemSeg should report the current values of seg.i_membrane. Could you please clarify why this does not work and why one must use PtrVector instead?

Thanks in advance.

Posted: **Mon Jun 16, 2014 6:07 am**

Could you please tell whether h.cvode.cache_efficient(1) works when using a constant time step?

Yes. In fact it probably has a more beneficial effect for the fixed step than the variable step. You can experiment with both methods and see what the difference is with and without the
invocation of cache_efficient. There won't be much difference for small models, but models that do not fit into high speed memory cache on the cpu can get a lot faster with it.

The problem with only calling

Code: Select all

self.iMemSeg[jseg] = seg.i_membrane

once during setup instead of every time step is that the statement does not define the two variables to be the same location (or the first is not a mirror of the second; or the first is not
an alias for the second).
Thus when seg.i_membrane changes, self.iMemSeg[jseg] does not get changed. When the statement is executed , self.iMemSeg[jseg] just gets assigned the present value of seg.i_membrane.

Posted: **Mon Jun 16, 2014 4:09 pm**

Thanks for clarification. I realize that this is what is happening. My confusion is with why this happens given that objects in python are passed by assignment not by value. So if the element of an array (or a list) were assigned to objects once, then the updates to the values stored in those objects should be callable from that array (or a list), right?

www.neuron.yale.edu

Creating an array of membrane currents

Creating an array of membrane currents

Re: Creating an array of membrane currents

Re: Creating an array of membrane currents

Re: Creating an array of membrane currents

Re: Creating an array of membrane currents

Re: Creating an array of membrane currents

Re: Creating an array of membrane currents

Re: Creating an array of membrane currents

Re: Creating an array of membrane currents