## Step Functions in NEURON

Anything that doesn't fit elsewhere.
delcin
Posts: 5
Joined: Tue Jul 26, 2016 5:08 pm

### Step Functions in NEURON

Hello,

I am working with the code in Pashut et al 2011 to make a realistic 2-D model of TMS on a layer V pyramidal neuron. Right now, I am having trouble understanding how NEURON is calculating step functions in this loop:

Code: Select all

``````for ii = 0, nseg {
xr = range.x[ii]
xrf= range.x[ii+1]

//x_xtra(xr:xrf)=xint.x[ii]:xint.x[ii+1] //Removed x and y_xtra, they aren't being used

//x_xtra(xr:xrf)=int(xint.x[ii]):xint.x[ii+1]     //[micrometer]

//y_xtra(xr:xrf)=yint.x[ii]:yint.x[ii+1]

XX_xtra(xr:xrf)=int(xint.x[ii]+XShift):int(xint.x[ii+1]+XShift) //[micrometer]
YY_xtra(xr:xrf)=int(yint.x[ii]+YShift):int(yint.x[ii+1]+YShift)

DX_xtra(xr:xrf)=xint.x[ii+1]-xint.x[ii]:xint.x[ii+1]-xint.x[ii] //[micrometer]
DY_xtra(xr:xrf)=yint.x[ii+1]-yint.x[ii]:yint.x[ii+1]-yint.x[ii]

//the value of the spatial part of the electric field in the direction of the cable
Exs_xtra(xr:xrf)=Ex.x[XX_xtra(xr)][YY_xtra(xr)]:Ex.x[XX_xtra(xrf)][YY_xtra(xrf)]
Eys_xtra(xr:xrf)=Ey.x[XX_xtra(xr)][YY_xtra(xr)]:Ey.x[XX_xtra(xrf)][YY_xtra(xrf)]

//The spatial derivative of the above value.  This is required
//in order to transform the axial current I=E*X/ri to the membrane
//current we will inject via the MOD file since, acording to cable
//theory i_m=-dI/dx.  Units are [V ms/m A microm]. This is
//the derivative per unit length.

BB=((Ex.x[int(xint.x[ii+1]+XShift)][int(yint.x[ii+1]+YShift)] - Ex.x[int(xint.x[ii]+XShift)][int(yint.x[ii]+YShift)])*DX_xtra(xr) + (Ey.x[int(xint.x[ii+1]+XShift)][int(yint.x[ii+1]+YShift)] - Ey.x[int(xint.x[ii]+XShift)][int(yint.x[ii]+YShift)])*DY_xtra(xr))/(DX_xtra(xr)^2+DY_xtra(xr)^2)

DEDA_xtra(xr:xrf)=BB:BB
``````
To see how this loop works, I chose a dendrite with 4 specified points and therefore 3 segments. The variable range is defined as:

Code: Select all

``````range = new Vector(nseg+2)
range.indgen(1/nseg)
range.sub(1/(2*nseg))
range.x[0]=0
range.x[nseg+1]=1``````
My understanding is that the authors meant for range to be the midpoints of sections, which makes sense as that is where NEURON defines its values of rx, diam etc. In a case of a dendrite with three segments, range is [0,.1667,.5,.8333,1].

However, when I printed out the values of the parameters defined in the loop, even though the first loop is meant to define values from x = 0 to x = .166667, in reality it defines values from x = 0 to x = .333, the end of first section. So instead of midpoints, the step function for all these parameters change values at points where segments change. Am I missing something? How can I fix this loop to define values for midpoints?

Thank you,

Dehan
ted
Posts: 5787
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

### Re: Step Functions in NEURON

That particular bit of programming wins this week's "obfuscated code contest."
Using "tapered dendrite syntax" to specify section geometry is not a particularly good way to do anything--maybe that's why almost nobody uses it. The problem is that it is such a great way to write code that confuses the user. Why? Because nseg must be specified before the taper is specified, and if nseg is changed afterward, the taper is destroyed. Consider this example:

Code: Select all

``````proc report() {
for (x) print x, x*L, diam(x)
}

create dend
access dend
L = 100
nseg = 5
diam(0:1) = 2:1
report()``````
NEURON prints
0 0 1.9
0.1 10 1.9
0.3 30 1.7
0.5 50 1.5
0.7 70 1.3
0.9 90 1.1
1 100 1.1
as expected.
Now do

Code: Select all

``````nseg = 3
report()``````
and NEURON prints

0 0 1.9
0.16666667 16.666667 1.9
0.5 50 1.5
0.83333333 83.333333 1.1
1 100 1.1

so taper is too shallow at the proximal and distal ends, and too steep in the middle.
Try to fix this by going back to nseg of 5 and how does diameter change with distance from the 0 end?
0 0 1.9
0.1 10 1.9
0.3 30 1.9
0.5 50 1.5
0.7 70 1.1
0.9 90 1.1
1 100 1.1
which isn't an improvement.

Guess what happens if you ever set dend's nseg to 1? Best to try it yourself and see what happens.

(assuming you have now done that experiment)--You like that? I don't.

pt3d syntax is far easier to use. The syntax is clear and the results make sense because NEURON remembers the 3d point data and reuses those if you ever change nseg. Consider this example:

Code: Select all

``````create axon
access axon
pt3dclear()

nseg = 5
report()``````
which generates this output
0 0 1.9
0.1 10 1.9
0.3 30 1.7
0.5 50 1.5
0.7 70 1.3
0.9 90 1.1
1 100 1.1
as it should.

Now do

Code: Select all

``````nseg = 3
report()``````
and you get
0 0 1.8333333
0.16666667 16.666667 1.8333333
0.5 50 1.5
0.83333333 83.333333 1.1666667
1 100 1.1666667
which is correct.

And if you do this

Code: Select all

``````create soma, dend
connect dend(0), soma(1)
soma {
pt3dclear()
}
dend {
pt3dclear()
}``````
what do you think happens to the coordinates of dend's 0 and 1 ends? Well, if you create a Shape plot or merely call
define_shape()
you'll find that
dend print x3d(0)
prints 10, and
dend print x3d(1)
prints 110--exactly what they should be, given the fact that soma is the root section and the proximal end of dend is attached to soma's 1 end (which is located at (10,0,0)). In other words, the xyz coords used to define the geometry of each non-root section are translated to where they should be so that the 0 end of the non-root section is at the same location as the point of attachment on its parent section.

The spatial coordinates of segment centers are quickly calculated by treating successive pt3d data points along a section as breakpoints in piecewise linear functions of range along the section and then evaluating those functions at the ranges that correspond to segment centers. For an example of this, examine interpxyz.hoc in extracellular_stim_and_rec.zip, which you can download from this link viewtopic.php?f=28&t=168 in the Hot tips area of the Forum.

What are the step functions mentioned in your post?
delcin
Posts: 5
Joined: Tue Jul 26, 2016 5:08 pm

### Re: Step Functions in NEURON

Thank you for the fast response Ted!

The code I am working with does use p3d syntax to specify section geometry. The loop I shared in my original post is meant to specify certain variables that correspond to midpoints of each segment in a section.

Looking at interpxyz.hoc it is clear that they based their code on it for interpolation. Range is defined exactly the way it is defined in interpxyz.hoc, to correspond to midpoints of the segments of the currently accessed section.

Then they use range to define step functions that change values in each of the midpoint of segments.If the currently accessed section was a dendrite with 3 specified points using p3d(), range would be (0,.16,.5,.81,1), then the authors use range to define parameters. For example,

Code: Select all

``````for ii = 0, nseg {
xr = range.x[ii]
xrf= range.x[ii+1]
XX_xtra(xr:xrf)=int(xint.x[ii]+XShift):int(xint.x[ii+1]+XShift)``````

xint, is also defined in the same way as it is in interpxyz.hoc, it is the interpolated x-coordinates of the specified points of the segment. In this case XX_xtra is meant to correspond to the interpolated x coordinate + Shift, which defines the shift in the whole cell's position. So overall XX_xtra is the x-coordinate of the segment in the extracellular environment.

The first time the loop above runs, xr is 0, whereas xrf is .16. Although I am not fully familiar with the syntax, I am assuming they are trying to set values for XX_xtra that stay the same from 0 to .16, then from ,16 to .5, etc. However, when I print XX_xtra, I see that it is changing values at (0,.33,.66,.99,1) meaning at the ends of segments. I am trying to change this loop so that XX_xtra changes values at the midpoints of segments.

Dehan
ted
Posts: 5787
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

### Re: Step Functions in NEURON

delcin wrote:The code I am working with does use p3d syntax to specify section geometry.
Good. After model setup execute define_shape() just to make sure that the xyz coords of the 0 end of each child section agrees with the xyz coordinate of the location on its parent to which it is attached. If you need to move the entire cell to a particular location in space, just use pt3dchange to change the xyz coordinates of each 3d point that belongs to the soma, then execute define_shape() again, and all xyz coordinates will be changed to reflect the cell's new position. This procedure will change the position of the soma's 3d points:

Code: Select all

``````x = 0
y = 0
z = 0
proc position() { local i
soma for i = 0, n3d()-1 {
pt3dchange(i, \$1-x+x3d(i), \$2-y+y3d(i), \$3-z+z3d(i), diam3d(i))
}
x = \$1  y = \$2  z = \$3
}``````
(note that x, y, and z must be declared outside of proc position()). If your model cell is an instance of a class defined by a template, that template will already contain a proc position() if (1) it was exported from a CellBuilder or Network Builder or (2) the template was written by someone who knew to include a proc position().
The loop I shared in my original post is meant to specify certain variables that correspond to midpoints of each segment in a section.
If you merely need to know the range values that correspond to segment centers (i.e. relative positions along the length of a section at which the segment centers are located), just iterate over the section with
for (x) statement
or
for (x,0) statement
The first version includes the nodes at the 0 and 1 ends, and the second omits those points.

If you need to know the xyz coordinates of segment centers, that's what proc grindaway in interpxyz.hoc does. It stores those values in the xtra mechanism's x, y, and z range variables.
they use range to define step functions that change values in each of the midpoint of segments.If the currently accessed section was a dendrite with 3 specified points using p3d(), range would be (0,.16,.5,.81,1), then the authors use range to define parameters. For example,

Code: Select all

``````for ii = 0, nseg {
xr = range.x[ii]
xrf= range.x[ii+1]
XX_xtra(xr:xrf)=int(xint.x[ii]+XShift):int(xint.x[ii+1]+XShift)``````
There's a couple of big problems in this for loop.

First
for ii = 0, nseg
tries to iterate nseg+1 points. That's one more point than there are segments. The final pass through this loop must produce a garbage result, even if passes 0..nseg-1 worked just fine (which seems doubtful).

Second, the right hand side of the statement that begins with XX_xtra forces the values on both sides of the
a:b
syntax to be integers. That guarantees roundoff errors unless one is very careful about jointly constraining the values of nseg and section lengths.

If you absolutely must reuse that code, bugs and all, go ahead. But expect your results to be challenged by anyone who ever reads it.
I am trying to change this loop so that XX_xtra changes values at the midpoints of segments.
What values do you want it to have? Are you simply trying to determine the locations of segment centers, e.g. for the purpose of simulating extracellular recording or stimulation?
delcin
Posts: 5
Joined: Tue Jul 26, 2016 5:08 pm

### Re: Step Functions in NEURON

I wrote out XX_xtra as an example to show how the loop was meant to work, but it makes much more sense to get the mid-points of the segments using for(x), so thank you for that!

The main function of this part of the code is to assign values for the spatial derivative of the electric field corresponding to the place in the field the segment is in. This is calculated by:

Code: Select all

``````                     for ii = 0, nseg {
xr = range.x[ii]
xrf= range.x[ii+1]
BB=((Ex.x[int(xint.x[ii+1]+XShift)][int(yint.x[ii+1]+YShift)] - Ex.x[int(xint.x[ii]+XShift)][int(yint.x[ii]+YShift)])*DX_xtra(xr) + (Ey.x[int(xint.x[ii+1]+XShift)][int(yint.x[ii+1]+YShift)] - Ey.x[int(xint.x[ii]+XShift)][int(yint.x[ii]+YShift)])*DY_xtra(xr))/(DX_xtra(xr)^2+DY_xtra(xr)^2)

DEDA_xtra(xr:xrf)=BB:BB

``````

Ex and Ey are 4000x4000 matrices corresponding to the magnitude of the electric field at the specified points, DX_xtra and DY_xtra are the distances from end-midpoint,midpoint-end for the beginning and end of the segment and, from midpoint-midpoint for the middle of the segment. I want to change this loop and get it to work correctly, meaning I want get DEDA correspond to the spatial derivative of the electric field between midpoints.

Thank you so much for your time,

Dehan
ted
Posts: 5787
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

### Re: Step Functions in NEURON

I downloaded that particular model entry and examined the 1 dimensional case. The code works, after a fashion, for the specific parameters used by the authors, but contains many crude hacks and deliberate inaccuracies that might escape notice, but it is brittle and bound to cause problems for anyone who attempts to reuse it and trust it outside of the domain in which those inaccuracies are small.

First some coments, and then the big problems.

RB.hoc
Default section is not the root section (default section is dend[0], root is dend[1]).
Each section's nseg is 100 (should be an odd number).
Why is dend[0] 8e4 um long but dend[1] is 1 um shorter?

xtra.mod
This is absolutely not the same as the xtra.mod that I developed and distributed.
Problem: it calculates a variable called i_m in the BREAKPOINT block.
i_m is what they're using to stimulate the cell.
i_m should be calculated before each cy' = f(y,t) setup,
which means it should be calculated in a BEFORE BREAKPOINT block.

The most serious problem: the use of int() in proc grindaway() inserts roundoff errors that will be noticeable as soon as the code is used with reasonable values for nseg and neurites of reasonable length. Example: with nseg = 3, dend[0] 10 um long, dend[1] 100 um long, it calculates dend[0]'s segment centers as being at x locations of -1, -5, and -8 um, and dend[1]'s segment centers as being at 16, 50, and 83 um. The original proc grindaway() that I developed calculates segment centers correctly and introduces no such errors.
delcin
Posts: 5
Joined: Tue Jul 26, 2016 5:08 pm

### Re: Step Functions in NEURON

ted wrote:RB.hoc
Default section is not the root section (default section is dend[0], root is dend[1]).
Each section's nseg is 100 (should be an odd number).
Why is dend[0] 8e4 um long but dend[1] is 1 um shorter?
I have no idea why these are the way they are but I have changed them since.
ted wrote:xtra.mod
This is absolutely not the same as the xtra.mod that I developed and distributed.
Problem: it calculates a variable called i_m in the BREAKPOINT block.
i_m is what they're using to stimulate the cell.
i_m should be calculated before each cy' = f(y,t) setup,
which means it should be calculated in a BEFORE BREAKPOINT block.
Is there any documentation I can study to understand how .mod files function? I don't think I understand why this is an issue. Thank you for pointing it out anyways though.

I believe the authors use int() because they need integer numbers to be able to call values from the electric field that is defined as a matrix. Would the rounding error matter much in a case which the size of the modeled electric field is very large compared to the neuron cell?

Also what exactly is the difference between

Code: Select all

``YY_xtra(xr:xrf)=yint.x[ii]+YShift:yint.x[ii+1]+YShift``
and

Code: Select all

``YY_xtra(xr)=yint.x[ii]+YShift``
? Don't these two achieve essentially the same thing?
ted
Posts: 5787
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

### Re: Step Functions in NEURON

ted wrote:xtra.mod
This is absolutely not the same as the xtra.mod that I developed and distributed.
Problem: it calculates a variable called i_m in the BREAKPOINT block.
i_m is what they're using to stimulate the cell.
i_m should be calculated before each cy' = f(y,t) setup,
which means it should be calculated in a BEFORE BREAKPOINT block.
Is there any documentation I can study to understand how .mod files function? I don't think I understand why this is an issue.
Well, you could read the original xtra.mod, which contains a comment that says why BEFORE BREAKPOINT is where the calculation of the driving function should be performed. Ah, but the authors of the model that you are now trying to reuse didn't bother to include those comments, did they? Did they at least say where they got the xtra.mod that they then modified?

In brief, if forcing functions are involved, their values must be calculated before they are used to advance the solution. That's the purpose of the BEFORE BREAKPOINT block--it's where to put statements that must be executed before the solution is advanced from one time to the next. Most models don't involve calculation of forcing functions in a mod file, which is why you won't see many mod files that have BEFORE BREAKPOINT blocks.

Anyway, your forcing function is i_m, and if the formula for i_m is in the BREAKPOINT block, i_m is calculated after the model's system equations are solved. Not good.
I believe the authors use int() because they need integer numbers to be able to call values from the electric field that is defined as a matrix.
There are principled methods for interpolation of values at points that lie off the grid.
Would the rounding error matter much in a case which the size of the modeled electric field is very large compared to the neuron cell?
The relevant comparison is the resolution of the grid compared to the length of a segment. The authors deliberately chose parameters that corresponded to very long segments, so a 1 um grid interval introduces only a small error. A model of a real cell (especially a cerebellar Purkinje cell) is likely to have many segments that are short enough that a 1 um grid is rather coarse. I'm all in favor of back of the envelope estimates and expedient hacks, but it's important to be explicit about such approximations and their limitations--yet the authors' code contains almost no comments.
Also what exactly is the difference between

Code: Select all

``YY_xtra(xr:xrf)=yint.x[ii]+YShift:yint.x[ii+1]+YShift``
and

Code: Select all

``YY_xtra(xr)=yint.x[ii]+YShift``
? Don't these two achieve essentially the same thing?
Do they? The're compact, but they're also as opaque as a bottle of ink. Can you explain, in a combination of plain english and standard mathematical notation, with maybe a numerical example or two, what these statements are supposed to do, and how they work? To simplify the task, you might want to assume a small value for nseg, like maybe 3 or 5.
delcin
Posts: 5
Joined: Tue Jul 26, 2016 5:08 pm

### Re: Step Functions in NEURON

Thank you for that explanation. It makes more sense now.

The statements I wrote would be in a loop that goes like this:

Code: Select all

``````forall {

for ii = 0, nseg {
xr = range.x[ii]
xrf= range.x[ii+1]
YY_xtra(xr:xrf)=yint.x[ii]+YShift:yint.x[ii+1]+YShift
}
}
``````
If we assume that the currently accessed section is a dendrite with 3 segments, range = [0,.1666,.4999,.8333,1]. yint is the interpolated y coordinates of the segment. YShift determines the position of the neuron in the electric field. YY_xtra is a range variable that corresponds to the interpolated y coordinates of a segment in the electric field.

When this loop runs, first time through it sets YY_xtra(0:.166), second time (.166:.499) etc. This makes me think that the range variable is set to change at the midpoints of the segments. However, when I print out YY_xtra, I see that setting YY_xtra(0:.166), sets YY_xtra from 0 to .33. I am essentially trying to find a way to modify range variables so that I can get them to change values at the midpoints and not at the endpoints as this is where rx() and area() are calculated.

YY_xtra here is just an example range variable that I don't really intend to use. As you pointed out there are better ways to do that. However I would like to be able to define range variables that change at midpoints to use when defining the electric field.
ted
Posts: 5787
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

### Re: Step Functions in NEURON

OK, let's forget about hoc's "taper over a range" syntax for the time being.
I would like to be able to define range variables that change at midpoints to use when defining the electric field.
Time to review some principles. The spread of electrical signals in nerve cells is described by a family of cable equations. There is one cable equation for each unbranched neurite in the cell (neurite = soma, axon, dendrite, spine neck, spine head, presynaptic terminal). Each neurite's cable equation is a partial differential equation that relates axial and radial (transmembrane) currents to membrane potential. At the anatomical and temporal scale of neuronal structures, these variables are all continuous functions of time and space.

To convert the cable equations into a form that is suitable for approximate numerical solution by a digital computer, it is first necessary to discretize both time and space. This means that the numerical solution will approximate the time course of electrical (and chemical) signals over a finite number of points in space, at a finite number of points in time. These points are called grid points.

Spatial discretization in NEURON involves approximating each neurite ("section") by a set of compartment of equal length. The number of compartments (segments) in a section, which is specified by the section parameter nseg, should be large enough that the transmembrane current densities within any compartment are nearly uniform over the surface of that compartment. The second spatial derivative of membrane potential is approximated by the central difference method. This has two important consequences.

1. The numerical error resulting from spatial discretization is proportional to (1/nseg)^2. This means that membrane potential will be second order accurate in space at the points that are located at segment centers.

2. Second order approximate values for the membrane potential at spatial locations between adjacent grid points can be calculated by simple linear interpolation. In other words, the membrane potential at grid points can be treated as a nodes (locations of possible slope change) in a piecewise linear approximation to membrane potential along the entire model.

"But I thought that it was assumed that the transmembrane potential over the surface of a segment was uniform."

No, the assumption is that the spatial grid is sufficiently fine that the transmembrane current densities within a compartment are nearly uniform over the surface of the compartment.

If you need to simulate the effects of an electrical or magnetic field on a neuron
AND
the spatial discretization of the neuron model is already sufficiently fine to capture the structure of that field
THEN
you can represent the effects of the field by perturbing currents or potentials that are applied to the model at locations that correspond to segment centers.

However, if the field has finer structure that you think may have a significant effect on your model, you may need to increase the resolution of your model neuron's spatial grid accordingly.