Page 1 of 1

Accessing segment diameter

Posted: Wed Mar 12, 2014 8:51 pm
by budd
Hi,

What is the correct way of accessing segment diameter from python? In a python version of 'cone.hoc' (see 'cone.py' below) seg.diam did not return the same values compared to a hoc call (see output below).

cone.py

Code: Select all

import neuron
import math

h = neuron.h
a = h.Section()       
a.Ra = 100
a.nseg = 10

a.push()  
h.pt3dclear()
for i in range(0,31):   
    x = math.pi*i/30.0
    h.pt3dadd(200*math.sin(x), 200*math.cos(x), 0, 100*math.sin(4*x))

for seg in a.allseg():
    print seg.x, seg.diam

h.psection()    
    
h("for (x) print x, diam(x)")  
output

Code: Select all

0.0 500.0
0.05 500.0
0.15 500.0
0.25 500.0
0.35 500.0
0.45 500.0
0.55 500.0
0.65 500.0
0.75 500.0
0.85 500.0
0.95 500.0
1.0 500.0
PySec_b72bc938 { nseg=10  L=628.031  Ra=100
	/*location 0 attached to cell 0*/
	/* First segment only */
	insert capacitance { cm=1}
	insert morphology { diam=500}
}
0 54.180325 
0.05 54.180325 
0.15 87.665604 
0.25 33.45362 
0.35 87.665604 
0.45 54.180325 
0.55 54.180325 
0.65 87.665604 
0.75 33.45362 
0.85 87.665604 
0.95 54.180325 
1 54.180325 

Re: Accessing segment diameter

Posted: Thu Mar 13, 2014 11:42 am
by ted
Interesting. I'd use

Code: Select all

for seg in a:
  print seg.x, seg.diam
just because it's easier to type--but it produces the same result you got.

The problem is not in how you are trying to access segment diameter. it's that the pt3d data have not yet been used to compute the segment diameters. However, if you simply execute

Code: Select all

for seg in a:
  print seg.x, seg.diam
AFTER the
h("for (x) print x, diam(x)")
statement, you'll see that Python reports diameters that are based on the pt3d data.

What's happening: the pt3dadd statements really are setting up the pt3d data structures, but translation of those structures into segment geometry is deferred until it is "absolutely necessary," that is, until a statement is executed that requires segment properties to be consistent with the pt3d data. This makes sense because models with detailed 3D anatomy may have geometry specifications that require hundreds or thousands of pt3dadd statements, and it would waste time if every pt3dadd statement forced updating of everything that depends on cell geometry.

Examples of statements that force updating of segment geometry include things that are absolutely essential for simulations to generate correct results, e.g. h.finitialize() which is used to initialize a model, and h.area(x) which calculates the area of the segment that contains location x where x lies in the interval 0..1 (sometimes used in model setup code when specifying biophysical properties). And in hoc, even diam(x) has the power to trigger segment geometry updates, as you discovered. But accessing segment diameter directly from Python doesn't--yet.

If this is fixable, the fix will probably soon appear in the latest development code. Whether or not it is fixable, a practical workaround is just to execute h.area(x) or h.finitialize() after the last pt3dadd statement has been executed.


Two final comments:
1. You really don't want pt3d data to include extremely small diam3d values, because those will result in extremely high series resistance between adjacent compartments. Examine your diam3d values and you'll find three.
2. It's a good idea for nseg to be odd. Read about this, and about effective strategies for choosing nseg, in the FAQ list
https://www.neuron.yale.edu/neuron/faq/ ... -questions

Re: Accessing segment diameter

Posted: Thu Mar 13, 2014 8:17 pm
by budd
Ted, thanks for your detailed reply explaining this python curiosity. From it I have gained a better understanding of the translation process of pt3d data into segment geometry.

Using the python statement h.finitialize() after the last pt3dadd statement instead of the hoc call h("for (x) print x, diam(x)") does indeed produce the same output.

Code: Select all

import neuron
import math

h = neuron.h
a = h.Section()       
a.Ra = 100
a.nseg = 10

a.push()  
h.pt3dclear()
for i in range(0,31):   # made 31 not 30 because of range() limits.
    x = math.pi*i/30.0
    h.pt3dadd(200*math.sin(x), 200*math.cos(x), 0, 100*math.sin(4*x))

h.finitialize()

for seg in a.allseg():
    print seg.x, seg.diam
output

Code: Select all

0.0 54.1803251938
0.05 54.1803251938
0.15 87.6656042459
0.25 33.4536201131
0.35 87.6656042459
0.45 54.1803251938
0.55 54.1803251938
0.65 87.6656042459
0.75 33.4536201131
0.85 87.6656042459
0.95 54.1803251938
1.0 54.1803251938
Regarding your two final comments, I used a python implementation of 'cone.hoc' to simply illustrate the basic problem not for simulation. I guessed from the use of the sine function segment diameter might reach zero, which is confirmed by printing out pt3d diameter data using h.diam3d() python function. I would point out that the original hoc code (http://www.neuron.yale.edu/neuron/stati ... n/cone.hoc) also has an even number of segments (nseg =10) :-) Nevertheless, I appreciate the intention.

Re: Accessing segment diameter

Posted: Thu Mar 13, 2014 9:37 pm
by ted
A lot of what ends up in my answers is there to convey max info to the widest audience while trying to prevent future problems.

Re: Accessing segment diameter

Posted: Fri Mar 14, 2014 10:52 am
by hines
The Python interface was not checking the Section recalc_area_ flag when evaluating diam. I pushed a fix to
http://www.neuron.yale.edu/hg/neuron/nr ... bc03849e21