Basic simulation and search for activation threshold

This tutorial will use the fiber model created in the fiber creation tutorial to run a basic simulation of fiber stimulation. Then we will demonstrate running a bisection search for the fiber’s activation threshold (i.e., the minimum stimulation amplitude needed to generate an action potential).

Create the fiber

As in the fiber creation tutorial, we use build_fiber() to create a 10 µm diameter MRG fiber.

from pyfibers import build_fiber, FiberModel

# create fiber model
n_nodes = 25
fiber = build_fiber(
    FiberModel.MRG_INTERPOLATION, diameter=10, n_nodes=n_nodes
)  # diameter in µm

import pyfibers

# Enable INFO level logging to see simulation progress
pyfibers.enable_logging()

Simulation setup

Before we can run a simulation, we need to create a stimulation waveform (i.e., I(t), the time‐course of the extracellular stimulation). We use a biphasic rectangular pulse in this tutorial. We also specify our simulation parameters as constants.

See stimulation module or the Custom simulation code for more information on creating different stimulation waveforms.

from scipy.interpolate import interp1d

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

Plot the waveform to see what it looks like.

import matplotlib.pyplot as plt
import numpy as np

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()

We also need extracellular potentials along the length of the fiber. We will use an extracellular point source for this tutorial, 250 µm from the fiber, positioned over its center. Our fiber has 265 sections, so we need 265 potentials (i.e., one for the middle of each section). Typically, extracellular potentials are generated with a “unit” current source (i.e., 1 mA).

In many cases, users will want to use electrical potentials generated using an outside source (e.g., a finite element model). See our documentation on extracellular potentials for more information.

fiber.potentials = fiber.point_source_potentials(
    0, 250, fiber.length / 2, 1, 10
)  # x, y, z (µm), i0 (mA), sigma (S/m)

plt.plot(fiber.longitudinal_coordinates, fiber.potentials[0])
plt.xlabel('Distance along fiber (µm)')
plt.ylabel('Electrical potential (mV)')
plt.title('Extracellular potentials')
plt.show()

Running a simulation

To run a simulation, we will use an instance of the ScaledStim class. This class is used to run simulations of a fiber’s response to extracellular stimulation.

See also

stimulation module contains information on ScaledStim parameters.

from pyfibers import ScaledStim

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

ScaledStim: 20 ms (dt=0.001 ms) (t_init_ss=-200 ms, dt_init_ss=5 ms)

We can use the run_sim() method of the ScaledStim class to run a simulation. By default, this method monitors for action potentials at the distal end of the fiber. It returns the number of detected action potentials and the time of the last action potential. Since we used a unit current source (1 mA), our stimamp (stimulation amplitude) here represents the current of the stimulus (in mA).

stimamp = -1.5  # technically unitless, but scales the unit (1 mA) stimulus to 1.5 mA
ap, time = stimulation.run_sim(stimamp, fiber)
print(f'Number of action potentials detected: {ap}')
print(f'Time of last action potential detection: {time} ms')

INFO:pyfibers.stimulation:Running amplitude(s): -1.5
INFO:pyfibers.stimulation:N aps: 1, time 0.379
Number of action potentials detected: 1.0
Time of last action potential detection: 0.379 ms

Run search for activation threshold

If we want to determine the activation threshold, we can use the find_threshold() method. This method returns the stimulation amplitude at which the fiber begins to activate and the number of generated action potentials. The threshold amplitude is calculated using a bisection search (see Algorithms in PyFibers).

amp, ap = stimulation.find_threshold(fiber)
print(f'Activation threshold: {amp} mA')

