trouble with e_extracellular and play

Anything that doesn't fit elsewhere.
Post Reply
rgiraud
Posts: 2
Joined: Thu Jul 16, 2020 10:04 am

trouble with e_extracellular and play

Post by rgiraud »

Hello
I'm having issues trying to use extracellular stimulation on a single unmyelinated axon. For the moment I'm computing the electric field along the axon induced by a stimulating electrode delivering a pulse stimulation, in a near future I will use fem simulation to get this field.

I want to set these computed values into e_extracellular. I have seen some examples about this method but when I plot the values set into e_extracellular I have something incoherent (all values are zeros).
Is there something I missed with the play method?
I use the following code

Code: Select all

def set_extra_stim(self,xe,ye,ze,sigma,tpulse,ttot):
        #xe,ye,ze are the position of the stimulation electrode
        #tpulse is the high time of the pulse 
        #ttot the total time of the stimulation
        deltat = 0.01 #timestep 
        tdelay=5 #pulse start at 5 ms
        i=0
        v_ext=[]
        v_e=[]
        t_ext = neuron.h.Vector(np.arange(ttot/deltat)*deltat)
        i0=[]
        time=0
        
        #create a single pulse as stimulus
        for k in range(int(ttot/deltat)):
            if time<tdelay or time>tdelay+tpulse:
                time=time+deltat
                i0.append(0)
            else:
                i0.append(-10)
                time=time+deltat
        
        for sec in self.allseclist: #list of all the sections
            for seg in sec:
                V0=[]
                r=math.sqrt((self.x-xe)**2+(self.y-ye)**2+(self.z[i]-ze)**2) #compute distance from the electrode
                for j in range(len(i0)):
                    V0.append(i0[j]/(4*math.pi*sigma*r)) #calculate the potential along the axon

                v_ext.append(h.Vector(V0))
                vstim=h.Vector(V0)
                v_e.append(V0)
                i=i+1
                vstim.play(seg._ref_e_extracellular, t_ext) 

Thank you for your help
Kind regards
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: trouble with e_extracellular and play

Post by ted »

When using an unfamiliar feature, work through a series of toy exercises before addressing the task that is giving you difficulty.

1. Use Vector.play to drive the value of a scalar.
2. Create a single section model with nseg==1 and use Vector.play to drive the segment's v.
3. Insert extracellular into the single section, single segment model and use Vector.play to drive the segment's e_extracellular.
4. Increase the section's length to 500 um, increase nseg to 5, and use Vector.play to drive each segment's e_extracellular.

I bet you discover the answer to your question before you even get to step 4.
rgiraud
Posts: 2
Joined: Thu Jul 16, 2020 10:04 am

Re: trouble with e_extracellular and play

Post by rgiraud »

Hello ted
Thank you very much for you reply, I tried to do your exercises here is the code for the step 4

Code: Select all

import neuron
from neuron import h
import neuron.gui
import matplotlib.pyplot as plt
import numpy as np
h.celsius = 33  # set temperature in celsius
h.tstop = 2e1  # set simulation duration (ms)
h.dt = 0.001  # 0.0005 # set time step (ms)
h.finitialize(-80)
h.v_init=-80

axon1 = h.Section(name = 'axon')
axon1.nseg = 5
axon1.L = 10000
axon1.diam = 1
axon1.insert('hh')
axon1.insert('extracellular')

t = h.Vector(np.linspace(0, 2 * np.pi, 50))
for seg in axon1:
    y = h.Vector(np.sin(t))
    y.play(seg._ref_e_extracellular, t)

vmem = neuron.h.List()
memvrec = neuron.h.Vector()
for seg in axon1:
    memvrec.record(seg._ref_e_extracellular, h.dt)
    vmem.append(memvrec)

h.run()
plt.figure()

for i in range(5):
	plt.plot(vmem[i])
plt.show()
I tried to apply a sinus waveform to each segment of the Axon using play and I recorded the e_extracellular potential to check if the sinus appeared on each segment. That's what I'm observing so I guess my code is working.

