Implementing extracellular stimulation

[Jgod]

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

my step of the project is to have that extracellular stimulation as a grid superposed on the model

Not sure what you mean by this. Do you mean that you know the values of extracellular potential at a set of points in space (what you call a "grid")? If so, your first task is to use that information to calculate the extracelular potential at locations that correspond to the anatomy of the cell you are stimulating. If you prefer to use the "line source approximation" approach, those locations would be straight line segments that either (1) approximate the trajectory of the centroids of the segments of your discretized model cell or (2) are simply the straight lines that run from the 0 to the 1 ends of the segments of your discretized model cell. If you prefer the "point source approximation", those locations could either be the centers of the straight lines that run from the 0 to 1 ends of the segments of your discretized model cell. or the midpoints of the 3D centroid paths of the segments.

I would like to then go in and change a single compartment of the grid to a specific value so that I can see the behavior it would have on the model.

"change a single compartment of the grid to a specific value"? This is not doable in a conductive medium (and it would violate Maxwell's equations). "Ye canna change the laws of physics." Perhaps you mean "change the extracellular potential for a single compartment of a model cell"--which is conceivably doable (and indeed membrane biophysicists have done this at single nodes of Ranvier).

Let me know which of these you have in mind, and then our discussion can become more specific.

[ted]

Yow, instead of quoting your message and replying to it, in a hurry I made the mistake of editing your message. At least I captured your questions so we can continue the discussion.

[Jgod]

I think by grid I meant that I have a matrix of voltage values that I already know from a specific simulation that I run. The values come out of multiple electrodes and seeing what is the minimum voltage for activation is in that specific area therefore with all the possible outputs forming some sort of grid. later on in the project i want to superpose that grid over M1 layer iv PCs. But right now just trying to implement that grid as a whole. I guess changing the specific value isn't needed as I know from the get go all my values. So I am assuming I should proceed with that first part of your reply?

[ted]

I am assuming I should proceed with that first part of your reply?

Please do.

[Jgod]

Hi,

I am working on a code where I extract small compartments of a model cell and put in an electrode (could be an IClamp? unsure about it) in each of the compartment. i am struggling with the implementing aspect of it (doing so with an arbitary value to get syntax before putting realistic values in). Here is my code:

Code: Select all

import os
import numpy as n
from neuron import h, gui  # use this line in linux to open up the GUI

# Load cell model
# Specify the directory path
directory_path = "/Users/jeremigodbout/branchedcells/"

# Load your cell model from the specified directory
file_name = "bcell_comp.hoc"
full_path = os.path.join(directory_path, file_name)
h.load_file(full_path)

# Access the soma section
soma = h.soma
# Access the axon section
axon = h.axon
# Access the apical section
apicals = h.ap
# Access the basal section
basal = h.bas

# Create a list to store midpoints for each section
midpoints_soma = []
midpoints_axon = []
midpoints_apicals0 = []
midpoints_apicals1 = []
midpoints_apicals2 = []
midpoints_basal = []

# Use this print statement to check if length is pulled correctly
print(axon.L)


# function to select segment resolution

# max length of compartment for HH = 50*sqrt(radius)
def nseg(section):
    nseg_value = (int(section.L/(n.sqrt(section.diam/2))))
            # int((section / (0.1 * h.lambda_f(100)) + .999) / 2) * 2 + 1)  # from neuron?
    return nseg_value


def calculate_midpoints(section, midpoints_list):
    nseg_value = nseg(section)
    for i in range(nseg_value):
        segments_l = (section.L / nseg_value)
        midpoints_list.append(segments_l * i + (segments_l / 2))
    return midpoints_list


# Calculate midpoints for each section
calculate_midpoints(soma, midpoints_soma)
calculate_midpoints(axon, midpoints_axon)
calculate_midpoints(apicals[0], midpoints_apicals0)
calculate_midpoints(apicals[1], midpoints_apicals1)
calculate_midpoints(apicals[2], midpoints_apicals2)
calculate_midpoints(basal, midpoints_basal)

# Print or use the midpoints for each section as needed
print("Soma Midpoints:", midpoints_soma)
print("Axon Midpoints:", midpoints_axon)
print("Apical Midpoints:", midpoints_apicals0)
print("Basal Midpoints:", midpoints_basal)

# Insert extra cell voltage capacities into the branches
h.nlayer_extracellular(1) # number of layers
axon.insert(h.extracellular)
apicals[0].insert(h.extracellular)
apicals[1].insert(h.extracellular)
apicals[2].insert(h.extracellular)
basal.insert(h.extracellular)
soma.insert(h.extracellular)

# Insert extracellular values into the cells
# should design field depending on values but choose hard ones for now
# set electrode positions with midpoints arrays
electrodes = []
def electrode_insert(section,midpoints_list):
    electrode = h.electrode(midpoints_list, sec=section)
    electrode.position(midpoints_list, sec=section)
    electrode.voltage = 1.0  # Set your desired extracellular voltage for each electrode
    electrodes.append(electrode)
    return electrode

# Insert electrodes
electrode_insert(soma, midpoints_soma)
electrode_insert(axon, midpoints_axon)
electrode_insert(apicals[0], midpoints_apicals0)
electrode_insert(apicals[1], midpoints_apicals1)
electrode_insert(apicals[2], midpoints_apicals2)
electrode_insert(basal, midpoints_basal)

It fails to properly give the compartments and electrode behavior. Any help is appreciated!
Thank you for your time,
Jérémi

[ted]

I am working on a code where I extract small compartments of a model cell and put in an electrode (could be an IClamp? unsure about it) in each of the compartment.

What do you mean by "put in an electrode"? Are you trying to inject current into the model cell, or are you trying to apply an extracellular stimulus?

get syntax before putting realistic values in

It is a good idea to make sure that syntax is correct. Two suggestions about verifying syntax. First, keep your "syntax testing code" as simple as possible. For example, if the test must involve a model cell, use only a single section with nseg = 1, unless the syntax in question absolutely requires nseg > 1 and/or more than one section. Second, read and re-read the relevant parts of the documentation.

Code: Select all

# Create a list to store midpoints for each section

Why? Given a section, it is very easy to discover the range values (normalized locations) of its segments' centers, e.g. for seg in soma: print(seg.x)

Code: Select all

# max length of compartment for HH = 50*sqrt(radius)

Why? Can you cite a reference that used this spatial discretization strategy?

Code: Select all

# Insert extra cell voltage capacities into the branches
h.nlayer_extracellular(1) # number of layers

Why 1?

Code: Select all

axon.insert(h.extracellular)
apicals[0].insert(h.extracellular)
. . .

Why not iterate over all the sections with one statement? for sec in soma.wholetree(): sec.insert('extracellular')

What do you think electrode_insert() does, why do you want to do that, and how can you discover if it is working properly?

Code: Select all

    electrode = h.electrode(midpoints_list, sec=section)

What is h.electrode?

[Jgod]

Apologies for the double post, tried to edit and messed it up.

Thank you for your response. By "put in an electrode", I meant apply an extracellular stimulus. The reason for trying to get a different segment division of the model due to the amount of voltage values I would be inputing into the cell, therefore it varies over very short distances. The spacial discretization strategy was discussed to be a rule of thumb for HH models by some teachers in a previous class I took although now that you mention it I cannot find much resources on such a thing (would need to look deeper/ask for the sources from that teacher). It did give me very small segments which I thought could be useful for my multiple voltage changes? For the layers being 1 that was a typo from my old code my layers is set as 2 (as I believe it is the default for extracellular stimulus in the documentation).

For the electrode_insert, I believe I was confusing its role with what sec.insert('extracellular') did. I was in a sense adding another "behavior" of my stimuli where it was already in place. And h.electrode was a result of my confusion.

Thank you for your time

[ted]

Apologies for the double post, tried to edit and messed it up.

OK, double post is now gone. BTW I'm an expert on messing up edits of posts.

By "put in an electrode", I meant apply an extracellular stimulus.

Got it.

The reason for trying to get a different segment division of the model due to the amount of voltage values I would be inputing into the cell, therefore it varies over very short distances.

What varies over short distances? Extracellular potential? How does creating a list of segment midpoints help with that? First, the range values that correspond to segment midpoints are available as a simple function call at any time--no need to calculate and store those values in advance. Second, to set up extracellular stimulation, you'll need the xyz coordinates of segment midpoints.

The spacial discretization strategy was discussed to be a rule of thumb for HH models

Spatial discretization is necessary for numerical solution of the cable equation, regardless of whether the cable in question is an axon, or a dendrite, or an electrical cable. "some teachers in a previous class I took"--well, ask for your money back. The d_lambda rule for spatial discretization has a sound, objective basis, verified by tests on many models.

For the layers being 1 that was a typo from my old code my layers is set as 2 (as I believe it is the default for extracellular stimulus in the documentation).

2 layers is the default for NEURON's extracellular "mechanism," regardless of whether or not extracelluar stimulation or recording is used.

For the electrode_insert . . .

OK no need for h.electrode.

With all that resolved, it's time to focus on to your actual task. Is your model cell's anatomy based on detailed morphometry (each section's shape is specified by pt3d statements), or is it a "stylized" model (where each section's geometry is specified in terms of diameter and length)? And what can you say about the extracellular stimulus itself--how are you calculating the extracellular potential as a function of xyz location in the volume conductor?

[Jgod]

I talked to my PI and the lambda rule is what we'll progress forward with. Makes more sense to use something already widely accepted.

the range values that correspond to segment midpoints are available as a simple function call at any time--no need to calculate and store those values in advance. Second, to set up extracellular stimulation, you'll need the xyz coordinates of segment midpoints.

My extracellular voltage would vary over distance and I was trying to create that list of midpoints so that I could then assign each midpoints their specific extracellular voltage. How would I go ahead and do that? My model is a stylized model and the extracellular is field that we created with a pipeline where we take multiple electrodes (from patient images) in the brain that when meshed together outputs a grid mesh of voltage activation values, which is what i want to use in a future cell model. Therefore, I have multiple varying values of voltage of a 3D mesh that I want to use for a specific simulation.

[ted]

Your "3D mesh" defines a "stimulated volume" in space. Your first task is to decide the location and orientation of your model cell in that volume. That's something that you and your PI must do.

Next, for each segment in your model, you have to do two things: (1) determine the coordinates of its center, and (2) from that determine the coordinates of the corresponding point in the "stimulated volume." Part (1) is simple--think of each section's x3d, y3d, and z3d data as the breakpoints of piecewise linear functions of arc length ("range", i.e. normalized path distance from its 0 end), and sample those functions at range values that correspond to segment centers. Part (2), the mapping between model coordinates and stimulated volume coordinates, is much easier than it sounds (it's just a translation and a couple of rotations, all doable with matrix algebra).

