.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples/models/plot_horsager2009.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_examples_models_plot_horsager2009.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_models_plot_horsager2009.py:


===============================================================================
Horsager et al. (2009): Predicting temporal sensitivity
===============================================================================

This example shows how to use the
:py:class:`~pulse2percept.models.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 :py:mod:`~pulse2percept.datasets`
subpackage.

The dataset can be loaded as follows:

.. GENERATED FROM PYTHON SOURCE LINES 20-25

.. code-block:: Python


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





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    (608, 21)



.. GENERATED FROM PYTHON SOURCE LINES 27-36

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:

.. GENERATED FROM PYTHON SOURCE LINES 36-41

.. code-block:: Python


    single_pulse = load_horsager2009(subjects='S05', electrodes='C3',
                                     stim_types='single_pulse')
    single_pulse






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>subject</th>
          <th>implant</th>
          <th>electrode</th>
          <th>task</th>
          <th>stim_type</th>
          <th>stim_dur</th>
          <th>stim_freq</th>
          <th>stim_amp</th>
          <th>pulse_type</th>
          <th>pulse_dur</th>
          <th>pulse_num</th>
          <th>interphase_dur</th>
          <th>delay_dur</th>
          <th>ref_stim_type</th>
          <th>ref_freq</th>
          <th>ref_amp</th>
          <th>ref_amp_factor</th>
          <th>ref_pulse_dur</th>
          <th>ref_interphase_dur</th>
          <th>theta</th>
          <th>source</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>S05</td>
          <td>ArgusI</td>
          <td>C3</td>
          <td>threshold</td>
          <td>single_pulse</td>
          <td>200.0</td>
          <td>NaN</td>
          <td>179.792755</td>
          <td>cathodic_first</td>
          <td>0.075</td>
          <td>1.0</td>
          <td>0.075</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>110.3</td>
          <td>fig_3</td>
        </tr>
        <tr>
          <th>1</th>
          <td>S05</td>
          <td>ArgusI</td>
          <td>C3</td>
          <td>threshold</td>
          <td>single_pulse</td>
          <td>200.0</td>
          <td>NaN</td>
          <td>97.391689</td>
          <td>cathodic_first</td>
          <td>0.150</td>
          <td>1.0</td>
          <td>0.150</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>110.3</td>
          <td>fig_3</td>
        </tr>
        <tr>
          <th>2</th>
          <td>S05</td>
          <td>ArgusI</td>
          <td>C3</td>
          <td>threshold</td>
          <td>single_pulse</td>
          <td>200.0</td>
          <td>NaN</td>
          <td>64.690443</td>
          <td>cathodic_first</td>
          <td>0.220</td>
          <td>1.0</td>
          <td>0.220</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>110.3</td>
          <td>fig_3</td>
        </tr>
        <tr>
          <th>3</th>
          <td>S05</td>
          <td>ArgusI</td>
          <td>C3</td>
          <td>threshold</td>
          <td>single_pulse</td>
          <td>200.0</td>
          <td>NaN</td>
          <td>33.039171</td>
          <td>cathodic_first</td>
          <td>0.530</td>
          <td>1.0</td>
          <td>0.530</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>110.3</td>
          <td>fig_3</td>
        </tr>
        <tr>
          <th>4</th>
          <td>S05</td>
          <td>ArgusI</td>
          <td>C3</td>
          <td>threshold</td>
          <td>single_pulse</td>
          <td>200.0</td>
          <td>NaN</td>
          <td>29.333355</td>
          <td>cathodic_first</td>
          <td>0.750</td>
          <td>1.0</td>
          <td>0.750</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>110.3</td>
          <td>fig_3</td>
        </tr>
        <tr>
          <th>5</th>
          <td>S05</td>
          <td>ArgusI</td>
          <td>C3</td>
          <td>threshold</td>
          <td>single_pulse</td>
          <td>200.0</td>
          <td>NaN</td>
          <td>24.721303</td>
          <td>cathodic_first</td>
          <td>0.950</td>
          <td>1.0</td>
          <td>0.950</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>110.3</td>
          <td>fig_3</td>
        </tr>
        <tr>
          <th>6</th>
          <td>S05</td>
          <td>ArgusI</td>
          <td>C3</td>
          <td>threshold</td>
          <td>single_pulse</td>
          <td>200.0</td>
          <td>NaN</td>
          <td>18.344511</td>
          <td>cathodic_first</td>
          <td>2.000</td>
          <td>1.0</td>
          <td>2.000</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>110.3</td>
          <td>fig_3</td>
        </tr>
        <tr>
          <th>7</th>
          <td>S05</td>
          <td>ArgusI</td>
          <td>C3</td>
          <td>threshold</td>
          <td>single_pulse</td>
          <td>200.0</td>
          <td>NaN</td>
          <td>14.649325</td>
          <td>cathodic_first</td>
          <td>4.000</td>
          <td>1.0</td>
          <td>4.000</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>NaN</td>
          <td>110.3</td>
          <td>fig_3</td>
        </tr>
      </tbody>
    </table>
    </div>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 42-50

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:

.. GENERATED FROM PYTHON SOURCE LINES 50-59

.. code-block:: Python


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




.. image-sg:: /examples/models/images/sphx_glr_plot_horsager2009_001.png
   :alt: plot horsager2009
   :srcset: /examples/models/images/sphx_glr_plot_horsager2009_001.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <Axes: ylabel='0'>



.. GENERATED FROM PYTHON SOURCE LINES 60-64

Simulating the model response
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

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

.. GENERATED FROM PYTHON SOURCE LINES 64-87

.. code-block:: Python


    from pulse2percept.models import Horsager2009Temporal
    model = Horsager2009Temporal()
    model.build()

    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')
    fig.tight_layout()




.. image-sg:: /examples/models/images/sphx_glr_plot_horsager2009_002.png
   :alt: plot horsager2009
   :srcset: /examples/models/images/sphx_glr_plot_horsager2009_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 88-107

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 :math:`\\theta` over time.

The process of finding the stimulus amplitude needed to achieve model output
:math:`\\theta` can be automated with the help of the
:py:meth:`~pulse2percept.models.Horsager2009Temporal.find_threshold` method.

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

.. GENERATED FROM PYTHON SOURCE LINES 107-132

.. code-block:: Python


    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'],
                             interphase_dur=row['interphase_dur'],
                             stim_dur=row['stim_dur'],
                             cathodic_first=True)
        # 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,
                                           bright_tol=0.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.legend()
    plt.title('Fig. 3B: S05 (C3)')




.. image-sg:: /examples/models/images/sphx_glr_plot_horsager2009_003.png
   :alt: Fig. 3B: S05 (C3)
   :srcset: /examples/models/images/sphx_glr_plot_horsager2009_003.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    Text(0.5, 1.0, 'Fig. 3B: S05 (C3)')



.. GENERATED FROM PYTHON SOURCE LINES 133-139

Fixed-duration pulse train thresholds
-------------------------------------

The same procedure can be repeated for
:py:class:`~pulse2percept.stimuli.BiphasicPulseTrain` stimuli to reproduce
Fig. 4.

.. GENERATED FROM PYTHON SOURCE LINES 139-167

.. code-block:: Python


    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'],
                                  interphase_dur=row['interphase_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,
                                           bright_tol=0.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.legend()
    plt.title('Fig. 4B: S05 (C3), 0.075 ms pulse width')





.. image-sg:: /examples/models/images/sphx_glr_plot_horsager2009_004.png
   :alt: Fig. 4B: S05 (C3), 0.075 ms pulse width
   :srcset: /examples/models/images/sphx_glr_plot_horsager2009_004.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    Text(0.5, 1.0, 'Fig. 4B: S05 (C3), 0.075 ms pulse width')



.. GENERATED FROM PYTHON SOURCE LINES 168-185

Other stimuli
-------------

Bursting pulse triplets
^^^^^^^^^^^^^^^^^^^^^^^

"Bursting pulse triplets" as shown in Fig. 7 are readily supported via the
:py:class:`~pulse2percept.stimuli.BiphasicTripletTrain` class.

Variable-duration pulse trains
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

A "variable-duration" pulse train is essentially
:py:class:`~pulse2percept.stimuli.BiphasicPulseTrain` cut to the length of
N pulses.

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

.. GENERATED FROM PYTHON SOURCE LINES 185-198

.. code-block:: Python


    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),
                                                200))
    pt.plot()




.. image-sg:: /examples/models/images/sphx_glr_plot_horsager2009_005.png
   :alt: plot horsager2009
   :srcset: /examples/models/images/sphx_glr_plot_horsager2009_005.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <Axes: ylabel='0'>



.. GENERATED FROM PYTHON SOURCE LINES 199-210

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:

.. GENERATED FROM PYTHON SOURCE LINES 210-222

.. code-block:: Python


    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)








.. GENERATED FROM PYTHON SOURCE LINES 223-226

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.

.. GENERATED FROM PYTHON SOURCE LINES 226-236

.. code-block:: Python


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








.. GENERATED FROM PYTHON SOURCE LINES 237-238

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

.. GENERATED FROM PYTHON SOURCE LINES 238-243

.. code-block:: Python


    anod_stand = MonophasicPulse(0.5 * amp_th, phase_dur, delay_dur=20)

    anod_test = MonophasicPulse(amp_test, phase_dur, delay_dur=delay_dur)








.. GENERATED FROM PYTHON SOURCE LINES 244-245

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

.. GENERATED FROM PYTHON SOURCE LINES 245-250

.. code-block:: Python


    from pulse2percept.stimuli import Stimulus

    latent_add = cath_stand.append(cath_test).append(anod_stand).append(anod_test)
    latent_add.plot()



.. image-sg:: /examples/models/images/sphx_glr_plot_horsager2009_006.png
   :alt: plot horsager2009
   :srcset: /examples/models/images/sphx_glr_plot_horsager2009_006.png
   :class: sphx-glr-single-img


.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <Axes: ylabel='0'>




.. rst-class:: sphx-glr-timing

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


.. _sphx_glr_download_examples_models_plot_horsager2009.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_horsager2009.ipynb <plot_horsager2009.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_horsager2009.py <plot_horsager2009.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_horsager2009.zip <plot_horsager2009.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_