First let me apologize for this delayed reply.
I've zipped my file (for the first two figures)and sent to your email
I received your email message, but it had no attachment. Try again?
The model used by Waxman et al. 1978 is based on what they used for
Brill, Waxman, Moore, and Joyner (1977)
Conduction velocity and spike configuration in myelinated fibres: computed dependence on internode distance.
J. Neurology, Neurosurgery, and Psychiatry. 40: 769-774
ModelDB has a re-implementation, by Michael Hines for NEURON, of the Brill et al. model; its accession number is 9848. I downloaded that and found that, while it reproduced the qualitative results in Figure 1 of the Brill paper, the reproduction was not quantitatively exact.
The same is true for Figure 2. The README file of ModelDB entry 9848 mentions reasons why Figure 2 might be difficult to reproduce exactly with NEURON, but what accounts for the inability to exactly reproduce the results in Figure 1?
The source code for ModelDB entry 9848 is rather complex and very clever, since it was designed to demonstrate qualitative reproduction of the results in Figure 1 and 2. To avoid that complexity and make it easier to examine the differences between Figure 1 and the NEURON model's results, I decided to re-implement the model myself.
Immediately I discovered that the model specification in Brill et al. was incomplete. For one thing, the authors changed the actual structure of the model, depending on internode length L. Quoting from the bottom of page 770 to the top of 771
"For an internodal length greater than 1000 um, a steady velocity of propagation obtained within two nodes of either end, and for these computations we used only 12 nodes. As L was decreased from this value, the number of nodes was increased in each case (up to 50 for an L of 25 um) so as to obtain a sufficient region of uniform velocity of propagation."
But they don't say how many nodes (and internodes) were actually used for models with L of 50, 100, or 200.
Brill et al. don't describe the stimulus they used. Waxman and Brill 1978 used "a twice-threshold current for 200 us" so maybe Brill et al. did the same. But internode length affects the amount of current required to trigger a spike. Did Brill et al. readjust the stimulus current amplitude every time they changed internode length? They don't say.
Brill et al. used a 5 us integration time step, so the precision of their measurements of ICT were limited to 5 us. Internode conduction time decreases with decreasing internode length, and for the model they
used, ICT is only about 5-6 us for a 25 um internode and about 20-25 us for a 200 um internode. There is no way they could have obtained the results shown in Fig. 1B unless they estimated threshold crossing times by applying linear interpolation to plots of v vs. t. That would be OK from a technical standpoint, because using the Crank-Nicholson method to compute v produces results that have second order accuracy in time, which means that values of v at times that lie between multiples of dt can be estimated by linear interpolation, and such interpolated values themselves have second order accuracy in time. BUT Brill et al. say nothing about using interpolation!
I should also mention the parameter that Brill et al. (and Waxman and Brill 1978) call "axoplasmic resistance per unit axon length" and represent as ra (see the Table on page 770). In a footnote on page 770, Brill et al. say it was "Calculated from specific axoplasmic resistance of 100 ohm-cm." You might think that this meant you should set Ra in your NEURON model to 100. But no, that would be incorrect. Why?
ra = 4*Ra / (PI * diam^2)
where ra is "resistance per unit length of axon" and Ra is "cytoplasmic resistivity." Solving for Ra:
Ra = rax * PI * diam^2 / 4
For rax = 1.26e8 ohm/cm, and diam = 10 um, the actual equivalent value of Ra to five decimal places is 98.960 ohm cm. If your own NEURON implementation set Ra to 100 ohm cm, that might contribute to differences between your model's results and the results published by Waxman and Brill. It wouldn't account for differences between Fig. 1 of Brill et al. and the results produced by Michael Hines's re-implementation of that work, because Hines set Ra to 98.96... ohm cm.
What about the fact that Brill et al. used a "modified" Crank-Nicholson integration method? It seems unlikely that the algorithm was identical to what NEURON uses when secondorder == 2. And did it make the same assumptions about boundary conditions between adjacent nodes and internodes that NEURON does?
I see that in Michael Hines's re-implementation for NEURON, the parameters used to reproduce Fig. 1 qualitatively are stimulus 10 nA pulse for 0.1 ms, dt == 0.01 (not 0.005), and ICT is calculated from interpolated threshold crossing times. Also, his code uses 50 nodes and 49 internodes, regardless of internode length. Even so, the results are qualitatively very similar to the published Fig. 1. My own
re-implementation, which also uses 50 nodes and 49 internodes, automatically adjusts stimulus amplitude to 2 x threshold every time internode length is changed. Also, I measured the ICT for the internode halfway down the length of the model axon, but didn't bother to extrapolate the exact time
at which the spike arrived at the two ends of that internode (just used a couple of APCount to detect the spike). For internode lengths > 200 um, the results I got for ICT and conduction velocity were very close to those produced by Michael's implementation.
So there are a lot of questions about what Brill et al. actually did, and neither Michael Hines nor I were able to reproduce their results exactly. But does it matter? Remember what Richard Hamming wrote:
The purpose of computing is insight, not numbers. Frankly, I think it's pretty amazing that the results generated by the code that Brill et al. used are so close to those that are generated by the NEURON re-implementations of their model. In particular, the NEURON re-implementations produce results that verify the qualitative features that underlie the conclusions reported by Brill et al..
If I were launching a new study based on what Waxman and Brill did, here's what I'd do. I'd keep the basic outline of the Brill et al. and Waxman and Brill models--axon diameter, node lengths, internode lengths, membrane properties of the nodes and of the myelinated and demyelinated regions. I'd definitely set Ra to 100 ohm cm, because the footnote on page 770 proves that's the value they intended to use. And instead of setting internode nseg to 10, I'd use an odd value--because you're using NEURON, not whatever FORTRAN code that Brill et al. were using--and base nseg on the d_lambda rule.
Also, I'd abandon Crank-Nicholson and use adaptive integration (cvode) instead. Why? Because adaptive integration makes it much easier to get high quality estimates of internode conduction time. In my own re-implementation, which used Crank-Nicholson with dt = 5 us, I used APCounts to detect the times of "threshold crossings" at the 0 and 1 ends of an internode. That was easy, but precision was limited by dt. Michael Hines's implementation also used Crank-Nicholson, but used Vector recording to capture the actual time course of v at the 0 and 1 ends, then used linear interpolation to refine the threshold crossing times. That produced much more accurate internode conduction times, but complicated his program.
The smart way to detect threshold crossing times is to use a pair of NetCons--one at the 0 end of the internode, and the other at the 1 end--and then use the NetCon class's record() method to detect the threshold crossing times. Read about NetCon.record() here
https://www.neuron.yale.edu/neuron/stat ... Con.record
"How does that help? Cvode doesn't force the simulator to find the exact time at which v crosses the 'spike time' threshold."
True. By default, detection of threshold crossings happens at time step boundaries. However, if you call the Cvode class's
condition_order() method with an argument of 2, threshold crossing times will be found by linear interpolation within the integration step in which the crossing happened. So just call
before you run the first simulation, and every run after that will use linear interpolation to produce refined estimates of threshold crossing times. Assuming that you are using adaptive integration, of course. You'll find documentation of cvode's condition_order method in the Programmer's Reference documentation of cvode.