INFO:pyfibers.stimulation:Beginning bounds search. Checking initial bounds...
INFO:pyfibers.stimulation:Bounds iter.   0: [   -0.0100,    -1.0000]
INFO:pyfibers.stimulation:Running amplitude(s): -1
INFO:pyfibers.stimulation:N aps: 1, time 0.437
INFO:pyfibers.stimulation:Running amplitude(s): -0.01
INFO:pyfibers.stimulation:N aps: 0, time None
INFO:pyfibers.stimulation:Found AP at 0.437 ms, subsequent runs will exit at 5.437 ms. Change 'exit_t_shift' to modify this.
INFO:pyfibers.stimulation:Bounds located after 0 iter.: [   -0.0100,    -1.0000]
INFO:pyfibers.stimulation:Beginning bisection search...
INFO:pyfibers.stimulation:Bisect. iter.   1: [   -0.0100,    -1.0000]  tol=0.9900/0.0100  next=   -0.5050
INFO:pyfibers.stimulation:Running amplitude(s): -0.505
INFO:pyfibers.stimulation:N aps: 0, time None
INFO:pyfibers.stimulation:Bisect. iter.   2: [   -0.5050,    -1.0000]  tol=0.4950/0.0100  next=   -0.7525
INFO:pyfibers.stimulation:Running amplitude(s): -0.7525
INFO:pyfibers.stimulation:N aps: 0, time None
INFO:pyfibers.stimulation:Bisect. iter.   3: [   -0.7525,    -1.0000]  tol=0.2475/0.0100  next=   -0.8762
INFO:pyfibers.stimulation:Running amplitude(s): -0.8762
INFO:pyfibers.stimulation:N aps: 0, time None
INFO:pyfibers.stimulation:Bisect. iter.   4: [   -0.8762,    -1.0000]  tol=0.1238/0.0100  next=   -0.9381
INFO:pyfibers.stimulation:Running amplitude(s): -0.9381
INFO:pyfibers.stimulation:N aps: 0, time None
INFO:pyfibers.stimulation:Bisect. iter.   5: [   -0.9381,    -1.0000]  tol=0.0619/0.0100  next=   -0.9691
INFO:pyfibers.stimulation:Running amplitude(s): -0.9691
INFO:pyfibers.stimulation:N aps: 1, time 0.477
INFO:pyfibers.stimulation:Bisect. iter.   6: [   -0.9381,    -0.9691]  tol=0.0319/0.0100  next=   -0.9536
INFO:pyfibers.stimulation:Running amplitude(s): -0.9536
INFO:pyfibers.stimulation:N aps: 1, time 0.548
INFO:pyfibers.stimulation:Bisect. iter.   7: [   -0.9381,    -0.9536]  tol=0.0162/0.0100  next=   -0.9459
INFO:pyfibers.stimulation:Running amplitude(s): -0.9459
INFO:pyfibers.stimulation:N aps: 0, time None
INFO:pyfibers.stimulation:Threshold found at stimamp = -0.9536 after 7 bisection iterations
INFO:pyfibers.stimulation:Validating threshold...
INFO:pyfibers.stimulation:Running amplitude(s): -0.9536
INFO:pyfibers.stimulation:N aps: 1, time 0.548
Activation threshold: -0.95359375 mA

Intracellular stimulation

We can also use intracellular stimulation by using the IntraStim class.

from pyfibers import IntraStim

# reset fiber potentials to zero to eliminate extracellular stimulation
fiber.potentials = np.zeros(len(fiber.coordinates))

# Intracellular stimulation arguments
clamp_kws = {
    'delay': 1.0,  # milliseconds
    'pw': 0.1,  # milliseconds, pulse width per pulse
    'dur': 10,  # milliseconds, total duration of stimulation
    'freq': 1000,  # pulse repetition frequency in Hz
    'amp': 1,  # nA, intracellular current amplitude. We'll set to 1 and scale with run_sim
}
loc = 0.5  # location of the intracellular electrode along the fiber

stimulation = IntraStim(
    istim_loc=loc, dt=time_step, tstop=time_stop, clamp_kws=clamp_kws
)

# Run simulation with intracellular stimulation
stimamp = 2  # technically unitless, but scales the unit (1 nA) stimulus to 2 nA
ap, time = stimulation.run_sim(stimamp, fiber)
print(f'Number of action potentials detected: {ap}')
print(f'Time of last action potential detection: {time} ms')

