stimulation module

Defines classes for running simulations using model fibers.

This module provides classes and functionalities to manage stimulation of model fibers. It includes enumerations for different threshold, termination, and bounds search modes, as well as a base Stimulation class and the more specialized ScaledStim class for extracellular stimulation and IntraStim class for intracellular stimulation.

enum pyfibers.stimulation.ThresholdCondition(value)

Different threshold search conditions.

• ThresholdCondition.ACTIVATION:

  Search for the minimal stimulus amplitude required to generate an action potential.

• ThresholdCondition.BLOCK:

  Search for the stimulus amplitude required to block propagation (i.e., prevent action potentials).

Member Type:

str

Valid values are as follows:

ACTIVATION = <ThresholdCondition.ACTIVATION: 'activation'>

BLOCK = <ThresholdCondition.BLOCK: 'block'>

enum pyfibers.stimulation.BoundsSearchMode(value)

Modes for adjusting bounds in the bounds search phase of finding threshold.

• BoundsSearchMode.PERCENT_INCREMENT: Adjust bounds by multiplying/dividing by a percentage factor.

• BoundsSearchMode.ABSOLUTE_INCREMENT: Adjust bounds by adding/subtracting a fixed increment.

Member Type:

str

Valid values are as follows:

PERCENT_INCREMENT = <BoundsSearchMode.PERCENT_INCREMENT: 'percent'>

ABSOLUTE_INCREMENT = <BoundsSearchMode.ABSOLUTE_INCREMENT: 'absolute'>

enum pyfibers.stimulation.TerminationMode(value)

Modes for determining when to terminate bisection search phase of finding threshold.

• TerminationMode.PERCENT_DIFFERENCE: Convergence is based on percentage difference between bounds.

• TerminationMode.ABSOLUTE_DIFFERENCE: Convergence is based on the absolute difference between bounds.

Member Type:

str

Valid values are as follows:

PERCENT_DIFFERENCE = <TerminationMode.PERCENT_DIFFERENCE: 'percent'>

ABSOLUTE_DIFFERENCE = <TerminationMode.ABSOLUTE_DIFFERENCE: 'absolute'>

enum pyfibers.stimulation.BisectionMean(value)

Mean type used during bisection search phase of finding threshold.

• BisectionMean.GEOMETRIC: Use geometric mean (sqrt(bottom * top)).

• BisectionMean.ARITHMETIC: Use arithmetic mean ((bottom + top) / 2).

Member Type:

str

Valid values are as follows:

GEOMETRIC = <BisectionMean.GEOMETRIC: 'geometric'>

ARITHMETIC = <BisectionMean.ARITHMETIC: 'arithmetic'>

class pyfibers.stimulation.Stimulation(dt=0.001, tstop=50, t_init_ss=-200, dt_init_ss=5, custom_run_sim=None)

Manage stimulation of NEURON simulations.

Provides methods to configure time stepping, run simulations, and perform threshold search. This class is not meant to be used as is; subclasses should override the Stimulation.run_sim() method. For example, the ScaledStim subclass provides extracellular stimulation capabilities.

For more details on using this class to write custom stimulation routines, see Custom Simulation Code.

Initialize the Stimulation class.

Parameters:

• dt (float) – Time step for the simulation (ms).

• tstop (float) – Total duration of the simulation (ms).

• t_init_ss (float) – Start time (<=0) used to let the system reach steady-state before t=0.

• dt_init_ss (float) – Larger time step used during the steady-state period (ms).

• custom_run_sim (Callable[..., tuple[int, float | None]] | None) – Custom simulation function provided by the user; otherwise, the subclass must override Stimulation.run_sim().

Raises:

RuntimeError – If NEURON sections/fibers have not been created yet.

Variables:

• dt – Time step for the simulation (ms).

• tstop – Total duration of the simulation (ms).

• t_init_ss – Start time (<=0) used to let the system reach steady-state before t=0 (ms).

• dt_init_ss – Larger time step used during the steady-state period (ms).

• custom_run_sim – Custom simulation function provided by the user.

• time – NEURON Vector recording the global simulation time.

run_sim(*args, **kwargs)

