Optimal d_lambda for measuring propagation speed

Anything that doesn't fit elsewhere.
Post Reply
davidsterratt
Posts: 28
Joined: Tue Jan 17, 2006 11:36 am

Optimal d_lambda for measuring propagation speed

Post by davidsterratt »

I am trying to measure the propagation speed in the myelinated axon model of McIntyre et al. (2002, http://senselab.med.yale.edu/ModelDb/sh ... model=3810).

I'm measuring the propagation speed by detecting the time the membrane potential passes -40mV using an APcount point process (code can be provided). I am using CVODE with cvode.atol = 0.001. I am looking at nodes away from the end of the axon to elimate end-effects.

I wanted to check that lambda was set to a sufficiently small value. Using the fixnseg.hoc script at http://www.neuron.yale.edu/neuron/stati ... ambda.html, I have changed d_lambda, and observed the propagation speed computed by dividing the distance between 8 nodes by the time taken for the AP to travel between them. Altering d_lambda, I get the following result:

Code: Select all

d_lambda        V (m/s)
0.002           43.45
0.005           49.31
0.01            54.25
0.02            58.81
0.05            55.69
0.1             57.39
0.2             57.39
I had expected that as I decreased d_lambda, the answer would start to converge, but it seems to be diverging. Is this behaviour to be expected? Admittedly, with d_lambda small, the number of compartments is rather large and impractical. However, it would be good to know if there is a reason known for this behaviour. Is there a lower limit on d_lambda?
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Optimal d_lambda for measuring propagation speed

Post by ted »

Strange results indeed. Do you get similar results if you simply
forall nseg*=3
a few times?

The code from ModelDB seems to have a relatively straightforward model setup section, but are you sure that the biophysical parameters are unaffected by changing nseg?

Code: Select all

I'm measuring the propagation speed by detecting the time the membrane potential passes -40mV using an APcount point process
Are you measuring the latencies at two different points, then calculating (distance between points)/(latency difference)? If not, that's the first thing I'd suggest doing. For a constant stimulus amplitude, changing discretization will change spike initiation latency, and this will affect conduction velocity estimates based on spike latency measured at a single point.

The actual location of spike initiation will depend somewhat on stimulus intensity relative to spike threshold, and in compartmental models the spike threshold is affected by spatial discretization. Use too large a stimulus, and it is possible for excess depolarizing current to spread farther into the axon, pushing the spike initiation zone farther from the stimulating electrode. A standard practice in experimental electrophysiology is to adjust stimulus intensity to twice the spike threshold. "Well, should I fiddle with stimulus intensity after every change in d_lambda?" I'd say no, just find what it should be with d_lambda = 0.1 and leave it at that while you try other values for d_lambda. That should be good enough to get reasonable qualitative results about the effect of d_lambda on conduction velocity.
I am looking at nodes away from the end of the axon to elimate end-effects.
Presumably you're stimulating at one end. Are the measurements being made far enough away from _both_ ends? One way to tell is with a space plot, using Keep Lines every so often to capture a few profiles of v. It would be nice to see at least one "complete spike waveform" between the two ends, and a nice long stretch in the middle where spike peak amplitude is the same from one node of Ranvier to the next.

An interesting phenomenon that has cropped up in some models is a more or less simultaneous spike intitiation along a long piece of axon. In unmyelinated axon this looks like a stretch of axon along which v is relatively flat for tens or hundreds of microns. In myelinated axon, it manifests as multiple adjacent nodes of Ranvier at which the spike waveform is nearly synchronous. This happens with particular combinations of parameters and stimulus conditions that favor a more or less uniform excitation along a long length of axon at nearly the same time. It's very easy to see in a space plot when you step through a simulation at 0.1 or 0.2 ms intervals.

Similar phenomena are seen in experimental observations of orthodromic spike intiation and propagation with voltage sensitive dyes of the axon initial segment and the first few hundred microns of the axon per se. In such cases the apparent conduction velocity is infinite. Of course, a slight asynchrony between nodes would produce a merely very large conduction velocity, and anything that impairs synchrony--whether it be a change of the stimulus amplitude, or a change of a model parameter (e.g. a conductance density), or a change of a simulation parameter (dt or spatial discretization)--would reduce the conduction velocity in the initial part of the axon.

Once you get away from the end of the axon where the stimulus is applied, ordinary propagation takes over and experimentally determined conduction velocity falls to "reasonable" levels.

So are you sure the axon is long enough to avoid "spike initiation zone effects"?
Is there a lower limit on d_lambda?
Whatever the lower limit on double precision floating point variable is.
davidsterratt
Posts: 28
Joined: Tue Jan 17, 2006 11:36 am

Re: Optimal d_lambda for measuring propagation speed

Post by davidsterratt »

Thanks very much for your very helpful reply Ted. I've tried responding to all your queries, but I've run out of time today, so the later ones aren't addressed very well.
Strange results indeed. Do you get similar results if you simply
forall nseg*=3 a few times?
It does vary a bit, though it looks more as though its converging. I'm using the file below, propspeed.hoc:

Code: Select all

load_file("MRGaxon.hoc")

// Switch on CVODE for greater accuracy
cvode_active(1)

// Record from nodes from injection site
objref apc[axonnodes]
for n=0,axonnodes-1 {
    node[n] apc[n] = new APCount(0.5)
    apc[n].thresh = -40
}

// Placeholders for vectors to store distance and time
objref s, time

// Specify which indicies count as the ends of the cable
objref endinds
endinds = new Vector(2, 1)
endinds.x(1) = axonnodes - 10 - 3

// Placeholders for vectors to store end distances and times
objref S, Time

// Graph to plot distance (s) against time (t)
objref g_st
g_st = new Graph(0)
g_st.view(0, 2, 0, deltax, 420, 400, 300, 300)

// Function to measure propagation speed
// Returns propagation speed in m/s
func measure_propspeed() {
    // Run simulation
    run()
    
    // Copy time and distance to vectors
    s = new Vector(axonnodes-10)
    time = new Vector(axonnodes-10)
    
    for i=0,axonnodes-10-1 {
        time.x(i) = apc[i+10].time
        s.x(i) = i*deltax
    }
    // Plot
//    g_st.erase()
    s.mark(g_st, time)
    g_st.exec_menu("View = plot")
    
    // Store start and end times
    S =    s.ind(endinds)
    Time = time.ind(endinds)
    
    // Plot fit
    S.line(g_st, Time)
    
    // Compute velocity
    dS = S.x(1) - S.x(0)
    dTime = Time.x(1) - Time.x(0)
    V = dS/dTime/1000
    
    return(V)
}

// User interface
V=0
xpanel("Propagation speed")
xbutton("Measure speed", "V = measure_propspeed()")
xfixedvalue("V")
xpanel(100,400)
If I just try what you suggested, I get the following error:

Code: Select all

oc>forall nseg*=3
oc>measure_propspeed()
IDA initialization failure, weighted norm of residual=1.96227
err=-6
/disk/scratch/sterratt//x86_64/bin/nrniv: variable step integrator error
 near line 5
 hoc_ac_ = istim
                ^
        fadvance()
      advance()
    step()
  continuerun(10)
and others
However, if I set the number of segments in the nodes to 1, things are OK:

Code: Select all

oc>measure_propspeed()
	57.148199 
oc>forall nseg*=3
oc>forsec "node" nseg=1
oc>measure_propspeed()
	50.086892 
oc>forall nseg*=3
oc>forsec "node" nseg=1
oc>measure_propspeed()
	48.607682 
oc>forall nseg*=3
oc>forsec "node" nseg=1
oc>measure_propspeed()
	48.160906 
oc>forall nseg*=3
oc>forsec "node" nseg=1
oc>measure_propspeed()
	48.318407 
oc>forall nseg*=3
oc>forall nseg*=3
oc>forsec "node" nseg=1
oc>measure_propspeed()
	48.35622 
By 9 segments per non-node segment, things seem to be converging.
The code from ModelDB seems to have a relatively straightforward model setup section, but are you sure that the biophysical parameters are unaffected by changing nseg?
I can't see from the code why they would be. None of the parameters depend on nseg, and when I change nseg, none of the parameters of a MYSA section change - I haven't checked others.
Are you measuring the latencies at two different points, then calculating (distance between points)/(latency difference)? If not, that's the first thing I'd suggest doing.
Yes, that's what I'm doing. I'm using the stimulus provided with the modeldb files, in which the midway node is stimulated, leading to an action potential propagating away in both directions. I'm measuring the latency between 1 node away and 8 nodes away from the injection node. A distance-time plot suggests that there are no initiation or end effects within this range.
Presumably you're stimulating at one end. Are the measurements being made far enough away from _both_ ends? One way to tell is with a space plot, using Keep Lines every so often to capture a few profiles of v. It would be nice to see at least one "complete spike waveform" between the two ends, and a nice long stretch in the middle where spike peak amplitude is the same from one node of Ranvier to the next.
If I do the space plot with keep lines, I see that there are six nodes between the spike initiation node and the end of the axon that have more-or-less the same peak amplitude. These nodes are the ones over which I'm measruing the speed.
An interesting phenomenon that has cropped up in some models is a more or less simultaneous spike intitiation along a long piece of axon. In unmyelinated axon this looks like a stretch of axon along which v is relatively flat for tens or hundreds of microns. In myelinated axon, it manifests as multiple adjacent nodes of Ranvier at which the spike waveform is nearly synchronous.
Time plots of the voltage at neigbouring nodes show that they aren't synchronous.
So are you sure the axon is long enough to avoid "spike initiation zone effects"?
I need to think more about this - I have to go now!
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Optimal d_lambda for measuring propagation speed

Post by ted »

There are many potential pitfalls in determining conduction velocity.

1. Spike initiation artifacts.
The stimulus used by the model authors is actually about 2*threshold, which is good for robust initialization. However, it is also large enough to dominate the rising phase of the spike at the stimulated node of Ranvier, and it spreads far along the axon, making a significant contribution to the early rising phase of spikes as far as three or four nodes of Ranvier away from the site of current injection. The peak amplitudes of the nodal spikes also change along the axon. A test with nseg = 27 throughout and dt = 0.001, with stimulus current injected at the 0 end of node[0] (and stimulus magnitude cut in half) reveals that spike peak amplitudes at nodes 0 - 3 are from ~0.2 to ~9 mV larger than in the middle of the axon.

2. Axon termination artifacts.
Simulation with nseg 27 and dt 0.001 reveals that sealed end effect is first noticable four nodes away from the terminal node.

So it's best to avoid the first four nodes and the last four nodes when trying to determine conduction velocity in this model. If the stimulus is applied to node[10], that doesn't leave a lot of nodes to work with.

But there is another concern as well--

3. Time measurement artifacts.
Solutions are computed at discrete points in time. Even if spike waveforms and amplitudes are identical, there is no guarantee that a v computed at any location will ever exactly equal the level that one might choose as the threshold for detecting spike passage time. Most measurements will be in error by an amount that could be almost as large as the time step. I see that you have chosen to set the threshold at -40 mV, which is a good choice since it is near the max slope of the rising phase. However, the spike passage times will only be first order correct--and the error will show up as slight scatter in conduction velocity measurements. Your spike passage times will be more accurate if you use the NetCon class's record() method, and specify cvode.condition_order(2). This will give you spike times that are interpolated, i.e. spike times with second order error. You'll want to read about NetCon record() and CVode condition_order() for more details.

So in summary I'd suggest moving the IClamp to node[0] (and cutting its peak current in half), avoiding nodes 0-4 and 16-20, and using NetCon record with second order threshold detection.

Meanwhile I'm going to run a few tests of my own to see whether the d_lambda strategy gives sufficient accuracy when those confounds have been eliminated.
davidsterratt
Posts: 28
Joined: Tue Jan 17, 2006 11:36 am

Re: Optimal d_lambda for measuring propagation speed

Post by davidsterratt »

Dear Ted,

apologies for not having replied sooner.

I have now made the changes you suggest, viz:
  • Stimulus now placed at node[0]
  • Stimulus strength now half of original
  • Measuring distance and time between node 5 and 15
  • Second order spike time detection using cvode.record()
Here is the code:

Code: Select all

load_file("MRGaxon.hoc")
//load_file("morphology.hoc")
load_file("propspeed.ses")

// Switch on CVODE for greater accuracy
cvode_active(1)

// Second order detection of thresholds
cvode.condition_order(2)

// Reduce stimulus current by half
istim = 0.5

// Record from nodes from injection site
objref null
objref nc[axonnodes], apt[axonnodes]
for i=0,axonnodes-1 {
    node[i] nc[i] = new NetCon(&v(0.5), null, -40, 0, 0)
    apt[i] = new Vector()               // Action potential times
    nc[i].record(apt[i])
}

// Placeholders for vectors to store distance and time
objref s, time

// Specify which indicies count as the ends of the cable
objref lim
lim = new Vector(2, 5)
lim.x(1) = axonnodes - 6

// Placeholders for vectors to store end distances and times
objref S, Time

// Graph to plot distance (s) against time (t)
objref g_st
g_st = new Graph(0)
g_st.view(0, 2, 0, deltax, 420, 400, 300, 300)

// Function to measure propagation speed
// Returns propagation speed in m/s
func measure_propspeed() {
    // Insert stimulus (needed if nseg in node has changed)
    stimul()
    
    // Move stimulus to node[0]
    node[0] stim.loc(0.5)

    // Run simulation
    run()
    
    // Copy time and distance to vectors
    s =    new Vector(axonnodes)
    time = new Vector(axonnodes)
    
    for i=0,axonnodes-1 {
        time.x(i) = apt[i].x(0)
        s.x(i) = i*deltax
    }
    // Plot
    //    g_st.erase()
    s.mark(g_st, time)
    g_st.exec_menu("View = plot")
    
    // Store start and end times
    S =    s.ind(lim)
    Time = time.ind(lim)
    
    // Plot fit
    S.line(g_st, Time)
    
    // Compute velocity
    dS = S.x(1) - S.x(0)
    dTime = Time.x(1) - Time.x(0)
    V = dS/dTime/1000
    
    return(V)
}

// User interface
V=0
xpanel("Propagation speed")
xbutton("Measure speed", "V = measure_propspeed()")
xfixedvalue("V")
xpanel(100,400) 

When changing nseg, I get the following results:

Code: Select all

oc>measure_propspeed()
	56.23173 
oc>forall nseg*=3
oc>forsec "node" nseg=1
oc>measure_propspeed()
	49.659936 
oc>forall nseg*=3
oc>forsec "node" nseg=1
oc>measure_propspeed()
	48.06019 
oc>forall nseg*=3
oc>forsec "node" nseg=1
oc>measure_propspeed()
	47.703338 
oc>forall nseg*=3
oc>forsec "node" nseg=1
oc>measure_propspeed()
	47.654984 
The numbers are up to 1 m/s smaller than before, but the pattern is still the same.

By reinserting the IClamp after nseg has been changed, I was also able to explore the effect of changing nseg uniformly, a bit:

Code: Select all

oc>measure_propspeed()
	56.30297 
oc>forall nseg*=3
oc>measure_propspeed()
	45.25735 
oc>forall nseg*=3
oc>measure_propspeed()
IDA initialization failure, weighted norm of residual=7.96376
err=-6
/disk/scratch/sterratt//x86_64/bin/nrniv: variable step integrator error
 near line 5
 hoc_ac_ = istim
                ^
        fadvance()
      advance()
    step()
  continuerun(10)
and others
There is quite a big difference between having nseg=1 in the nodes (49.7 m/s) and nseg=3 in the nodes (45.2 m/s). NEURON doesn't seem to like more than 3 segments in a node.

Here are the results when changing d_lambda:

Code: Select all

oc>load_file("fixnseg.hoc")
	1 
oc>measure_propspeed()
	56.30297 
oc>d_lambda=0.1
oc>geom_nseg()
oc>measure_propspeed()
	56.299369 
oc>d_lambda=0.05
oc>geom_nseg()
oc>measure_propspeed()
	54.724368 
oc>d_lambda=0.02
oc>geom_nseg()
oc>measure_propspeed()
	54.064961 
oc>d_lambda=0.01
oc>geom_nseg()
oc>measure_propspeed()
	53.552325 
oc>d_lambda=0.005
oc>geom_nseg()
oc>measure_propspeed()
	48.489984 
oc>print node[0].nseg
1 
oc>d_lambda=0.002
oc>geom_nseg()
oc>measure_propspeed()
/disk/scratch/sterratt//x86_64/bin/nrniv: subscript out of range x
 near line 21
 hoc_ac_ = istim
                ^
        measure_propspeed()
I had some previous runs in which I didn't have 2nd order spike detection, in which I was able to get a result with d_lambda=0.0002. The velocity then was 43.2 m/s, a considerable drop from 48.4 m/s. The number of segments in each node as 3 with this value of d_lambda.

In summary, it seems that the number of segments in each node has quite a bit effect on how fast the propagation is. I suppose this isn't that surprising, as propagation through a single segment is, almost by definition, instantaneous. However, I'm still not clear whether the nodes should have one segment or three.
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Optimal d_lambda for measuring propagation speed

Post by ted »

Well, David, this is one pain in the neck model.
In summary, it seems that the number of segments in each node has quite a bit effect on how fast the propagation is. I suppose this isn't that surprising
It surprises me. Just to be sure there isn't something strange in your code that I'm missing, I trimmed away all but the "specification of biological properties" from the ModelDB entry, wrote some simulation control and data analysis code of my own, and did a few tests. I did this partly because the model has nonuniform membrane properties, and although resting potential is nearly uniform, it takes about 50 ms of model time with dt = 0.025 to reach steady state.

After a proper initialization, v at the nodes is slightly more depolarized than in the internodes. This difference varies along the length of the axon, and is most prominent at the two ends (with the spatial grid specified by d_lambda = 0.1, v ~ -79.955 at node 0 and 1, ~ -79.96 at node 10). I thought this might make a difference in a model that has a persistent sodium current, as this one does. I used a custom proc init that allows the model to reach steady state before starting the actual run. Also, in addition to detecting spike times at nodes 5 and 15, the instrumentation includes spike detection at nodes 8 and 12, between which resting potential is closer to being uniform.

But no, the qualitative finding is much the same as yours: discretization of the nodes of Ranvier has a big effect on simulation results. That said, discretization of time has a big effect too. Here are the conduction velocities, in mm/second, from three sets of simulations (using second order spike time detection of course):
--one set with fixed dt ranging from 0.01 to 1e-5 ms
--one set with adaptive integration with atol ranging from 1e-3 to 1e-5 (using the Atol Scale Tool to automatically adjust the scaling of individual state tolerances)
--one set with adaptive integration (atol 1e-3) using d_lambda 0.1, 0.01, and 0.003

Code: Select all

Fixed time step
        nodes            8-12
dt       5-15    8-12   delta
1e-2    47917   46000
1e-3    55288   54762    8762
1e-4    56152   56166    1404
1e-5    56237   56235      69

Adaptive integration
        nodes
atol     5-15    8-12
1e-3    56249   56247
1e-4    56248   56245
1e-5    56249   56245

Back to atol 1e-3
          total   nodes            8-12
d_lambda   nseg    5-15    8-12   elta
0.1         461   56249   56247
0.01       3901   53525   53524    2723
0.003     12541   48375   48377    5147
Note: up to this point, nseg was 1 for all nodes of Ranvier.
0.003     12583   43642   43642    4735  node nseg = 3
0.003     12625   42649   43648    4741  node nseg = 5
And, as you found, setting node nseg = 9 generated an error message.
I wasn't surprised that making dt smaller sped up the spike, or that making the grid finer in the internodes reduced the conduction velocity. And it was reassuring to see that adaptive integration with atol = 1e-3 was sufficient for temporal accuracy. But the effect of changing nseg in the nodes was very surprising.

I thought that the persistent sodium current might be the culprit (affecting dynamics in the subthreshold region where other currents are very small), but reducing its conductance density to 0 in all nodes had only a small quantitative effect and absolutely no significant qualitative effect.
I'm still not clear whether the nodes should have one segment or three.
I'd say they need more than 3. I'd increase nseg as much as necessary to see velocity converge to some value. IDA initialization failure prevents use of nseg = 7 or larger for the nodes of Ranvier with adaptive integration, but that doesn't preclude the use of fixed time steps. Fixed dt = 1e-5 is almost as good as adaptive integration, and works with node nseg = 9 and d_lambda = 0.003 elsewhere--gives conduction velocities of 41974 (nodes 5-15) and 41975 (nodes 8-12), but takes 1085 s to complete a 0.9 ms simulation* on my PC, not counting time to initialize to steady state (but that's done with dt = 0.025 so is much quicker). Is this close enough, will a further reduction of dt or increase of node nseg affect conduction velocity significantly? Left as an exercise to the reader . . .

