The model was originally number 3810 from Model DB, but I made changes so that the axons are created in a template and the Stim.hoc code applies the stimulus values in the text file StimOne.txt.
There are numerous possible sources of problems. It is difficult to write templates properly. Your stimulus is expressed as a voltage, but McIntyre's original code used current. All of that could be perfectly fine, and the problem could lie elsewhere--inherent in the original source code but never exposed if stimuli weren't sufficiently large.
Speaking of which, have you tried large stimuli with the original model? I just did, and it produces the same error message when the stimulus current is increased by a factor of 1e3 (i.e. from the original 2 nA to 2000 nA). So there's the first thing to investigate and resolve. Looking at AXNODE.mod I see a lot of strange vtrap functions, and some attempts to avoid numerical
underflow (??). Plotting the voltage dependence of mp_inf, m_inf, h_inf, s_inf, tau_mp, tau_m, tau_h, and tau_s, I see that m_inf has a strange discontinuity in the vicinity of -34 mV, at which the original formula for m_inf would have a 0 in the denominator. There's probably an error in the corresponding vtrap function.
For that matter, I'm not sure there aren't errors in _all_ of the vtrap functions, but tracking them down and eliminating them* is going to be a pain in the neck. But you won't know if any of them are causing your "out of range" error messages until you have eliminated all of them. Also, you'll want to check at every step, by running a simulation with McIntyre's original model, to make sure that you haven't broken anything--you'll want to verify that using your revised mod file produces results that are identical (or at least not significantly different from) those produced when the original mod file was used.
*--The obvious one with m_inf came about through the use of multiple vtrap functions and a proliferation of parameter names. This is a minor development and debugging nightmare. There should really be only _one_ vtrap function, to which you pass parameters as necessary in order to calculate the value of a or b in procedure evaluate_fct().
To help get you started, here's a little background on vtrap and some hints:
HH-style model rate equations often contain expressions that are equivalent to
x / (exp(x/y) - 1)
As x/y -> 0, the numerator and denominator -> 0, which can cause problems.
This can be avoided by noting that, as x/y -> 0,
exp(x/y) - 1 = 1 + x/y + ((x/y)^2)/2 + . . . - 1
approaches
x/y + (x^2)/((y^2)*2)
so
x / (exp(x/y) - 1)
approaches
x / (x/y + (x^2)/((y^2)*2)) = y / (1 + x/(2*y))
which is closely approximated by
y * (1 - x/(2*y))
when x/(2*y) << 1.
In NMODL, this can be implemented as a function
Code: Select all
FUNCTION vtrap(x,y) {
if (fabs(x/y) > 1e-6) {
vtrap = x/(exp(x/y) - 1)
} else {
vtrap = y*(1 - x/y/2) : from Taylor's series
}
}
You can use this function to replace McIntyre's original vtrap, which as far as I can tell is not used anywhere in this mod file.
You'll want to replace each
vtrapN(v) where N = 1, 2, etc.
with something of the form
c*vtrap(x, y) where c is a constant, x is an expression involving v, and y is another constant.
For example, examining the body of function vtrap1 I see that the original formula was
ampA*(v+ampB)/(1 - exp(-(v+ampB)/ampC))
so
a = q10_1*vtrap1(v)
is really
a = q10_1*ampA*(v+ampB)/(1 - exp(-(v+ampB)/ampC))
which means that
c = q10_1*ampA
and you have to get this
(v+ampB)/(1 - exp(-(v+ampB)/ampC))
into the same form as
x/(exp(x/y) - 1)
This is accomplished by the substitutions x = -(v+ampB) and y = ampC.
Then the statement
a = q10_1*vtrap1(v)
becomes
a = q10_1*ampA*vtrap(-(v+ampB), ampC)
Subject of course to verification, which comes by compiling the mod file, running a new simulation, and verifying that the results are the same as with the original mod file.
Repeat this process for vtrap2, 6, 7, and 8, and if you do it right, your tau and inf values will not have strange discontinuities in them. At that point you will be able to reexamine what happens to the original McIntyre model with large stimulus intensities, and also to try the revised mod file with your own hoc code.