Reducing complex models to equivalent cables

[mattmc88]

Is it possible to use NEURON to convert a model (which consists of passive biophysics distributed over a realistic morphology imported from neurolucida via cvapp) into a single cable of varying diameter, using a morphoelectrotonic transform? See, for example (Burke RE, Journal of Computational Neuroscience 9, 31-47, 2000) for a variety of transforms and their relative (dis)advantages. I would like to do this in order to construct small compartmental models in hardware which mimic the modeled cell's response.

[ted]

It's done all the time, but there's no gold standard for how to do it;
as you have noted, there are lots of approaches to desiging reduced
compartmental models, e.g.
Bush, P.C. and Sejnowski, T.J. Reduced compartmental models of
neocortical pyramidal cells. Journal of Neuroscience Methods 46:159-166,
1993.
to cite just one. NEURON has no built-in tool that does the job for you,
but it has everything you need to do the task yourself.

[mattmc88]

I have been trying to modify the attshape.hoc file in order to output log(A) for every point, without much luck. What I really would like is the diameter as a funtion of electrotonic distance for each branch. Has anyone out there done something like this before, or have any suggestions for how to do it?

[ted]

I have been trying to modify the attshape.hoc file

I hope you were working on a duplicate copy of attshape.hoc, in some directory
other than where attshape.hoc is located. It's not a good idea to do anything to
the installed files or directories themselves; too easy to break stuff or--much
worse--introduce a bug that lets the program run but produce incorrect results.

What I really would like is the diameter as a funtion of electrotonic distance
for each branch.

It's pretty straightforward to use the Impedance class to discover all sorts of things
about a model. The following demo shows how.

Code: Select all

load_file("nrngui.hoc")

// toy model cell
create soma, dend, axon
access soma
connect axon(0), soma(0)
connect dend(0), soma(1)

soma {
  diam = 30  L = 30
  insert pas
  g_pas /= 20
  e_pas = -65  // same as HH resting potential
}

dend {
  diam(0:1) = 2:0.5  L = 600
  insert pas
  g_pas /= 20
  e_pas = -65  // same as HH resting potential
}

axon {
  diam = 1  L = 1000
  insert hh
}

// specify reasonable values for nseg--see
//   http://www.neuron.yale.edu/neuron/static/docs/d_lambda/d_lambda.html
freq = 100      // Hz, frequency at which AC length constant will be computed
d_lambda = 0.1
func lambda_f() { // currently accessed section, $1 == frequency
  return 1e5*sqrt(diam/(4*PI*$1*Ra*cm))
}

proc geom_nseg() {
  soma area(0.5) // make sure diam reflects 3d points
  forall { nseg = int((L/(d_lambda*lambda_f(freq))+0.9)/2)*2 + 1  }
}

geom_nseg()

// prepare to use Impedance class

// always a good idea to finitialize before computing impedance
v_init=-65
finitialize(v_init)

// demonstrate use of impedance class
objref zz
zz = new Impedance()
FREQ = 0 // Hz

WHERE = 0.5 // location in the soma that is the reference point
soma distance(0, WHERE)  // sets origin for distance calculations

proc calcZ() {
  soma zz.loc(WHERE)  // sets origin for impedance calculations
  /*
    zz.compute() assumes conductances depend only on v, i.e.
    ignores the impedance contributions of gating state differential equations
  */
  zz.compute(FREQ, 1) // takes the impedance contributions of 
                      // gating state differential equations into account
                      // but requires mechanisms to be compatible with CVODE

  print "x distance(x) input(x) input_phase(x) transfer(x) transfer_phase(x) ratio(x)"
  forall {
    print secname()
    for (x) print x, distance(x), zz.input(x), zz.input_phase(x), zz.transfer(x), zz.transfer_phase(x), zz.ratio(x)
  }
}

calcZ()

[yml23]

I'm working with Matt on this equivalent cable problem and have been playing with the code you posted.

1) I assume that imp.ratio(freq) is the attenuation value we need to calculate L = log(A), correct?

2) I don't understand where the "x" in the "ratio(x)" of your code comes from. According to the procedure declaration in impedanc.hoc, the input argument is freq. Where does x come from?

