Kilohertz (kHz) Frequency Stimulation Block

This tutorial will use the fiber model created in the simulation tutorial to run a bisection search for the fiber’s block threshold (i.e., the minimum stimulation amplitude needed to stop the progression of a propagating action potential).

Create the fiber and set up simulation

As before, we create a fiber and electrical potentials.

import pyfibers

# Enable logging to see block threshold search progress
pyfibers.enable_logging()

import numpy as np
import matplotlib.pyplot as plt

import scipy.signal as sg
from pyfibers import build_fiber, FiberModel, ScaledStim

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

# define unscaled extracellular potentials (in response to a unitary current)
fiber.potentials = fiber.point_source_potentials(0, 250, fiber.length / 2, 1, 10)
plt.plot(fiber.longitudinal_coordinates, fiber.potentials)
plt.xlabel('Distance along fiber (μm)')
plt.ylabel('Electrical potential (mV)')
plt.title('Extracellular potentials')

# turn on saving of transmembrane potential (Vm) so that we can plot it later
fiber.record_vm()

[Vector[1],
 Vector[2],
 Vector[3],
 Vector[4],
 Vector[5],
 Vector[6],
 Vector[7],
 Vector[8],
 Vector[9],
 Vector[10],
 Vector[11],
 Vector[12],
 Vector[13],
 Vector[14],
 Vector[15],
 Vector[16],
 Vector[17],
 Vector[18],
 Vector[19],
 Vector[20],
 Vector[21],
 Vector[22],
 Vector[23],
 Vector[24],
 Vector[25]]

We’ll create a 20 kHz square wave, using a longer time step than in previous tutorials so we have time to observe the fiber’s behavior before, during, and after block. We will set the block signal to only be on from t=50 to t=100 ms. Then we’ll add our waveform to a new ScaledStim instance.

# Setup for simulation
time_step = 0.001  # ms
time_stop = 150  # ms


# create time vector to make waveform with scipy.signal
# blank out first and last thirds of simulation so we can see block turn on and off
def waveform(t):  # noqa: D
    on, off = 50, 100
    frequency = 20  # kHz (because our time units are ms)
    if t > on and t < off:
        return sg.square(2 * np.pi * frequency * t)
    return 0


# Create stimulation object
blockstim = ScaledStim(waveform=waveform, dt=time_step, tstop=time_stop)

Let’s plot the whole waveform, as well as zoom in to show where it turns on and off.

# whole waveform (cannot see individual pulses)
time_steps = np.arange(0, time_stop, time_step)
plt.plot(time_steps, [waveform(t) for t in time_steps])
plt.xlim([0, 150])
plt.xlabel('Time (ms)')
plt.ylabel('Unscaled Stimulation Amplitude')
plt.gcf().set_size_inches((10, 5))

# create figure to show kilohertz waveform details (turning on and off)
fig, ax = plt.subplots(1, 2, figsize=(10, 5), sharey=True)

# zoom in: kHz turning on
ax[0].plot(time_steps, [waveform(t) for t in time_steps])
ax[0].set_xlim([49.8, 50.2])
ax[0].set_xlabel('Time (ms)\n\ntruncated x-limits to show waveform')
ax[0].set_ylabel('Unscaled Stimulation Amplitude')

# zoom in: kHz turning off
ax[1].plot(time_steps, [waveform(t) for t in time_steps])
ax[1].set_xlim([99.8, 100.2])
ax[1].set_xlabel('Time (ms)\n\ntruncated x-limits to show waveform')

Text(0.5, 0, 'Time (ms)\n\ntruncated x-limits to show waveform')

Let’s try stimulating our fiber. We will use a stimulation amplitude low enough that the membrane voltage appears unaffected. This uses the run_sim() method.

ap, time = blockstim.run_sim(-0.5, fiber)

INFO:pyfibers.stimulation:Running: -0.5
INFO:pyfibers.stimulation:N aps: 0, time None

Let’s create a function to plot our response to the block waveform, which we’ll use going forward. Call it now to plot the fiber response.