Run a simulation using either a custom_run_sim method or a subclass override.

Return type:

tuple[int, float | None]

pre_run_setup(fiber, ap_detect_threshold=-30)

Prepare the simulation environment before running.

This method sets the temperature, initializes the membrane potential, configures AP detection, and optionally balances the membrane currents for certain fiber models (e.g., Tigerholm). It also applies any intracellular stimulation parameters if provided.

Parameters:

• fiber (Fiber) – The Fiber object for which the simulation will be configured.

• ap_detect_threshold (float) – The voltage threshold for detecting action potentials (mV).

Return type:

None

static ap_checker(fiber, ap_detect_location=0.9, precision=3, check_all_apc=True)

Check how many action potentials occurred at a specified node location.

Parameters:

• fiber (Fiber) – The Fiber object to evaluate for APs.

• ap_detect_location (float) – Normalized location along the fiber in [0,1] to check for APs.

• precision (int) – Decimal places to round the detected AP time.

• check_all_apc (bool) – If True, raise a warning if APs occur elsewhere but not at the detect location.

Return type:

tuple[int, float | None]

Returns:

A tuple (num_aps, last_ap_time). If no APs are detected, last_ap_time is None.

Raises:

AssertionError – If the detected AP time is non-positive.

static threshold_checker(fiber, block=False, ap_detect_location=0.9, block_delay=0, thresh_num_aps=1, check_all_apc=True)

Determine whether a stimulation was above or below threshold, for activation or block.

Parameters:

• fiber (Fiber) – The Fiber object to evaluate.

• ap_detect_location (float) – Normalized location in [0,1] where APs are detected.

• block (bool) – If True, check for block threshold; otherwise, check for activation threshold.

• block_delay (float) – Time after simulation start to check for block (ms).

• thresh_num_aps (int) – Number of APs that constitutes a suprathreshold response.

• check_all_apc (bool) – Passed to Stimulation.ap_checker() for additional warnings.

Return type:

bool

Returns:

True if stimulation is suprathreshold; False if subthreshold.

Raises:

• NotImplementedError – If block is True and thresh_num_aps != 1.

• RuntimeError – If no APs are detected at all in a block threshold search.

find_threshold(fiber, condition=ThresholdCondition.ACTIVATION, bounds_search_mode=BoundsSearchMode.PERCENT_INCREMENT, bounds_search_step=10, termination_mode=TerminationMode.PERCENT_DIFFERENCE, termination_tolerance=1, stimamp_top=-1, stimamp_bottom=-0.01, max_iterations=50, exit_t_shift=5, bisection_mean=BisectionMean.ARITHMETIC, block_delay=0, thresh_num_aps=1, silent=False, **kwargs)

Perform a bisection search to find the threshold stimulus amplitude.

This method first locates bounds where one amplitude is subthreshold and

another is suprathreshold (bounds search phase). Then, it repeatedly narrows the bounds until they converge based on the specified termination mode and tolerance (bisection search phase). Ideally, the initial bounds should be chosen such that stimamp_top is supra-threshold and stimamp_bottom is sub-threshold.

Note that enums (ThresholdCondition, BoundsSearchMode,

TerminationMode, BisectionMean) can be provided as an enum member (e.g. ThresholdCondition.ACTIVATION) or as the member’s string value (e.g. “activation”).

See also

For more details on the threshold search process, see Algorithms in PyFibers.

Parameters:

• fiber (Fiber) – Instance of pyfibers.fiber.Fiber to apply stimulation to.

• condition (ThresholdCondition) – The threshold condition (ThresholdCondition.ACTIVATION or ThresholdCondition.BLOCK).

• bounds_search_mode (BoundsSearchMode) – The bounds search mode (BoundsSearchMode.PERCENT_INCREMENT or BoundsSearchMode.ABSOLUTE_INCREMENT).

• bounds_search_step (float) – The iterative increase/decrease of the upper/lower bound during bounds search - if bounds_search_mode is “percent” this is the percentage increase/decrease - if bounds_search_mode is “absolute” this is the absolute increase/decrease

• termination_mode (TerminationMode) – The termination mode (TerminationMode.PERCENT_DIFFERENCE or TerminationMode.ABSOLUTE_DIFFERENCE).

