nseg odd?

Extending NEURON to handle reaction-diffusion problems.

Moderators: hines, wwlytton, ramcdougal

Post Reply
gartland

nseg odd?

Post by gartland »

Hi,
I noticed in one of your Hot Tips posts (http://www.neuron.yale.edu/phpBB/viewto ... f=28&t=562) you wrote:
Note that chemical signals tend to have a much shorter effective length constant than
electircal signals do. Therefore nseg will probably have to be much larger than if you were
interested only in electrical signals. As always, choose odd values for nseg, and test
for adequate spatial accuracy by following this protocol:
It make sense to me why one would need a finer spatial scale to study diffusion since it's length constant will most likely be on a fine spatial scale, however why does the number of segments need to be odd?

I have had some issues with modeling diffusion and I have convinced myself (hopefully incorrectly) that the problem is due to adding too many segments (or shells in radial diffusion). Without getting in to specifics (yet :) ), I was thinking that the model might be constrained by the fact that there can only be flux between two adjacent compartments. If the segments get too thin then you'd actually want to have molecules diffusing through more than one segment in less than a single time step. I guess I was thinking that perhaps the model wasn't appropriately shrinking the time-steps to deal with the thin segments. Now that I write that, I bet it does do that appropriately, so maybe I'm just not interpreting my results correctly? Still, what about the case where a diffusion coefficient is very large (simulating instantaneous equilibrium), will the time-step shrink enough so that an ion could equilibrate fast enough in all the compartments? For my own sanity I would love to know that NEURON will perform correctly in these cases (motivation below).

Let's consider the radial diffusion of calcium with influx and eflux happening only in the outer shell:

1) [Ca] transients should be the same in the case where nshells = 1 and nshells = 100 with diffusion coefficient = very large

2) As nshells approaches a large number, the results should approach a single solution

Are these tests reasonable? Does the calcium diffusion model in the NEURON book pass these tests? I can't put my finger on the code I used to test these right now, but I seem to remember that these tests do not pass. I'll run them again if you're interested though...

Now I'm running into discrepancies in similar controls for longitudinal diffusion and I'm wondering if I'm just confused, or if there really could be problems with these "control experiments" that I'm trying. Anyone have thoughts about these questions? Any insight is welcomed. Thanks!
ted
Site Admin
Posts: 5810
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: nseg odd?

Post by ted »

gartland wrote:why does the number of segments need to be odd?
nseg doesn't absolutely have to be odd. But it should be. See
Why should I use an odd value for nseg?
in the FAQ list http://www.neuron.yale.edu/neuron/faq/general-questions
I have had some issues with modeling diffusion and I have convinced myself (hopefully incorrectly) that the problem is due to adding too many segments (or shells in radial diffusion).
Maybe. In order to describe concentration with a continuous variable, and avoid worrying about whether or not the movement of individual particles must be tracked, the spatial grid for diffusion should be coarse enough that each compartment contains at least a few hundred ions/molecules/whatever it is that is diffusing. For smaller numbers of diffusing entities, concentration can still be described with a continuous variable but its value requires a statistical interpretation (very low concentrations represent "averages" of what you would observe if you could perform repeated experiments with all other variables remaining the same).
I guess I was thinking that perhaps the model wasn't appropriately shrinking the time-steps to deal with the thin segments.
NEURON does nothing to time steps that you don't tell it to. If you are using adaptive integration, NEURON chooses dt according to the rules you specify.
what about the case where a diffusion coefficient is very large (simulating instantaneous equilibrium), will the time-step shrink enough so that an ion could equilibrate fast enough in all the compartments?
The adaptive integrators are quite fastidious about changing dt (and integration order) as needed to satisfy the user-specified error tolerances. Use the GUI's VariableStepControl, and be sure to click on the "Atol Scale Tool" button and perform and Analysis Run, then click on the Rescale button so that tolerances are automatically rescaled for each state variable. Then run a series of simulations in which you make the error criterion ("Absolute Tolerance") progressively smaller from run to run, and compare results. Stop when you get to Absolute Tolerances such that you don't see a significant difference from run to run.
Does the calcium diffusion model in the NEURON book pass these tests? I can't put my finger on the code I used to test these right now, but I seem to remember that these tests do not pass. I'll run them again if you're interested though...
If you're interested in using NEURON in your work, maybe you'd like to resolve the question you just raised before investing further effort with NEURON in other tasks.
Now I'm running into discrepancies in similar controls for longitudinal diffusion and I'm wondering if I'm just confused, or if there really could be problems with these "control experiments" that I'm trying.
I'm not much of a mind reader. The discussion will have to be more specific before I can say any more.
gartland

