example_10_extracellular_current_stimulation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "# Example: Extracellular current stimulation #\n",
    "This is an example of **``LFPy``** running in an **``IPython notebook``**. To run through this example code and produce output, press **``<shift-Enter>``** in each code block below.\n",
    "\n",
    "First step is to **import ``LFPy``** and other packages for analysis and plotting:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sys\n",
    "if sys.version < '3':\n",
    "    from urllib2 import urlopen\n",
    "else:    \n",
    "    from urllib.request import urlopen\n",
    "import ssl\n",
    "from warnings import warn\n",
    "import zipfile\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.collections import LineCollection\n",
    "import LFPy\n",
    "import neuron"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fetch Hay et al. 2011 model files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.isfile('L5bPCmodelsEH/morphologies/cell1.asc'):\n",
    "    #get the model files:\n",
    "    u = urlopen('http://senselab.med.yale.edu/ModelDB/eavBinDown.asp?o=139653&a=23&mime=application/zip',\n",
    "                context=ssl._create_unverified_context())\n",
    "    localFile = open('L5bPCmodelsEH.zip', 'wb')\n",
    "    localFile.write(u.read())\n",
    "    localFile.close()\n",
    "    #unzip:\n",
    "    myzip = zipfile.ZipFile('L5bPCmodelsEH.zip', 'r')\n",
    "    myzip.extractall('.')\n",
    "    myzip.close()\n",
    "\n",
    "#compile mod files every time, because of incompatibility with Mainen96 files:\n",
    "if \"win32\" in sys.platform:\n",
    "    pth = \"L5bPCmodelsEH/mod/\"\n",
    "    warn(\"no autompile of NMODL (.mod) files on Windows.\\n\" \n",
    "         + \"Run mknrndll from NEURON bash in the folder L5bPCmodelsEH/mod and rerun example script\")\n",
    "    if not pth in neuron.nrn_dll_loaded:\n",
    "        neuron.h.nrn_load_dll(pth+\"nrnmech.dll\")\n",
    "    neuron.nrn_dll_loaded.append(pth)\n",
    "else:\n",
    "    os.system('''\n",
    "              cd L5bPCmodelsEH/mod/\n",
    "              nrnivmodl\n",
    "              ''')\n",
    "    neuron.load_mechanisms('L5bPCmodelsEH/mod/')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Function to calculate external potential from point current source"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_extracellular_stimuli(dbs_params, cell):\n",
    "    \"\"\" Function to calculate and apply external potential \"\"\"\n",
    "    x0, y0, z0 = dbs_params['position']\n",
    "    sigma = 0.3\n",
    "    ext_field = np.vectorize(lambda x,y,z: 1 / (4 * np.pi* sigma * np.sqrt((x0 - x)**2 + (y0 - y)**2 + (z0 - z)**2)))\n",
    "    \n",
    "    ### MAKING THE EXTERNAL FIELD\n",
    "    n_tsteps = int(cell.tstop / cell.dt + 1)\n",
    "    t = np.arange(n_tsteps) * cell.dt\n",
    "    pulse = np.zeros(n_tsteps)\n",
    "    start_time = dbs_params['start_time']\n",
    "    end_time = dbs_params['end_time']\n",
    "    start_idx = np.argmin(np.abs(t - start_time))\n",
    "    end_idx = np.argmin(np.abs(t - end_time))\n",
    "    pulse[start_idx:end_idx] = dbs_params['amp'] * 1000.\n",
    "\n",
    "    v_cell_ext = np.zeros((cell.totnsegs, n_tsteps))\n",
    "    v_cell_ext[:, :] = ext_field(cell.x.mean(axis=-1), cell.y.mean(axis=-1), cell.z.mean(axis=-1)\n",
    "                                ).reshape(cell.totnsegs, 1) * pulse.reshape(1, n_tsteps)\n",
    "    cell.insert_v_ext(v_cell_ext, t)\n",
    "    return ext_field, pulse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Defining plotting function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotting_function(dbs_params, cell, ext_field, pulse):\n",
    "    \n",
    "    x0, y0, z0 = dbs_params['position']\n",
    "    amp = dbs_params['amp'] * 1000.\n",
    "    \n",
    "    ### PLOTTING THE RESULTS\n",
    "    cell_plot_idxs = [0,  cell.get_closest_idx(x=-150., y=0, z=750.)]    \n",
    "    cell_plot_colors = {cell_plot_idxs[idx]: plt.cm.Greens_r(1./(len(cell_plot_idxs) + 1) * idx + 0.1) \n",
    "                        for idx in range(len(cell_plot_idxs))}\n",
    "\n",
    "    #print cell.xmid[cell_plot_idxs[1]], cell.ymid[cell_plot_idxs[1]], cell.zmid[cell_plot_idxs[1    ]]\n",
    "\n",
    "    # Plotting the morphology\n",
    "    plt.figure(figsize=(16,9))\n",
    "\n",
    "    v_field_ext = np.zeros((50, 200))\n",
    "    xf = np.linspace(np.min(cell.x[:, 1]), np.max(cell.x[:, 1]), 50)\n",
    "    zf = np.linspace(np.min(cell.z[:, 1]), np.max(cell.z[:, 1]), 200)\n",
    "    for xidx, x in enumerate(xf):\n",
    "        for zidx, z in enumerate(zf):\n",
    "            v_field_ext[xidx, zidx] = ext_field(x, 0, z) * amp\n",
    "    plt.subplots_adjust(hspace=0.5)\n",
    "    plt.subplot(121, aspect='equal', xlabel='x [$\\mu m$]', ylabel='z [$\\mu m$]', \n",
    "                xlim=[-400, 400], xticks=[-400, 0, 400], title='Green dots: Measurement points')\n",
    "    plt.imshow(v_field_ext.T, extent=[np.min(cell.x[:, 1]), np.max(cell.x[:, 1]), \n",
    "                                      np.min(cell.z[:, 1]), np.max(cell.z[:, 1])], \n",
    "               origin='lower', interpolation='nearest', cmap=plt.cm.bwr_r, vmin=-150, vmax=150)\n",
    "\n",
    "    plt.plot(cell.x.T, cell.z.T, 'k')\n",
    "    \n",
    "    \n",
    "    plt.colorbar(label='mV')\n",
    "    [plt.plot([cell.x[idx, 0], cell.x[idx, 0]], [cell.z[idx, 0], cell.z[idx, 0]], c='gray', zorder=1) \n",
    "     for idx in range(cell.totnsegs)]\n",
    "    [plt.plot(cell.x.mean(axis=-1)[idx], cell.z.mean(axis=-1)[idx], 'o', c=cell_plot_colors[idx], ms=12) \n",
    "     for idx in cell_plot_idxs]\n",
    "    \n",
    "    l, = plt.plot(x0, z0, 'y*', ms=12)\n",
    "    plt.legend([l], [\"point current source\"], frameon=False)\n",
    "\n",
    "    # Plotting the membrane potentials\n",
    "    plt.subplot(222, title='Membrane potential', xlabel='Time [ms]', ylabel='mV', ylim=[-80, 20])\n",
    "    [plt.plot(cell.tvec, cell.vmem[idx, :], c=cell_plot_colors[idx], lw=2) for idx in cell_plot_idxs]\n",
    "\n",
    "    # Plotting the input currents\n",
    "    ax1 = plt.subplot(224, ylim=[-2*np.max(np.abs(pulse / 1000)), 2*np.max(np.abs(pulse / 1000))], \n",
    "                      ylabel='$\\mu$A', title='Injected current')\n",
    "    ax1.plot(cell.tvec, pulse / 1000)\n",
    "\n",
    "    #plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define cell parameters used as input to cell-class\n",
    "cellParameters = {\n",
    "    'morphology'    : 'L5bPCmodelsEH/morphologies/cell1.asc',\n",
    "    'templatefile'  : ['L5bPCmodelsEH/models/L5PCbiophys3.hoc',\n",
    "                       'L5bPCmodelsEH/models/L5PCtemplate.hoc'],\n",
    "    'templatename'  : 'L5PCtemplate',\n",
    "    'templateargs'  : 'L5bPCmodelsEH/morphologies/cell1.asc',\n",
    "    'passive' : False,\n",
    "    'nsegs_method' : None,\n",
    "    'dt' : 2**-6,\n",
    "    'tstart' : -159,\n",
    "    'tstop' : 30,\n",
    "    'v_init' : -60,\n",
    "    'celsius': 34,\n",
    "    'pt3d' : True,\n",
    "    'extracellular': True,\n",
    "}\n",
    "\n",
    "# delete old sections from NEURON namespace\n",
    "LFPy.cell.neuron.h(\"forall delete_section()\")\n",
    "\n",
    "# Initialize cell instance, using the LFPy.Cell class\n",
    "cell = LFPy.TemplateCell(**cellParameters)\n",
    "cell.set_rotation(x=4.729, y=-3.166)\n",
    "\n",
    "# The parameters of the extracellular point current source \n",
    "dbs_params = {'position': [-50., 0, 0.],\n",
    "              'amp': -15., # uA,\n",
    "              'start_time': 10.,\n",
    "              'end_time': 12.,\n",
    "              }\n",
    "\n",
    "ext_field, pulse = make_extracellular_stimuli(dbs_params, cell)\n",
    "\n",
    "cell.simulate(rec_vmem=True)\n",
    "\n",
    "plotting_function(dbs_params, cell, ext_field, pulse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}