However I couldn't see where was mistake in the previous code I sent and while trying to debug it I found something very unusual. My code is now working but I have no idea why. The following code is the exact same as the one I sent before but I added a plot(vstim) at the end of the code and it's now working.

Code: Select all

def set_extra_stim(self,xe,ye,ze,sigma,tpulse,ttot):
        #xe,ye,ze are the position of the stimulation electrode
        #tpulse is the high time of the pulse 
        #ttot the total time of the stimulation
        deltat = 0.01 #timestep 
        tdelay=5 #pulse start at 5 ms
        i=0
        v_ext=[]
        v_e=[]
        t_ext = neuron.h.Vector(np.arange(ttot/deltat)*deltat)
        i0=[]
        V0=[]
        time=0
        
        #create a single pulse as stimulus
        for k in range(int(ttot/deltat)):
            if time<tdelay or time>tdelay+tpulse:
                time=time+deltat
                i0.append(0)
            else:
                i0.append(-10)
                time=time+deltat
        print('io',len(i0))

        
        plt.figure()
        for sec in self.allseclist: #list of all the sections
            for seg in sec:
                V0=[]
                r=math.sqrt((self.x-xe)**2+(self.y-ye)**2+(self.z[i]-ze)**2) #compute distance from the electrode

                for j in range(len(i0)):
                    V0.append(i0[j]/(4*math.pi*sigma*r*10**-6)) #calculate the potential along the axon

                v_ext.append(h.Vector(V0))
                vstim=h.Vector(V0)
                
                v_e.append(V0)
                i=i+1
                vstim.play(seg._ref_e_extracellular, t_ext) 
                plot(vstim)

I tried to apply some delay instead or to plot another variable but only this specific configuration works. If I try to decrement the plot from the for seg loop to the for sec loop only the first segment of each section is computed properly. I have no idea why this is happening do you have an explanation ?
Thank you for your help
Best regards
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: trouble with e_extracellular and play

Post by ted »

Actually the "step 4" code isn't doing what you think it is. Examine the time course of each segment's e_extracellular and you will see that only the last segment's e_extracellular is being driven by the sine wave. Why? The first pass through the loop

Code: Select all

for seg in axon1:
    y = h.Vector(np.sin(t))
    y.play(seg._ref_e_extracellular, t)
creates a new Vector and sets up Vector.play into the first segment's e_extracellular. Good so far. But every pass after that destroys the Vector (and Vector.play) that was built by the previous pass. Why? An object persists only if one or more variables reference it. The statement
y = h.Vector(np.sin(t))
makes y reference a new instance of the h.Vector class. That decreases the reference count for any object that y previously referred to. If that object's reference class drops to 0, poof! no more object!

So only the last pass through the loop does something that persists. How do you verify this? Plot each segment's e_extracellular vs. t and you'll see that only the last one is driven by Vector.play. Besides, I have made this exact mistake myself, so I recognize it. The cure? Append each new Vector to a list before starting the next pass. You could use a Python list, like this:

Code: Select all

yvecs = []
t = h.Vector(np.linspace(0, 2 * np.pi, 50))
for seg in axon1:
    y = h.Vector(np.sin(t))
    y.play(seg._ref_e_extracellular, t)
    yvecs.append(y)
To verify that this works, execute

Code: Select all

for foo in yvecs:
  print(foo)
to confirm that yvecs contains the right type and number of objects, and that they differ from each other.

Or you could use a hoc list (h.List) as you did in the second loop.

Code: Select all

vmem = neuron.h.List()
memvrec = neuron.h.Vector()
for seg in axon1:
    memvrec.record(seg._ref_e_extracellular, h.dt)
    vmem.append(memvrec)
But even this loop has a problem. What does

Code: Select all

for foo in vmem:
  print(foo)
tell you? What's the fix? The statement that creates memvrec
memvrec = neuron.h.Vector()
shouldn't be outside of the loop--it should be the first statement inside the loop.
Post Reply