Multi compartmental Cerebellar Granule cell optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Readme\n",
    "\n",
    "\n",
    "* This notebook can be either run with \"Cell\" - \"Run All\", or with \"Cell\" - \"Run Cells\" (from the above menu)\n",
    "  * With the first option (\"Run All\"), the entire script will be automatically executed till the end.\n",
    "  * With the second option (\"Run Cells\"), only the code of a single cell will be executed. Clicking THIS it is possibile to procede, cell by cell, till the end of the code.\n",
    "  \n",
    " ![skipforward](https://raw.githubusercontent.com/antonelepfl/testvue/master/skip_next.png \"size=10\")\n",
    "\n",
    "\n",
    "* Each time a notebook is executed, a kernel status symbol is displayed in the top right corner of the notebook\n",
    "  * When there is a circle bullet •, it means that the kernel is running and the user have to be sure to not interfere with the code execution\n",
    "  * When there is an empty circle bullet ○, it means that the kernel is idle and the user can interfere with the notebook.\n",
    "  \n",
    "See Guidebook at https://collab.humanbrainproject.eu/#/collab/1655/nav/18580"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hint for this notebook\n",
    "\n",
    "Since the simulations are computational intensive, this optimization can be run on the Collab platform with a maximum of 6 individuals for only 6 generations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preparation\n",
    "#### Define the Simulation parameters\n",
    "\n",
    "Define the number of individuals and the number of generations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install NEURON-nightly --upgrade"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!rm -rf /opt/app-root/lib/nrn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NUMBER_INDIVIDUALS=6 # Number of individuals in offspring\n",
    "NUMBER_GENERATIONS=6 #Maximum number of generations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#installation of BluePyOpt\n",
    "!pip3 install bluepyopt --upgrade --no-dependencies >/dev/null"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This custom class define the parameters and section to create a custom built axon.\n",
    "It is used to define the temperature, Vinit and the tables used in the MOD files too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import of BluePyOpt\n",
    "import bluepyopt as bpop\n",
    "import bluepyopt.ephys as ephys\t\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "\n",
    "class GrCMorphologymulti(ephys.morphologies.NrnFileMorphology):\n",
    "  def __init__(self, morphology_path, do_replace_axon):\n",
    "\n",
    "        super(GrCMorphologymulti, self).__init__(morphology_path, do_replace_axon)\n",
    "\t\n",
    "  @staticmethod\n",
    "  def replace_axon(sim=None, icell=None):\n",
    "        \"\"\"Replace axon\"\"\"\n",
    "\n",
    "        # Define origin of distance function\n",
    "        sim.neuron.h.distance(sec=icell.soma[0])\n",
    "\n",
    "        for section in icell.axonal:\n",
    "            sim.neuron.h.delete_section(sec=section)\n",
    "\n",
    "        sim.neuron.h.execute('create axon[2]', icell)\n",
    "\n",
    "#Hillock\n",
    "        icell.axon[0].nseg = 1\n",
    "        icell.axon[0].L = 2.5\n",
    "        icell.axon[0].diam = 1.5\n",
    "        icell.axon[0].Ra = 100\n",
    " \n",
    "#axon \n",
    "        icell.axon[1].nseg = 3\n",
    "        icell.axon[1].L = 70\n",
    "        icell.axon[1].diam = 0.3\n",
    "        icell.axon[1].Ra = 100\n",
    "\n",
    "\n",
    "#Connections of the axon to the hillock and to the soma      \n",
    "        icell.axon[0].connect(icell.soma[0], 1.0, 0.0)\n",
    "        icell.axon[1].connect(icell.axon[0], 1.0, 0.0)\n",
    "        \n",
    "        icell.Hilock.append(sec = icell.axon[0])\n",
    "        icell.AIS.append(sec = icell.axon[1]) \n",
    "    \n",
    "#Temperature and v_init\n",
    "        sim.neuron.h.celsius = 37\n",
    "        sim.neuron.h.v_init = -70\n",
    "        \n",
    "#Ionic channels tables\n",
    "        sim.neuron.h.usetable_GRC_KA = 0\n",
    "        sim.neuron.h.usetable_GRC_KIR  = 0\n",
    "        sim.neuron.h.usetable_GRC_KM = 0\n",
    "        sim.neuron.h.usetable_GRC_KCA = 0\n",
    "        sim.neuron.h.usetable_GRC_CA = 0\n",
    "        sim.neuron.h.usetable_GRC_KV = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To run the optimization, it is necessary to download specific files, from the collab storage space, to the place where the python notebook will run. \n",
    "The files downloaded are: the neuron morphology (asc - neurolucida format) and the ionic channels (mod files)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This command will compile the MOD files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "!unzip -o granule_multi.zip"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!nrnivmodl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First step in the definition of the BluePyOpt Granular cell template.\n",
    "- The morphology, which can be a neurolucida file or an swc file.\n",
    "- The locations for the ionic channels.\n",
    "The basic sectionlists are somatic, axonal, apical and basal. \n",
    "It is possibile to define custom sections and sectionslists"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Morphology location\n",
    "morph = GrCMorphologymulti('GrC2017.asc', do_replace_axon = True)\n",
    "\n",
    "#Locations\n",
    "dend_loc = ephys.locations.NrnSeclistLocation('dend', seclist_name='dend')\n",
    "\n",
    "somatic_loc = ephys.locations.NrnSeclistLocation('somatic', seclist_name='somatic')\n",
    "\n",
    "axon0 = ephys.locations.NrnSeclistLocation('Hilock', seclist_name='Hilock')\n",
    "axon1 = ephys.locations.NrnSeclistLocation('AIS', seclist_name='AIS')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This cell contains all the information to place the passive properties, with specific parameters, on the morphological locations previously defined. For example: the membrane capacitance on all sections, will have a name cm_all, a param name, taken from the mod file, the value, the location and the fact that is a fixed parameter (frozen = True)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#DENDRITES.\n",
    "\n",
    "#Passive properties \n",
    "cm_all = ephys.parameters.NrnSectionParameter(\n",
    "        name='cm_all', \n",
    "        param_name='cm',\n",
    "        value=1,\n",
    "        locations=[dend_loc,somatic_loc, axon0, axon1],\n",
    "        frozen=True)\n",
    "\n",
    "ena = ephys.parameters.NrnSectionParameter(\n",
    "        name='ena_all', \n",
    "        param_name='ena',\n",
    "        value=87.39,\n",
    "        locations=[axon0, axon1],\n",
    "        frozen=True)\n",
    "\n",
    "ek = ephys.parameters.NrnSectionParameter(\n",
    "        name='ek_all', \n",
    "        param_name='ek',\n",
    "        value=-84.69,\n",
    "        locations=[dend_loc, somatic_loc, axon0, axon1],\n",
    "        frozen=True)\n",
    "\n",
    "eleak = ephys.parameters.NrnSectionParameter(\n",
    "        name='eleak_all', \n",
    "        param_name='el_GRC_LKG1',\n",
    "        value = -16.5,\n",
    "        locations=[dend_loc, somatic_loc, axon0, axon1],\n",
    "        frozen=True)\n",
    "\n",
    "eca = ephys.parameters.NrnSectionParameter(\n",
    "        name='eca_all', \n",
    "        param_name='eca',\n",
    "        value = 129.33 ,\n",
    "        locations=[dend_loc, somatic_loc],\n",
    "        frozen=True)\n",
    "\n",
    "Ra = ephys.parameters.NrnSectionParameter(\n",
    "        name='ra_all', \n",
    "        param_name='Ra',\n",
    "        value = 100,\n",
    "        locations=[dend_loc, somatic_loc, axon0, axon1],\n",
    "        frozen=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ionic mechanism are loaded and inserted in the respective sections. The name must be unique instead the prefix is the same as the prefix contained in the MOD files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#CHANNELS\n",
    "\n",
    "#Dendrites\n",
    "GRC_KCA_dend = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='GRC_KCA_dend',\n",
    "        prefix='GRC_KCA',\n",
    "        locations=[dend_loc])\n",
    "\n",
    "GRC_CA_dend = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='GRC_CA_dend',\n",
    "        prefix='GRC_CA',\n",
    "        locations=[dend_loc])\n",
    "\n",
    "GrC_Leak_gaba = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='GRC_LKG2_dend',\n",
    "        prefix='GRC_LKG2',\n",
    "        locations=[dend_loc])\n",
    "\n",
    "#Soma\n",
    "GRC_KA_soma = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='GRC_KA',\n",
    "        prefix='GRC_KA',\n",
    "        locations=[somatic_loc])\n",
    "\n",
    "GRC_KIR_soma = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='GRC_KIR',\n",
    "        prefix='GRC_KIR',\n",
    "        locations=[somatic_loc])\n",
    "\n",
    "GRC_KM_soma = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='GRC_KM',\n",
    "        prefix='GRC_KM',\n",
    "        locations=[somatic_loc])\n",
    "\n",
    "#Mixed soma, dendrites and axon.\n",
    "\n",
    "GrC_Leak = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Leak_all',\n",
    "        prefix='GRC_LKG1',\n",
    "        locations=[dend_loc, somatic_loc, axon0, axon1])\n",
    "\n",
    "GRC_NA = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='GRC_NA',\n",
    "        prefix='GRC_NA',\n",
    "        locations=[axon0, axon1])\n",
    "\n",
    "GRC_KV = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='GRC_KV',\n",
    "        prefix='GRC_KV',\n",
    "        locations=[axon0, axon1])\n",
    "\n",
    "GRC_Calc = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='GRC_CALC',\n",
    "        prefix='GRC_CALC',\n",
    "        locations=[dend_loc, somatic_loc])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the section in which, each single ionic channel, receives its conductance range. The name is to be unique, the param_name is, again, taken from each MOD file. The bounds are the ranges of conductances and, since are a range, are not frozen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#CHANNELS PARAM\n",
    "\n",
    "#Dendrites\n",
    "GrC_Leak_param_dend = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gl_GRC_LKG1_dend',\n",
    "        param_name='gl_GRC_LKG1',\n",
    "        locations=[dend_loc],\n",
    "        bounds=[1.6908225e-06, 0.0002479873],\n",
    "        frozen=False)\n",
    "\n",
    "GrC_Leak_gaba_param_dend = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='ggaba_GRC_LKG2_dend',\n",
    "        param_name='ggaba_GRC_LKG2',\n",
    "        locations=[dend_loc],\n",
    "        bounds=[3.57216e-06, 0.0005239168],\n",
    "        frozen=False)\n",
    "\n",
    "GrC_GRC_KCA_param_dend = ephys.parameters.NrnSectionParameter(\n",
    "        name='gkbar_GRC_KCA_dend',\n",
    "        param_name='gkbar_GRC_KCA',\n",
    "        bounds=[0.000285, 0.0476],\n",
    "        locations=[dend_loc],\n",
    "        frozen=False)\n",
    "\n",
    "GrC_GRC_CA_param_dend = ephys.parameters.NrnSectionParameter(\n",
    "        name='gcabar_GRC_CA_dend',\n",
    "        param_name='gcabar_GRC_CA',\n",
    "        bounds=[0.000438, 0.143],\n",
    "        locations=[dend_loc],\n",
    "        frozen=False)\n",
    "\n",
    "GrC_Calc_param_dend = ephys.parameters.NrnSectionParameter(\n",
    "        name='beta_GRC_CALC_dend',\n",
    "        param_name='beta_GRC_CALC',\n",
    "        bounds=[0.6, 0.61],\n",
    "        locations=[dend_loc],\n",
    "        frozen=False)\n",
    "\n",
    "GrC_cai0_Calc_param_dend = ephys.parameters.NrnSectionParameter(\n",
    "        name='cai0_GRC_CALC_dend',\n",
    "        param_name='cai0_GRC_CALC',\n",
    "        bounds=[0.002, 0.0025],\n",
    "        locations=[dend_loc],\n",
    "        frozen=False)\n",
    "\n",
    "#Soma\n",
    "GRC_KIR_param_soma = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gkbar_GRC_KIR_soma',\n",
    "        param_name='gkbar_GRC_KIR',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[0.000191, 0.0318],\n",
    "        frozen=False)\n",
    "\n",
    "GRC_KA_param_soma = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gkbar_GRC_KA_soma',\n",
    "        param_name='gkbar_GRC_KA',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[0.0004, 0.1],\n",
    "        frozen=False) \n",
    "\n",
    "GrC_Leak_param_soma = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gl_GRC_LKG1_soma',\n",
    "        param_name='gl_GRC_LKG1',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[8.04e-06, 0.00134],\n",
    "        frozen=False)\n",
    "\n",
    "GRC_KM_param_soma = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gkbar_GRC_KM_soma',\n",
    "        param_name='gkbar_GRC_KM',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[1.875e-05, 0.003125],\n",
    "        frozen=False) \n",
    "\n",
    "GrC_Calc_param_soma = ephys.parameters.NrnSectionParameter(\n",
    "        name='beta_GRC_CALC_soma',\n",
    "        param_name='beta_GRC_CALC',\n",
    "        bounds=[1.49, 1.5],\n",
    "        locations=[somatic_loc],\n",
    "        frozen=False)\n",
    "\n",
    "GrC_cai0_Calc_param_soma = ephys.parameters.NrnSectionParameter(\n",
    "        name='cai0_GRC_CALC_soma',\n",
    "        param_name='cai0_GRC_CALC',\n",
    "        bounds=[0.0001, 0.00013],\n",
    "        locations=[somatic_loc],\n",
    "        frozen=False)\n",
    "\n",
    "#Axon\n",
    "#Hilock\n",
    "\n",
    "GRC_NA_param_axon0 = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gnabar_GRC_NA_axon0',\n",
    "        param_name='gnabar_GRC_NA',\n",
    "        locations=[axon0],\n",
    "        bounds=[0.01791, 3.886],\n",
    "        frozen=False)  \n",
    "\n",
    "GRC_KV_param_axon0 = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gkbar_GRC_KV_axon0',\n",
    "        param_name='gkbar_GRC_KV',\n",
    "        locations=[axon0],\n",
    "        bounds=[0.00267, 0.445],\n",
    "        frozen=False) \n",
    "\n",
    "GrC_Leak_param_axon0 = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gl_GRC_LKG1_axon0',\n",
    "        param_name='gl_GRC_LKG1',\n",
    "        locations=[axon0],\n",
    "        bounds=[7.214e-06, 0.0025],\n",
    "        frozen=False)\n",
    "\n",
    "#axon\n",
    "GRC_NA_param_axon1 = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gnabar_GRC_NA_axon1',\n",
    "        param_name='gnabar_GRC_NA',\n",
    "        locations=[axon1],\n",
    "        bounds=[0.000174, 0.029],\n",
    "        frozen=False)  \n",
    "\n",
    "GRC_KV_param_axon1 = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gkbar_GRC_KV_axon1',\n",
    "        param_name='gkbar_GRC_KV',\n",
    "        locations=[axon1],\n",
    "        bounds=[0.00033, 0.0559],\n",
    "        frozen=False) \n",
    "\n",
    "GrC_Leak_param_axon1 = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gl_GRC_LKG1_axon1',\n",
    "        param_name='gl_GRC_LKG1',\n",
    "        locations=[axon1],\n",
    "        bounds=[6.4319e-07, 0.0001071],\n",
    "        frozen=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creating the template\n",
    "\n",
    "To create the cell template, we pass all the objects, previously prepared, to the constructor of the template.\n",
    "The section lists for the axon, are declared here and used in the morphology class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GrCmulti = ephys.models.CellModel(\n",
    "        name='SMGrCmulti',\n",
    "        morph=morph,\n",
    "        mechs=[GRC_CA_dend,\n",
    "               GRC_KCA_dend,\n",
    "               GrC_Leak,\n",
    "               GrC_Leak_gaba,\n",
    "               GRC_Calc,\n",
    "               GRC_KIR_soma, \n",
    "               GRC_KA_soma, \n",
    "               GRC_KM_soma,\n",
    "               GRC_NA,\n",
    "               GRC_KV],\n",
    "        \n",
    "        params=[cm_all,\n",
    "            ena,\n",
    "            ek,\n",
    "            eleak,\n",
    "            eca,\n",
    "            Ra,\n",
    "            GrC_GRC_CA_param_dend, \n",
    "            GrC_GRC_KCA_param_dend, \n",
    "            GrC_Leak_param_dend,\n",
    "            GrC_Leak_gaba_param_dend,\n",
    "            GrC_Calc_param_dend, \n",
    "            GrC_cai0_Calc_param_dend,\n",
    "            GRC_KIR_param_soma,\n",
    "            GRC_KA_param_soma,\n",
    "            GrC_Leak_param_soma,\n",
    "            GRC_KM_param_soma,\n",
    "            GrC_Calc_param_soma,\n",
    "            GrC_cai0_Calc_param_soma,\n",
    "            GRC_NA_param_axon0,\n",
    "            GRC_KV_param_axon0,\n",
    "            GrC_Leak_param_axon0,\n",
    "            GRC_NA_param_axon1,\n",
    "            GRC_KV_param_axon1,\n",
    "            GrC_Leak_param_axon1])  \n",
    "\n",
    "#custom sectionlists for the axon\n",
    "GrCmulti.seclist_names.append('Hilock')\n",
    "GrCmulti.seclist_names.append('AIS')\n",
    "\n",
    "print(GrCmulti)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Location of the electrode for the current injections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soma_loc = ephys.locations.NrnSeclistCompLocation(\n",
    "        name='soma',\n",
    "        seclist_name='somatic',\n",
    "        sec_index=0,\n",
    "        comp_x=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Stimuli, recordings location and protocols. For each protocol there is a recording location and a stimulus, in this case are both placed on the the soma.\n",
    "\n",
    "The delay, duration and total duration can be changed but, more simulation time requires more computational time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sweep_protocols = []\n",
    "for protocol_name, amplitude in [('step1', 0.01), ('step2', 0.016), ('step3', 0.022)]:#\n",
    "    stim = ephys.stimuli.NrnSquarePulse(\n",
    "                step_amplitude=amplitude,\n",
    "                step_delay=100,\n",
    "                step_duration=2000,\n",
    "                location=soma_loc,\n",
    "                total_duration=2200)\n",
    "    rec = ephys.recordings.CompRecording(\n",
    "            name='%s.soma.v' % protocol_name,\n",
    "            location=soma_loc,\n",
    "            variable='v')\n",
    "    protocol = ephys.protocols.SweepProtocol(protocol_name, [stim], [rec])\n",
    "    sweep_protocols.append(protocol)\n",
    "threestep_protocol = ephys.protocols.SequenceProtocol('twostep', protocols=sweep_protocols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Loading of NEURON as the simulator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nrn = ephys.simulators.NrnSimulator()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These conductances are based on the original non optimized model to give an idea of how the optimized cells have to responde to be declared validated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "default_params = {'gcabar_GRC_CA_dend':0.012771,\n",
    "              'gkbar_GRC_KCA_dend':0.003411,\n",
    "              'gl_GRC_LKG1_dend': 0.000025,\n",
    "              'ggaba_GRC_LKG2_dend': 0.000049,\n",
    "              'beta_GRC_CALC_dend': 0.609785, \n",
    "              'cai0_GRC_CALC_dend': 0.0025,\n",
    "              'gkbar_GRC_KIR_soma': 0.002689,\n",
    "              'gkbar_GRC_KA_soma':0.009488,\n",
    "              'gl_GRC_LKG1_soma':0.000127,\n",
    "              'gkbar_GRC_KM_soma':0.000209,\n",
    "              'beta_GRC_CALC_soma':1.496075,\n",
    "              'cai0_GRC_CALC_soma':0.000120,\n",
    "              'gnabar_GRC_NA_axon0': 0.386557,\n",
    "              'gkbar_GRC_KV_axon0':0.027823,\n",
    "              'gl_GRC_LKG1_axon0':0.000120,\n",
    "              'gnabar_GRC_NA_axon1':0.002765,\n",
    "              'gkbar_GRC_KV_axon1':0.003350,\n",
    "              'gl_GRC_LKG1_axon1':0.000010}\n",
    "\n",
    "responses = threestep_protocol.run(cell_model=GrCmulti, param_values=default_params, sim=nrn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the response traces:\n",
    "\n",
    "3 traces are plotted in a graph, each trace generated by a different current injection step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_responses(responses):\n",
    "    plt.subplot(3,1,1)\n",
    "    plt.plot(responses['step1.soma.v']['time'], responses['step1.soma.v']['voltage'], label='step1')\n",
    "    plt.legend()\n",
    "    plt.subplot(3,1,2)\n",
    "    plt.plot(responses['step2.soma.v']['time'], responses['step2.soma.v']['voltage'], label='step2')\n",
    "    plt.legend()\n",
    "    plt.subplot(3,1,3)\n",
    "    plt.plot(responses['step3.soma.v']['time'], responses['step3.soma.v']['voltage'], label='step3')\n",
    "    plt.legend()\n",
    "    plt.tight_layout()\n",
    "\n",
    "plot_responses(responses)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "List of all the features to be used for each step of current injection.\n",
    "All these information were taken from in vitro traces obtain from experiments on mice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "efel_feature_means = {'step1': {'AP_height': 20.93,\n",
    "                    'ISI_CV':0.261,\n",
    "                    'AHP_depth_abs_slow':-52.69,\n",
    "                    'AP_width':0.665,\n",
    "                    'voltage_base':-68.5,\n",
    "                    'AHP_depth_abs':-59.21,\n",
    "                    'time_to_first_spike':31.9,\n",
    "                    'adaptation_index2':0.1062,\n",
    "                    'mean_frequency': 25}, \n",
    "\n",
    "                    'step2': {'AP_height': 19.255,\n",
    "                    'ISI_CV':0.14,\n",
    "                    'AHP_depth_abs_slow':-48.935,\n",
    "                    'AP_width':0.695,\n",
    "                    'voltage_base':-68.77,\n",
    "                    'AHP_depth_abs':-58.3,\n",
    "                    'time_to_first_spike':19.0,\n",
    "                    'adaptation_index2':0.034,\n",
    "                    'mean_frequency': 40}, \n",
    "                  \n",
    "                    'step3': {'AP_height': 17.645,\n",
    "                    'ISI_CV':0.148,\n",
    "                    'AHP_depth_abs_slow':-32.67,\n",
    "                    'AP_width':0.7135,\n",
    "                    'voltage_base':-69.125,\n",
    "                    'AHP_depth_abs':-57.191,\n",
    "                    'time_to_first_spike':14.65,\n",
    "                    'adaptation_index2':0.029,\n",
    "                    'mean_frequency': 50}}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Objectives. The function that try to match what was defined above, for example the speed, with the actual result in the individual obtained during each generation.\n",
    "\n",
    "As per BluePyOpt default, each feature will be an objective."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objectives = []\n",
    "\n",
    "for protocol in sweep_protocols:\n",
    "    stim_start = protocol.stimuli[0].step_delay\n",
    "    stim_end = stim_start + protocol.stimuli[0].step_duration\n",
    "    for efel_feature_name, mean in efel_feature_means[protocol.name].items():\n",
    "        feature_name = '%s.%s' % (protocol.name, efel_feature_name)\n",
    "        feature = ephys.efeatures.eFELFeature(\n",
    "                    feature_name,\n",
    "                    efel_feature_name=efel_feature_name,\n",
    "                    recording_names={'': '%s.soma.v' % protocol.name},\n",
    "                    stim_start=stim_start,\n",
    "                    stim_end=stim_end,\n",
    "                    exp_mean=mean,\n",
    "                    exp_std=0.1 * mean)\n",
    "        objective = ephys.objectives.SingletonObjective(feature_name, feature)\n",
    "        objectives.append(objective)\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Objective evaluator. Which calculates the quality of each results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "score_calc = ephys.objectivescalculators.ObjectivesCalculator(objectives) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cell Evaluator. The conductances that will be evaluated during the optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_evaluator = ephys.evaluators.CellEvaluator(\n",
    "        cell_model=GrCmulti,\n",
    "        param_names=['gcabar_GRC_CA_dend', \n",
    "              'gkbar_GRC_KCA_dend',     \n",
    "              'gl_GRC_LKG1_dend',\n",
    "              'ggaba_GRC_LKG2_dend',\n",
    "              'beta_GRC_CALC_dend', \n",
    "              'cai0_GRC_CALC_dend',\n",
    "              'gkbar_GRC_KIR_soma',\n",
    "              'gkbar_GRC_KA_soma',\n",
    "              'gl_GRC_LKG1_soma',\n",
    "              'gkbar_GRC_KM_soma',\n",
    "              'beta_GRC_CALC_soma',\n",
    "              'cai0_GRC_CALC_soma',\n",
    "              'gnabar_GRC_NA_axon0',\n",
    "              'gkbar_GRC_KV_axon0',\n",
    "              'gl_GRC_LKG1_axon0',\n",
    "              'gnabar_GRC_NA_axon1',\n",
    "              'gkbar_GRC_KV_axon1',\n",
    "              'gl_GRC_LKG1_axon1'],\n",
    "    \n",
    "        fitness_protocols={threestep_protocol.name: threestep_protocol},\n",
    "        fitness_calculator=score_calc,\n",
    "        sim=nrn)\n",
    "\n",
    "print(cell_evaluator.evaluate_with_dicts(default_params))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optimization. Number of offsprings and max generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "optimisation = bpop.optimisations.DEAPOptimisation(\n",
    "        evaluator=cell_evaluator,\n",
    "        offspring_size = NUMBER_INDIVIDUALS,\n",
    "        seed=os.getenv('BLUEPYOPT_SEED'))\n",
    "\n",
    "\n",
    "final_pop, hall_of_fame, logs, hist = optimisation.run(max_ngen=NUMBER_GENERATIONS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluation of the best individual based on its simulations, which are performed with the code below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#best individual\n",
    "best_ind = hall_of_fame[0]\n",
    "print( 'Best individual: ', best_ind)\n",
    "print( 'Fitness values: ', best_ind.fitness.values)\n",
    "\n",
    "best_ind_dict = cell_evaluator.param_dict(best_ind)\n",
    "print( cell_evaluator.evaluate_with_dicts(best_ind_dict))\n",
    "\n",
    "plot_responses(threestep_protocol.run(cell_model=GrCmulti, param_values=best_ind_dict, sim=nrn))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot of the best individual and the min fitness obtained during the optimization "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import numpy\n",
    "gen_numbers = logs.select('gen')\n",
    "min_fitness = logs.select('min')\n",
    "max_fitness = logs.select('max')\n",
    "plt.plot(gen_numbers, min_fitness, label='min fitness')\n",
    "plt.xlabel('generation #')\n",
    "plt.ylabel('score (# std)')\n",
    "plt.legend()\n",
    "plt.xlim(min(gen_numbers) - 1, max(gen_numbers) + 1) \n",
    "plt.ylim(0.9*min(min_fitness), 1.1 * max(min_fitness)) \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}