Problem Mapping segments to 3D: Possible Range variable bug?

[Keivan]

I was trying to use map_segments_to_3d() procedure originally writen by Terrence Brannon, to assign each nseg a unique 3D address. this is the code:

Code: Select all

: 3-D mapping of model geometry 
NEURON {
	SUFFIX d3
	RANGE x, y, z
}
ASSIGNED {x(micron)  y(micron) z(micron)}

Code: Select all

load_file("morphology.nrn")

forall {insert d3}

proc endpt() {
  P=(n3d()-1)*$1

  x_d3($1)=x3d(P)
  y_d3($1)=y3d(P)
  z_d3($1)=z3d(P)
}

proc fracpt() { local posn, A, i
  A=$1
  posn=$2
  i=$3

  x_d3(posn)=x3d(i-1) + (x3d(i) - x3d(i-1))*A
  y_d3(posn)=y3d(i-1) + (y3d(i) - y3d(i-1))*A
  z_d3(posn)=z3d(i-1) + (z3d(i) - z3d(i-1))*A
}

// make 3-d mapping of cell sections
proc map_segments_to_3d() {local i,x,D,alpha
    forall {
		i=0
		endpt(0)
		for (x) if (x > 0 && x < 1) {
			while (arc3d(i) / L < x) {
				i += 1
			}
			D=arc3d(i) - arc3d(i-1)
			if (D <= 0) {
				printf("\t\t * %s had a D < 0\n", secname())
			}
			alpha = (x * L - arc3d(i-1))/D
			fracpt(alpha,x,i)
		}
		endpt(1)
	}
}
map_segments_to_3d()

but when I test it with the following code

Code: Select all

dendA5_011111111111 {
	for (x) {print "x=",x," x_d3(x)=",x_d3(x)," y_d3(x)=",y_d3(x)," z_d3(x)=",z_d3(x)}

	//this part show that the algorithm is correct and working
	i=0
	endpt(0)
	print "x_d3(0)=",x_d3(0)," y_d3(0)=",y_d3(0)," z_d3(0)=",z_d3(0)
	for (x) if (x > 0 && x < 1) {
		while (arc3d(i) / L < x) {
			i += 1
		}
		print "x_d3(x)=",x_d3(x)," y_d3(x)=",y_d3(x)," z_d3(x)=",z_d3(x)
	}
	endpt(1)
	print "x_d3(1)=",x_d3(1)," y_d3(1)=",y_d3(1)," z_d3(1)=",z_d3(1)
}

the result show that although the algorithm seems correct but the range variables x,y,z defined in d3 mechanism, behave buggy.
results of the test code:

Code: Select all

x=0  x_d3(x)=43.779999  y_d3(x)=219.50999  z_d3(x)=2.97 
x=0.5  x_d3(x)=43.779999  y_d3(x)=219.50999  z_d3(x)=2.97 
x=1  x_d3(x)=43.779999  y_d3(x)=219.50999  z_d3(x)=2.97 
//notice: value of range variables does not change above

// but the algorithm seems correct
x_d3(0)=36.290001  y_d3(0)=206.33  z_d3(0)=1.3200001 
x_d3(x)=36.290001  y_d3(x)=206.33  z_d3(x)=1.3200001 
x_d3(1)=43.779999  y_d3(1)=219.50999  z_d3(1)=2.97 

[ted]

A good question, and it's a good thing that you tested that code before relying on it.

The algorithm is indeed clever, but it is doomed to fail because of the way range variables are implemented. Recall that NEURON achieves 2nd order accuracy in space by using the central difference method for spatial discretization. The central difference method assumes that the spatial grid is sufficiently fine that the value of a range variable over the length of a compartment (segment) is well approximated by its value in the middle of the segment.

The implication of all this is that foo(x) refers to the value of the range variable foo at the middle of the segment that contains the point that corresponds to x. This fact has two important consequences.

1. If x == 0, foo(x) refers to the value of foo at the middle of the first segment, i.e. at x == 1/(2*nseg). Likewise, if x == 1, foo(x) refers to the value of foo at the middle of the last segment, i.e. at x == 1 - 1/(2*nseg).

So if nseg = 1, the "for loop"
for (x) foo(x)=f(x)
results in foo(0.5) == f(1)
and of course foo(0) and foo(1) will also == f(1) because they merely return the value of foo(0.5).

And you can see that
for (x) foo(x)=f(x)
will _always_ fail to assign the correct value to foo in the last segment, regardless of the value of nseg--because foo(1-1/(2*nseg)) will _always_ be assigned the value of f(1) instead of f(1-1/(2*nseg)). This error afflicts the complex model that was published by Poirazi et al., and it is likely to be a common problem in many models that use "for loops" to assign values to range variables.

The way to avoid this error is simple: never assign a value to a range variable at the 0 or 1 end of a section. Instead, use the
for (x,0) foo(x)=f(x)
syntax, which omits the 0 and 1 ends (see documentation of "for" at
http://www.neuron.yale.edu/neuron/stati ... r.html#for
)

2. If x lies on the boundary between two segments, foo(x) will refer to the value of foo in one of the segments, but which one? The decision will be settled by roundoff error. This is one reason why it is best to use an odd value for nseg (see the item
Why should I use an odd value for nseg?
in the FAQ list http://www.neuron.yale.edu/neuron/faq/general-questions ).


Anyway, back to the code you posted: the fix, of course, is to omit the assignment at x == 1.