*--since the model is initialized to steady state by jumping to -100 ms and advancing to -50 ms before resetting to 0 ms, I can have the IClamp start at 0.1 ms (to provide a short pre-stimulus baseline) and terminate the simulation after only 0.9 ms yet still capture the entire rising phase and the first 1/3 of the falling phase of the spike at node 20.
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Optimal d_lambda for measuring propagation speed

Post by ted »

An addendum: there are a few minor bugs in the NMODL code for the axnode mechanism. The first is in FUNCTION vtrap2, whose original source code is

Code: Select all

FUNCTION vtrap2(x) {
        if (fabs((x+bmpB)/bmpC) < 1e-6) {
                vtrap2 = -bmpA*bmpC
        }else{
                vtrap2 = (bmpA*(-(x+bmpB))) / (1 - Exp((x+bmpB)/bmpC))
        }
}
This has a discontinuity at the critical point, i.e. when the argument equals -bmpB. The assignment statement

Code: Select all

                vtrap2 = -bmpA*bmpC
should be changed to

Code: Select all

                vtrap2 = bmpA*bmpC
i.e. the minus sign should be removed.

This error caused reversal of the sign of return value of vtrap2(v) at a single membrane potential, producing an isolated discontinuity in the voltage dependence of mp_inf, i.e. at v = -bmpB.