Re: nseg odd?

Post by gartland »

Your post was very helpful, thank you. Naively, I hadn't realized that I should be monitoring/adjusting the tolerances as I wrecklessly made compartments extremely small and made diffusion extremely fast. Theoretically, the tests I referred to above should pass for a model of multi-compartment diffusion, but practically they are impossible to test unless I change the tolerances in extreme ways as well.
ted
Site Admin
Posts: 5810
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: nseg odd?

Post by ted »

Error tolerances, whether absolute or relative, are not affected by compartment size.
gartland

Re: nseg odd?

Post by gartland »

Thanks. Yes that makes sense to me now.
Is it possible that my problem was that I wasn't running my model with an adaptive integrator like CVODE, so that the time steps were not shrinking to meet the needs of the tolerances I have specified? I noticed today that although one MOD file uses a kinetic block and reports no errors when it is compiled, the MOD file that I was using to generate spiking currents in the cell reports that it is incompatible with CVODE when I compile it. I guess this means that the whole simulation must be running with the backward Euler method then? I just went back through the chapter about this in the NEURON book a bit which was very helpful, but now I'd like to be able to check what method NEURON is using with my model. Is there a way to get that information from the hoc command line (I'm running NEURON through python for better or for worse so i don't have access to the GUI)?

Also is it possible to see what the actual time step widths used over the course of the simulation as in the figures in the book?

Previously when my model was running with the fixed step integrator, was the time step set by me specifying "dt" before initialization? I remember reading at some point that changing "dt" should only change the vectors that i get when i use the "record" method, but is this only true when using an adaptive integrator?

Feel free just to point me towards a section of the book or the web site. I really appreciate your help though!
ted
Site Admin
Posts: 5810
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: nseg odd?

Post by ted »

gartland wrote:Is it possible that my problem was that I wasn't running my model with an adaptive integrator like CVODE, so that the time steps were not shrinking to meet the needs of the tolerances I have specified?
Maybe. Pitfalls lurk everywhere.
the MOD file that I was using to generate spiking currents in the cell reports that it is incompatible with CVODE
Sounds like you fell victim to some ancient NMODL code. There's a lot of it floating around. Still works (which is good news for the old hoc files that rely on it) but should be avoided in new program development.
I guess this means that the whole simulation must be running with the backward Euler method then?
I guess. You can easily find out by watching what happens to dt in the course of a simulation.