• termination_tolerance (float) – Difference between upper and lower bounds that indicates convergence - absolute difference if termination_mode is “absolute” - percentage difference if termination_mode is “percent”

• stimamp_top (float) – Initial upper-bound stimulus amplitude to test.

• stimamp_bottom (float) – Initial lower-bound stimulus amplitude to test.

• max_iterations (int) – Maximum attempts to find bounding amplitudes before bisection.

• exit_t_shift (float) – Extra time (ms) after an AP is detected, beyond which the simulation can be cut short.

• bisection_mean (BisectionMean) – The bisection mean type (BisectionMean.ARITHMETIC or BisectionMean.GEOMETRIC).

• block_delay (float) – Time (ms) after start to check for a blocked AP, used in block searches.

• thresh_num_aps (int) – Number of action potentials for threshold search - if activation, suprathreshold requires detected aps >= thresh_num_aps - if block, suprathreshold requires detected aps < thresh_num_aps

• silent (bool) – If True, suppress print statements for the search process.

• kwargs – Additional arguments passed to the run_sim method.

Return type:

tuple[float, tuple[int, float | None]]

Returns:

A tuple (threshold_amplitude, (num_detected_aps, last_detected_ap_time)).

Raises:

RuntimeError – If contradictory bounding conditions occur or if the search fails to converge.

supra_exit(fiber, ap_detect_location, thresh_num_aps=1)

Determine if simulation can be exited early, for activation threshold searches only.

Parameters:

• fiber (Fiber) – The Fiber object to check for an action potential.

• ap_detect_location (float) – Normalized location to check in [0,1].

• thresh_num_aps (int) – Number of APs required to consider it suprathreshold.

Return type:

bool

Returns:

True if the specified number of APs has occurred at or before this time.

static end_excitation_checker(fiber, multi_site_check=True, fail_on_end_excitation=True)

Check for end-excitation.

Determines activation sites by finding local minima in each node’s AP time. If an AP is detected near in or adjacent to the passive end nodes, raise an error or issue a warning based on fail_on_end_excitation.

Parameters:

• fiber (Fiber) – The Fiber object to check.

• multi_site_check (bool) – If True, warn if multiple activation sites are detected.

• fail_on_end_excitation (bool | None) – Controls handling of end-excitation: - True: Raise RuntimeError if end-excitation is detected. - False: Only warn if end-excitation is detected. - None: Skip the check entirely.

Return type:

bool

Returns:

True if end excitation is detected, False otherwise.

Raises:

RuntimeError – If end excitation is detected and fail_on_end_excitation is True.

threshsim(stimamp, fiber, condition=ThresholdCondition.ACTIVATION, block_delay=0, thresh_num_aps=1, **kwargs)

Run a single stimulation trial at a given amplitude and check for threshold.

Parameters:

• stimamp (float) – Stimulus amplitude to apply.

• fiber (Fiber) – The Fiber object to stimulate.

• condition (ThresholdCondition) – Threshold condition (ThresholdCondition.ACTIVATION or ThresholdCondition.BLOCK).

• block_delay (float) – If condition=BLOCK, time after which AP detection is considered blocked (ms).

• thresh_num_aps (int) – Number of APs required to be considered suprathreshold.

• kwargs – Additional arguments for the run_sim method.

Return type:

tuple[bool, tuple[int, float | None]]

Returns:

A tuple (is_suprathreshold, (num_aps, last_ap_time)).

class pyfibers.stimulation.ScaledStim(waveform, dt=0.001, tstop=50, t_init_ss=-200, dt_init_ss=5, pad_waveform=True, truncate_waveform=True)

Manage extracellular stimulation of model fibers.

This class takes one or more waveforms, each of which is expected to match one corresponding set of fiber potentials (i.e., one source) in the fiber being stimulated. Therefore, if you have N potential sets on the fiber, you must provide N waveforms, each describing the time course of stimulation for that source.

The waveform can be either:

• A callable that accepts a single float argument for the stimulation time (in ms) and returns

the waveform value at that time.

• A list of such callables with length N if there are N sources

See also