Similar errors are present in FUNCTIONs vtrap7 and vtrap8. These have similar effects (an isolated discontinuity at one membrane potential) and can be fixed in a similar way: eliminating the - sign in the first clause of the FUNCTION's conditional statement.

This is not a fatal problem and would not cause a serious qualitative difference in any simulation in which membrane potential is free to change. However, I will notify the authors and the curator of ModelDB.

The results reported in my previous post were generated with a revised version of this mechanism.
davidsterratt
Posts: 28
Joined: Tue Jan 17, 2006 11:36 am

Re: Optimal d_lambda for measuring propagation speed

Post by davidsterratt »

Thanks very much for your help Ted - I'm glad it's not just me who's finding the model tricky! I've modified the AXNODES mod file as you suggest.

I've set up the initialisation that you suggest, and have now just finished running a simulation with a fixed time step of dt=1e-4, d_lambda=0.003 and nseg=9 in nodes. I found that I needed a stimulation amplitude greater than 1nA, so used the 2nA in the original model. With this setup, I got a propagation velocity of 41910 mm/s (nodes 5-15), which is not far off what you obtain with dt=1e-5 (41974 mm/s). The simulation time was 494.29s, including initialisation.

I still need to do simulations with dt=1e-5. I tried these earlier today, but they took a long time, and I wasn't putting in enough current, so I didn't get any spikes. I will run one later.
davidsterratt
Posts: 28
Joined: Tue Jan 17, 2006 11:36 am