def block_plot(on=None):
    """Plot the response of a fiber to kHz stimulation."""  # noqa: DAR101
    plt.figure(figsize=(10, 5))
    ax = plt.gca()

    on = on or [50, 100]  # Set default

    ax.plot(
        blockstim.time,
        fiber.vm[fiber.loc_index(0.9)],
        label=r"$V_m(t)$ at 90% fiber length",
    )

    ax.set_xlabel('Time (ms)')
    ax.set_ylabel('Transmembrane Potential (mV)')

    ax.axvspan(on[0], on[1], alpha=0.5, color='red', label='block on')

    # label the synapse activation with vertical black lines
    if fiber.stim:
        for istim_time in np.arange(15, 150, 10):
            label = 'Intracellular stim' if istim_time == 15 else '__'
            ax.axvline(
                x=istim_time, color='black', alpha=0.5, label=label, linestyle='--'
            )

    ax.legend(loc='upper right')
    plt.xlim(-5, 155)


block_plot()
plt.ylim([-82, 20])

(-82.0, 20.0)

kHz Excitation (below block threshold)

At higher stimulation amplitudes (but still below block threshold), the high frequency stimulation can cause repeated spiking. Let’s turn up the stimulus and watch this happen.

ap, time = blockstim.run_sim(
    -1.5, fiber
)  # NOTE: could happen either above or below threshold
block_plot()

INFO:pyfibers.stimulation:Running: -1.5
INFO:pyfibers.stimulation:N aps: 27, time 99.386

Transmission (with onset response to kHz)

We wanted to stop action potentials, but we are creating many of them! Let’s try turning up the stimulation. You will see an onset response of spiking caused by the stimulation which fades after some time.

ap, time = blockstim.run_sim(-2.5, fiber)
block_plot()

INFO:pyfibers.stimulation:Running: -2.5
INFO:pyfibers.stimulation:N aps: 2, time 52.566

We need to know whether we are truly blocking. We have been recording from one end of the fiber (0.9). Let’s add spiking activity to the other end of the fiber (0.1). At t=15 ms, the activity turns on and a spike occurs every 10 ms. Note that, as we saw in the plot of our extracellular potentials, stimulation is occurring in the center of the fiber (0.5).

We add intrinsic activity using add_intrinsic_activity().

# add intrinsic activity, run with same kHz amplitude as above plot
# loc = 0.1 to avoid erroneous end excitation
fiber.add_intrinsic_activity(loc=0.1, start_time=15, avg_interval=10, num_stims=14)

# plot potentials along fiber, intrinsic activity location, and action potential monitoring location
plt.figure(figsize=(8, 4))
plt.plot(
    fiber.longitudinal_coordinates,
    fiber.potentials[0],
    label='Extracellular potentials',
)
intrinsic_activity_loc = fiber.longitudinal_coordinates[
    fiber.loc_index(0.1, target='sections')
]
monitoring_loc = fiber.longitudinal_coordinates[fiber.loc_index(0.9, target='sections')]
plt.axvline(
    intrinsic_activity_loc,
    color='r',
    linestyle='--',
    label='Intrinsic activity location',
)
plt.axvline(monitoring_loc, color='b', linestyle='--', label='AP monitoring location')
plt.legend()
plt.xlabel('Distance along fiber (μm)')
plt.ylabel('Electrical potential (mV)')
plt.title('Fiber potentials and activity locations')

Text(0.5, 1.0, 'Fiber potentials and activity locations')

We will stimulate using the same amplitude as before.

ap, time = blockstim.run_sim(-2.5, fiber)
block_plot()

INFO:pyfibers.stimulation:Running: -2.5
INFO:pyfibers.stimulation:N aps: 14, time 145.456

This is called “transmission”, where the action potentials continue unhindered down the axon. The action potential takes time to propagate from one end of the fiber to the other, hence the offset between intracellular stimulation pulses and recorded action potentials.

kHz Block

Finally, if we stimulate at an even higher level, we will observe a complete cessation of spike transmission after the block onset response fades.

ap, time = blockstim.run_sim(-3, fiber)
block_plot()

INFO:pyfibers.stimulation:Running: -3
INFO:pyfibers.stimulation:N aps: 10, time 145.455

Search for kHz block threshold

