Horsager et al. (2009): Predicting temporal sensitivity

This example shows how to use the Horsager2009Model.

The model introduced in [Horsager2009] assumes that electrical stimulation leads to percepts that quickly increase in brightness (over the time course of ~100ms) and then slowly fade away (over the time course of seconds).

The model was fit to perceptual sensitivity data for a number of different pulse trains, which are available in the datasets subpackage.

The dataset can be loaded as follows:

from pulse2percept.datasets import load_horsager2009
data = load_horsager2009()
(608, 21)

Single-pulse thresholds

Loading the data

The data includes a number of thresholds measured on single-pulse stimuli. We can load a subset of these data; for example, for subject S05 and Electrode C3:

single_pulse = load_horsager2009(subjects='S05', electrodes='C3',
subject implant electrode task stim_type stim_dur stim_freq stim_amp pulse_type pulse_dur pulse_num interphase_dur delay_dur ref_stim_type ref_freq ref_amp ref_amp_factor ref_pulse_dur ref_interphase_dur theta source
0 S05 ArgusI C3 threshold single_pulse 200.0 NaN 179.792755 cathodic_first 0.075 1.0 0.075 NaN NaN NaN NaN NaN NaN NaN 110.3 fig_3
1 S05 ArgusI C3 threshold single_pulse 200.0 NaN 97.391689 cathodic_first 0.150 1.0 0.150 NaN NaN NaN NaN NaN NaN NaN 110.3 fig_3
2 S05 ArgusI C3 threshold single_pulse 200.0 NaN 64.690443 cathodic_first 0.220 1.0 0.220 NaN NaN NaN NaN NaN NaN NaN 110.3 fig_3
3 S05 ArgusI C3 threshold single_pulse 200.0 NaN 33.039171 cathodic_first 0.530 1.0 0.530 NaN NaN NaN NaN NaN NaN NaN 110.3 fig_3
4 S05 ArgusI C3 threshold single_pulse 200.0 NaN 29.333355 cathodic_first 0.750 1.0 0.750 NaN NaN NaN NaN NaN NaN NaN 110.3 fig_3
5 S05 ArgusI C3 threshold single_pulse 200.0 NaN 24.721303 cathodic_first 0.950 1.0 0.950 NaN NaN NaN NaN NaN NaN NaN 110.3 fig_3
6 S05 ArgusI C3 threshold single_pulse 200.0 NaN 18.344511 cathodic_first 2.000 1.0 2.000 NaN NaN NaN NaN NaN NaN NaN 110.3 fig_3
7 S05 ArgusI C3 threshold single_pulse 200.0 NaN 14.649325 cathodic_first 4.000 1.0 4.000 NaN NaN NaN NaN NaN NaN NaN 110.3 fig_3

Creating the stimulus

To recreate Fig. 3 in the paper, where the model fit to single-pulse stimuli is shown, we first need to recreate the stimulus used in the figure.

For example, we can create a stimulus from a single biphasic pulse (0.075 ms phase duration) with amplitude 180 uA, lasting 200 ms in total:

import numpy as np
from pulse2percept.stimuli import BiphasicPulse
phase_dur = 0.075
stim_dur = 200
pulse = BiphasicPulse(180, phase_dur, interphase_dur=phase_dur,
                      stim_dur=stim_dur, cathodic_first=True)
pulse.plot(time=np.linspace(0, 10, num=10000))
plot horsager2009
<Axes: ylabel='0'>

Simulating the model response

The model’s response to this stimulus can be visualized as follows:

from pulse2percept.models import Horsager2009Temporal
model = Horsager2009Temporal()

percept = model.predict_percept(pulse, t_percept=np.arange(stim_dur))

max_bright = percept.data.max()

import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(pulse.time, -20 + 10 * pulse.data[0, :] / pulse.data.max(),
        linewidth=3, label='pulse')
ax.plot(percept.time, percept.data[0, 0, :], linewidth=3, label='percept')
ax.plot([0, stim_dur], [max_bright, max_bright], 'k--', label='max brightness')
ax.plot([0, stim_dur], [0, 0], 'k')

ax.set_xlabel('Time (s)')
ax.set_ylabel('Predicted brightness (a.u.)')
ax.set_xlim(0, stim_dur)
fig.legend(loc='center right')
plot horsager2009

Finding the threshold current

Finally, we need to find the “threshold current” to ultimately reproduce Fig. 3. In the real world, the threshold current is defined as the stimulus amplitude needed to elicit a detectable phosphene (e.g.) 50% of the time. This threshold current typically differs for every stimulus, stimulated electrode, and patient.

In the model, there is no notion of “seeing something 50% of the time”. Instead, the model was assumed to reach threshold if the model response exceeded some constant \(\\theta\) over time.

The process of finding the stimulus amplitude needed to achieve model output \(\\theta\) can be automated with the help of the find_threshold method.

We will run this method on every data point from the ones selected above:

amp_th = []
for _, row in single_pulse.iterrows():
        # Set up a biphasic pulse with amplitude 1uA - the amplitude will be
        # up-and-down regulated by find_threshold until the output matches
        # theta:
    stim = BiphasicPulse(1, row['pulse_dur'],
    # Find the current that gives model output theta. Search amplitudes in the
    # range [0, 300] uA. Stop the search once the candidate amplitudes are
    # within 1 uA, or the model output is within 0.1 of theta:
    amp_th.append(model.find_threshold(stim, row['theta'],
                                       amp_range=(0, 300), amp_tol=1,

plt.semilogx(single_pulse.pulse_dur, single_pulse.stim_amp, 's', label='data')
plt.semilogx(single_pulse.pulse_dur, amp_th, 'k-', linewidth=3, label='model')
plt.xticks([0.1, 1, 4])
plt.xlabel('pulse duration (ms)')
plt.ylabel('threshold current (uA)')
plt.title('Fig. 3B: S05 (C3)')
Fig. 3B: S05 (C3)
Text(0.5, 1.0, 'Fig. 3B: S05 (C3)')

Fixed-duration pulse train thresholds

The same procedure can be repeated for BiphasicPulseTrain stimuli to reproduce Fig. 4.

from pulse2percept.stimuli import BiphasicPulseTrain

# Load the data:
fixed_dur = data[(data.stim_type == 'fixed_duration') &
                 (data.subject == 'S05') &
                 (data.electrode == 'C3') &
                 (data.pulse_dur == 0.075)]

# Find the threshold:
amp_th = []
for _, row in fixed_dur.iterrows():
    stim = BiphasicPulseTrain(row['stim_freq'], 1, row['pulse_dur'],
                              stim_dur=row['stim_dur'], cathodic_first=True)
    amp_th.append(model.find_threshold(stim, row['theta'],
                                       amp_range=(0, 300), amp_tol=1,

plt.semilogx(fixed_dur.stim_freq, fixed_dur.stim_amp, 's', label='data')
plt.semilogx(fixed_dur.stim_freq, amp_th, 'k-', linewidth=3, label='model')
plt.xticks([5, 15, 75, 225])
plt.xlabel('frequency (Hz)')
plt.ylabel('threshold current (uA)')
plt.title('Fig. 4B: S05 (C3), 0.075 ms pulse width')
Fig. 4B: S05 (C3), 0.075 ms pulse width
Text(0.5, 1.0, 'Fig. 4B: S05 (C3), 0.075 ms pulse width')

Other stimuli

Bursting pulse triplets

“Bursting pulse triplets” as shown in Fig. 7 are readily supported via the BiphasicTripletTrain class.

Variable-duration pulse trains

A “variable-duration” pulse train is essentially BiphasicPulseTrain cut to the length of N pulses.

For example, the following recreates a pulse train used in Fig. 5B:

from pulse2percept.stimuli import BiphasicPulseTrain

n_pulses = 2
freq = 3
amp = 180
phase_dur = 0.075
pt = BiphasicPulseTrain(freq, amp, phase_dur, interphase_dur=phase_dur,
                        n_pulses=n_pulses, cathodic_first=True,
                        stim_dur=np.maximum(np.ceil(n_pulses * 1000.0 / freq),
plot horsager2009
<Axes: ylabel='0'>

Latent addition

“Latent addition” stimuli only show up in the supplementary materials (see Fig. S2.2).

They are pseudo-monophasic pulse pairs, where the anodic phases were presented 20 ms after the end of the second cathodic pulse.

The initial cathodic pulse always has a fixed amplitude of 50% of the single pulse threshold:

from pulse2percept.stimuli import MonophasicPulse

# Phase duration:
phase_dur = 0.075

# Single-pulse threshold determines this current:
amp_th = 20

# Cathodic phase of the standard pulse::
cath_stand = MonophasicPulse(-0.5 * amp_th, phase_dur)

The delay between the start of the conditioning pulse and the start of the test pulse was varied systematically (between 0.15 and 12 ms). The amplitude of the second pulse was varied to determine thresholds.

# Delay was varied between 0.15 and 12 ms:
delay_dur = 12

# Vary this current to determine threshold:
amp_test = 45

# Cathodic phase of the test pulse (delivered after a delay):
cath_test = MonophasicPulse(-amp_test, phase_dur, delay_dur=delay_dur)

The anodic phase were always presented 20 ms after the second cathodic phase:

The last step is to concatenate all the pulses into a single stimulus:

from pulse2percept.stimuli import Stimulus

latent_add = cath_stand.append(cath_test).append(anod_stand).append(anod_test)
plot horsager2009
<Axes: ylabel='0'>

Total running time of the script: ( 0 minutes 1.368 seconds)

Gallery generated by Sphinx-Gallery