Re: Optimal d_lambda for measuring propagation speed

Post by davidsterratt »

After a long time, I've had reason to look at this again. I've now got graphs (and a spreadsheet) of the propagation speed at various values of dt and d_lambda which I could share. However, the picture is much as described earlier in this thread. There is a difference in conduction velocity of about 2.5% between the dt=0.0001 and dt=0.003, but a difference of roughly 25% between d_lambda=0.0003 and d_lambda=0.1. The effects of changing dt and d_lambda are roughly additive, and below d_lambda=0.01, the conduction velocities appear to be converging.

I've also been looking at the Moore & al (1978) model of myelinated axon (http://senselab.med.yale.edu/modeldb/sh ... model=9851), which doesn't use the extracellular mechanism. I've found that the velocity depends very weakly (if at all) on the number of segments in the nodes:

Code: Select all

oc>load_file("fig2.hoc")
oc>forall nodes nseg=1
oc>velocity()  
	23.786949 
oc>forall nodes nseg=3
oc>velocity()
	22.782776 
oc>forall nodes nseg=9
oc>velocity()
	22.658735 
oc>forall nodes nseg=27
oc>velocity()
	22.648105 
oc>forall nodes nseg=81
oc>velocity()
	22.646899 
I'm therefore wondering if the need to go to very small compartment sizes in the McIntyre & al model is a function of the extracellular
mechanism. If so, is this a fundamental property of the ODEs that have to be solved, or could it possibly be a bug in how the extracellular mechanism is used in this model or how it's implemented in NEURON? I was prompted to think about bugs by reading the documentation of the extracellular mechanism in the Programmers' reference!

