{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Horsager et al. (2009): Predicting temporal sensitivity\n\nThis example shows how to use the\n:py:class:`~pulse2percept.models.Horsager2009Model`.\n\nThe model introduced in [Horsager2009]_ assumes that electrical stimulation\nleads to percepts that quickly increase in brightness (over the time course\nof ~100ms) and then slowly fade away (over the time course of seconds).\n\nThe model was fit to perceptual sensitivity data for a number of different\npulse trains, which are available in the :py:mod:`~pulse2percept.datasets`\nsubpackage.\n\nThe dataset can be loaded as follows:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from pulse2percept.datasets import load_horsager2009\ndata = load_horsager2009()\ndata.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Single-pulse thresholds\n\n### Loading the data\n\nThe data includes a number of thresholds measured on single-pulse stimuli.\nWe can load a subset of these data; for example, for subject S05 and\nElectrode C3:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "single_pulse = load_horsager2009(subjects='S05', electrodes='C3',\n stim_types='single_pulse')\nsingle_pulse" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Creating the stimulus\n\nTo recreate Fig. 3 in the paper, where the model fit to single-pulse stimuli\nis shown, we first need to recreate the stimulus used in the figure.\n\nFor example, we can create a stimulus from a single biphasic pulse\n(0.075 ms phase duration) with amplitude 180 uA, lasting 200 ms in total:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nfrom pulse2percept.stimuli import BiphasicPulse\nphase_dur = 0.075\nstim_dur = 200\npulse = BiphasicPulse(180, phase_dur, interphase_dur=phase_dur,\n stim_dur=stim_dur, cathodic_first=True)\npulse.plot(time=np.linspace(0, 10, num=10000))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Simulating the model response\n\nThe model's response to this stimulus can be visualized as follows:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from pulse2percept.models import Horsager2009Temporal\nmodel = Horsager2009Temporal()\nmodel.build()\n\npercept = model.predict_percept(pulse, t_percept=np.arange(stim_dur))\n\nmax_bright = percept.data.max()\n\nimport matplotlib.pyplot as plt\nfig, ax = plt.subplots(figsize=(12, 5))\nax.plot(pulse.time, -20 + 10 * pulse.data[0, :] / pulse.data.max(),\n linewidth=3, label='pulse')\nax.plot(percept.time, percept.data[0, 0, :], linewidth=3, label='percept')\nax.plot([0, stim_dur], [max_bright, max_bright], 'k--', label='max brightness')\nax.plot([0, stim_dur], [0, 0], 'k')\n\nax.set_xlabel('Time (s)')\nax.set_ylabel('Predicted brightness (a.u.)')\nax.set_xlim(0, stim_dur)\nfig.legend(loc='center right')\nfig.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Finding the threshold current\n\nFinally, we need to find the \"threshold current\" to ultimately reproduce\nFig. 3.\nIn the real world, the threshold current is defined as the stimulus amplitude\nneeded to elicit a detectable phosphene (e.g.) 50% of the time.\nThis threshold current typically differs for every stimulus, stimulated\nelectrode, and patient.\n\nIn the model, there is no notion of \"seeing something 50% of the time\".\nInstead, the model was assumed to reach threshold if the model response\nexceeded some constant $\\\\theta$ over time.\n\nThe process of finding the stimulus amplitude needed to achieve model output\n$\\\\theta$ can be automated with the help of the\n:py:meth:`~pulse2percept.models.Horsager2009Temporal.find_threshold` method.\n\nWe will run this method on every data point from the ones selected above:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "amp_th = []\nfor _, row in single_pulse.iterrows():\n # Set up a biphasic pulse with amplitude 1uA - the amplitude will be\n # up-and-down regulated by find_threshold until the output matches\n # theta:\n stim = BiphasicPulse(1, row['pulse_dur'],\n interphase_dur=row['interphase_dur'],\n stim_dur=row['stim_dur'],\n cathodic_first=True)\n # Find the current that gives model output theta. Search amplitudes in the\n # range [0, 300] uA. Stop the search once the candidate amplitudes are\n # within 1 uA, or the model output is within 0.1 of theta:\n amp_th.append(model.find_threshold(stim, row['theta'],\n amp_range=(0, 300), amp_tol=1,\n bright_tol=0.1))\n\nplt.semilogx(single_pulse.pulse_dur, single_pulse.stim_amp, 's', label='data')\nplt.semilogx(single_pulse.pulse_dur, amp_th, 'k-', linewidth=3, label='model')\nplt.xticks([0.1, 1, 4])\nplt.xlabel('pulse duration (ms)')\nplt.ylabel('threshold current (uA)')\nplt.legend()\nplt.title('Fig. 3B: S05 (C3)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fixed-duration pulse train thresholds\n\nThe same procedure can be repeated for\n:py:class:`~pulse2percept.stimuli.BiphasicPulseTrain` stimuli to reproduce\nFig. 4.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from pulse2percept.stimuli import BiphasicPulseTrain\n\n# Load the data:\nfixed_dur = data[(data.stim_type == 'fixed_duration') &\n (data.subject == 'S05') &\n (data.electrode == 'C3') &\n (data.pulse_dur == 0.075)]\n\n# Find the threshold:\namp_th = []\nfor _, row in fixed_dur.iterrows():\n stim = BiphasicPulseTrain(row['stim_freq'], 1, row['pulse_dur'],\n interphase_dur=row['interphase_dur'],\n stim_dur=row['stim_dur'], cathodic_first=True)\n amp_th.append(model.find_threshold(stim, row['theta'],\n amp_range=(0, 300), amp_tol=1,\n bright_tol=0.1))\n\nplt.semilogx(fixed_dur.stim_freq, fixed_dur.stim_amp, 's', label='data')\nplt.semilogx(fixed_dur.stim_freq, amp_th, 'k-', linewidth=3, label='model')\nplt.xticks([5, 15, 75, 225])\nplt.xlabel('frequency (Hz)')\nplt.ylabel('threshold current (uA)')\nplt.legend()\nplt.title('Fig. 4B: S05 (C3), 0.075 ms pulse width')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Other stimuli\n\n### Bursting pulse triplets\n\n\"Bursting pulse triplets\" as shown in Fig. 7 are readily supported via the\n:py:class:`~pulse2percept.stimuli.BiphasicTripletTrain` class.\n\n### Variable-duration pulse trains\n\nA \"variable-duration\" pulse train is essentially\n:py:class:`~pulse2percept.stimuli.BiphasicPulseTrain` cut to the length of\nN pulses.\n\nFor example, the following recreates a pulse train used in Fig. 5B:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from pulse2percept.stimuli import BiphasicPulseTrain\n\nn_pulses = 2\nfreq = 3\namp = 180\nphase_dur = 0.075\npt = BiphasicPulseTrain(freq, amp, phase_dur, interphase_dur=phase_dur,\n n_pulses=n_pulses, cathodic_first=True,\n stim_dur=np.maximum(np.ceil(n_pulses * 1000.0 / freq),\n 200))\npt.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Latent addition\n\n\"Latent addition\" stimuli only show up in the supplementary materials\n(see Fig. S2.2).\n\nThey are pseudo-monophasic pulse pairs, where the anodic phases were\npresented 20 ms after the end of the second cathodic pulse.\n\nThe initial cathodic pulse always has a fixed amplitude of 50% of the single\npulse threshold:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from pulse2percept.stimuli import MonophasicPulse\n\n# Phase duration:\nphase_dur = 0.075\n\n# Single-pulse threshold determines this current:\namp_th = 20\n\n# Cathodic phase of the standard pulse::\ncath_stand = MonophasicPulse(-0.5 * amp_th, phase_dur)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The delay between the start of the conditioning pulse and the start of the\ntest pulse was varied systematically (between 0.15 and 12 ms).\nThe amplitude of the second pulse was varied to determine thresholds.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Delay was varied between 0.15 and 12 ms:\ndelay_dur = 12\n\n# Vary this current to determine threshold:\namp_test = 45\n\n# Cathodic phase of the test pulse (delivered after a delay):\ncath_test = MonophasicPulse(-amp_test, phase_dur, delay_dur=delay_dur)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The anodic phase were always presented 20 ms after the second cathodic phase:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "anod_stand = MonophasicPulse(0.5 * amp_th, phase_dur, delay_dur=20)\n\nanod_test = MonophasicPulse(amp_test, phase_dur, delay_dur=delay_dur)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The last step is to concatenate all the pulses into a single stimulus:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from pulse2percept.stimuli import Stimulus\n\nlatent_add = cath_stand.append(cath_test).append(anod_stand).append(anod_test)\nlatent_add.plot()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.9" } }, "nbformat": 4, "nbformat_minor": 0 }