3) I see how this demo, and the attshape.hoc code, calculate the attenuation over the entire dendrite. How do I calculate it for individual segments? It's pretty useless over a large dendrite. Can you give a simple explanation or example of how to access data in the individual segments of an existing imported morphology (instead of creating new segments)? Do I just defint nseg small enough that the new segments are always smaller than the imported segments?

Thanks,
Yaniv Lazimy

[ted]

1) I assume that imp.ratio(freq) is the attenuation value we need to calculate L = log(A), correct?

It's the voltage transfer ratio, defined as "voltage downstream/voltage upstream",
which is necessarily <= 1. Attenuation, however, is the inverse of this. You will find
the Impedance class described here
http://www.neuron.yale.edu/neuron/stati ... edanc.html

2) I don't understand where the "x" in the "ratio(x)" of your code comes from. According to the procedure declaration in impedanc.hoc, the input argument is freq. Where does x come from?

In this context, x is the location of a node (location in space at which NEURON
computes a solution), in units of arc length (normalized distance) along the currently
accesed section. Read about
for (var) stmt
here
http://www.neuron.yale.edu/neuron/stati ... r.html#for
You should probably also read chapters 5 and 6 of The NEURON Book.
http://ftp.neuron.yale.edu/ftp/ted/book ... xedref.pdf
http://ftp.neuron.yale.edu/ftp/ted/book ... xedref.pdf

3) I see how this demo, and the attshape.hoc code, calculate the attenuation over the entire dendrite. How do I calculate it for individual segments? It's pretty useless over a large dendrite. Can you give a simple explanation or example of how to access data in the individual segments of an existing imported morphology (instead of creating new segments)? Do I just defint nseg small enough that the new segments are always smaller than the imported segments?