Assuming that the very small d_lamnbda and dt values do give the correct answer, I suspect it will be possible to recover the original conduction velocities by scaling up the conductances in the node of Ranvier by a uniform factor. However, I haven't looked at this systematically.
hines
Site Admin
Posts: 1682
Joined: Wed May 18, 2005 3:32 pm

Re: Optimal d_lambda for measuring propagation speed

Post by hines »

Although the concept of d_lambda is sound, the implementation is useless in this model because of the neglect of xraxial.

With
cvode.dae_init_dteps(1e-14)
one gains enough initialzation power with IDA to allow forall nseg = 81 and cvode.atol(1e-4).

Looking at the trajectory for node[2].v(0.5) it looks like somthing between linear and quadratic convergence in velocity is obtained with
forall nseg = 9, 27, 81

It is a bit of a puzzle to me why so many segments are required for convergence, why it is not quadratic, and whether only certain regions of abrupt changes in xraxial need a lot of segments.

I found it helpful to use a tstop of 1.5 and look at spaceplots of v and vext in the neighborhood of a node.
As I get time. I'll continue to think about this.
davidsterratt
Posts: 28
Joined: Tue Jan 17, 2006 11:36 am

Re: Optimal d_lambda for measuring propagation speed

Post by davidsterratt »

Many thanks for the reply Mike.

