the NMODL reference book

Where does one find that?

Here's my take on implementing extracellular stimulation, and using at_time in NMODL.

1. In extracellular stimulation, the potential adjacent to the extracellular surface of the cell is some function f of space and time. The most direct way to implement a computational representation of extracellular stimulation is to control the time course of the potential at the external surface of the model cell.

The alternative is the "activating function" approach, in which the stimulus is represented by injecting current into the model cell. There was a time when the activating function approach might have had a slight computational advantage over the direct approach, but that was decades ago when computers were slow and had limited memory.

Few people who use the activating function approach seem to be aware of the assumptions and limitations of that approach, which is probably understandable because its advocates don't seem to write about them. Here are a couple of things to think about:

How do those assumptions hold up near the ends of a fiber of finite length?

What does the activating function approach predict would happen to a fiber in a field whose second spatial derivative is 0? (and it's easy to create such a field: just apply a voltage across a parallel pair of plates in a conductive medium)

It's so easy to eliminate such concerns. Just use the direct approach to extracellular stimulation, as implemented in xtra.

2. It is easiest to implement the direct approach if one separates f into its spatial and temporal components i.e. f(x, t) = g(x) h(t) where x is understood to be a 1D, 2D, or 3D vector. [One might ask "What if f's time course is also a function of position, as would happen if there are electrodes at two or more locations, and the stimulus current delivered by at least one electrode has a different waveform than what is applied via the others?" This situation is easy to deal with, but first let's focus on the simpler case where the normalized waveform of the current delivered by every electrode is identical to the normalized waveform delivered by every other electrode.]

xtra takes care of the variation in space. It's a "density" (aka "distributed") mechanism, so it will exist in all segments of all sections into which it has been inserted. But its "stimulus" variable is is GLOBAL, so is only needs to be driven in one segment, and all segments will get the same waveform.

What remains is generating the stimulus time course. at_time is deprecated because it is difficult to work with. Deprecated means "don't use this in new code." The way to go is either Vector.play, or to combine algebraic functions and state machines that use events; the latter can be implemented in hoc, Python, or (if absolutely necessary) NMODL (but not using at_time).

And if you use adaptive integration, always check to make sure that the integrator doesn't miss parameter discontinuities--be sure to read the Python or hoc documentation for Vector.play.

3. What if f's time course is a function of position? For example, suppose there are two stimulating electrodes a and b, placed at different locations xa and xb, and each delivers a current with a different normalized waveform e.g. ia and ib (which might be sine waves with slightly different phase or frequency). To represent this, copy xtra.mod to xtra2.mod, change its SUFFIX to xtra2, and revise the code to implement the following substitutions:

rx with rxa, rxb

x, y, z with xa, ya, za and xb, yb, zb

is with isa and isb

In the INITIAL and BEFORE BREAKPOINT blocks

ex is then assigned the value (isa*rxa + isb*rxb)*(1e6)

and at the time of model setup be sure that your hoc or Python code calculates and assigns the appropriate values for rxa and rxb.