Threshold data from Horsager et al. (2009)

This example shows how to use the Horsager et al. (2009) dataset.

[Horsager2009] used a set of psychophysical detection tasks to determine perceptual thresholds in two Argus I users.

Important

You will need to install Pandas (pip install pandas) for this dataset.

A typical task involved turning on a single electrode in the array and asking the subject if they saw something. If the subject did not detect a visual percept (i.e., a phosphene), the stimulus amplitude was increased. Using a staircase procedure, the researchers were able to determine the current amplitude at which subjects were able to detect a phosphene 50% of the time. This amplitude is called the threshold current, and it is stored in the “stim_amp” column of the dataset.

Note

Mean threshold currents were extracted from the main publication and its supplementary materials using Webplot Digitizer. Therefore, these data are provided without any guarantee of correctness or completeness.

Loading the dataset

The dataset can be loaded as a Pandas DataFrame:

from pulse2percept.datasets import load_horsager2009
data = load_horsager2009()
print(data)

    subject implant electrode  ... ref_interphase_dur  theta    source
0       S05  ArgusI        C3  ...                NaN  110.3     fig_3
1       S05  ArgusI        C3  ...                NaN  110.3     fig_3
2       S05  ArgusI        C3  ...                NaN  110.3     fig_3
3       S05  ArgusI        C3  ...                NaN  110.3     fig_3
4       S05  ArgusI        C3  ...                NaN  110.3     fig_3
..      ...     ...       ...  ...                ...    ...       ...
603     S06  ArgusI        D1  ...              0.075  147.0  fig_s3.2
604     S06  ArgusI        D1  ...              0.075  147.0  fig_s3.2
605     S06  ArgusI        D1  ...              0.075  147.0  fig_s3.2
606     S06  ArgusI        D1  ...              0.075  147.0  fig_s3.2
607     S06  ArgusI        D1  ...              0.075  147.0  fig_s3.2

[608 rows x 21 columns]

Inspecting the DataFrame tells us that there are 552 threshold measurements (the rows) each with 21 different attributes (the columns).

These attributes include specifiers such as “subject”, “electrode”, and “stim_type”. We can print all column names using:

data.columns

Index(['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'],
      dtype='object')

Note

The meaning of all column names is explained in the docstring of the load_horsager2009 function.

For example, “stim_type” corresponds to the different stimulus types that were used in the paper:

data.stim_type.unique()

array(['single_pulse', 'fixed_duration', 'variable_duration',
       'fixed_duration_supra', 'bursting_triplets',
       'bursting_triplets_supra', 'latent_addition'], dtype=object)

The first entry, “single_pulse”, corresponds to the single biphasic pulse used to produce Figure 3. We can verify which figure a data point came from by inspecting the “source” column.

To select all the “single_pulse” rows, we can index into the DataFrame as follows:

print(data[data.stim_type == 'single_pulse'])

   subject implant electrode  ... ref_interphase_dur  theta    source
0      S05  ArgusI        C3  ...                NaN  110.3     fig_3
1      S05  ArgusI        C3  ...                NaN  110.3     fig_3
2      S05  ArgusI        C3  ...                NaN  110.3     fig_3
3      S05  ArgusI        C3  ...                NaN  110.3     fig_3
4      S05  ArgusI        C3  ...                NaN  110.3     fig_3
..     ...     ...       ...  ...                ...    ...       ...
75     S06  ArgusI        D1  ...                NaN  132.0  fig_s3.1
76     S06  ArgusI        D1  ...                NaN  132.0  fig_s3.1
77     S06  ArgusI        D1  ...                NaN  132.0  fig_s3.1
78     S06  ArgusI        D1  ...                NaN  132.0  fig_s3.1
79     S06  ArgusI        D1  ...                NaN  132.0  fig_s3.1

[80 rows x 21 columns]

This leaves us with 80 rows, some of which come from Figure 3, others from Figure S3.1 in the supplementary material.

An alternative to indexing into the DataFrame is to load only a subset of the data:

print(load_horsager2009(stim_types='single_pulse'))

   subject implant electrode  ... ref_interphase_dur  theta    source
0      S05  ArgusI        C3  ...                NaN  110.3     fig_3
1      S05  ArgusI        C3  ...                NaN  110.3     fig_3
2      S05  ArgusI        C3  ...                NaN  110.3     fig_3
3      S05  ArgusI        C3  ...                NaN  110.3     fig_3
4      S05  ArgusI        C3  ...                NaN  110.3     fig_3
..     ...     ...       ...  ...                ...    ...       ...
75     S06  ArgusI        D1  ...                NaN  132.0  fig_s3.1
76     S06  ArgusI        D1  ...                NaN  132.0  fig_s3.1
77     S06  ArgusI        D1  ...                NaN  132.0  fig_s3.1
78     S06  ArgusI        D1  ...                NaN  132.0  fig_s3.1
79     S06  ArgusI        D1  ...                NaN  132.0  fig_s3.1

[80 rows x 21 columns]

Plotting the data

Arguably the most important column is “stim_amp”. This is the current amplitude of the different stimuli (single pulse, pulse trains, etc.) used at threshold.

We might be interested in seeing how threshold amplitude varies as a function of pulse duration. We could either use Matplotlib to generate a scatter plot or use pulse2percept’s own visualization function:

from pulse2percept.viz import scatter_correlation
scatter_correlation(data.pulse_dur, data.stim_amp)

<Axes: xlabel='pulse_dur', ylabel='stim_amp'>

scatter_correlation above generates a scatter plot of the stimulus amplitude as a function of pulse duration, and performs linear regression to calculate a correlation $r$ and a $p$ value. As expected from the literature, now it becomes evident that stimulus amplitude is negatively correlated with pulse duration (no matter the exact stimulus used).

Recreating the stimuli

To recreate the stimulus used to obtain a specific data point, we need to use the values specified in the different columns of a particular row.

For example, the first row of the dataset specifies the stimulus used to obtain the threshold current for subject S05 on Electrode C3 using a single biphasic pulse:

row = data.loc[0, :]
print(row)

subject                          S05
implant                       ArgusI
electrode                         C3
task                       threshold
stim_type               single_pulse
stim_dur                       200.0
stim_freq                        NaN
stim_amp                  179.792755
pulse_type            cathodic_first
pulse_dur                      0.075
pulse_num                        1.0
interphase_dur                 0.075
delay_dur                        NaN
ref_stim_type                    NaN
ref_freq                         NaN
ref_amp                          NaN
ref_amp_factor                   NaN
ref_pulse_dur                    NaN
ref_interphase_dur               NaN
theta                          110.3
source                         fig_3
Name: 0, dtype: object

Using BiphasicPulse, we can create the pulse to specifications:

from pulse2percept.stimuli import BiphasicPulse
stim = BiphasicPulse(row['stim_amp'], row['pulse_dur'],
                     interphase_dur=row['interphase_dur'],
                     stim_dur=row['stim_dur'], electrode=row['electrode'])

We can zoom in on the first 10 ms to see the stimulus waveform:

import numpy as np
stim.plot(time=np.arange(0, 10, 0.01))

<Axes: ylabel='C3'>