INFO:pyfibers.stimulation:Running amplitude(s): 2
INFO:pyfibers.stimulation:Num. APs: 5, last AP time: 9.396
Number of action potentials detected: 5.0
Time of last action potential detection: 9.396 ms

As before, we can run a search for activation threshold.

amp, ap = stimulation.find_threshold(fiber, stimamp_top=1, stimamp_bottom=0)
print(f'Activation threshold: {amp} nA')

INFO:pyfibers.stimulation:Beginning bounds search. Checking initial bounds...
INFO:pyfibers.stimulation:Bounds iter.   0: [   +0.0000,    +1.0000]
INFO:pyfibers.stimulation:Running amplitude(s): 1
INFO:pyfibers.stimulation:Num. APs: 1, last AP time: 1.424
INFO:pyfibers.stimulation:Running amplitude(s): 0
INFO:pyfibers.stimulation:Num. APs: 0, last AP time: None
INFO:pyfibers.stimulation:Found AP at 1.424 ms, subsequent runs will exit at 6.4239999999999995 ms. Change 'exit_t_shift' to modify this.
INFO:pyfibers.stimulation:Bounds located after 0 iter.: [   +0.0000,    +1.0000]
INFO:pyfibers.stimulation:Beginning bisection search...
INFO:pyfibers.stimulation:Bisect. iter.   1: [   +0.0000,    +1.0000]  tol=1.0000/0.0100  next=   +0.5000
INFO:pyfibers.stimulation:Running amplitude(s): 0.5
INFO:pyfibers.stimulation:Num. APs: 0, last AP time: None
INFO:pyfibers.stimulation:Bisect. iter.   2: [   +0.5000,    +1.0000]  tol=0.5000/0.0100  next=   +0.7500
INFO:pyfibers.stimulation:Running amplitude(s): 0.75
INFO:pyfibers.stimulation:Num. APs: 0, last AP time: None
INFO:pyfibers.stimulation:Bisect. iter.   3: [   +0.7500,    +1.0000]  tol=0.2500/0.0100  next=   +0.8750
INFO:pyfibers.stimulation:Running amplitude(s): 0.875
INFO:pyfibers.stimulation:Num. APs: 0, last AP time: None
INFO:pyfibers.stimulation:Bisect. iter.   4: [   +0.8750,    +1.0000]  tol=0.1250/0.0100  next=   +0.9375
INFO:pyfibers.stimulation:Running amplitude(s): 0.9375
INFO:pyfibers.stimulation:Num. APs: 0, last AP time: None
INFO:pyfibers.stimulation:Bisect. iter.   5: [   +0.9375,    +1.0000]  tol=0.0625/0.0100  next=   +0.9688
INFO:pyfibers.stimulation:Running amplitude(s): 0.9688
INFO:pyfibers.stimulation:Num. APs: 0, last AP time: None
INFO:pyfibers.stimulation:Bisect. iter.   6: [   +0.9688,    +1.0000]  tol=0.0312/0.0100  next=   +0.9844
INFO:pyfibers.stimulation:Running amplitude(s): 0.9844
INFO:pyfibers.stimulation:Num. APs: 1, last AP time: 1.485
INFO:pyfibers.stimulation:Bisect. iter.   7: [   +0.9688,    +0.9844]  tol=0.0159/0.0100  next=   +0.9766
INFO:pyfibers.stimulation:Running amplitude(s): 0.9766
INFO:pyfibers.stimulation:Num. APs: 1, last AP time: 1.562
INFO:pyfibers.stimulation:Threshold found at stimamp = 0.9766 after 7 bisection iterations
INFO:pyfibers.stimulation:Validating threshold...
INFO:pyfibers.stimulation:Running amplitude(s): 0.9766
INFO:pyfibers.stimulation:Num. APs: 3, last AP time: 7.5
Activation threshold: 0.9765625 nA

See also

Analyze simulation results shows how to analyze simulation outputs (Vm, gating, currents).