Implementing extracellular potential from sinusoidal sources

[drunja40]

I am only few weeks into the subject and started with literature review, in particular this paper https://www.sciencedirect.com/science/a ... 1220303720 , because it is very relevant to the topic of Interferential current stimulation and also has source code, which is great to start learning on practice.

The authors model extracellular field as a summation of two point sources injecting current at certain frequencies (2000 Hz and 2005 Hz).
About the implementation they say that:
the potential field for a monopole current source I0 located at (xs, ys, zs)in the anisotropic volume as:

and then
To calculate the overall potential field under quasistatic conditions, we summed the individual potential fields generated by each electrode pair:


As I see they reuse spatial field from Formula1 in Formula2, which I read as "in known homogenous medium each current source generates independent and predictably oscillating voltage potential".
Now I try to look at their code (I only watched few tutorials on NEURON yet) and understand how is this done.

Their full model code is on the bottom of this post, but here is how they compute Formula1 voltages and assign them to each section of the cell:

Code: Select all

    VS1=Voltages(nrn_xyz[:,0],nrn_xyz[:,1],nrn_xyz[:,2],-5e-3,+5e-3,+0e-3,+I,sigmaXX,sigmaYY,sigmaZZ,'V')['V']
    VS2=Voltages(nrn_xyz[:,0],nrn_xyz[:,1],nrn_xyz[:,2],+5e-3,+5e-3,+0e-3,-I,sigmaXX,sigmaYY,sigmaZZ,'V')['V']
    VS3=Voltages(nrn_xyz[:,0],nrn_xyz[:,1],nrn_xyz[:,2],-5e-3,-5e-3,+0e-3,+I,sigmaXX,sigmaYY,sigmaZZ,'V')['V']
    VS4=Voltages(nrn_xyz[:,0],nrn_xyz[:,1],nrn_xyz[:,2],+5e-3,-5e-3,+0e-3,-I,sigmaXX,sigmaYY,sigmaZZ,'V')['V']

    V1 = VS1+VS2  # voltage of first pair
    V2 = VS3+VS4  # voltage of second pair

    recordingDict = dict()
    location = 0.5
    variable = {'v'}
    iSec=0
    # MyPoints = list()
    for Sec in MyCell.get_secs():
        Sec.rx1_xtra=V1[iSec]
        Sec.rx2_xtra=V2[iSec]
        iSec=iSec+1
        if ("node" in Sec.name()):
            for vars in variable:
                rec_dict(recordingDict, Sec, location, vars)

Right after that they assign sinusoidal time-series to some variables, but these are nowhere connected to previously calculated potentials. I need some help in understanding this.

Code: Select all

   
    h.load_file("stdrun.hoc")
    h.init()
    h.finitialize(h.v_init)
    h.fcurrent()

    h.tstop = Duration
    Times=np.arange(0,h.tstop+h.dt,h.dt)
    Sim1=Amps*np.sin(2*np.pi*Freq/1000.0*Times)
    Sim2=Amps*np.sin(2*np.pi*(Freq+DFreq)/1000.0*Times)

    t = h.Vector(Times)
    Stim1=h.Vector(Sim1*1e-3)
    Stim2=h.Vector(Sim2*1e-3)
    Stim1.play(h._ref_is1_xtra, t, 0)
    Stim2.play(h._ref_is2_xtra, t, 0)

Can someone please explain

• How these two simple formulas are implemented in given model.

• Why can't the total voltage be assigned as in

  Code: Select all

  for Sec in MyCell.get_secs():
          Sec.rx1_xtra=V1[iSec]
          Sec.rx2_xtra=V2[iSec]
          iSec=iSec+1
  

• What exactly is achieved by

  Code: Select all

      
      h.tstop = Duration
      Times=np.arange(0,h.tstop+h.dt,h.dt)
      Sim1=Amps*np.sin(2*np.pi*Freq/1000.0*Times)
      Sim2=Amps*np.sin(2*np.pi*(Freq+DFreq)/1000.0*Times)
  
      t = h.Vector(Times)
      Stim1=h.Vector(Sim1*1e-3)
      Stim2=h.Vector(Sim2*1e-3)
      Stim1.play(h._ref_is1_xtra, t, 0)
      Stim2.play(h._ref_is2_xtra, t, 0)
      

I understand that some logic of the model is hidden in .hoc files (which, it seems, authors also reused from someone). Is this the missing piece then and in order to grasp what is happening I also need to read into several .hoc files?

Sorry for pretty beginners questions, I was just trying to get fast-tracked to NEURON via reverse-engineering a well explained source-code, but I guess I just need to spend more time reading this forum, the documentation and watching youtube hands-on sessions.


