
.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples\models\plot_granley2021_biphasic.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_granley2021_biphasic.py>`
        to download the full example code.

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

.. _sphx_glr_examples_models_plot_granley2021_biphasic.py:


=========================================================================================
Granley et al. (2021): Effects of Biphasic Pulse Parameters with the BiphasicAxonMapModel
=========================================================================================

This example shows how to use the
:py:class:`~pulse2percept.models.BiphasicAxonMapModel` to model the effects of 
biphasic pulse train parameters phosphene appearance in an epiretinal
implant such as :py:class:`~pulse2percept.implants.ArgusII`. 

Biphasic pulse trains are a commonly used type of stimulus in visual prostheses. 
This model enhances the :py:class:`~pulse2percept.models.AxonMapModel` to reflect
the effects of the amplitude, frequency, and pulse duration on threshold,
phosphene size, brightness, and streak length, according to previous
psychophysical and electrophysiological studies.

The :py:class:`~pulse2percept.models.BiphasicAxonMapModel` shares the same underlying 
assumptions as the axon map model. Namely, an axon's sensitivity to electrical stimulation
is assumed to decay exponentially with...

*  distance along the axon from the soma (:math:`d_s`), with spatial decay
   constant :math:`\lambda`,
*  distance from the stimulated electrode (:math:`d_e`), with spatial decay 
   constant :math:`\rho`.

In the biphasic model, the radial decay rate :math:`\rho` is scaled by :math:`F_{size}`,
the axonal decay rate :math:`\lambda` is scaled by :math:`F_{streak}`, and the brightness 
contribution from each electrode is scaled by :math:`F_{bright}`. These 3 equations are called
effect models. The final equation for the brightness intensity for a pixel located at polar 
coordinates :math:`(r, \theta)` is given by:

.. math::

    I =  \max_{axon}\sum_{elecs}F_\mathrm{bright} \exp\left(\frac{-d_{e}^2}{2\rho^2 F_\mathrm{size} } + 
            \frac{-d_{s}^2}{2\lambda^2 F_\mathrm{streak} }\right).


Basic Model Usage
-----------------
The biphasic axon map model can be instantiated and ran similarly to other models,
with the exception that all stimuli are required to be :py:class:`~pulse2percept.stimuli.BiphasicPulseTrain`

.. GENERATED FROM PYTHON SOURCE LINES 44-52

.. code-block:: Python


    import matplotlib.pyplot as plt
    import numpy as np
    from pulse2percept.implants import ArgusII
    from pulse2percept.models import BiphasicAxonMapModel
    from pulse2percept.stimuli import BiphasicPulseTrain
    model = BiphasicAxonMapModel(rho=200, lam=800)








.. GENERATED FROM PYTHON SOURCE LINES 54-56

Parameters you don't specify will take on default values. You can inspect
all current model parameters as follows:

.. GENERATED FROM PYTHON SOURCE LINES 56-59

.. code-block:: Python


    print(model)





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

 .. code-block:: none

    BiphasicAxonMapModel(ax_segments_range=(0, 50), 
                         axon_pickle='axons.pickle', 
                         axons_range=(-180, 180), 
                         bright_model=DefaultBrightModel, 
                         eye='RE', grid_type='rectangular', 
                         ignore_pickle=False, lam=800, 
                         loc_od=(15.5, 1.5), 
                         min_ax_sensitivity=0.001, 
                         min_current_spread=1e-08, 
                         n_ax_segments=500, n_axons=1000, 
                         n_gray=None, n_jobs=14, n_threads=14, 
                         ndim=[2], noise=None, rho=200, 
                         size_model=DefaultSizeModel, 
                         spatial=BiphasicAxonMapSpatial, 
                         step=0.25, 
                         streak_model=DefaultStreakModel, 
                         temporal=None, thresh_percept=0, 
                         verbose=True, 
                         vfmap=Watson2014Map(ndim=2), 
                         xrange=(-15, 15), yrange=(-15, 15))




.. GENERATED FROM PYTHON SOURCE LINES 60-74

The most important parameters are ``rho`` and ``lam``, which control the 
radial and axonal current spread, respectively. The parameters ``a0``-``a9`` are 
coefficients for the size, streak, and bright models, which will be discussed
later in this example.

The rest of the parameters are shared with 
:py:class:`~pulse2percept.models.AxonMapModel`. For full details on these 
parameters, see the Axon Map Tutorial


Next, build the model to perform expensive, one time calculations,
and specify a visual prosthesis from the
:py:mod:`~pulse2percept.implants` module. Models with an axon map are well 
suited for epiretinal implants, such as Argus II.

.. GENERATED FROM PYTHON SOURCE LINES 74-77

.. code-block:: Python

    model.build()
    implant = ArgusII()








.. GENERATED FROM PYTHON SOURCE LINES 78-90

.. important ::

    You need to build a model only once. After that, you can apply any number
    of stimuli -- or even apply the model to different implants -- without
    having to rebuild (which takes time).

    However, if you change model parameters
    (e.g., by directly setting ``model.a5 = 2``), you will have to
    call ``model.build()`` again for your changes to take effect.


You can visualize the location of the implant and the axon map

.. GENERATED FROM PYTHON SOURCE LINES 90-96

.. code-block:: Python


    model.plot()
    implant.plot()
    plt.show()





.. image-sg:: /examples/models/images/sphx_glr_plot_granley2021_biphasic_001.png
   :alt: plot granley2021 biphasic
   :srcset: /examples/models/images/sphx_glr_plot_granley2021_biphasic_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 97-107

As mentioned above, the Biphasic Axon Map Model only accepts 
:py:class:`~pulse2percept.stimuli.BiphasicPulseTrain`
stimuli with no :py:attr:`~pulse2percept.stimuli.BiphasicPulseTrain.delay_dur`. 
The amplitude given to the BiphasicPulseTrain
is interpreted as amplitude as a factor of threshold (i.e. an amp of 1 means 
1xTh)

You can easily assign BiphasicPulseTrains to electrodes with a dictionary
The following creates a train with 20Hz frequency, 1xTh amplitude, and 0.45ms
pulse / phase duration.

.. GENERATED FROM PYTHON SOURCE LINES 107-111

.. code-block:: Python


    implant.stim = {'A4' : BiphasicPulseTrain(20, 1, 0.45)}
    implant.stim.plot()




.. image-sg:: /examples/models/images/sphx_glr_plot_granley2021_biphasic_002.png
   :alt: plot granley2021 biphasic
   :srcset: /examples/models/images/sphx_glr_plot_granley2021_biphasic_002.png
   :class: sphx-glr-single-img


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

 .. code-block:: none


    <Axes: ylabel='A4'>



.. GENERATED FROM PYTHON SOURCE LINES 112-113

Finally, you can predict the percept resulting from stimulation

.. GENERATED FROM PYTHON SOURCE LINES 113-118

.. code-block:: Python


    percept = model.predict_percept(implant)
    ax = percept.plot()
    ax.set_title('Predicted percept')
    plt.show()



.. image-sg:: /examples/models/images/sphx_glr_plot_granley2021_biphasic_003.png
   :alt: Predicted percept
   :srcset: /examples/models/images/sphx_glr_plot_granley2021_biphasic_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 119-120

Increasing the frequency will make phosphenes brighter

.. GENERATED FROM PYTHON SOURCE LINES 120-128

.. code-block:: Python

    fig, axes = plt.subplots(1, 2, sharex=True, sharey=True)
    implant.stim = {'A4' : BiphasicPulseTrain(50, 1, 0.45)}
    new_percept = model.predict_percept(implant)
    new_percept.plot(ax=axes[1])
    percept.plot(ax=axes[0], vmax=new_percept.max())
    axes[0].set_title("20 Hz")
    axes[1].set_title("40 Hz")
    plt.show()



.. image-sg:: /examples/models/images/sphx_glr_plot_granley2021_biphasic_004.png
   :alt: 20 Hz, 40 Hz
   :srcset: /examples/models/images/sphx_glr_plot_granley2021_biphasic_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 129-131

Note that without setting vmax, matplotlib automatically rescales images to
have the same max brightness and the difference isn't visible

.. GENERATED FROM PYTHON SOURCE LINES 133-134

Increasing amplitude increases both size and brightness

.. GENERATED FROM PYTHON SOURCE LINES 134-143

.. code-block:: Python

    fig, axes = plt.subplots(1, 2, sharex=True, sharey=True)
    implant.stim = {'A4' : BiphasicPulseTrain(20, 3, 0.45)}
    new_percept = model.predict_percept(implant)
    new_percept.plot(ax=axes[1])
    percept.plot(ax=axes[0], vmax=new_percept.max())
    axes[0].set_title("1xTh")
    axes[1].set_title("3xTh")
    plt.show()




.. image-sg:: /examples/models/images/sphx_glr_plot_granley2021_biphasic_005.png
   :alt: 1xTh, 3xTh
   :srcset: /examples/models/images/sphx_glr_plot_granley2021_biphasic_005.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 144-146

Increasing pulse duration decreases threshold, thus indirectly causing an 
increase in size and brightness (amp factor is increased)

.. GENERATED FROM PYTHON SOURCE LINES 146-155

.. code-block:: Python

    fig, axes = plt.subplots(1, 2, sharex=True, sharey=True)
    implant.stim = {'A4' : BiphasicPulseTrain(20, 1, 4)}
    new_percept = model.predict_percept(implant)
    new_percept.plot(ax=axes[1])
    percept.plot(ax=axes[0], vmax=new_percept.max())
    axes[0].set_title("0.45ms")
    axes[1].set_title("4ms")
    plt.show()




.. image-sg:: /examples/models/images/sphx_glr_plot_granley2021_biphasic_006.png
   :alt: 0.45ms, 4ms
   :srcset: /examples/models/images/sphx_glr_plot_granley2021_biphasic_006.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 156-158

If you account for the change in threshold by decreasing amplitude, then 
the only affect of increasing pulse duration is the streak length decreasing

.. GENERATED FROM PYTHON SOURCE LINES 158-167

.. code-block:: Python

    fig, axes = plt.subplots(1, 2, sharex=True, sharey=True)
    implant.stim = {'A4' : BiphasicPulseTrain(20, 0.023835, 20)}
    new_percept = model.predict_percept(implant)
    new_percept.plot(ax=axes[1])
    percept.plot(ax=axes[0], vmax=new_percept.max())
    axes[0].set_title("0.45ms")
    axes[1].set_title("20ms, 0.02xTh")
    plt.show()




.. image-sg:: /examples/models/images/sphx_glr_plot_granley2021_biphasic_007.png
   :alt: 0.45ms, 20ms, 0.02xTh
   :srcset: /examples/models/images/sphx_glr_plot_granley2021_biphasic_007.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 168-172

This illustrates another important point: The amplitude used for the Biphasic
model is relative to the threshold current at 0.45ms pulse duration. Since larger 
pulse durations have been shown to reduce the threshold amplitude needed, the 
0.02xTh amplitude used in the previous plot still is able to produce a phosphene.

.. GENERATED FROM PYTHON SOURCE LINES 174-192

Changing Effect Models
----------------------
All of the 'effects' plotted above (e.g. size increasing with amplitude)
are controlled by the effect models :math:`F_{bright}`, :math:`F_{size}`, and
:math:`F_{streak}`. The variables 
``bright_model``, ``size_model``, and ``streak_model`` encode the 
effects models.

These default to :py:class:`~pulse2percept.models.granley2021.DefaultBrightModel`,
:py:class:`~pulse2percept.models.granley2021.DefaultSizeModel`, and 
:py:class:`~pulse2percept.models.granley2021.DefaultStreakModel` respectively, which
implement the simple scaling functions described in `Granley et al. (2021) <[Granley2021]>`_.


The coefficients ``a0``-``a9`` parametrize these effect models. While the default values
are likely to work for most cases, they can be customized to be patient specific. 
Notice how we only have to change the value given to the `BiphasicAxonMapModel`, 
and it is automatically passed down to the effect models.

.. GENERATED FROM PYTHON SOURCE LINES 194-197

.. code-block:: Python

    model.a5 = 0
    print(model.size_model.a5)





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

 .. code-block:: none

    0




.. GENERATED FROM PYTHON SOURCE LINES 198-203

For example, ``a0`` and ``a1`` control how threshold changed with pulse duration: 
:math:`amp = (A_0*pdur + A_1)^{-1}*amp`. Thus, pulse duration threshold 
scaling can easily be disabled by setting ``a0`` to 0 and ``a1`` to 1. If we increase 
pulse duration like we did previously, we will now see that only streak length decreases, 
and we no longer have to change amplitude to account for change in threshold

.. GENERATED FROM PYTHON SOURCE LINES 203-218

.. code-block:: Python

    model = BiphasicAxonMapModel(rho=200, lam=800)
    model.a0 = 0
    model.a1 = 1
    model.build()
    fig, axes = plt.subplots(1, 2, sharex=True, sharey=True)
    implant.stim = {'A4' : BiphasicPulseTrain(20, 1, 0.45)}
    percept = model.predict_percept(implant)
    implant.stim = {'A4' : BiphasicPulseTrain(20, 1, 20)}
    new_percept = model.predict_percept(implant)
    new_percept.plot(ax=axes[1])
    percept.plot(ax=axes[0], vmax=new_percept.max())
    axes[0].set_title("0.45ms")
    axes[1].set_title("20ms")
    plt.show()




.. image-sg:: /examples/models/images/sphx_glr_plot_granley2021_biphasic_008.png
   :alt: 0.45ms, 20ms
   :srcset: /examples/models/images/sphx_glr_plot_granley2021_biphasic_008.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 219-232

Similarly, ``a2``-``a4`` control brightness scaling; ``a5``-``a6`` control size scaling, and
``a7``-``a9`` control streak length scaling. For more details on these parameters,
see the effect models documentation, or [Granley2021]_ 

Advanced Usage
----------------------

Custom Effect Models
=====================
For most cases, using the provided, default implementation of the effect models
will probably be enough. However, the effect models are completely modular, and 
can be replaced by any python callable with the parameters frequency, amplitude, 
and pulse duration. For example, we can easily change the model to no longer scale size

.. GENERATED FROM PYTHON SOURCE LINES 232-248

.. code-block:: Python

    model = BiphasicAxonMapModel(rho=200, lam=800)
    def size_modulation(freq, amp, pdur):
        return 1
    model.size_model = size_modulation
    model.build()

    fig, axes = plt.subplots(1, 2, sharex=True, sharey=True)
    implant.stim = {'A4' : BiphasicPulseTrain(20, 1, 0.45)}
    percept = model.predict_percept(implant)
    implant.stim = {'A4' : BiphasicPulseTrain(20, 3, 0.45)}
    new_percept = model.predict_percept(implant)
    new_percept.plot(ax=axes[1])
    percept.plot(ax=axes[0], vmax=new_percept.max())
    axes[0].set_title("1xTh")
    axes[1].set_title("3xTh")
    plt.show()



.. image-sg:: /examples/models/images/sphx_glr_plot_granley2021_biphasic_009.png
   :alt: 1xTh, 3xTh
   :srcset: /examples/models/images/sphx_glr_plot_granley2021_biphasic_009.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 249-257

The stimuli with larger amplitude created a brighter, but equally-sized phosphene


The effect models can even be a class, and can have its own parameters, 
which can be shared with the overarching BiphasicAxonMapModel itself (e.g. an effect 
model can depend on ``rho``, and if ``model.rho`` is changed, ``rho`` will also be changed in
the effect model). For an example of this, 
see :py:class:`~pulse2percept.models.granley2021.DefaultSizeModel`


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

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


.. _sphx_glr_download_examples_models_plot_granley2021_biphasic.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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