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
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
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.