For more information on defining potentials, see Supplying Extracellular Potentials.

** Example Usage **

# Assuming you have already created a fiber "my_fiber"
# Add potentials to the fiber
my_fiber.potentials = electrical_potentials_array


# Create function which takes in a time and returns a waveform value
def my_square_pulse(t: float) -> float:
    if t > 0.1 and t <= 0.2:  # if time is during the pulse
        return 1  # on
    else:  # if time is outside the pulse
        return 0  # off


# Same waveform as above but using scipy.interpolate.interp1d
my_square_pulse = interp1d(
    [0, 0.1, 0.2, 50],  # start, on, off, stop, in ms
    [0, 1, 0, 0],  # waveform values at those times
    kind="previous",
)

# Create a :class:`ScaledStim` object
dt, tstop = 0.001, 50
stim = ScaledStim(waveform=my_square_pulse, dt=dt, tstop=tstop)

# Run the simulation
stim.run_sim(-1, my_fiber)

# Calculate threshold
stim.find_threshold(my_fiber, condition="activation")

Initialize the ScaledStim class.

Parameters:

• waveform (ndarray | list[Callable[[float], float]] | Callable[[float], float]) – Callable or list of callables (e.g., function). Each callable should accept a single float argument for the simulation time (in ms) and return the waveform value at that time. It is recommended that you provide unit waveforms (maximum absolute value of 1). Support for waveform as a numpy array is retained for backwards compatibility, and will be removed in a future release.

• dt (float) – Main simulation time step (ms).

• tstop (float) – Total simulation duration (ms).

• t_init_ss (float) – Time (<=0) to reach steady-state prior to the main simulation.

• dt_init_ss (float) – Large time step used during the steady-state period (ms).

• pad_waveform (bool) – If True, extend the waveform with zeros to match simulation time.

• truncate_waveform (bool) – If True, truncate the waveform if it exceeds the simulation time.

Variables:

• waveform – The waveform(s) to be applied to the fiber potentials. See ScaledStim for more info.

• pad – If True, extend the waveform with zeros to match simulation time.

• truncate – If True, truncate the waveform if it exceeds the simulation time.

run_sim(stimamp, fiber, ap_detect_location=0.9, exit_func=<function ScaledStim.<lambda>>, exit_func_interval=100, exit_func_kws=None, use_exit_t=False, fail_on_end_excitation=True, ap_detect_threshold=-30)

Run a simulation with a given stimulus amplitude(s).

See also

For more information on the underlying math, see the Algorithms docs page.

Parameters:

• stimamp (float | list[float]) – Amplitude to scale the product of extracellular potentials and waveform. - Should be a single float for one source - If stimamp is a single float and there are multiple sources, the same stimamp is applied to all sources - If stimamp is a list of floats, each float is applied to the corresponding source

• fiber (Fiber) – The Fiber object to stimulate.

• ap_detect_location (float) – Normalized location in [0,1] to check for APs.

• exit_func (Callable[..., bool]) – Callback to check if the simulation can be ended early (e.g., upon detection of an AP).

• exit_func_interval (int) – How often (in time steps) to call exit_func.

• exit_func_kws (dict | None) – Additional arguments for exit_func.

• use_exit_t (bool) – If True, simulation will stop after self._exit_t (if set).

• fail_on_end_excitation (bool | None) – Behavior for end excitation detection - If True, raise error if end excitation is detected - If False, continue simulation if end excitation is detected - If None, do not check for end excitation

• ap_detect_threshold (float) – Threshold for detecting action potentials (default: -30 mV)

Raises:

RuntimeError – If NaNs are detected in membrane potentials or if required setup (e.g., istim) is missing.

Return type:

tuple[int, float]

Returns:

Tuple (number_of_APs, time_of_last_AP).

class pyfibers.stimulation.IntraStim(dt=0.001, tstop=50, t_init_ss=-200, dt_init_ss=5, istim_ind=None, istim_loc=None, clamp_kws=None)

Manage intracellular stimulation of model fibers.

This class extends the Stimulation class to provide intracellular stimulation capabilities. The intracellular stimulation is managed via a custom h.trainIClamp mechanism. This mechanism allows for repeated square pulses of current to be injected into a fiber. Its arguments are provided as clamp_kws when creating an instance of this class.