Next, for each segment, you'll have to determine the waveform and amplitude of extracellular potential at its center. This task is simplified by the following fact:
If your conductive medium is purely resistive and linear, and it contains two or more electrodes (one or more of these could be a ground) that are stationary in space, and the current delivered by each electrode is the product
aj f(t)
where aj is an electrode-specific amplitude and f(t) is a normalized waveform that is the same over all electrodes,
then
the extracellular potential v(t) at any point x,y,z in the conductive medium will be the product
V(x,y,z) f(t)
where V(x,y,z) is a scalar that depends only on location, and f(t) is the normalized stimulus waveform.

This fact means that you need to do two things: determine the normalized stimulus waveform, and determine the value of the scalar V at the center of each segment.

It's easy to determine the normalized stimulus waveform. Just normalize the waveform of extracellular potential at any of the points in your 3D mesh. As a sanity check, and to verify that your extracellular stimulus was calculated properly in the first place, you should probably compare normalized waveforms from several (dozens or hundreds of) different locations in your 3D mesh. A significant difference between them would be hard to reconcile with "linearity" and "purely resistive".

If you already have stimulus waveforms calculated over a 3D grid, then you have what you need to get the values of f(t) and the values of V(x,y,z) (for the latter, just find the time tx at which |f(t)| is largest, and then
for each point in the grid
V(x,y,z) at that point = amplitude of the stimulus waveform at that x,y,z

Most likely, none of your model's segment centers will correspond directly to any of the points in your 3D mesh. This means you'll have to use interpolation to discover the V(x,y,z) at each of the segment centers. There are many interpolation strategies, so you and your PI will need to choose one. Who knows, maybe the software that calculated your 3D mesh can do that for you.

"OK, suppose I have f(t) and I know the V(x,y,z) for each segment center. What next?"

That's the easiest part. I'll address it in a subsequent post.

[Jgod]

From further discussions and coding, we have concluded that we want (for now) a singular point where the voltage is in each segment as our wave form as a bolus (in the future a time varying voltage will be implemented). We have values that we want to apply to the middle of x,y,z coordinates (i.e. center of segment, center of disk & length). I also know my V(x,y,z) for each segment center. How can I proceed to implement that. I have worked through your previous suggestions and have written a code that follows your suggestions and runs.

"OK, suppose I have f(t) and I know the V(x,y,z) for each segment center. What next?"

That's the easiest part. I'll address it in a subsequent post.

So knowing what I stated above, how could I proceed?

Also thank you so much for all the help its been very insightful and helpful!! :)

