Adding extracellular stimulus to network of neurons

[EleBern]

Hello,

I am trying to add an extracellular stimulus to a network of single compartment cells. However, the results running a simulation without any extracellular stimulus and with an unrealistically strong stimulus generate the exact same voltage traces.

I checked the value of e_extracellular and vext[0], and they both are as expected in the range of +100 mV. The membrane potential, v, stays at about -70 mV. In the documentation I read that “When the extracellular mechanism is present, v refers to the membrane potential and vext (i.e. vext[0]) refers to the extracellular potential just next to the membrane. Thus the internal potential is v+vext (but see BUGS).”

I was not able to find any reference to these bugs. My code is in hoc, the neuron version I’m using is 8.1.0. Could you please help figure out this issue?

Thank you,
Kind regards

[ted]

The difference between what you expected to happen, and what actually happened, has nothing to do with software, let alone software bugs. You have rediscovered a consequence of Maxwell's equations and Gauss's law. Suggested reading in this Forum:
viewtopic.php?p=10782#p10782

Does this mean you have to abandon your network of single compartment model neurons? No. You could just attach a current source to each model cell, and make that current clamp deliver a current with appropriate time course and amplitude. What is "appropriate" for any particular cell depends on the cell class (which implies a particular geometry of soma, axon, and dendrites--yes, your models are single compartments, but extracellular stimulation should have a stronger effect on a pyramidal cell than on a granule cell because of their different shapes and sizes), and the presumed anatomical location and orientation of any particular model cell. And you'll want to do this in a way that is plausible and will convince reviewers.

[EleBern]

Of course, thank you for pointing it out Ted!

[ted]

Given a point neuron, extracellular stimulation can be represented by a current source attached to the model cell, but the question then is: how to determine the amplitude and time course of that current? Amplitude is clearly a function of assumed location of the model cell in the conductive medium, but should also take into account the cell's intended shape and orientation relative to the potential gradient in the medium. For example, if a point neuron model is used to represent a biological neuron with extended geometry and "anatomical polarity" (i.e. the spike trigger zone in the biological neuron is closer to one end of the cell than the other), it should be more sensitive to a properly oriented extracellular potential gradient than, say, a point neuron representation of a small, relatively symmetric shape like a small inhibitory neuron. I'd start looking at other network models that involved point neurons to see how their authors dealt with this issue--e.g. models implemented with Brian or NEST.

[EleBern]

Hi Ted,

I've had a look at models and papers and I've found a few examples of how extracellular stimulation is added to single compartment neuron models.

Thank you for your help!

[ted]

Haven't been much help yet--you've done all the work so far.

Implementing stimulation will be easy if you only need a single rectangular pulse--the IClamp class will do for that.

But what if your stimulus current should have a more complex time course, e.g. two or more monophasic rectangular pulses, one or more biphasic pulses possibly separated by "silent" intervals, a swept sine wave ("zap" function)? The most convenient way to implement those is with a custom mod file. Let me know and I can point you in the right direction.

[EleBern]

At the moment I am using a sinuisodal wave as a stimulus, just to test if I'm applying it to the neurons correctly.
The goal is to then use more complex stimuli (likely multiple monophasic rectangular pulses).
I plan to calculate the stimulus current separately from the network simulation and save the current to be injected in a neuron at each time step in a vector.

[ted]

I asked because DBS protocols use biphasic, charge-balanced waveforms.

Vector play works for arbitrary waveforms. If the conductive medium is nondispersive ("purely resistive"), then the time course of extracelluar potential vx(x,y,z) will be of the form A(x,y,z)*f(t) where f is a function of time with maximum amplitude 1 and A is a function of position. In other words, vx has the same "shape" ("waveform") everywhere; the only difference between vx at different locations is its magnitude. This makes it easy to spawn a bunch of vectors that can be used to drive extracellular potential (or the current delivered by an IClamp* to a single compartment model cell), e.g. in pseudocode

create a "normalized waveform" vector f
for each location of interest in your model
make a copy of f
multiply the copy with the appropriate A(x,y,z)
use the scaled copy to drive e_extracellular at (x,y,z) or IClamp.i attached to the cell that is located at (x,y,z)

The Vector class's methods for copying and scaling a vector operate at machine language speed.

*--An aside about the waveform of the current delivered by an IClamp: if the extracellular stimulus is much shorter than membrane time constant (e.g. 1/10th or less), the detailed time course of the current waveform won't matter much, because of the low pass filtering effect of the cell membrane. What will matter is the net charge delivered by the current injection.

A hint about construction of the stimulus waveform:

If the waveform is a sequence of one or more pulses, using Vector play with interpolation (set the "continuous" argument to 1 (hoc) or True (Python)) can greatly simplify your specification of f. Just provide the values of f at t==0 and on both sides of each abrupt change of f. Example:
for a single pulse that starts at 5 ms and lasts 1 ms, the t and corresponding f values are
0, 0
5, 0
5, 1
6, 1
6, 0
7, 0 // this last "point" ensures that f remains 0 forever after t==6 -- read the documentation Vector play
Yes, just 6 value pairs. And if you change dt, you don't have to change these values (assuming that dt is a whole number fraction of 1).

<< NTC revised this after initial posting >>

[EleBern]

Hi Ted,

In the past month I've been playing around with the possibility of switching to a multi-compartment model.
While parameter tuning the new multi-compartment model in my network, I started getting a series of warning messages as below:

Code: Select all

exp(700.58) out of range, returning exp(700)
No more errno warnings during this execution

I checked the mod files and I can't find any obvious discontinuities in any gating variable.
The warning doesn't suggest where the issue may be. Do you have any suggestions on how I may pinpoint the issue?

Thank you very much

[ted]

I've never seen this msg signify a discontinuity. It's usually caused by NMODL code with rate constants that involve exponential functions. Typically the exp is in a denominator where limiting the value returned by exp(x) to exp(700) will do no harm. Nothing will blow up numerically, and your simulations should produce mathematically reasonable results.