Re: Extracellular stimulation and recording

[darkleo]

(In response to
Extracellular stimulation and recording
viewtopic.php?t=168)

Dear colleagues,
I'm working on a direct current model of transcutaneous electrical stimulation of the spinal cord. The volume conductor model was thought to be purely resistive, with multiple layers of tissues (skin, muscles, bones, etc.). COMSOL was used to simulate the electric field induced in the spinal cord. In this way, I hope to simulate how the field affects the excitability of motoneurons (particularly somata and dendrites). My question is, how can I apply DC stimulation using the "extracellular" and "xtra" mechanisms? Will the distribution of the electrical potential calculated with the FEM over the conductive volume (gray matter) in this case be the input for the model's extracellular electrical potential? Or should I rely on the electric field? Should I calculate rx using V (Calculated Electric Potential) / I (input current), which is defined as "transfer resistance between the stimulus electrode and the local node"?

[ted]

direct current model of transcutaneous electrical stimulation of the spinal cord. The volume conductor model was thought to be purely resistive, with multiple layers of tissues (skin, muscles, bones, etc.). COMSOL was used to simulate the electric field induced in the spinal cord.
. . . how can I apply DC stimulation using the "extracellular" and "xtra" mechanisms? Will the distribution of the electrical potential calculated with the FEM over the conductive volume (gray matter) in this case be the input for the model's extracellular electrical potential?

Yes.

Should I calculate rx using V (Calculated Electric Potential) / I (input current), which is defined as "transfer resistance between the stimulus electrode and the local node"?

You could, if that's convenient. In this case, you would have told COMSOL the shape and location of the stimulating electrode, and the time course of the applied current, and used COMSOL to calculate the time course of extracellular potential in the conductive medium, then interpolated the extracellular potentials at the locations that correspond to the centers of your model's segments). Then, the transfer resistance between any given location (x,y,z) and the stimulating electrode would be rx(x, y, z) = v(x, y, z, t)/i(t). You want v and i to be in units such that rx will be in megohms, e.g. mV and nA, or V and uA. If the tissue is nondispersive (i.e. purely resistive), then the values of rx can be calculated using a very simple, brief waveform--for example, a rectangular pulse starting at t = 1 mS that lasts 1 mS.

You might find it useful to read this post viewtopic.php?p=20130#p20130.

You might also want to look at this article
Reilly JP. Survey of numerical electrostimulation models. Phys Med Biol. 2016 Jun 21;61(12):4346-63. doi: 10.1088/0031-9155/61/12/4346. Epub 2016 May 25. PMID: 27223870.
Let me know if you have difficulty getting a copy of this article.
Source code for the NEURON model is available from https://modeldb.science/239006
This illustrates use of xstim (similar to xtra but just a bit easier to use because it omits the stuff that calculates local field potential generated by neuronal activity).
Also shows how to use fzap.mod and fsquare.mod as function generators to drive extracellular potential.

--Ted

[darkleo]

Thank you Ted!

I implemented the program in Python using the suggestions you provided. In this version, I evaluated the electric potential, exported from the data obtained in COMSOL (V). Now, I would like to check the magnitude and the individual components of the Electric Field (Ex, Ey, and Ez) in the morphologically realistic motoneuron (Soma, Dendrites, AHill, and IS). In this case, is it enough to consider the distance between the nodes? Something like what was done in "stimsphere"?

Code: Select all

def _discretize(self):
    for sec in self.all:
        nseg = int( (sec.L/(0.1*self.lambda_f(sec, 100))+0.9)/2 )*2 + 1            
        sec.nseg = nseg

def set_extracellular(self):
    for sec in self.all:
        sec.insert('extracellular')
        
        # Default values
        sec.xg[0] = self._xg0 #1e9  #Internodes: 1e-9 (Reilly 2016)
        sec.xg[1] = self._xg1 #1e9  
        sec.xc[0] = self._xc0 #0.0
        sec.xc[1] = self._xc1 #0.0
        sec.xraxial[0] = self._xraxial0 #1e9
        sec.xraxial[1] = self._xraxial1 #1e9

        sec.insert('xstim')
        # link extracellular and xstim
        for seg in sec:
            seg.xstim._ref_ex = seg._ref_e_extracellular
            seg.x_xstim, seg.y_xstim, seg.z_xstim = self.midpoint(sec, seg.x)
                