[ted]

What do you mean by "a singular point"?
What does "our wave form as a bolus" mean?

[Jgod]

Here is a sketch of what I meant by single point of stimuli in a segment:

https://imgur.com/tU9v8sB

I believe what we meant as a "bolus" is that single point of voltages is on/off for a specific period of time as if it was clamped (VClamp at that middle point) (I think this is the right terminology for what I am trying to achieve?)

[ted]

Bolus is a time honored term for qualitative description of the time course of injecting a drug, but not so for electrical stimulation, which has its own well established vocabulary. For the sake of clarity and mutual understanding, let's stick with the latter: one speaks of pulses or other waveforms. What you want is a brief pulse--and if the amplitude is to remain constant while the pulse is on, it's a brief rectangular pulse.

Assuming that you know the extracellular potential during the stimulus pulse at each location that corresponds to a segment center, you could in principle use the Vector class's play() method to force e_extracellular at each segment center to follow the desired time course: start at 0 when t==0, jump to the appropriate value at the start of the pulse, and fall back to 0 at the end of the pulse. Conceptually straightforward, and easy to implement. In pythonish pseudocode

Code: Select all

create a Vector to hold important stimulus time values
fill that Vector with these values: 0, ton, ton, toff, toff, 1e9
# where ton and toff are the times at which the pulse starts and stops, respectively

