Resampling high-resolution potentials

Sometimes you may want to sample potentials agnostic to the final fiber compartment coordinates (e.g., model several fibers along the same trajectory with different diameters, ultrastructure, or membrane mechanisms). This tutorial details how to resample these potentials to match the coordinates of a specific fiber.

Note

This tutorial will use a lot of code explained in the simulation tutorial, so it is recommended to review that before proceeding.

Generate high-resolution potentials

This tutorial provides an example of repurposing electrical potentials that were sampled at high spatial resolution. Users may use external softwares to calculate extracellular potentials (e.g., COMSOL, ANSYS). In this example we will use a gaussian distribution with 5 um spacing between coordinates.

Note

The spacing does not have to be uniform, but the distance between consective points must be suffieciently small to not affect your simulation results. Your coordinates must be one-dimensional arc-lengths along the length of the fiber. If your coordinates are three dimensional, you can use a function such as scipy.spatial.distance.euclidean() to calculate the arc-length between each coordinate, or use a 3D fiber path.

from scipy.interpolate import interp1d
from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt

n_coords = 50000

supersampled_potentials = norm.pdf(np.linspace(-1, 1, n_coords), 0, 0.05) * 10
coords = np.cumsum([1] * n_coords)

plt.scatter(coords, supersampled_potentials)
plt.title('Extracellular potentials')
plt.xlabel('Position along fiber (\u03bcm)')
plt.ylabel('Potential (mV)')
plt.show()

Create a fiber

For this tutorial, we will create a model Fiber using the MRG model as in the fiber creation tutorial. Instead of specifying the number of coordinates, we will specify the length of our fiber as the length of our super-sampled fiber coordinates.

from pyfibers import build_fiber, FiberModel

fiber_length = np.amax(coords) - np.amin(coords)

fiber = build_fiber(FiberModel.MRG_INTERPOLATION, diameter=10, length=fiber_length)
print(fiber)

MRG_INTERPOLATION fiber of diameter 10 µm and length 49382.20 µm 
	node count: 45, section count: 485. 
	Fiber is not 3d.

To obtain potential values at the center of each fiber compartment, we must resample our high-resolution “super samples” potentials. We can use the resample_potentials() method of the fiber object to do this.

fiber.potentials = fiber.resample_potentials(supersampled_potentials, coords)

plt.scatter(fiber.coordinates[:, 2], fiber.potentials)
plt.title('Extracellular potentials')
plt.xlabel('Position along fiber (um)')
plt.ylabel('Potential (mV)')
plt.show()

Our potentials are not centered on the fiber. By default, the resampling occurs with the zero point of the supersampled potentials and the fiber aligned. We can center the fiber along the supersampled path by setting center=True. We can also have the potentials replace the current ones automatically with inplace=True.

fiber.resample_potentials(supersampled_potentials, coords, center=True, inplace=True)

plt.scatter(fiber.coordinates[:, 2], fiber.potentials)
plt.title('Extracellular potentials')
plt.xlabel('Position along fiber (um)')
plt.ylabel('Potential (mV)')
plt.show()

Simulation

As before, we will create a monophasic stimulation waveform.

time_step = 0.001  # milliseconds
time_stop = 20  # milliseconds

start, on, off = 0, 0.1, 0.2  # milliseconds
waveform = interp1d(
    [start, on, off, time_stop], [0, 1, 0, 0], kind="previous"
)  # monophasic rectangular pulse

time_steps = np.arange(0, time_stop, time_step)
plt.plot(time_steps, waveform(time_steps))
plt.xlim(0, 1)
plt.title('Stimulation waveform')
plt.xlabel('Time (ms)')
plt.ylabel('Amplitude')
plt.show()

and then run a fixed amplitude simulation:

from pyfibers import ScaledStim

# Create instance of Stimulation class
stimulation = ScaledStim(waveform=waveform, dt=time_step, tstop=time_stop)

ap, time = stimulation.run_sim(-1.5, fiber)
print(f'Number of action potentials detected: {ap}')
print(f'Time of last action potential detection: {time} ms')

Running: -1.5
	N aps: 1, time 0.529
Number of action potentials detected: 1.0
Time of last action potential detection: 0.529 ms