Stimulus Activation/Deactivation in NMODL

NMODL and the Channel Builder.
Post Reply
ajpc6d
Posts: 33
Joined: Mon Oct 24, 2022 6:33 pm

Stimulus Activation/Deactivation in NMODL

Post by ajpc6d »

I made the following NMODL file to deliver a transmembrane current injection into a cell. The RANGE variable amp receives a played-in vector from a Python script. Below, I show the original method (which doesn't work quite right) and the second method (which does).
The difference between them is only evident for nanosecond-scale (or GHz frequency) pulses; in essence, the excitation threshold found in NEURON for this time/frequency regime is well below expected values (with method #1).
My question is why is this a necessary step? Does the divergence of the two methods indicate some clumsiness in my variable time-step method (in the 'parent' Python file)? Again, the NMODL file now does what I need, so I'm only curious if there is a general rule of usage I could/should learn from this episode.

Code: Select all

NEURON {
    SUFFIX direct_inject
    NONSPECIFIC_CURRENT i
    RANGE i, amp, delay, pw
}

UNITS {
    (mV) = (millivolt)
    (uA) = (microamp)
    (mA) = (milliamp)
}

PARAMETER {
    :amperage = -2e-2 (uA/cm2)
    
   }

ASSIGNED {
    i (mA/cm2)
    amp (uA/cm2)
    delay (ms)
    pw (ms)

}

INITIAL {
    i = 0}

BREAKPOINT {
    : FIRST, I TRIED THIS LINE. THE FUNCTION DID NOT PERFORM AS EXPECTED. 
    : i = -amp*(1e-3)
    
    : WHEN THE ABOVE FAILED, I IMPLEMENTED THE CODE BELOW
    LOCAL X
    at_time(delay) 		: THESE TWO at_time() LINES WERE TRIED WITH THE FIRST METHOD, BUT DID NOT SUFFICE
    at_time(delay+pw)
    if (t > delay && t < delay+pw) {
        X = 1				: THIS "GATING VARIABLE" APPEARS TO BE THE KEY
    }
    else {
        X = 0
    }
    i = -amp*(1e-3)*X 		: WITH THIS IMPLEMENTATION, THE FUNCTION WORKS AS EXPECTED
}
ted
Site Admin
Posts: 6314
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Stimulus Activation/Deactivation in NMODL

Post by ted »

Well, of course you don't need to write your own NMODL code to implement a stimulus that consists of a single pulse; IClamp will handle that. You're probably looking for a way to deliver a series of stimuli.

Your first approach would work if you were using fixed time step integration with a dt <= stimulus duration. It wouldn't work with adaptive integration (nicknamed "CVODE") because the integrator would very likely choose a time step that completely skips over the brief stimulus.

Your second approach works because it notifies NEURON's adaptive integrators about the start and end of the stimulus. However, at_time is deprecated (one problem with it is the difficulty of implementing anything but the most trivial usage).

NMODL with events is the best way to implement waveforms that are described by algebraic functions (e.g. amplitude or frequency modulated sine wave, sum of two or more functions). Why events? Events provide the most convenient way to control stimulus start, duration, number of repetitions, sudden parameter changes, etc..

How to learn more about using events to generate stimulus waveforms? Or maybe even something that already exists and will do what you want without any new programming on your part? Search the Forum for mentions of pulsedistrib in 2022 and izap in 2023.

And of course Vector.play (no NMODL code at all) is the easiest (if not always most efficient of memory usage) way to implement stimuli with arbitrary waveforms.
ajpc6d
Posts: 33
Joined: Mon Oct 24, 2022 6:33 pm

Re: Stimulus Activation/Deactivation in NMODL

Post by ajpc6d »

From https://nrn.readthedocs.io/en/latest/py ... tml#NetCon, I see
The target must be a POINT_PROCESS or ARTIFICIAL_CELL that defines a NET_RECEIVE procedure.
From the NMODL reference book*, I see
An onset event, generated by the system when the connecting NetCon’s source passed threshold t - delay ago
From these two pieces of information, can I conclude that NET_RECEIVE is incompatible with a SUFFIX? And if so, is there a standardized procedure to use in parallel with Vector.play to guarantee appropriate time steps? I had thought CVode.event would be best, but maybe there's another option I haven't encountered yet.

*https://www.neuron.yale.edu/neuron/stat ... /nmodl.htm
Last edited by ajpc6d on Sat Aug 12, 2023 5:02 pm, edited 1 time in total.
ted
Site Admin
Posts: 6314
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Stimulus Activation/Deactivation in NMODL

Post by ted »

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