create a list that will be used to hold the stimulus vectors
for each section:
  for each segment:
    get the value vx that you want e_extracellular to be for this particular segment
    create a new Vector to hold the local stimulus amplitude values
    fill this "stimulus Vector" with the following values: 0, 0, vx, vx, 0, 0
    make this "stimulus Vector" play its values into this segment's e_extracellular
    append the stimulus Vector to the stimulus vector list

# and now every time you run a simulation, e_extracellular at all those segment centers
# will follow the time course you wanted

By
"make this Vector play its values into this segment's e_extracellular"
I mean make it do "vector play with interpolation". Assuming the stimulus time vector is called stimtvec, the stimulus amplitude vector is called stimampvec, and the segment is called seg, then executing this statement
stimampvec.play(seg._ref_e_extracellular, stimtvec, TRUE)
will make e_extracellular at this location follow the desired time course during a simulation.

Be sure to read the documentation of Vector play at nrn.readthedocs.io

[Jgod]

Hi Ted, back with more questions. I got the stimuli code to work pretty well on a simple model and working to size it up. On top of that, I am trying to find positions of each sections relative to the soma so that I can know their relative (x,y,z) position to a point source with known distance (x,y,z) to the soma. What would be the best way to go about this? Also, the model that I am currently using is only on the x-axis right now. What would be the best approach to making such a model "twist" or "turn" to create the need of (x,y,z) coordinates to represent it (that may not be possible I am unsure).

Thank you for your time,
Jérémi