Of course the d_lambda won't work because of xraxial (I should have thought of this!). I suppose one might be able to work out an equivalent with a bit of algebra, but perhaps the best approach is just to increase nseg and the time resolution until convergence is reached.

I hadn't managed to get cvode working, but I hadn't been using the dae_init_eps() method to initialise.

Looking at space plots of v and vext is informative. You can see the voltage differences between either end of the node (~0.15mV during the rising phase of the action potential). There is no difference in vext. It takes order of 0.2ms for the peak of the AP to pass through the node.

In short, I'm now satisfied that there are arguments for using what seems like ridiculously large numbers of segments in nodes (and perahps elsewhere). To round things off, what I would like to check is that the cvode method and the dt methods converge on the same conduction velocity. I'll let you know how I get on.
davidsterratt
Posts: 28
Joined: Tue Jan 17, 2006 11:36 am

Re: Optimal d_lambda for measuring propagation speed

Post by davidsterratt »

I've now derived the conduction velocity from simulations using fixed dt and CVODE methods for various numbers of segments per section.

Fixed dt:

Code: Select all

      dt nsegall       V
1  1e-03       1 55.0239
2  1e-03       3 44.5736
3  1e-03       9 41.2186
4  1e-03      27 40.0697
5  1e-03      81 39.9306
6  1e-04       1 56.0976
7  1e-04       3 45.1512
8  1e-04       9 41.6818
9  1e-04      27 40.6792
10 1e-04      81 40.4218
11 1e-05       1 56.2072
12 1e-05       3 45.2168
13 1e-05       9 41.7347
14 1e-05      27 40.7382
15 1e-05      81 40.4901
CVODE:

