
.. 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): Temporal sensitivity of the epiretinal percept
===============================================================================

Detection threshold depends on how charge is delivered in time.
[Horsager2009]_ measured thresholds in Argus I users while varying pulse
duration and pulse-train frequency, and fit a cascade of linear filters and a
nonlinearity.

This example reproduces Figs. 3B and 4B: threshold current against pulse
duration and against frequency, for subject S05, electrode C3.

.. GENERATED FROM PYTHON SOURCE LINES 15-24

.. code-block:: Python


    import matplotlib.pyplot as plt
    import numpy as np
    from scipy.optimize import brentq

    from pulse2percept.datasets import load_horsager2009
    from pulse2percept.models.retina import Horsager2009Temporal
    from pulse2percept.stimuli import BiphasicPulse, BiphasicPulseTrain








.. GENERATED FROM PYTHON SOURCE LINES 26-34

What the model does
-------------------

The model is temporal only, so it needs no implant. Brightness rises over
about 100 ms and decays over seconds; this integration window makes
threshold depend on pulse timing.

A single cathodic-first biphasic pulse, 0.075 ms per phase, 180 uA:

.. GENERATED FROM PYTHON SOURCE LINES 34-56

.. code-block:: Python


    model = Horsager2009Temporal()
    model.build()

    stim_dur = 200
    pulse = BiphasicPulse(180, 0.075, interphase_dur=0.075, stim_dur=stim_dur,
                          cathodic_first=True)
    percept = model.predict_percept(pulse, t_percept=np.arange(stim_dur))

    fig, ax = plt.subplots(figsize=(10, 4))
    ax.plot(pulse.time, -20 + 10 * pulse.data[0, :] / pulse.data.max(),
            linewidth=2, label='pulse')
    ax.plot(percept.time, percept.data[0, 0, :], linewidth=2, label='percept')
    ax.axhline(percept.data.max(), color='k', linestyle='--',
               label='max brightness')
    ax.axhline(0, color='k')
    ax.set_xlabel('time (ms)')
    ax.set_ylabel('predicted brightness (a.u.)')
    ax.set_xlim(0, stim_dur)
    ax.legend(loc='center right')
    fig.tight_layout()




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





.. GENERATED FROM PYTHON SOURCE LINES 57-65

Defining threshold
------------------

Behavioral threshold is the amplitude detected on 50% of trials.
[Horsager2009]_ equates it with the amplitude at which the peak model
response reaches a constant :math:`\theta`, fit per subject and electrode
and included in the dataset. Threshold is then a 1D root search over
amplitude:

.. GENERATED FROM PYTHON SOURCE LINES 65-77

.. code-block:: Python



    def threshold_amp(make_stim, theta, amp_range=(0, 300)):
        """Amplitude (uA) whose predicted max brightness matches ``theta``"""
        def objective(amp):
            stim = make_stim(amp)
            percept = model.predict_percept(stim,
                                            t_percept=np.arange(stim.duration))
            return percept.data.max() - theta
        return brentq(objective, *amp_range)









.. GENERATED FROM PYTHON SOURCE LINES 78-83

Threshold vs pulse duration (Fig. 3B)
-------------------------------------

Longer pulses need less current. Each condition is simulated at its
published pulse and interphase duration:

.. GENERATED FROM PYTHON SOURCE LINES 83-106

.. code-block:: Python


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

    amp_th = []
    for _, row in single_pulse.iterrows():
        def pulse_at(amp, row=row):
            return BiphasicPulse(amp, row['pulse_dur'],
                                 interphase_dur=row['interphase_dur'],
                                 stim_dur=row['stim_dur'],
                                 cathodic_first=True)

        amp_th.append(threshold_amp(pulse_at, row['theta']))

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




.. image-sg:: /examples/models/images/sphx_glr_plot_horsager2009_002.png
   :alt: Fig. 3B: S05, electrode C3
   :srcset: /examples/models/images/sphx_glr_plot_horsager2009_002.png
   :class: sphx-glr-single-img


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

 .. code-block:: none


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



.. GENERATED FROM PYTHON SOURCE LINES 107-113

Threshold vs frequency (Fig. 4B)
--------------------------------

At fixed train duration, higher frequencies deliver more pulses into the
integration window, so less current per pulse is needed. The stimulus is now
a :py:class:`~pulse2percept.stimuli.BiphasicPulseTrain`:

.. GENERATED FROM PYTHON SOURCE LINES 113-137

.. code-block:: Python


    fixed_dur = load_horsager2009(subjects='S05', electrodes='C3',
                                  stim_types='fixed_duration')
    fixed_dur = fixed_dur[fixed_dur.pulse_dur == 0.075]

    amp_th = []
    for _, row in fixed_dur.iterrows():
        def train_at(amp, row=row):
            return BiphasicPulseTrain(row['stim_freq'], amp, row['pulse_dur'],
                                      interphase_dur=row['interphase_dur'],
                                      stim_dur=row['stim_dur'],
                                      cathodic_first=True)

        amp_th.append(threshold_amp(train_at, row['theta']))

    plt.figure()
    plt.semilogx(fixed_dur.stim_freq, fixed_dur.stim_amp, 's', label='data')
    plt.semilogx(fixed_dur.stim_freq, amp_th, 'k-', linewidth=2, label='model')
    plt.xticks([5, 15, 75, 225], ['5', '15', '75', '225'])
    plt.xlabel('frequency (Hz)')
    plt.ylabel('threshold current (uA)')
    plt.legend()
    plt.title('Fig. 4B: S05, electrode C3, 0.075 ms pulses')




.. image-sg:: /examples/models/images/sphx_glr_plot_horsager2009_003.png
   :alt: Fig. 4B: S05, electrode C3, 0.075 ms pulses
   :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. 4B: S05, electrode C3, 0.075 ms pulses')



.. GENERATED FROM PYTHON SOURCE LINES 138-153

Both curves follow the measured thresholds over about two orders of
magnitude, using the published :math:`\theta` without further fitting.

What this does not establish
----------------------------

* Temporal only: no phosphene location or shape.
* :math:`\theta` and the filter parameters were fit per subject and
  electrode; thresholds for a new electrode require comparable data.
* Measured on Argus I, cathodic-first pulses, few subjects. Brightness
  above threshold is not modeled here.
* The paper's other conditions use the same procedure: bursting triplets
  (:py:class:`~pulse2percept.stimuli.BiphasicTripletTrain`), variable
  duration trains (``n_pulses``), and latent addition (appended
  :py:class:`~pulse2percept.stimuli.MonophasicPulse` objects).


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

   **Total running time of the script:** (0 minutes 1.569 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>`_
