diam(0:1) = A:B in Python

When Python is the interpreter, what is a good
design for the interface to the basic NEURON
concepts.

Moderator: hines

Post Reply
vellamike

diam(0:1) = A:B in Python

Post by vellamike » Wed Dec 01, 2010 12:18 pm

I hope there isn't a better way to do this already but I wrote a python method to reproduce the output of diam(0:1)=zero_bound:one_bound in hoc. Could someone check this works correctly as I'm still not convinced I understand precisely how diam(0:1) works.

Code: Select all

####################################################
#    method *should* produce the same output as
#    sec.diam(0:1)=zero_bound:one_bound         
####################################################

def taper_diam(sec,zero_bound,one_bound):
       
    i=0 #segment number (1-based counting)
    x=0
    for seg in sec:
        i+=1
        x+=1/(sec.nseg+1)*i
        seg.diam=(one_bound-zero_bound)*x+zero_bound

ted
Site Admin
Posts: 5645
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: diam(0:1) = A:B in Python

Post by ted » Wed Dec 01, 2010 1:05 pm

How to check it out for yourself:
1. Read this and see if it doesn't help clear up any questions you might have
http://www.neuron.yale.edu/neuron/stati ... /nc4p2.htm
2. Fire up two NEURON sessions, one using Python and the other with hoc as the interpreter, and compare the results of your Python code against the results of the hoc code.

vellamike

Re: diam(0:1) = A:B in Python

Post by vellamike » Thu Dec 02, 2010 10:53 am

Thanks, I did as you suggested and actually had to modify my code quite a bit, here's the method that seems to work correctly:

Code: Select all

def taper_diam(sec,zero_bound,one_bound):
    dx=1.0/(sec.nseg)
    x=dx/2
    for seg in sec:
        seg.diam=(one_bound-zero_bound)*x+zero_bound
        x+=dx
Is there any way in the Python implementation of Neuron to do something along the lines of rangevar(xmin:xmax) without writing your own methods like this?

ted
Site Admin
Posts: 5645
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: diam(0:1) = A:B in Python

Post by ted » Thu Dec 02, 2010 11:34 am

Good question.
rangevar(xmin:xmax) = minval:maxval
seems quite nonpythonic; maybe somebody who is fluent in Python can propose a reasonable alternative that could be incorporated into the Python interface to NEURON.

Barring that, if you're dead set on specifying taper the way hoc does, you can always cheat.

Code: Select all

h("create dend")
h("dend nseg = 5")
h("dend diam(0:1) = 7:3")
Then

Code: Select all

h("dend for (x,0) print x, diam(x)")
or more pythonically

Code: Select all

for seg in h.dend:
  print seg.x, seg.diam
returns

Code: Select all

0.1 6.6 
0.3 5.8 
0.5 5 
0.7 4.2 
0.9 3.4 
In any case, I suspect that if there isn't a compact Python equivalent for hoc's "range variable taper" syntax, it won't matter much. Few existing models implemented in hoc use the "taper" syntax. For whatever reason, most modelers that use hoc have either assumed cylindrical neurites with constant diameter, or resorted to detailed morphometric data which makes the "taper" syntax irrelevant. On those occasions where a parameter must be varied according to some mathematical rule, it is easy enough to write a procedure that applies the rule to those neurites that need it.

ted
Site Admin
Posts: 5645
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: diam(0:1) = A:B in Python

Post by ted » Thu Dec 02, 2010 12:22 pm

Perhaps the best way to deal with algorithmic variation of range variables is to conceive of the problem not as one that requires magic syntax, but as a special instance of the general problem of iterating over members of a set and taking the appropriate action for each member. Stated in pseudocode

Code: Select all

for each internal node in a section
  rangevar at that node = function(node position)
Restated in Python

Code: Select all

for seg in sec:
  seg.param = f(seg.x)
which is understandable at a glance, unlike hoc's taper syntax.

mctavish
Posts: 74
Joined: Tue Mar 14, 2006 6:10 pm
Location: New Haven, CT

Re: diam(0:1) = A:B in Python

Post by mctavish » Thu Dec 02, 2010 2:19 pm

This type of linear mapping can be easily handled with numpy's linspace.

Code: Select all

import numpy
from neuron import h as nrn
from itertools import izip

def taper_diam(sec,zero_bound,one_bound):
    for (seg, d) in izip(sec, numpy.linspace(zero_bound, one_bound, sec.nseg)):
        seg.diam=d
        
# Test
dend = nrn.Section(name='dend')
dend.L = 100
dend.nseg = 5
taper_diam(dend, 10, 2)
for seg in dend:
    print seg.diam

ted
Site Admin
Posts: 5645
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: diam(0:1) = A:B in Python

Post by ted » Thu Dec 02, 2010 4:01 pm

Excellently pythonic, but it produces a result that differs from hoc's taper syntax.
dend diam(0:1) = 2:3
results in these diameters
2.1
2.3
2.5
2.7
2.9
With Tom's Python code,
taper_diam(dend, 2, 3)
produces these diameters
2.0
2.25
2.5
2.75
3.0
This result may seem cosmetically preferable to the result of the hoc syntax, but it's a poorer approximation to a neurite that tapers from 2 at the 0 end to 3 at the 1 end. The diameter is too small for range<0.5 and too large for range>0.5. NEURON's approach to spatial discretization is based on the central difference approximation of the second derivative in space, which requires that the diameter of any compartment be the mean of the actual anatomical diameter over the length of that compartment.

mctavish
Posts: 74
Joined: Tue Mar 14, 2006 6:10 pm
Location: New Haven, CT

Re: diam(0:1) = A:B in Python

Post by mctavish » Thu Dec 02, 2010 4:23 pm

Good eye, Ted! I was coding what I thought Mike wanted. Here is a hoc equivalent, which is quite comparable to what he wrote before:

Code: Select all

import numpy
from neuron import h as nrn
from itertools import izip

def taper_diam(sec,zero_bound,one_bound):
    dx = 1./sec.nseg
    for (seg, x) in izip(sec, numpy.arange(dx/2, 1, dx)):
        seg.diam=(one_bound-zero_bound)*x+zero_bound
       
# Test
dend = nrn.Section(name='dend')
dend.L = 100
dend.nseg = 5
taper_diam(dend, 2, 3)
for seg in dend:
    print seg.diam

ted
Site Admin
Posts: 5645
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: diam(0:1) = A:B in Python

Post by ted » Thu Dec 02, 2010 4:44 pm

That should do it!

Post Reply