Multi-electrode extracellular stimulation

Anything that doesn't fit elsewhere.
Post Reply
Posts: 3
Joined: Fri Sep 08, 2017 7:58 am

Multi-electrode extracellular stimulation

Post by ChrisK »

Hey NEURON-community,

I am using NEURON quite a while, but until now it has been limited to standard models provided.

Now, I am in the field of neuroprosthetics, where I want to model different stimulation patterns elicited by multiple extracellular stimulation electrodes.

As a start, I used the package "extracellular_stim_and_rec" provided by ted.
I modified the stim.hoc in order to generate a desired biphasic pulse.

I am now wondering: How could I (most easily) implement multiple (instead of only one) extracellular electrodes and how could I add these electrodes to the electrode-GUI provided?

A search did not give me helpful results.

I would be more than pleased, if you could help me with my project!

Thank you very much in advance.

Posts: 3
Joined: Fri Sep 08, 2017 7:58 am

Re: Multi-electrode extracellular stimulation

Post by ChrisK »

Hey there,

can anybody maybe help me with figuring out the aforementioned problem?

Thank you very much in advance!

Site Admin
Posts: 5810
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine

Re: Multi-electrode extracellular stimulation

Post by ted »

It doesn't get any simpler than this.

Step 1:
For each electrode j, solve the Poisson equation to find the closed form function fj(x,y,z)
that specifies the transfer resistance between electrode j
and any point (x,y,z) in the conductive medium.
Step 2:

Code: Select all

for each electrode j
  for each compartment k in a model
    plug the coordinates (xk, yk, zk) of that compartment into fj
    to calculate the transfer resistance between compartment k and electrode j
Step 3:
Revise xtra.mod so that it can handle the time course of each electrode's current ij and the transfer resistance between any given compartment and that electrode.

Clearly this becomes cumbersome if there are more than a very small number of electrodes, and each electrode delivers a current that has a different time course.

Alternatively, you could use a separate program such as COMSOL to calculate the detailed time course of extracellular potential at the center of each compartment k, and then use the Vector class to play those time courses into e_extracellular.

Implementation of either approach is greatly simplified if each electrode current ik can be described as the product ak*g(t) i.e. each electrode's current is the product of a (possibly unique) constant scale factor and a common time course. Even a situation in which there are two or three different time courses isn't too terrible (e.g. electrode pair 1&2 delivers a current with waveform g12(t) and electrode pair 2&3 delivers a current with waveform g23(t)).
Posts: 3
Joined: Fri Sep 08, 2017 7:58 am

Re: Multi-electrode extracellular stimulation

Post by ChrisK »

Hey ted, hey community,

thanks a lot for the help so far.
This really let me do the things I needed to do.

I added a 2nd extracellular electrode and adapted scripts accordingly.

I did some validation tests by comparing the modified version of extracellular_stim_and_rec with the original one.
By comparing values, for example of V(0.5) at the soma and soma_e_extracellular, results were the same.

I switched of one electrode or the other and compared it to the original version.
At a second attempt I superimposed the two electrodes at the same position where the sum of the two amplitudes is the same as the amplitude of the original file.

I only got different result for vrec.
I adjusted field.hoc according to another post in the forum (see code from above).
Unfortunately, results were different here as soon as a stimulus was applied.
The general waveform was the same, the value were higher / lower though.

Does anybody have an idea?

Especially for calcrx.hoc I am not sure whether I calculte rx1_xtra / rx2_xtra in the right way.

- calcrx.hoc

Code: Select all

// assume monopolar stimulation and recording
// electrode coordinates:
// for this test, default location is 50 microns horizontal from the cell's 0,0,0
XE1 = 50  // um
YE1 = 0
ZE1 = 0

XE2 = -50  // um
YE2 = 0
ZE2 = 0