Example Usage

# Assuming you have already created a fiber "my_fiber"
clamp_kws = {
    "delay": 1,  # ms
    "pw": 0.1,  # ms
    "dur": 50,  # ms
    "freq": 100,  # Hz
    "amp": 1,  # nA, recommended to set 1 for scaling purposes
}
istim_ind = 0  # Stimulate the first section
dt, tstop = 0.001, 50
istim = IntraStim(dt=dt, tstop=tstop, istim_ind=istim_ind, clamp_kws=clamp_kws)

# Run a simulation with a given stimulation amplitude
stimamp = 2
istim.run_sim(2, my_fiber)

# Run a threshold search
istim.find_threshold(my_fiber, condition="activation")

Initialize IntracellularStim class.

Parameters:

• dt (float) – time step for simulation [seconds]

• tstop (float) – time step for simulation [seconds]

• t_init_ss (float) – the time (<=0ms) for the system to reach steady state before starting the simulation [ms]

• dt_init_ss (float) – the time step used to reach steady state [ms]

• istim_ind (int) – the Fiber section index (unmyelinated) or node of Ranvier number (myelinated) receiving stimulation

• istim_loc (float) – node location along the Fiber (using NEURON style indexing)

• clamp_kws (dict) – keyword arguments for the h.trainIClamp. All optional, default given in parentheses. - ‘delay’: (0) the delay from the start of the simulation to the onset of the intracellular stimulation [ms] - ‘pw’: (1) the pulse duration of the intracellular stimulation [ms] - ‘dur’: (50) the duration from the start of the simulation to the end of the intracellular stimulation [ms] - ‘freq’: (100) the intracellular pulse repetition rate [Hz] - ‘amp’: (1) the intracellular stimulation amplitude [nA] Note that amp is scaled by the stimamp parameter in run_sim.

Variables:

• istim – the NEURON h.trainIClamp object for intracellular stimulation

• istim_record – the NEURON Vector recording the intracellular stimulation current

• istim_ind – the Fiber section index or node of Ranvier number receiving stimulation

• istim_loc – the node location along the Fiber receiving stimulation

• istim_params – the dictionary of intracellular stimulation parameters (see _add_istim())

_add_istim(fiber)

Create instance of h.trainIClamp for intracellular stimulation.

This method is not called by the user, see :class:IntraStim. Note that trainIClamp is a mod file included in this package. It is an extension of NEURON’s built-in IClamp that allows repeated square pulses.

Parameters:

fiber (Fiber) – The Fiber object to attach intracellular stimulation to.

Return type:

IntraStim

Returns:

The Stimulation instance (self).

run_sim(stimamp, fiber, ap_detect_location=0.9, exit_func=<function IntraStim.<lambda>>, exit_func_interval=100, exit_func_kws=None, use_exit_t=False, fail_on_end_excitation=True, ap_detect_threshold=-30)

Run a simulation for a single stimulation amplitude.

Parameters:

• stimamp (float) – Amplitude to be applied to extracellular stimulation - Should be a single float for one source - If stimamp is a single float and there are multiple sources, the same stimamp is applied to all sources - If stimamp is a list of floats, each float is applied to the corresponding source

• fiber (Fiber) – The Fiber to be stimulated.

• ap_detect_location (float) – Location to detect action potentials (percent along fiber)

• exit_func (Callable) – Function to call to check if simulation should be exited

• exit_func_interval (int) – Interval to call exit_func

• exit_func_kws (dict) – Keyword arguments to pass to exit_func

• use_exit_t (bool) – If True, use the time returned by exit_func as the simulation end time

• fail_on_end_excitation (bool) – Behavior for end excitation detection - if True, raise error if end excitation is detected - if False, continue simulation if end excitation is detected - if None, do not check for end excitation

• ap_detect_threshold (float) – Threshold for detecting action potentials (default: -30 mV)

Raises:

RuntimeError – If NaNs are detected in fiber potentials

Return type:

tuple[int, float]

Returns:

Number of detected APs and time of last detected AP.