Code: Select all

  cvodeatol nsegall       V
1     1e-04       1 56.1773
2     1e-04       3 45.2875
3     1e-04       9 41.7438
4     1e-04      27 40.7437
5     1e-04      81 40.6401
At all values of nseg the agreement between CVODE and fixed dt=1e-5 is within 1%. For both fixed-dt and CVODE, the velocity has more-or-less converged by nseg=27. I'm therefore happy to use the CVODE method, probably with nseg=81 for good measure.

The reason that CVODE hadn't been working for me was that I'd been using

Code: Select all

cvode.condition_order(2)
and it had been causing NEURON to crash.
nikkip
Posts: 28
Joined: Mon Jun 24, 2013 6:09 pm

Re: Optimal d_lambda for measuring propagation speed

Post by nikkip »

Setting cvode.condition_order(2) also causes my NEURON to crash. Did you ever resolve this issue?
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Optimal d_lambda for measuring propagation speed

Post by ted »

Comments on spatial discretization and the d_lambda rule

As Michael Hines noted, the problem may be related to the fact that the calculation of lambda ignores xraxial.

The d_lambda rule was originally conceived in the context of modeling unmyelinated neurites, and the implementation of it used in the CellBuilder and elsewhere is couched in terms of the factors that govern the length constant of a simple cable at frequencies well above the membrane time constant: neurite diameter, cytoplasmic resistivity, and specific capacitance of the cell membrane.

