Mono 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 4 individuals for only 4 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 --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 = 4\n",
    "NUMBER_GENERATIONS = 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Installation of BluePyOpt\n",
    "!pip install bluepyopt --upgrade >/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, V_init and the tables used in the MOD files too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import bluepyopt as bpop\n",
    "import bluepyopt.ephys as ephys\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "\n",
    "\n",
    "class GrcMonoMorphology(ephys.morphologies.NrnFileMorphology):\n",
    "  def __init__(self, morphology_path, do_replace_axon):\n",
    "    super(GrcMonoMorphology, self).__init__(morphology_path, do_replace_axon)\n",
    "    \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",
    "        #Temperature and v_init\n",
    "        sim.neuron.h.celsius = 30\n",
    "        sim.neuron.h.v_init = -80\n",
    "        \n",
    "        #Ionic channels tables\n",
    "        sim.neuron.h.usetable_GrG_Na = 0\n",
    "        sim.neuron.h.usetable_GrC_pNa = 0\n",
    "        sim.neuron.h.usetable_GrC_CaHVA = 0\n",
    "        sim.neuron.h.usetable_GrG_KV = 0\n",
    "        sim.neuron.h.usetable_GrC_KA = 0\n",
    "        sim.neuron.h.usetable_GrC_Kir = 0\n",
    "        sim.neuron.h.usetable_GrC_KCa = 0\n",
    "        sim.neuron.h.usetable_GrG_KM = 0\n",
    "        \n",
    "        #sim.neuron.h.nrn_load_dll('$/opt/app-root/src/tost/x86_64/')\n"
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "!unzip -o granule_mono.zip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This command will compile the `MOD` files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!nrnivmodl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Defining passive properties\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "morph = GrcMonoMorphology('GrCmorphmono.asc', do_replace_axon = True) \n",
    "somatic_loc = ephys.locations.NrnSeclistLocation(\\\n",
    "    'somatic', seclist_name='somatic')"
   ]
  },
  {
   "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. \n",
    "\n",
    "* `name`: Name of the property, which must be unique\n",
    "* `param_name`: Name as taken from the `mod` file\n",
    "* `value`: The value of the property\n",
    "* `locations`: The location of where this property is being used, for example soma (somatic_loc).\n",
    "* `frozen`: Indicating this property has a fixed value (`frozen=true`) or is subject to optimization (`frozen=false`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Passive properties \n",
    "cm_param_soma = ephys.parameters.NrnSectionParameter(\n",
    "        name='cm_soma',\n",
    "        param_name='cm',\n",
    "        value=1.0,\n",
    "        locations=[somatic_loc],\n",
    "        frozen=True)\n",
    "\n",
    "ena = ephys.parameters.NrnSectionParameter(\n",
    "        name='ena_soma',\n",
    "        param_name='ena',\n",
    "        value=87.39,\n",
    "        locations=[somatic_loc],\n",
    "        frozen=True)\n",
    "\n",
    "ek = ephys.parameters.NrnSectionParameter(\n",
    "        name='ek_soma',\n",
    "        param_name='ek',\n",
    "        value=-84.69,\n",
    "        locations=[somatic_loc],\n",
    "        frozen=True)\n",
    "\n",
    "el = ephys.parameters.NrnSectionParameter(                            \n",
    "        name='el_soma',\n",
    "        param_name='el_GrC_Lkg1',\n",
    "        value=-58,\n",
    "        locations=[somatic_loc],\n",
    "        frozen=True)\n",
    "\n",
    "egaba = ephys.parameters.NrnSectionParameter(                         \n",
    "        name='egaba_soma',\n",
    "        param_name='egaba_GrC_Lkg2',\n",
    "        value=-65,\n",
    "        locations=[somatic_loc],\n",
    "        frozen=True)\n",
    "\n",
    "eca = ephys.parameters.NrnSectionParameter(\n",
    "        name='eca_soma',\n",
    "        param_name='eca',\n",
    "        value=129.33,\n",
    "        locations=[somatic_loc],\n",
    "        frozen=True)\n",
    "\n",
    "Ra = ephys.parameters.NrnSectionParameter(\n",
    "        name='Ra_soma',\n",
    "        param_name='Ra',\n",
    "        value=100,\n",
    "        locations=[somatic_loc],\n",
    "        frozen=True)\n",
    "\n",
    "#ionic mechanisms\n",
    "Mn_Na = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Na_mono',\n",
    "        prefix='GrG_Na',\n",
    "        locations=[somatic_loc])\n",
    "\n",
    "Mn_Nar = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Nar_mono',\n",
    "        prefix='GrG_Nar',\n",
    "        locations=[somatic_loc])\n",
    "\n",
    "Mn_Nap = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Nap_mono',\n",
    "        prefix='GrC_pNa',\n",
    "        locations=[somatic_loc])\n",
    "\n",
    "Mn_leak1 = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Leak1_mono',\n",
    "        prefix='GrC_Lkg1',\n",
    "        locations=[somatic_loc])\n",
    "\n",
    "Mn_leak2 = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Leak2_mono',\n",
    "        prefix='GrC_Lkg2',\n",
    "        locations=[somatic_loc])\n",
    "\n",
    "Mn_Ca = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Ca_mono',\n",
    "        prefix='GrC_CaHVA',\n",
    "        locations=[somatic_loc])\n",
    "\n",
    "Mn_Calc = ephys.mechanisms.NrnMODMechanism(      \n",
    "        name='Calc_mono',\n",
    "        prefix='Calc',\n",
    "        locations=[somatic_loc])\n",
    "\n",
    "Mn_Kir = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Kir_mono',\n",
    "        prefix='GrC_Kir',\n",
    "        locations=[somatic_loc])\n",
    "\n",
    "Mn_Kv = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Kv_mono',\n",
    "        prefix='GrG_KV',\n",
    "        locations=[somatic_loc])\n",
    "\n",
    "Mn_Kca = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Kca_mono',\n",
    "        prefix='GrC_KCa',\n",
    "        locations=[somatic_loc])\n",
    "\n",
    "Mn_Ka = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Ka_mono',\n",
    "        prefix='GrC_KA',\n",
    "        locations=[somatic_loc])\n",
    "\n",
    "Mn_km = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='km_mono',\n",
    "        prefix='GrG_KM',\n",
    "        locations=[somatic_loc])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the section in which the conductance range for each ionic channel is defined. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#parameters\n",
    "\n",
    "Mn_Na_param = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gnabar_GrG_Na_mono',\n",
    "        param_name='gnabar_GrG_Na',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[0.0104, 0.0156],\n",
    "        frozen=False)\n",
    "\n",
    "Mn_Nar_param = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gnabar_GrG_Nar_mono',\n",
    "        param_name='gnabar_GrG_Nar',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[0.0004, 0.0006],\n",
    "        frozen=False)\n",
    "\n",
    "Mn_Nap_param = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gnabar_GrC_pNa_mono',\n",
    "        param_name='gnabar_GrC_pNa',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[1.60e-5, 2.40e-5],\n",
    "        frozen=False)\n",
    "\n",
    "Mn_leak1_param = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gl_GrC_Lkg1_mono',\n",
    "        param_name='gl_GrC_Lkg1',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[4.54e-5, 6.82e-5],\n",
    "        frozen=False)\n",
    "\n",
    "Mn_leak2_param = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='egaba_GrC_Lkg2_mono',\n",
    "        param_name='egaba_GrC_Lkg2',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[1.74e-5, 2.60e-5],\n",
    "        frozen=False)\n",
    "\n",
    "Mn_Ca_param = ephys.parameters.NrnSectionParameter(\n",
    "        name='gcabar_GrC_CaHVA_mono',\n",
    "        param_name='gcabar_GrC_CaHVA',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[0.000368, 0.000552],\n",
    "        frozen=False)\n",
    "\n",
    "Mn_Calc_param = ephys.parameters.NrnSectionParameter(           \n",
    "        name='beta_Calc_mono',\n",
    "        param_name='beta_Calc',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[1.5, 1.51],\n",
    "        frozen=False)\n",
    "\n",
    "Mn_Kir_param = ephys.parameters.NrnSectionParameter(\n",
    "        name='gkbar_GrC_Kir_mono',\n",
    "        param_name='gkbar_GrC_Kir',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[0.00072, 0.00108],\n",
    "        frozen=False)\n",
    "\n",
    "Mn_Kv_param = ephys.parameters.NrnSectionParameter(\n",
    "        name='gkbar_GrG_KV_mono',\n",
    "        param_name='gkbar_GrG_KV',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[0.0024, 0.0036],\n",
    "        frozen=False)\n",
    "\n",
    "Mn_Kca_param = ephys.parameters.NrnSectionParameter(\n",
    "        name='gkbar_GrC_KCa_mono',\n",
    "        param_name='gkbar_GrC_KCa',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[0.0032, 0.0048],\n",
    "        frozen=False)\n",
    "\n",
    "Mn_Ka_param = ephys.parameters.NrnSectionParameter(\n",
    "        name='gkbar_GrC_KA_mono',\n",
    "        param_name='gkbar_GrC_KA',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[0.0032, 0.0048],\n",
    "        frozen=False)\n",
    "\n",
    "Mn_km_param = ephys.parameters.NrnSectionParameter(\n",
    "        name='gkbar_GrG_KM_mono',\n",
    "        param_name='gkbar_GrG_KM',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[0.00028, 0.00042],\n",
    "        frozen=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Creating the template\n",
    "\n",
    "To create the cell template, we pass all the previously prepared objects 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": [
    "grcmono = ephys.models.CellModel(                     \n",
    "        name='grcmono',\n",
    "        morph=morph,\n",
    "        mechs=[Mn_Na,\n",
    "               Mn_Nar, \n",
    "               Mn_Nap, \n",
    "               Mn_leak1, \n",
    "               Mn_leak2, \n",
    "               Mn_Ca, \n",
    "               Mn_Calc, \n",
    "               Mn_Kir, \n",
    "               Mn_Kv, \n",
    "               Mn_Kca, \n",
    "               Mn_Ka, \n",
    "               Mn_km],\n",
    "\n",
    "        params=[cm_param_soma, \n",
    "                ena, \n",
    "                ek, \n",
    "                el,\n",
    "                egaba,   \n",
    "                eca, \n",
    "                Ra, \n",
    "                Mn_Na_param, \n",
    "                Mn_Nar_param, \n",
    "                Mn_Nap_param, \n",
    "                Mn_leak1_param , \n",
    "                Mn_leak2_param , \n",
    "                Mn_Ca_param, \n",
    "                Mn_Calc_param, \n",
    "                Mn_Kir_param, \n",
    "                Mn_Kv_param, \n",
    "                Mn_Kca_param, \n",
    "                Mn_Ka_param, \n",
    "                Mn_km_param]) \n",
    "        \n",
    "print(grcmono)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Placing 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": [
    "#### Defining the stimuli, recordings location and protocols. \n",
    "\n",
    "For each protocol there is a recording location and a stimulus; in this example both are placed on the soma.\n",
    "\n",
    "The delay, duration and total duration can be changed, but a longer simumation time also 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 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 respond to be declared valid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "default_params = {'gnabar_GrG_Na_mono': 0.013, \n",
    "                  'gnabar_GrG_Nar_mono': 0.0005, \n",
    "                  'gnabar_GrC_pNa_mono': 2e-5, \n",
    "                  'gl_GrC_Lkg1_mono': 5.68e-5,\n",
    "                  'egaba_GrC_Lkg2_mono': 2.17e-5,\n",
    "                  'gcabar_GrC_CaHVA_mono': 0.00046,\n",
    "                  'beta_Calc_mono': 1.5,\n",
    "                  'gkbar_GrC_Kir_mono': 0.0009,\n",
    "                  'gkbar_GrG_KV_mono': 0.003,\n",
    "                  'gkbar_GrC_KCa_mono': 0.004,\n",
    "                  'gkbar_GrC_KA_mono': 0.004,\n",
    "                  'gkbar_GrG_KM_mono': 0.00035}\n",
    "\n",
    "responses = threestep_protocol.run(cell_model=grcmono, param_values=default_params, sim=nrn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plotting the response traces\n",
    "\n",
    "Three 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()\n"
   ]
  },
  {
   "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",
    "                                'mean_frequency': 55.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",
    "                                'mean_frequency': 84.36}, \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",
    "                                'mean_frequency': 109.8}}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Objectives\n",
    "\n",
    "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.05 * mean)\n",
    "        objective = ephys.objectives.SingletonObjective(\n",
    "            feature_name,\n",
    "            feature)\n",
    "        objectives.append(objective)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Objective evaluator\n",
    "\n",
    "Calculating the quality of each results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "        \n",
    "score_calc = ephys.objectivescalculators.ObjectivesCalculator(objectives)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cell Evaluator\n",
    "\n",
    "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=grcmono,\n",
    "        param_names=['gnabar_GrG_Na_mono',\n",
    "                     'gnabar_GrG_Nar_mono',\n",
    "                     'gnabar_GrC_pNa_mono', \n",
    "                     'gl_GrC_Lkg1_mono',\n",
    "                     'egaba_GrC_Lkg2_mono',\n",
    "                     'gcabar_GrC_CaHVA_mono',\n",
    "                     'beta_Calc_mono',\n",
    "                     'gkbar_GrC_Kir_mono',\n",
    "                     'gkbar_GrG_KV_mono',\n",
    "                     'gkbar_GrC_KCa_mono',\n",
    "                     'gkbar_GrC_KA_mono',\n",
    "                     'gkbar_GrG_KM_mono'],\n",
    "        fitness_protocols={threestep_protocol.name: threestep_protocol},\n",
    "        fitness_calculator=score_calc,\n",
    "        sim=nrn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Optimization\n",
    "\n",
    "Number of offsprings and max generation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(cell_evaluator.evaluate_with_dicts(default_params))\n",
    "\n",
    "\n",
    "optimisation = bpop.optimisations.DEAPOptimisation(\n",
    "    evaluator=cell_evaluator,\n",
    "    offspring_size = NUMBER_INDIVIDUALS,\n",
    "    seed=os.getenv('BLUEPYOPT_SEED'))\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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Final population:', final_pop)\n",
    "import numpy as np\n",
    "\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=grcmono, param_values=best_ind_dict, sim=nrn))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot of the best individual and the minimum 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()\n",
    " "
   ]
  },
  {
   "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
}