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

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

.. _sphx_glr_examples_models_plot_dynaphos.py:


===============================================================================================================
van der Grinten, de Ruyter van Steveninck, Lozano et al. (2023): Phosphene simulation using cortical prostheses
===============================================================================================================

This example shows how to apply the
:py:class:`~pulse2percept.models.cortex.DynaphosModel` (original and official 
implementation available `here <https://github.com/neuralcodinglab/dynaphos>`_) to an
:py:class:`~pulse2percept.implants.cortex.Orion` implant.

The Dynaphos model assumes that all stimuli are applied as biphasic pulse
trains. Unlike other models in p2p, this model is not separated into spatial and temporal components, and
must be run as a single composite model.

The model can be instantiated and run in three steps.

Creating the model
------------------

The first step is to instantiate the
:py:class:`~pulse2percept.models.cortex.DynaphosModel` class by calling its
constructor method.

.. GENERATED FROM PYTHON SOURCE LINES 25-34

.. code-block:: Python


    import numpy as np
    import matplotlib as mpl
    import matplotlib.pyplot as plt
    from pulse2percept.stimuli import BiphasicPulseTrain
    from pulse2percept.implants.cortex import Orion
    from pulse2percept.models.cortex import DynaphosModel
    model = DynaphosModel()








.. GENERATED FROM PYTHON SOURCE LINES 36-38

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

.. GENERATED FROM PYTHON SOURCE LINES 38-41

.. code-block:: Python


    print(model)





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

 .. code-block:: none

    DynaphosModel(a50=1.057631326853325e-07, 
                  a_thr=9.141886000943878e-08, dt=20, 
                  excitability=675, freq=300, 
                  grid_type='rectangular', 
                  kappa_trace=13.95528162, n_gray=None, 
                  noise=None, p_dur=0.17, regions=['v1'], 
                  rheobase=23.9, sig_slope=19152642.500946816, 
                  step=0.25, tau_act=111.111111, 
                  tau_trace=1967655.20573, verbose=True, 
                  vfmap=Polimeni2006Map, xrange=(-5, 5), 
                  yrange=(-5, 5))




.. GENERATED FROM PYTHON SOURCE LINES 42-73

This reveals a number of other parameters to set, such as:

* ``xrange``, ``yrange``: The extent of the visual field to be simulated,
  specified as a range of x and y coordinates (in degrees of visual angle,
  or dva). For example, we are currently sampling x values between -20 dva
  and +20dva, and y values between -15 dva and +15 dva.
* ``step``: The resolution (in dva) at which to sample the visual field.
  For example, we are currently sampling at 0.25 dva in both x and y
  direction.
* ``dt``: The time-step of the model in ms. This determines the frame-rate of the
  outputted percept.
* ``regions``: The regions of the visual cortex to simulate. Currently, the
  Dynaphos model is only defined for the v1 region.
* ``freq``, ``p_dur``: The default frequency and pulse duration for the 
  stimulus. This is used to encode a non-pulse train stimulus as a biphasic 
  pulse train.
* ``excitability``:  The excitability constant which determines current 
  spread, in uA/mm^2.
* ``rheobase``: The rheobase current constant, in uA.
* ``tau_trace``: The trace decay constant, in ms.
* ``kappa_trace``: The stimulus effect modifier
* ``tau_act``: The activation decay constant, in ms.
* ``a_thr``: The tissue activation threshold, under which a phosphene is not
  generated. Note that this is a threshold for the internal activation value, 
  unlike other models' ``thresh_percept`` which is a threshold for the final
  percept brightness.
* ``sig_slope``, ``a50``: The slope of the sigmoidal brightness curve and
  activation value for which the brightness reaches half of its maximum.

To change parameter values, either pass them directly to the constructor
above or set them by hand, like this:

.. GENERATED FROM PYTHON SOURCE LINES 73-76

.. code-block:: Python


    model.step = 0.05