Code: Select all

forsec all { nseg = int((L/(0.1*lambda_f(100))+.999)/2)*2 + 1  }
where lambda_f(freq) is the length constant of a neurite in um at frequency freq (Hz) and is defined in nrn/lib/hoc/stdlib.hoc

However, some representations of myelinated axons are actually "nested" cables in which internodes have two or more parallel paths for axial current flow--in the simplest case, along the interior of the axon itself, and along a thin cylindrical gap that lies between the axon membrane and the myelin. The resistance/length for axial current through the latter path may be comparable to, or larger, than it is for axial current through the interior of the axon, and this will affect the eigenvalues that govern signal attenuation along the internode.

The question is how to take this into account. In a myelinated internode with such a gap, these three elements
interior of the axon
axon membrane
gap between axon membrane and myelin
can be viewed as constituting a cable with specific membrane capacitance cm (specific membrane capacitance of the axon itself), diameter diam (diameter of the axon itself), and total resistance/length (when one considers a closed current loop) equal to the sum of cytoplasmic longitudinal resistance/length and extracellular longitudinal resistance/length, that is
Ra*4/(PI*diam^2) + Rag/(w*PI*diam^2)
where
-- Ra and Rag are the resistivities of the cytoplasm and the fluid in the gap between axon membrane and myelin
-- diam is axon diameter (not including the thickness of the gap and myelin)
-- w is the ratio wgap/diam where wgap is the width of the gap beteween the axon's membrane and the myelin.

So the total resistance/length is
Ra*4/(PI*diam^2) + Rag/(w*PI*diam^2) = (1 + Rag/(4*w*Ra)) * Ra * 4/(PI*diam^2)
which is the same as the resistance/length for a simple cable with diameter diam
and cytoplasmic resistivity k*Ra where k = 1 + Rag/(4*w*Ra)

This suggests basing spatial discretization of an internode on the AC length constant calculated for an axon with diameter diam, specific membrane capacitance cm, and a "scaled" cytoplasmic resistivity that is equal to k*Ra. Here's a bit of code that implements this strategy:

Code: Select all

/* assumes:
--sralist is a SectionList ("sra" for "scaled Ra") whose elements are internode sections 
  that have a gap between axon and myelin
--each section in sralist is a cylinder
--wgap is the same for all */
proc chop_internodes() { local tmp, w
  forsec $o1 {
    tmp = Ra
    w = wgap/diam(0.5)
    Ra = Ra + Rag/(4*w)
    nseg = int((L/(0.1*lambda_f(100))+.999)/2)*2 + 1
    Ra = tmp // restore it for possible future use
  }
}
chop_internodes(sralist)
This can of course be recast in a form that uses xraxial instead of Rag, wgap, and diam, but for now I will leave that to the interested reader.
Post Reply