It may be a good idea to rethink this from a different perspective.
Read chapter 5, and then consider the following:
Suppose you have a computational model of a neuron, which has been discretized
with a grid that is sufficiently fine to give good numerical accuracy, so that the model
has good resolution in both space and time.
Also suppose that the computational model preserves anatomically relevant structures,
i.e. for each unbranched neurite in the biological original, there is a corresponding
structure in the computational model (a section, using NEURON's terminology).
This is what you get if your computational model is created from the anatomical data
by CVAPP or NEURON's Import3D tool, and you have applied the d_lambda rule with
a sufficient discretization criterion (0.1 length constant at 100 Hz is generally adequate,
but it's always a good idea to check by asserting
forall nseg*=3
and doing a test run to be sure).
Such a model is both computationally efficient and easily managed.

Can you give a simple explanation or example of how to access data in the individual segments of an existing imported morphology (instead of creating new segments)?

The use of "segments" in this sentence is ambiguous. Do you mean that you don't
want to treat individual neurites as whole entities, which are then discretized as
described in chapter 5? Do you want instead to construct a computational model
in which the spatial grid has a node for each adjacent pair of xyzdiam measurements?
Does each adjacent pair of xyzdiam measurements really define an anatomically
or functionally relevant subdivision of the cell? Most measurement points are chosen
somewhat idiosyncratically by the person who slaved for long hours over the
microscope.

Do I just defint nseg small enough that the new segments are always smaller than the imported segments?

Why bother if the spatial grid selected by the d_lambda rule is adequate for accuracy?

[yml23]

Thanks for the clarifying information form your last post, Ted.
I've been trying different values of nseg:

v_init=-65
finitialize(v_init)
objref imp
imp = new Impedance()
FREQ = 100

WHERE = 0.1
soma distance(0, WHERE) // sets origin for distance calculations

proc calcZ() {local d, r1, r2
soma imp.loc(WHERE)
imp.compute(FREQ, 1)
rmin = 1
(x) ratio(x)"
forall {
print secname()
nseg = 1
print nseg
print "x distance(x) log(A)"
for(x=0; x<=1; x=(x+1/nseg)) print x, distance(x), log(1/imp.ratio(x))
}
}

calcZ()

The above code runs, but whenever I use a value of nseg that is not equal to 1, it crashes. What am I missing here?

[ted]

The above code runs, but whenever I use a value of nseg that is not equal to 1, it crashes. What am I missing here?

Beats the heck out of me why it runs at all. This line
(x) ratio(x)"
is mangled and should have generated an error message, but maybe it's flashing by
too fast to be seen. Until you're sure about your code, it's always best to work with a
toy model (one that has only a couple of sections).

But the bigger problem is: this code treats nseg as if it were an ordinary constant that
can be changed without consequence, like imin and imax in for i = imin, imax . . .

nseg governs the structure of the model's equations. Each section has its own value
of nseg that specifies how many ODEs are used to represent the spread of electrical
signals along that section. Changing nseg for any section will change the family of
equations that describe the model, so that you have to reinitialize the model and
compute a new simulation (or, in this case, a new set of impedances).

If you just want to access the values of range variables at section midpoints, don't
change nseg to t--the proper syntax is rangevar(0.5). If you want to see how nseg
affects the computed impedances etc., you must first assign the new nseg values,
then call finitialize(), and only then call calcZ.

[yml23]

With regards to your post form 8/31 - yes, I do want to create a node at every xyzdiam point. I have seen much discussion of nseg and d_lambda, but as yet I have not found a way to assign a node to every data point. I understand your suggestion about the d_lambda rule. However, I made these reconstructions specifically with diameter resolution in mind, and the fidelity of our model depends entirely on the accuracy of our morphological data. This cable transform will be used to create a model cell we will use to calibrate our amplifier response, so I would like to introduce as few a priori assumptions as possible and obtain the maximum fidelity to the morphological data. For our purposes, computation is not an issue - we only need to run this once to generate an equivalent cable.

Seems counterintuitive that it would be easier to go through an entire interpolation process with nseg and d_lambda than to simply use the data points from the .hoc file.

[ted]

I do want to create a node at every xyzdiam point.

This is the computational equivalent of empty resolution, but it's your decision. You'll
need to convert every adjacent pair of xyzdiam measurements into a separate
section. You'll have to write your own code to do that translation--but is the result
worth the effort? Have you run a test to see how this differs from letting NEURON
handle all this for you? Specifically, pick a neurite with lots of diameter variation along
its length, and implement two models of it: one that you build your way, and another
that you build with a single series of pt3dadd statements. Insert pas, set Ra and Cm
to reasonable values, and use the d_lambda rule to discretize both. Then compare
the electrical properties of the two models: DC input resistance at either end, or
amplitude and time course of the impulse response at, say, 0.2 ms or longer.

I made these reconstructions specifically with diameter resolution in mind, and the fidelity of our model depends entirely on the accuracy of our morphological data.

Good, but have you considered all sources of error? Like fixation artifact,
which is often anisotropic and can introduce ~5% error into measurements of length
and diameter. Light microscopy has serious limitations, especially when cells are
HRP-filled. Remember, 0.1 um uncertainty about the diameter of a 1 um neurite
means 10% uncertainty about surface area, and 20% uncertainty about cross-
sectional area--and a lot of cell membrane that is "electrically seen" from the soma
is out on those fine dendrites. Also, reconstruction systems often impose their own
"quantization error" and lower limit of resolution on measurable diameter changes.
Then there's the fact that cross-sections generally aren't cylindrical anyway, which
wreaks havoc with trying to calculate axial resistance--who knows what % error this
introduces, but some of the profiles that are seen on EM look pretty scary. Not to
mention the quesiton of what values to use for cytoplasmic resistivity and specific
membrane capacitance--reported values for both span quite a range.

You might want to try this expeiriment: make a simple passive model of your cell
(let Ra = 80, cm = 1, and use the d_lambda rule to handle the spatial grid).
Then use NEURON's Impedance tools to examine its electrical properties (here's a tutorial
http://www.neuron.yale.edu/neuron/stati ... zclass.htm).
Then change diameter by 0.1 um throughout and see what happens
forall for (x) diam+=0.1
or
forall for (x) diam-=0.1
Pay special attention to voltage transfer from dendrites to soma (Vin), and also
to input impedance (Zin), particularly in the distal dendrites.
Also see what happens when you perturb Ra or cm by +- 10%
forall for (x) cm*=1.1
or
forall Ra*=1.1
(Ra is a section variable, so you don't have to specify it at each internal node).

Seems counterintuitive that it would be easier to go through an entire interpolation process with nseg and d_lambda than to simply use the data points from the .hoc file.

You're "simply using the data points from the hoc file" if you just let NEURON read all
the pt3dadd statements for each section and take care of computing area, etc., itself,
instead of forcing it to wrestle with thousands of little sections.