Full model code:

Code: Select all

def model(Duration,Freq,DFreq,Amps,YY):
    fs = 1000*50*4  # Hz, sampling rate

    h.cvode_active(0)  # turn off variable time step
    h.steps_per_ms = fs*1e-3  # steps per ms for NEURON simulation
    h.dt = 1/(fs*1e-3)  # ms, sample spacing
    h.celsius = h.celsius  # temperature at which to run simulation
    h.v_init = -80.  # mV, initial voltage to begin simulations
    ZZ=0
    MyTraj = np.array([[-50e-3,YY,ZZ], [+51e-3,YY,ZZ]])
    
    MyPath = os.getcwd() + '/'
    CELL_DIR = MyPath+'cells/'
    CELL_FILE_NAME = 'NewMRGaxon_for_Python.hoc'
    MyD = 8.7
    MyCell = MRG(axon_trajectory=MyTraj, fiberD=MyD,
                CELL_FILE_NAME=CELL_FILE_NAME, CELL_DIR=CELL_DIR)

    I = 1.  # current source amplitude in A

    # find xyz coordinates of center of each axon compartment
    nrn_xyz = list()
    for sec in MyCell.get_secs():
        nrn_xyz.append(MyCell.center_3Dcoords(sec))

    nrn_xyz = np.array(nrn_xyz)*1e-6  # m, convert from um to m

    sigmaXX=0.6
    sigmaYY=0.083
    sigmaZZ=0.083

    VS1=Voltages(nrn_xyz[:,0],nrn_xyz[:,1],nrn_xyz[:,2],-5e-3,+5e-3,+0e-3,+I,sigmaXX,sigmaYY,sigmaZZ,'V')['V']
    VS2=Voltages(nrn_xyz[:,0],nrn_xyz[:,1],nrn_xyz[:,2],+5e-3,+5e-3,+0e-3,-I,sigmaXX,sigmaYY,sigmaZZ,'V')['V']
    VS3=Voltages(nrn_xyz[:,0],nrn_xyz[:,1],nrn_xyz[:,2],-5e-3,-5e-3,+0e-3,+I,sigmaXX,sigmaYY,sigmaZZ,'V')['V']
    VS4=Voltages(nrn_xyz[:,0],nrn_xyz[:,1],nrn_xyz[:,2],+5e-3,-5e-3,+0e-3,-I,sigmaXX,sigmaYY,sigmaZZ,'V')['V']

    V1 = VS1+VS2  # voltage of first pair
    V2 = VS3+VS4  # voltage of second pair

    recordingDict = dict()
    location = 0.5
    variable = {'v'}
    iSec=0
    # MyPoints = list()
    for Sec in MyCell.get_secs():
        Sec.rx1_xtra=V1[iSec]
        Sec.rx2_xtra=V2[iSec]
        iSec=iSec+1
        if ("node" in Sec.name()):
            for vars in variable:
                rec_dict(recordingDict, Sec, location, vars)

    MyCell.record(recordingDict)

    h.load_file("stdrun.hoc")
    h.init()
    h.finitialize(h.v_init)
    h.fcurrent()

    h.tstop = Duration
    Times=np.arange(0,h.tstop+h.dt,h.dt)
    Sim1=Amps*np.sin(2*np.pi*Freq/1000.0*Times)
    Sim2=Amps*np.sin(2*np.pi*(Freq+DFreq)/1000.0*Times)

    t = h.Vector(Times)
    Stim1=h.Vector(Sim1*1e-3)
    Stim2=h.Vector(Sim2*1e-3)
    Stim1.play(h._ref_is1_xtra, t, 0)
    Stim2.play(h._ref_is2_xtra, t, 0)

    h.run()
    return MyCell
    

[ted]

Name of the file that contains "Full model code:" ?

[drunja40]

Name of the file that contains "Full model code:" ?

The file name is AxonActivity.py
download link https://github.com/mirzakhalili/Mirzakh ... stems-2020

[ted]

they assign sinusoidal time-series to some variables, but these are nowhere connected to previously calculated potentials

Presumably you want to understand how extracellular stimulation is implemented.

First read the Programmer's Reference documentation of extracellular. Note that (1) extracellular adds additional layers to the cable equation, (2) with appropriate values for the xc, xg, and xraxial parameters, a segment's vext will be essentially equal to the value of that segment's e_extracellular.

Next read the Programmer's Reference documentation of the Vector class's play() method to see that a precalculated waveform can be used as a forcing function to drive the value of a variable.

In principle, one could

Code: Select all