.. GENERATED FROM PYTHON SOURCE LINES 77-80

Then build the model. This is a necessary step before you can actually use
the model to predict a percept, as it performs a number of expensive setup
computations (e.g., building the grid):

.. GENERATED FROM PYTHON SOURCE LINES 80-83

.. code-block:: Python


    model.build()





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

 .. code-block:: none


    DynaphosModel(a50=1.057631326853325e-07, 
                  a_thr=9.141886000943878e-08, dt=20, 
                  excitability=675, freq=300, 
                  grid_type='rectangular', 
                  kappa_trace=13.95528162, n_gray=None, 
                  noise=None, p_dur=0.17, regions=['v1'], 
                  rheobase=23.9, sig_slope=19152642.500946816, 
                  step=0.05, tau_act=111.111111, 
                  tau_trace=1967655.20573, verbose=True, 
                  vfmap=Polimeni2006Map, xrange=(-5, 5), 
                  yrange=(-5, 5))



.. GENERATED FROM PYTHON SOURCE LINES 84-102

.. 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 important model parameters outside the constructor
    (e.g., by directly setting ``model.step = 0.25``), you will have to
    call ``model.build()`` again for your changes to take effect.

Assigning a stimulus
--------------------
The second step is to specify a visual prosthesis from the
:py:mod:`~pulse2percept.implants` module.

In the following, we will create an 
:py:class:`~pulse2percept.implants.cortex.Orion` implant 
at the default location.

.. GENERATED FROM PYTHON SOURCE LINES 102-105

.. code-block:: Python


    implant = Orion()








.. GENERATED FROM PYTHON SOURCE LINES 106-108

You can inspect the location of the implant with respect to the visual
cortex using the built-in plot methods:

.. GENERATED FROM PYTHON SOURCE LINES 108-112

.. code-block:: Python


    model.plot()
    implant.plot()




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


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

 .. code-block:: none


    <Axes: xlabel='x (mm)', ylabel='y (mm)'>



.. GENERATED FROM PYTHON SOURCE LINES 113-124

By default, the plots will be added to the current Axes object.
Alternatively, you can pass ``ax=`` to specify in which Axes to plot.

The easiest way to assign a stimulus to the implant is to pass a NumPy array
that specifies the current amplitude to be applied to every electrode in the
implant.

Note that all stimuli passed to Dynaphos must have a time component.

For example, the following sends 100 microamps to all electrodes of the
implant, at 300Hz with a phase duration of 0.17ms:

.. GENERATED FROM PYTHON SOURCE LINES 124-133

.. code-block:: Python


    stim_freq = 300  # stimulus frequency (Hz)
    phase_dur = 0.17  # duration of the cathodic/anodic phase (ms)
    stim_dur = 2000  # stimulus duration (ms)
    stim_amp = 100  # stimulus current (uA)
    implant.stim = {e: BiphasicPulseTrain(amp=stim_amp, freq=stim_freq, 
                                          phase_dur=phase_dur, stim_dur=stim_dur) 
                    for e in implant.electrode_names}








.. GENERATED FROM PYTHON SOURCE LINES 134-138

Predicting the percept
----------------------
The third step is to apply the model to predict the percept resulting from
the specified stimulus. Note that this may take some time on your machine:

.. GENERATED FROM PYTHON SOURCE LINES 138-141

.. code-block:: Python


    percept = model.predict_percept(implant)








.. GENERATED FROM PYTHON SOURCE LINES 142-147

The output of the model is a :py:class:`~pulse2percept.percepts.Percept`
object that contains a time series with the predicted brightness of the
visual percept at every time step.

We can view the brightest frame as follows:

.. GENERATED FROM PYTHON SOURCE LINES 147-150

.. code-block:: Python

    brightest_frame = percept.max(axis='frames')
    plt.imshow(brightest_frame, cmap='gray')




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


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

 .. code-block:: none


    <matplotlib.image.AxesImage object at 0x00000251C610AE10>