If you need to fix this (or are merely curious), look for a mod file with a BREAKPOINT block that contains a statement like
SOLVE procedurename
In that file you will also find a
PROCEDURE procedurename
block that will contain statements similar to
m = m + (minf-m)*exp(-dt/mtau)
What you need to do is either
(1) forget about CVODE and stick with fixed dt
(2) replace this mod file with a mod file that implements the identical voltage-gated current but in a way that is CVODE-compatible
(3) convert this mod file to a CVODE-compatible form (see http://www.neuron.yale.edu/neuron/stati ... l#Channels)--or get someone to do it for you
(4) if all you need is a pair of currents that generate spikes, find a different spike generating mechanism that is CVODE-compatible
I'd like to be able to check what method NEURON is using with my model. Is there a way to get that information from the hoc command line
Examine the Programmer's Reference documentation of the CVode class. Specifically, read about the method called current_method().
I'm running NEURON through python for better or for worse so i don't have access to the GUI
Really? What if you execute the hoc statement
load_file("nrngui.hoc")
from Python (properly wrapped, that is, like any other hoc statement you'd call from Python). Do it before the statements that create sections.
is it possible to see what the actual time step widths used over the course of the simulation as in the figures in the book?
Use the Vector class's record() method to capture t, then after the end of the simulation clone the Vector and apply the Vector class's deriv() method to the clone. In hoc this would read
// statements to execute prior to run()
objref tvec
tvec.record(&t)
. . .
// statements to execute after run()
objref dtvec
dtvec = new Vector()
dtvec = tvec.c
dtvec.deriv(1) // now dtvec.x(i) = tvec.x(i) - tvec.x(i-1) for i = 1..tvec.size()-1
Be sure to read about the Vector class and these methods in the Programmer's Reference. Or at least about the record() method, if you're more comfortable with vector manipulations in Python.
Previously when my model was running with the fixed step integrator, was the time step set by me specifying "dt" before initialization?
Yes, but you can easily verify this for yourself (record t to a Vector, etc.).
I remember reading at some point that changing "dt" should only change the vectors that i get when i use the "record" method
Never heard of this, don't even know quite what it means. If you can point me to where you found this, I'd like to see it for myself. Probably needs to be revised for clarity.
gartland

Re: nseg odd?

Post by gartland »

I've been busily trying many of the suggestions you had, but I am stuck again. My general line of thinking is this:

I have a working model with longitudinal diffusion using the backwards Euler method of integration. To check the model accuracy I increase nseg by nseg * 3 iteratively but the model does not approach any consistent result. I blamed this on the fact that I was using a fixed time step and that really I probably needed to be increasing time resolution with spatial resolution in order to have the model approach stable, accurate, results.

Rather than changing both nseg and dt together I figured it would be a good idea to let CVODE change dt for me according to reasonable tolerances set in my MOD file.

What does work for me: A single compartment model with a single calcium channel mechanism voltage stepped with a SEClamp run with either CVODE or backward Euler (confirmed by recording variable t-steps and calling current_method())

Also works: A single compartment (multi-segment) model with the calcium channel and a simple KINETIC based model without LONGITUDINAL_DIFFUSION with both CVODE and b Euler

Also works: Same as above but now add LONGITUDINAL_DIFFUSION. This ONLY works with the backward Euler method. When I turn on CVODE i get a "Segmentation fault" right after i call "neuron.run(tstop) " If i either comment out the LONGITUDINAL_DIFFUSION lines in the mechanism or call CVode().active(0) then the model runs again.

Any ideas what I am running into? Should I try to alter the run loop myself to debug and find what is causing the Seg fault? Thanks for your thoughts!
ted
Site Admin
Posts: 5810
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: nseg odd?

Post by ted »

Sounds like I need to be able to try your model myself. Zip up just enough hoc, ses, and mod files to reproduce your model and send the archive to me
ted dot carnevale at yale dot edu
and I'll tell you what I find out. If you have any question about what to send and what to omit, see
What to include in a zip file, and what to leave out
viewtopic.php?f=28&t=560
gartland

Re: nseg odd?

Post by gartland »

Thank you for your help resolving this issue! After switching NEURON to a variable step integrator (CVODE) the simulation seems to converge on a result with relatively few increases in segment number. I had problems getting CVODE to work correctly. You suggested that I upgrade to the newest version of NEURON because some bugs with longitudinal diffusion had been found (specifically to the python version i guess?). When I recompiled, the model ran fine with CVODE so other python users, take note! If anyone has any questions about these or similar issues feel free to message me! Thanks again!
ted
Site Admin
Posts: 5810
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: nseg odd?

Post by ted »

gartland wrote:because some bugs with longitudinal diffusion had been found (specifically to the python version i guess?)
Nope. The bugs were in the computation of longitudinal diffusion, quite apart from anything to do with python.
Post Reply