for each segment
  create a corresponding Vector
  fill that Vector with the desired time course of that segment's e_extracellular
  use Vector play to make that Vector drive that segment's e_extracellular

But if the conductive medium is linear and nondispersive, one can take advantage of the fact that extracellular potential at any point xyz is the weighted sum of the stimulus currents, where the weight for any stimulus electrode j is the transfer resistance between that electrode and xyz. This is implemented in cells/xtra.mod, which you will want to read. The "SUFFIX xtra" declaration in its NEURON block means that xtra is a "density" or "distributed" mechanism, so inserting it into a section creates a separate instance of xtra in each of that section's segments. My original version of xtra.mod assumed only one stimulus electrode, but the authors have added a second stimulus eletrode and made corresponding revisions to xtra.mod. e_extracellular for any particular segment will be the sum
rx1*is1 + rx2*is2
(omitting scale factors)
where is1 and is2 are the currents delivered by the two stimulating electrodes, and rx1 and rx2 are the corresponding transfer resistances between those electrodes and the segment. rx1 and rx2 are RANGE variables (can differ from segment to segment) because the coupling between each segment and the stimulus electrodes depends on the geometry of the electrodes and the position of the segment relative to them. is1 and is2 however are GLOBAL variables (same values in all segments) because they are independent of segment location.

The Vector play() statements force the values of is1 and is2 to follow the desired time course during a simulation.

"But I don't see e_extracellular in this mod file."

Right. You see the POINTER variable ex. That is linked to the segment's e_extracellular by setpointer statements in setpointers.hoc (one setpointer statement per segment, since each segment has its own e_extracellular). This makes the NMODL variable ex reference the same memory location used by the segment's e_extracellular. Of course that could be done via a Python statement, but the hoc code worked so the authors decided not to bother rewriting it (why waste time?).

One thing I don't understand is why the authors reverted to assigning values to ex and er in a BREAKPOINT block instead of using BEFORE BREAKPOINT and AFTER SOLVE. Their version of xtra makes ex and er lag behind the actual stimulus current at any given time. And of course PROCEDURE f() is unused (and therefore useless). You should simply remove PROCEDURE f(), and replace the BREAKPOINT block with

BEFORE BREAKPOINT { : before each cy' = f(y,t) setup
ex = is1*rx1*(1e6)+is2*rx2*(1e6)
}
AFTER SOLVE { : after each solution step
er = (10)*rx1*im*area+(10)*rx2*im*area
}

or perhaps the more readable

BEFORE BREAKPOINT { : before each cy' = f(y,t) setup
ex = (1e6)*(is1*rx1 + is2*rx2)
}
AFTER SOLVE { : after each solution step
er = (10)*area*(rx1*im + rx2*im)
}

[drunja40]

thanks a lot for such an extensive answer with links. It will take time to go step by step through every statement, but I already can follow without much details and, at least, trace how other files are used during the python script execution.

[drunja40]

Thanks again for your reply, after diving a bit more into the code and your answer I have better understanding.
I also made that script run on AzureML cloud, which speeds up the simulations by a lot.
Now I can more effectively experiment with adjusting parameters and altering the code.

Primarily I am interested in:

1) Customizing the extracellular field. I plan to create a COMSOL model (non-homogeneous, but still relatively simple) to generate time-dependent extracellular field in absence of neurons and import it into NEURON afterwards. Therefore V-field won't be analytically computable any more.
Can I then apply the strategy you advised?

for each segment
create a corresponding Vector
fill that Vector with the desired time course of that segment's e_extracellular
use Vector play to make that Vector drive that segment's e_extracellular

Then the underlying time-vector will be limited by COMSOL simulation duration, but that should be ok.
Are there any other things to consider here? or how would you do that?

2) Another hick-up I see is that unlike in the current version of code, I won't have "infinite" space any more,rather I will need to map electrode locations in COMSOL to simulated extracellular in NEURON, so that the two simulations' geometries are aligned. Can't see yet, how to do that efficiently, but there are few papers claiming that and that sounds like a technical problem.

3) Another future adjustment will be the number of stimulating electrodes in extracellular (COMSOL). It might not be bipolar, but multipolar.
Do you suggest gradually making changes to:
- the python code AxonActivity.py to increase the number of Vectors

Code: Select all

 
 for Sec in MyCell.get_secs():
        Sec.rx1_xtra = V1[iSec]
        Sec.rx2_xtra = V2[iSec]

- xtra.mod file to increase the number of 'er_xtra'
- etc. until the program does not crash and gives meaningful results?

Or import COMSOL-generated extracellular field as a single number per space\time, without considering every electrode-generated field separately? Then I would only need to revert back to original xtra.mod to one source, right?