.. GENERATED FROM PYTHON SOURCE LINES 151-154

Plotting the percept-over-time next to the applied stimulus reveals 
that the model predicts the perceived brightness to increase rapidly and
then drop off slowly (over the time course of seconds).

.. GENERATED FROM PYTHON SOURCE LINES 154-170

.. code-block:: Python


    fig, ax = plt.subplots(figsize=(10, 5))
    ax.plot(implant.stim.time,
            -0.02 + 0.01 * implant.stim.data[0, :] / implant.stim.data.max(),
            linewidth=3, label='pulse')
    ax.plot(percept.time, np.max(np.max(percept.data, axis=1), axis=0), linewidth=3, label='percept')
    ax.plot([0, stim_dur], [percept.max(), percept.max()], 'k--', label='max brightness')
    ax.plot([0, stim_dur], [0, 0], 'k')

    ax.set_xlabel('time (ms)')
    ax.set_ylabel('predicted brightness (a.u.)')
    ax.set_yticks(np.arange(0, 1.01, 0.1))
    ax.set_xlim(0, stim_dur)
    fig.legend(loc='center')
    fig.tight_layout()




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





.. GENERATED FROM PYTHON SOURCE LINES 171-176

Brightness as a function of amplitude
-------------------------------------

The paper reports that phosphene brightness is affected
by amplitude modulation.

.. GENERATED FROM PYTHON SOURCE LINES 176-199

.. code-block:: Python


    # Use the amplitude values from the paper:
    amps = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100]

    # Output brightness in 20ms time steps for 700ms:
    t_percept = np.arange(0, 700, 20)

    # Initialize an empty list that will contain the brightness-over-time
    # For each amplitude
    bright_over_time = []

    for amp in amps:
        # For each value in the `amps` vector, now stored as `amp`, do:
        # 1. Generate a pulse train with amplitude `amp`, frequency 300Hz,
        #    166ms duration, and pulse duration 0.17ms
        implant.stim = { e: BiphasicPulseTrain(freq=300, amp=amp, phase_dur=0.17,
                                          stim_dur=166)
                        for e in implant.electrode_names }
        # 2. Run the model:
        percept = model.predict_percept(implant, t_percept=t_percept)
        # 3. Save the brightness over time
        bright_over_time.append(np.max(np.max(percept.data, axis=1), axis=0))








.. GENERATED FROM PYTHON SOURCE LINES 200-201

This allows us to reproduce the brightness-over-time plot in Fig. 3: 

.. GENERATED FROM PYTHON SOURCE LINES 201-219

.. code-block:: Python


    fig, ax = plt.subplots(1, figsize=(7.5,3))

    cmap = mpl.cm.YlOrBr
    norm = mpl.colors.Normalize(vmin=0, vmax=100)
    fig.colorbar(None,ax=ax,cmap=cmap,norm=norm,shrink=0.5,orientation='vertical',ticks=[0,100],label="Stim. Amplitude (uA)")

    ax.set_yticks([0.2, 0.6, 1.0])
    ax.set_xticks([0.0, 0.2, 0.4, 0.6])
    ax.set_ylim([0,1.1])
    ax.set_xlim([0,0.7])
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.set(ylabel='Brightness', xlabel='Time [s]')

    for i, amp in enumerate(amps):
        ax.plot(t_percept / 1000, bright_over_time[i], color=cmap(amp/100), linewidth=5)




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





.. GENERATED FROM PYTHON SOURCE LINES 220-224

Note: For the above plot, note that we plot the generated percept
brightness (which is limited by phosphene size & the tissue activation threshold)
while [vanderGrinten2023]_ plots the internal brightness state. The dashed lines in 
[vanderGrinten2023]_ represent time points at which a percept is not generated, 
as the activation threshold was not reached.


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

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


.. _sphx_glr_download_examples_models_plot_dynaphos.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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