proc setrx() {  // now expects xyc coords as arguments
  forall {
    if (ismembrane("xtra")) {
// avoid nodes at 0 and 1 ends, so as not to override values at internal nodes
      for (x,0) {
        r = sqrt((x_xtra(x) - xe)^2 + (y_xtra(x) - ye)^2 + (z_xtra(x) - ze)^2)
//        r = sqrt((x_xtra(x) - $1)^2 + (y_xtra(x) - $2)^2 + (z_xtra(x) - $3)^2)
        // 0.01 converts rho's cm to um and ohm to megohm
        // if electrode is exactly at a node, r will be 0
        // this would be meaningless since the location would be inside the cell
        // so force r to be at least as big as local radius
        if (r==0) r = diam(x)/2
        if (1 == $4) rx1_xtra(x) = (rho / 4 / PI)*(1/r)*0.01
		if (2 == $4) rx2_xtra(x) = (rho / 4 / PI)*(1/r)*0.01

create sElecA
create sElecB 

proc put_sElecA() {
	sElecA {
		// make it 1 um long
		pt3dadd($1-0.5, $2, $3, 1)
		pt3dadd($1+0.5, $2, $3, 1)

proc put_sElecB() {
	sElecB {
		// make it 1 um long
		pt3dadd($1-0.5, $2, $3, 1)
		pt3dadd($1+0.5, $2, $3, 1)

put_sElecA(XE1, YE1, ZE1)
put_sElecB(XE2, YE2, ZE2)

objref gElec // will be a Shape that shows extracellular electrode location
gElec = new Shape(0)  // create it but don't map it to the screen yet
gElec.view(-245.413, -250, 520.827, 520, 629, 104, 201.6, 201.28)
objref pElecA
objref pElecB  // bogus PointProcess just to show stim location

sElecA pElecA = new PointProcessMark(0.5) // middle of sElec
sElecB pElecB = new PointProcessMark(0.5) // middle of sElec

gElec.point_mark(pElecA, 2) // make it red
gElec.point_mark(pElecB, 3) // make it red

proc setelec() {
	xe = $1
	ye = $2
	ze = $3
	electrode_num = $4
	setrx(xe, ye, ze, electrode_num)
	if(1 == $4)	put_sElecA(xe, ye, ze)
	if(2 == $4)	put_sElecB(xe, ye, ze)

// setrx(50, 0, 0)  // put stim electrode at (x, y, z)

setelec(XE1, YE1, ZE1,1)  // put stim electrode at (x, y, z)
setelec(XE2, YE2, ZE2,2) 

// print "Use setelec(x, y, z) to change location of extracellular recording electrode"

xpanel("Extracellular Electrode Location", 0)
  xlabel("xyz coords in um")
  xvalue("x1", "XE1", 1, "setelec(XE1,YE1,ZE1,1)", 0, 1)
  xvalue("y1", "YE1", 1, "setelec(XE1,YE1,ZE1,1)", 0, 1)
  xvalue("z1", "ZE1", 1, "setelec(XE1,YE1,ZE1,1)", 0, 1)
  xvalue("x2", "XE2", 1, "setelec(XE2,YE2,ZE2,2)", 0, 1)
  xvalue("y2", "YE2", 1, "setelec(XE2,YE2,ZE2,2)", 0, 1)
  xvalue("z2", "ZE2", 1, "setelec(XE2,YE2,ZE2,2)", 0, 1)
- field.hoc

Code: Select all

// $Id: field.hoc,v 1.2 2005/09/10 23:02:15 ted Exp $

//vrec = 0  // extracellularly recorded potential

proc fieldrec() {
	vrec = 0
		forall {
			if (ismembrane("xtra")) {
				// avoid nodes at 0 and 1 ends, which shouldn't make any contribution
				for (x,0) vrec += er1_xtra(x) + er2_xtra(x)

proc init() {

proc advance() {

For the attaching the stimuli, I modified proc attach_stim() (see above).
I am not very sure but the updates of is1_xtra and is2_xtra are correct, right?

- stim.hoc

Code: Select all

proc attach_stim() {
// since is_xtra is GLOBAL, we only need to specify
// for one instance of xtra, i.e. at just one internal node
// of only one section that contains xtra
  forall {  // check each section to find one that has xtra
    //if (ATTACHED__ == 0) {  // don't bother if stim is already attached to something
      if (ismembrane("xtra")) {, stim_timeA, 1) // "interpolated" play, stim_timeB, 1) // "interpolated" play
        //ATTACHED__ = 1

For automated reading and saving of some results, I wanted to create a script and used code I found here in the forum.

Code: Select all

objref vsoma, fok
proc postprocess() {
// proc does not work, but no idea why

  	vsoma = new Vector()
  	fok = new File()
	print soma.v(0.5)

proc go() {

xbutton("Go", "go()")
xpanel(200, 200)
My target file is always empty and I cannot figure out why. I proceeded the same way as others already did in the forum.
Maybe I was coding to long at a time and oversaw some details?

Creating my own simple neuron for further test purposes of my biphasic, adaptable stimulation pulses, I created an axon, a soma and a dendrite in the CellBuilder.
Then I added xtra to proc biophys().
It seems to work fine and hence, it should be correct, right?

Code: Select all

proc biophys() {
  forsec all {
    Ra = 35.4
    cm = 1
    insert pas
      g_pas = 0.001
      e_pas = -70
    insert hh
      gnabar_hh = 0.12
      gkbar_hh = 0.036
      gl_hh = 0.0003
      el_hh = -54.3
    insert extracellular
      xraxial[0] = 1e+009
      xraxial[1] = 1e+009
      xg[0] = 1e+009
      xg[1] = 1e+009
      xc[0] = 0
      xc[1] = 0
      e_extracellular = 0
	insert xtra
	  rx1_xtra = 1
	  rx2_xtra = 1
	  x_xtra = 0
	  y_xtra = 0
	  z_xtra = 0

Again, thank you very much for your help! I read and tried a lot in the last weeks in NEURON, but I just couldn't solve the issues mentioned above.
If you want, I can provide a link to GitHub for my adjusted code for the biphasic, multielectrode stimulation.

Thank you very much and have a nice day!

Post Reply