def fetch_extracellular_field(self):
    #Load from Matlab-Comsol/File (...)
    for sec in self.all:
        for seg in sec:
            self.e_extracellular_records.append(h.Vector().record(seg._ref_e_extracellular))
            seg.xstim.rx = self.data_fem[idx_count][3]/self.current_input_fem[0] #(mV/nA)

def extracellular_stim(self, delay=100, dur=50, amp = 1):
        self.square_stimulation = h.Fsquare(0.5)
        self.square_stimulation._ref_x = h._ref_is_xstim
        self.square_stimulation.amp1 = amp
        self.square_stimulation.amp2 = 0
        self.square_stimulation.delay = delay
        self.square_stimulation.dp = dur
        self.square_stimulation.num = 1

[ted]

The electric field at any location p with coordinates (x,y,z) is the gradient of electric potential v at that location. You need to discover the locations of the internal nodes (segment centers) of all sections in your model cell, and then use that information to get COMSOL to tell you what the field is at those points. Or you could export a huge table (actually 3 tables) that contains the field sampled at regular intervals in a volume that contains your model cell, and then interpolate the field at each segment center from that table.

Here are a couple of comments about the code examples you included:

Don't ever assign values to any range variable at the 0 and 1 nodes. This applies not only to cm, densities of ion channels or pump mechanisms, compartmental volumes--it also applies to xg, xc, xraxial, and e_extracellular.

"Why?"

Because when you access a range variable other than membrane potential at the 0 or 1 end of a section, you're actually accessing the value of that range variable at the adjacent internal node of that section. Suppose you want to assign values to range variable foo as a function of location, and you do that by iterating over all nodes in segment bah, including the nodes at 0 and 1. Your assignment to foo at the node at 0 actually assigns an incorrect value to foo at the first internal node (which is the value of foo in the first segment). No big deal, because iterating to the first internal node will overwrite your mistake with the correct value for that segment. When you get to the last segment in bah, you assign foo the correct value for that segment, but what happens when you get to the node at 1? Yes, your iteration assigns a value to foo at 1, but that actually overwrites the last segment's correct value with a value that is incorrect. By the way, there is a very widely re-used model implemented in hoc, and available from ModelDB, which actually does this. The error probably isn't fatal for any qualitative result of that particular modeling project, so I won't name any names right here in public (and by now the authors--and many others--know who they are) but just be careful.

"So, if I just iterate over segments, there's no problem."

Right. But the problem could also occur in any code that simply assigns values to range variables at 0 and 1--as you in your sample code.

Bottom line is use the nodes at 0 and 1 for only these two purposes:
as a place to attach a point process or another section
or
as a location where membrane potential may be of interest--and membrane potential at the 0 and 1 ends are calculated not by numerical integration, but algebraically as the weighted sum of the potentials at the centers of all immediately adjacent segments e.g. the potential at some node vj which is located at the 0 or 1 end of a section is

Code: Select all

vj = { SUMMA (vk / rjk) } / SUMMA (1/rjk)
         k                    k

where each vk is the membrane potential at the middle of an adjacent segment, and rjk is the longitudinal resistance between node j and node k.

[darkleo]

Again..., Thank you!

Time to review chapter 5.

I adjusted it to something like (and for all range variables):

Code: Select all

def set_extracellular(self):
        idx_count = 0
        for sec in self.all:
            sec.insert('extracellular')
            sec.insert('xstim')

            for node in sec.allseg():
                if 0 < node.x < 1:
                    node.xg[0] = self._xg0
                    node.xg[1] = self._xg1
                    node.xc[0] = self._xc0
                    node.xc[1] = self._xc1
                    node.xraxial[0] = self._xraxial0
                    node.xraxial[1] = self._xraxial1

                    # link extracellular and xstim
                    node.xstim._ref_ex = node._ref_e_extracellular
                    node.x_xstim, node.y_xstim, node.z_xstim = self.midpoint(sec, node.x)
                    node.Ex_xstim, node.Ey_xstim, node.Ez_xstim= self.data_fem[idx_count][5:8]
                    idx_count += 1