If we want to determine the block threshold, we can use the find_threshold() method. This method returns the stimulation amplitude at which the fiber begins to block transmission and the number of generated action potentials. We can plot the response of the membrane potential and gating variables at threshold stimulation.

This uses the same algorithm as a search for activation threshold, but here suprathreshold stimulation is considered the absence (rather than the presence) of detected action potentials. Thus, intracellular stimulation is key to determining subthreshold stimulation.

The block_delay parameter defines when the function will start looking for action potentials to evaluate whether transmission is occurring. This should be far enough after block turns on to avoid any false positives due to the onset response. We will also truncate the end of the simulation so that it ends when block turns off.

# create new waveform that starts at t=0 to reduce simulation time
def waveform_truncated(t):  # noqa: D
    on, off = 10, 60
    frequency = 20  # kHz (because our time units are ms)
    if t > on and t < off:
        return sg.square(2 * np.pi * frequency * t)
    return 0


# Assign new waveform
blockstim.waveform = waveform_truncated

# Stop simulation when block turns off
blockstim.tstop = 60


# stimamp_top and stimamp_bottom shown to be above and below block threshold, respectively, in previous cells
amp, _ = blockstim.find_threshold(
    fiber,
    condition="block",
    stimamp_top=-3,
    stimamp_bottom=-2.5,
    exit_t_shift=None,
    block_delay=25,
)

print(f'Block threshold: {amp} mA')

block_plot(on=[10, 60])
plt.xlim(0, 60)

INFO:pyfibers.stimulation:Running: -3
INFO:pyfibers.stimulation:N aps: 2, time 12.655
INFO:pyfibers.stimulation:Running: -2.5
INFO:pyfibers.stimulation:N aps: 6, time 55.745
INFO:pyfibers.stimulation:Search bounds: top=-3, bottom=-2.5
INFO:pyfibers.stimulation:Beginning bisection search
INFO:pyfibers.stimulation:Search bounds: top=-3, bottom=-2.5
INFO:pyfibers.stimulation:Running: -2.75
INFO:pyfibers.stimulation:N aps: 3, time 55.956
INFO:pyfibers.stimulation:Search bounds: top=-3, bottom=-2.75
INFO:pyfibers.stimulation:Running: -2.875
INFO:pyfibers.stimulation:N aps: 2, time 12.637
INFO:pyfibers.stimulation:Search bounds: top=-2.875, bottom=-2.75
INFO:pyfibers.stimulation:Running: -2.8125
INFO:pyfibers.stimulation:N aps: 2, time 12.627
INFO:pyfibers.stimulation:Search bounds: top=-2.8125, bottom=-2.75
INFO:pyfibers.stimulation:Running: -2.78125
INFO:pyfibers.stimulation:N aps: 3, time 56.076
INFO:pyfibers.stimulation:Search bounds: top=-2.8125, bottom=-2.78125
INFO:pyfibers.stimulation:Running: -2.796875
INFO:pyfibers.stimulation:N aps: 3, time 56.243
INFO:pyfibers.stimulation:Search bounds: top=-2.8125, bottom=-2.796875
INFO:pyfibers.stimulation:Running: -2.804688
INFO:pyfibers.stimulation:N aps: 2, time 12.626
INFO:pyfibers.stimulation:Threshold found at stimamp = -2.804688
INFO:pyfibers.stimulation:Validating threshold...
INFO:pyfibers.stimulation:Running: -2.804688
INFO:pyfibers.stimulation:N aps: 2, time 12.626
Block threshold: -2.8046875 mA

(0.0, 60.0)

Re-excitation

It is important to be careful with the upper bound for your block threshold search, as “re-excitation” can occur at suprathreshold stimulation amplitudes. Thus, it is usually safer to set your initial search bounds subthreshold, to avoid any risk of starting with a top bound that causes re-excitation.

blockstim.tstop = 150
blockstim.waveform = waveform  # Set back to original waveform.
ap, time = blockstim.run_sim(
    -250, fiber
)  # NOTE: could happen either above or below threshold
block_plot()

INFO:pyfibers.stimulation:Running: -250
INFO:pyfibers.stimulation:N aps: 22, time 145.471