Purkinje cell simulation and validation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mouse Purkinje cell model\n",
    "\n",
    "This use case shows the steps to select and validate the results obtained after an otimization of a multicomparmental model of mouse Purkinje cell.\n",
    "\n",
    "- The morphology was taken from neuromorpho.org \n",
    "\n",
    ">  (http://neuromorpho.org/neuron_info.jsp?neuron_name=Purkinje-slice-ageP43-6)\n",
    "\n",
    "\n",
    "- The passive and active properties were taken from the published Purkinje cell model which can be found on modeldb\n",
    "\n",
    ">  (https://senselab.med.yale.edu/modeldb/ShowModel.cshtml?model=229585)\n",
    "\n",
    "> Masoli, S., Solinas, S., & D’Angelo, E. (2015). Action potential processing in a detailed Purkinje cell model reveals a critical role for axonal compartmentalization. Frontiers in Cellular Neuroscience, 9(February), 1–22. http://doi.org/10.3389/fncel.2015.00047\n",
    "\n",
    "- The range of the conductances, to define the range for the optimization, were adapted from:\n",
    "\n",
    "> Masoli, S., Solinas, S., & D’Angelo, E. (2015). Action potential processing in a detailed Purkinje cell model reveals a critical role for axonal compartmentalization. Frontiers in Cellular Neuroscience, 9(February), 1–22. http://doi.org/10.3389/fncel.2015.00047\n",
    "\n",
    "> Masoli, S., & D’Angelo, E. (2017). Synaptic Activation of a Detailed Purkinje Cell Model Predicts Voltage-Dependent Control of Burst-Pause Responses in Active Dendrites. Frontiers in Cellular Neuroscience, 11(September), 1–18. http://doi.org/10.3389/fncel.2017.00278\n",
    "\n",
    "- The validation is based on mouse cell experimental data provided by inlab (SP1) and literature:\n",
    "\n",
    "> Ransdell, J. L., Dranoff, E., Lau, B., Lo, W. L., Donermeyer, D. L., Allen, P. M., & Nerbonne, J. M. (2017). Erratum: Loss of Navβ4-Mediated Regulation of Sodium Currents in Adult Purkinje Neurons Disrupts Firing and Impairs Motor Coordination and Balance (Cell Reports (2017) 19(3) (532–544) (S2211124717304291) (10.1016/j.celrep.2017.03.068)). Cell Reports, 20(6), 1502. http://doi.org/10.1016/j.celrep.2017.07.072"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Technical Details\n",
    "\n",
    "The optimization procedure requires a lot of computational power and the available supercomputer (Julia - Julich), accesible from the Unicore system, is not strong enough to do the job in a reasonable amount of time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to run\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 notebook will be executed till the end.\n",
    "- With the second option (\"Run Cells\"), only the code of a single cell will be executed. Clicking the following symbol ![skip-next](https://raw.githubusercontent.com/antonelepfl/testvue/master/skip_next.png) will execute the notebook cell-by-cell, giving the possibility to interfere with the configuration in each cell.\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",
    "More details can be taken from the [Guidebook](https://lbologna.github.io/hbp-sp6-guidebook/)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparation\n",
    "The following cells will setup the access to the Collab Storage and download the required morphologies and ion channels."
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "This cell is used to prevent interference from previous runs of this notebooks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "!rm *step* \n",
    "!rm *indiv*\n",
    "!rm -r valid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Downloading morphology, ionic channels, previosly optimized results to the location where this notebooks runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2021-04-21 10:01:58--  https://object.cscs.ch/v1/AUTH_c0a333ecf7c045809321ce9d9ecdfdea/cerebellum_circuits/Purkinje_usecase/Purkinje_res.zip\n",
      "Resolving object.cscs.ch (object.cscs.ch)... 148.187.25.200, 148.187.25.201, 148.187.25.202\n",
      "Connecting to object.cscs.ch (object.cscs.ch)|148.187.25.200|:443... connected.\n",
      "HTTP request sent, awaiting response... 200 OK\n",
      "Length: 108516299 (103M) [application/zip]\n",
      "Saving to: ‘Purkinje_res.zip’\n",
      "\n",
      "100%[======================================>] 108,516,299 36.1MB/s   in 2.9s   \n",
      "\n",
      "2021-04-21 10:02:01 (36.1 MB/s) - ‘Purkinje_res.zip’ saved [108516299/108516299]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "!wget https://object.cscs.ch:443/v1/AUTH_c0a333ecf7c045809321ce9d9ecdfdea/cerebellum_circuits/Purkinje_usecase/Purkinje_res.zip"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "!unzip -o Purkinje_multi.zip\n",
    "!unzip -o Purkinje_res.zip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Modules used for the extraction of the features from voltage time traces and to plot them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import collections\n",
    "import numpy as np\n",
    "import efel\n",
    "import matplotlib.pyplot as plt\n",
    "import subprocess\n",
    "\n",
    "feature_values = collections.defaultdict(dict)\n",
    "conductvalues = np.genfromtxt(\"final_pop.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot of the traces belonging to a single individual\n",
    "\n",
    "After executing the next cell you can select the number of the sim to be plotted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ipywidgets as widgets\n",
    "import numpy as np\n",
    "conductvalues = np.genfromtxt(\"final_pop.txt\")\n",
    "\n",
    "\n",
    "sim_number = widgets.BoundedIntText(\n",
    "    value=0,\n",
    "    min=0,\n",
    "    max=len(conductvalues),\n",
    "    step=1,\n",
    "    description='Sim num:',\n",
    "    disabled=False\n",
    ")\n",
    "display(sim_number)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the selected sim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "# Location of the recording. For the Purkinje cell there is only one.\n",
    "location = '' \n",
    "\n",
    "steps = {\n",
    "    1: \"spontaneuos firing\",\n",
    "    2: \"0.5nA current injection\",\n",
    "    3: \"1nA current injection\",\n",
    "}\n",
    "\n",
    "for current_inj_step in range (1,4): #Number of the current step used in the optimization. 1, first step, 2 second step and so on.\n",
    "    filetxt = str(sim_number.value) + '_sim_step' + str(current_inj_step) + location +'.txt'\n",
    "    if os.path.isfile(filetxt) == False:\n",
    "        print( 'file not found')\n",
    "        pass\n",
    "    \n",
    "    elif os.path.isfile(filetxt):    \n",
    "        trace = np.genfromtxt(filetxt)\n",
    "        \n",
    "        time, voltage = zip(*trace)\n",
    "        \n",
    "        fig, ax = plt.subplots()\n",
    "        ax.plot(time,voltage, 'b', label='spikes, %s' % (steps[current_inj_step]))\n",
    "\n",
    "        legend = ax.legend(loc='upper right', shadow=True)\n",
    "\n",
    "        frame = legend.get_frame()\n",
    "        frame.set_facecolor('0.90')\n",
    "\n",
    "        plt.xlabel(\"Time (ms)\")\n",
    "        plt.ylabel(\"Membrane voltage (mv) \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining some functions to automatically select, analyze and save the results from the voltage traces.\n",
    "The selection process is based on the frequency, calculated with the ISI values, for three steps:\n",
    "\n",
    "- Spontaneuos firing\n",
    "- 0.5nA positive current injection\n",
    "- 1nA positive current injection\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trace_features(trace_init, trace_end, indiv, stepindex, location): \n",
    "    \n",
    "    filename = str(indiv) + '_sim_step' + str(stepindex) + location + '.txt'\n",
    "    print( \"trace number\", indiv)\n",
    "    \n",
    "    print( filename)\n",
    "    fileeps = str(indiv) + '_sim_step' + str(stepindex) + '_window' + location +'.eps'\n",
    "    if os.path.isfile(filename) and os.path.isfile(fileeps):\n",
    "        print( filename, 'trovato')\n",
    "        pass\n",
    "    \n",
    "    elif os.path.isfile(filename):\n",
    "        trace = np.genfromtxt(filename)\n",
    "\n",
    "        time = []\n",
    "        voltage = []\n",
    "\n",
    "        for x in range(len(trace)):\n",
    "            time.append(trace[x,0])\n",
    "            voltage.append(trace[x,1])\n",
    "        \n",
    "        fig, ax = plt.subplots()\n",
    "        ax.plot(time, voltage, 'b', label='spikes')\n",
    "\n",
    "        legend = ax.legend(loc='upper right', shadow=True)\n",
    "\n",
    "        frame = legend.get_frame()\n",
    "        frame.set_facecolor('0.90')\n",
    "\n",
    "        plt.xlabel(\"time (ms)\")\n",
    "        plt.ylabel(\"membrane voltage (mv) \")\n",
    "\n",
    "        plt.savefig('feat_indiv_' + str(indiv) + '_sim_step' + str(stepindex) +  '_start_'+ str(trace_init) + '_end_'+ str(trace_end) + '_window' + location + '.eps')\n",
    "        plt.close()\n",
    "        \n",
    "        traces = []\n",
    "\n",
    "        trace = {}  \n",
    "        trace['T'] = time\n",
    "        trace['V'] = voltage\n",
    "        trace['stim_start'] = [trace_init]\n",
    "        trace['stim_end'] = [trace_end]\n",
    "        traces = [trace]  \n",
    "\n",
    "        result = efel.getFeatureValues(traces, ['mean_frequency',\n",
    "                                                'Spikecount',\n",
    "                                                'ISI_CV'])\n",
    "    \n",
    "        for feature_name, feature_list in result[0].items():\n",
    "            if feature_list is not None and len(feature_list) > 0:\n",
    "                feature_values[feature_name] = np.mean(feature_list)\n",
    "            else:\n",
    "                feature_values[feature_name] = 0\n",
    "\n",
    "                \n",
    "        filename = 'feat_indiv_' + str(indiv) + '_sim_step' + str(stepindex) + '_start_'+ str(trace_init) + '_end_'+ str(trace_end) + '_window' + location + '.csv'\n",
    "        fo = open(filename, \"w\")\n",
    "        for e,v in feature_values.items():\n",
    "            fo.write(e+ ' ' + str(v)+\"\\n\")    \n",
    "        fo.close() \n",
    "        \n",
    "def save_feat_name(trace_init, trace_end, indiv, stepindex, location):\n",
    "    global final_name_matrix\n",
    "    filename = 'feat_indiv_' + str(indiv) + '_sim_step' + str(stepindex) + '_start_'+ str(trace_init) + '_end_'+ str(trace_end) + '_window' + location + '.csv'\n",
    "    trace = np.genfromtxt(filename, dtype = [('feat_name','S20'),('value','f8')])\n",
    "\n",
    "    name = []\n",
    "    for x in trace:\n",
    "        name.append(x['feat_name'])\n",
    "\n",
    "    trace_name = np.matrix(name)\n",
    "    add_number_gen = ['Indiv_' + str(indiv)]\n",
    "    trace_indiv = np.matrix(add_number_gen)\n",
    "    x_size,y_size = trace_name.shape\n",
    "    trace_name_mod = np.reshape(trace_name,(y_size ,1)) \n",
    "    final_name_matrix = np.append(trace_indiv, trace_name_mod, axis=0)  \n",
    "\n",
    "def matrix_builder(filename, stepindex,location):\n",
    "    global final_value_matrix  \n",
    "    trace = np.genfromtxt(filename, dtype = [('feat_name','S20'),('value','f8')])\n",
    "    \n",
    "    name = []\n",
    "    value = []\n",
    "    for x in trace:\n",
    "        value.append(x['value'])\n",
    "\n",
    "#value of each feature\n",
    "    trace_value = np.matrix(value)\n",
    "    \n",
    "    if location == '':\n",
    "        location = '_soma'\n",
    "    else:\n",
    "        pass\n",
    "    \n",
    "    stepindex_ext = [location + '_step_' +  str(stepindex)]\n",
    "    \n",
    "    trace_step = np.matrix(stepindex_ext)\n",
    "    \n",
    "    x_size,y_size = trace_value.shape\n",
    "    trace_value_mod = np.reshape(trace_value,(y_size ,1))\n",
    "\n",
    "    final_value_matrix = np.append(trace_step, trace_value_mod, axis=0)    \n",
    "\n",
    "\n",
    "def trace_load(trace_init, trace_end, indiv, stepindex, location): \n",
    "    global base_matrix_value\n",
    "    filename = 'feat_indiv_' + str(indiv) + '_sim_step' + str(stepindex) + '_start_'+ str(trace_init) + '_end_'+ str(trace_end) + '_window' + location + '.csv'\n",
    "    if os.path.isfile(filename) == True:\n",
    "        matrix_builder(filename,stepindex, location)\n",
    "        base_matrix_value = np.append(base_matrix_value, final_value_matrix, axis=1)\n",
    "        base_matrix_value[0,0] =  'Indiv_' + str(indiv)\n",
    "        print( 'done')\n",
    "        np.savetxt('Feature'+ '_indiv_' + str(indiv) + '_start_'+ str(trace_init) + '_end_'+ str(trace_end) + '.csv', base_matrix_value, delimiter = ',', fmt=\"%10s\")\n",
    "    \n",
    "    if os.path.isfile(filename) == False:\n",
    "        print( filename, 'not found')\n",
    "        pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def selector(indiv):\n",
    "    threshold = {\n",
    "        1: 0.05,\n",
    "        2: 0.3,\n",
    "        3: 0.36\n",
    "    }\n",
    "    \n",
    "    \n",
    "    if base_matrix_value.shape <= (3,3):\n",
    "        print( 'matrix not valid')\n",
    "    else:     \n",
    "            \n",
    "        realy_valid = []\n",
    "        realy_bad = []\n",
    "                         \n",
    "        for i in range(1,4):   \n",
    "            \n",
    "            ISI = float(base_matrix_value[3,i])\n",
    "            \n",
    "            if ISI <= threshold[i]:\n",
    "                realy_valid.append(i)\n",
    "            if ISI >= threshold[i]:\n",
    "                realy_bad.append(i)\n",
    "                \n",
    "\n",
    "        import os    \n",
    "       \n",
    "        if len(realy_valid) == 3:\n",
    "            print( 'Physiological Purkinje cell')\n",
    "            print( 'indiv', indiv)\n",
    "            \n",
    "            if os.path.isdir(\"valid/\") == False:\n",
    "                os.makedirs(\"valid/\")\n",
    "                \n",
    "                actualdir = \"valid/\"\n",
    "                indiv_number = str(indiv) + '_'\n",
    "                indiv_number_feat_single = \"feat_indiv_\" + str(indiv) + '_'\n",
    "                indiv_number_feat_matrix = \"Feature_indiv_\" + str(indiv)  + '_'\n",
    "                \n",
    "                if indiv == 0:\n",
    "                    pass\n",
    "                else:\n",
    "                    subprocess.call(\"mv \" + indiv_number + \"* \" + actualdir, shell=True)\n",
    "                    subprocess.call(\"mv \" + indiv_number_feat_single + \"* \" + actualdir, shell=True)\n",
    "                    subprocess.call(\"mv \" + indiv_number_feat_matrix + \"* \" + actualdir, shell=True)\n",
    "                \n",
    "            elif os.path.isdir(\"valid/\") == True:\n",
    "                print( 'dir exist!', \"valid/\")\n",
    "                print( 'indiv', indiv)\n",
    "                \n",
    "                actualdir = \"valid/\"\n",
    "                indiv_number = str(indiv) + '_'\n",
    "                indiv_number_feat_single = \"feat_indiv_\" + str(indiv) + '_'\n",
    "                indiv_number_feat_matrix = \"Feature_indiv_\" + str(indiv) + '_'\n",
    "                                \n",
    "                if indiv == 0:\n",
    "                    pass\n",
    "                else:\n",
    "                    subprocess.call(\"mv \" + indiv_number + \"* \" + actualdir, shell=True)\n",
    "                    subprocess.call(\"mv \" + indiv_number_feat_single + \"* \" + actualdir, shell=True)\n",
    "                    subprocess.call(\"mv \" + indiv_number_feat_matrix + \"* \" + actualdir, shell=True)     \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mkdir -p valid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell starts the process to:\n",
    "- extract of the features from the voltage-traces generated from the last iteration of the optimization.\n",
    "- Each individual properties are save together in a series of files.\n",
    "- Each individual is validated against the ISI values obtained from SP1 data and literature (Ransdell and Narbonne, 2017).\n",
    "- The individual that are within the parameters are moved to the \"good\" directory. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soma = ''\n",
    "\n",
    "for first_step in range(1,4):\n",
    "    trace_features(0,2100,0,first_step,soma)\n",
    "    \n",
    "save_feat_name(0,2100,0,1,soma)\n",
    "\n",
    "base_matrix_value = final_name_matrix\n",
    "\n",
    "for indiv in range(0,len(conductvalues)):\n",
    "    for step in range(1,4): #number of steps used during the optimization procedure.\n",
    "            trace_features(0,2100,indiv,step,soma)  \n",
    "            trace_load(0,2100,indiv,step,soma)  \n",
    "    selector(indiv)\n",
    "    base_matrix_value = final_name_matrix        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This cell reads the valid results and show them as a series of numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "valid_dir = \"valid/\"\n",
    "good_results = []\n",
    "for i in range(len(conductvalues)):\n",
    "    filename_final = valid_dir + str(i) + '_sim_step1' + '.txt'\n",
    "    fileeps_final = valid_dir + str(i) + '_sim_step1' +'.eps'\n",
    "    if os.path.isfile(filename_final) == True:\n",
    "        good_results.append(i)\n",
    "        \n",
    "print( good_results        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compiling the mod files for each ionic channel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!nrnivmodl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following two cells contain the same code, used during the optimization, modified to simulate the valid results.\n",
    "\n",
    "The following class (evoked by \"do_replace_axon = true\" in the following cell) overwrite the one provided with BluePyOpt and permits to define the properties and quantity of the sections composing the axon. This is due the fact that Purkinje cells have a myelinated axon composed by Ranvier Nodes and myelinated sections."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import bluepyopt as bpop\n",
    "import bluepyopt.ephys as ephys\t\n",
    "import math\n",
    "\n",
    "class nrnPCmorphology(ephys.morphologies.NrnFileMorphology):\n",
    "  def __init__(self, morphology_path, do_replace_axon):\n",
    "    super(nrnPCmorphology, self).__init__(morphology_path, do_replace_axon)\n",
    "\n",
    "  @staticmethod\n",
    "  def replace_axon(sim=None, icell=None):\n",
    "        \"\"\"Replace axon\"\"\"\n",
    "        ais_diams = [icell.axon[0].diam, icell.axon[0].diam]\n",
    "\n",
    "        # Define origin of distance function\n",
    "        sim.neuron.h.distance(sec=icell.soma[0])\n",
    "        for section in icell.axonal:\n",
    "            # If distance to soma is larger than 60, store diameter\n",
    "            if sim.neuron.h.distance(0.5, sec=section) > 60:\n",
    "                ais_diams[1] = section.diam\n",
    "                break\n",
    "\n",
    "        for section in icell.axonal:\n",
    "            sim.neuron.h.delete_section(sec=section)\n",
    "\n",
    "        # Create new axon array\n",
    "        sim.neuron.h.execute('create axon[11]', icell)\n",
    "\n",
    "        icell.axon[0].nseg = 1\n",
    "        icell.axon[0].L = 17\n",
    "        icell.axon[0].diam = 0.97\n",
    "        \n",
    "        icell.axon[1].nseg = 1\n",
    "        icell.axon[1].L = 4\n",
    "        icell.axon[1].diam = 0.97\n",
    "        \n",
    "        icell.axon[2].nseg = 1\n",
    "        icell.axon[2].L = 100\n",
    "        icell.axon[2].diam = 0.73\n",
    "        \n",
    "        icell.axon[3].nseg = 1\n",
    "        icell.axon[3].L = 4\n",
    "        icell.axon[3].diam = 0.73\n",
    "        \n",
    "        icell.axon[4].nseg = 1\n",
    "        icell.axon[4].L = 100\n",
    "        icell.axon[4].diam = 0.73\n",
    "        \n",
    "        icell.axon[5].nseg = 1\n",
    "        icell.axon[5].L = 4\n",
    "        icell.axon[5].diam = 0.73\n",
    "        \n",
    "        icell.axon[6].nseg = 1\n",
    "        icell.axon[6].L = 100\n",
    "        icell.axon[6].diam = 0.73\n",
    "        \n",
    "        icell.axon[7].nseg = 1\n",
    "        icell.axon[7].L = 4\n",
    "        icell.axon[7].diam = 0.73\n",
    "        \n",
    "        icell.axon[8].nseg = 1\n",
    "        icell.axon[8].L = 100\n",
    "        icell.axon[8].diam = 0.73\n",
    "        \n",
    "        for i in range(0,11):\n",
    "            icell.axonal.append(sec = icell.axon[i])\n",
    "            icell.all.append(sec = icell.axon[i])\n",
    "\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",
    "        icell.axon[2].connect(icell.axon[1], 1.0, 0.0)\n",
    "        icell.axon[3].connect(icell.axon[2], 1.0, 0.0)\n",
    "        icell.axon[4].connect(icell.axon[3], 1.0, 0.0)\n",
    "        icell.axon[5].connect(icell.axon[4], 1.0, 0.0)\n",
    "        icell.axon[6].connect(icell.axon[5], 1.0, 0.0)\n",
    "        icell.axon[7].connect(icell.axon[6], 1.0, 0.0)\n",
    "        icell.axon[8].connect(icell.axon[7], 1.0, 0.0)\n",
    "\n",
    "\n",
    "        icell.AIS.append(sec = icell.axon[0])\n",
    "        icell.AISK.append(sec = icell.axon[1]) \n",
    "        icell.myelin0.append(sec = icell.axon[2]) \n",
    "        icell.rn0.append(sec = icell.axon[3]) \n",
    "        icell.myelin1.append(sec = icell.axon[4]) \n",
    "        icell.rn1.append(sec = icell.axon[5]) \n",
    "        icell.myelin2.append(sec = icell.axon[6]) \n",
    "        icell.rn2.append(sec = icell.axon[7]) \n",
    "        icell.myelin3.append(sec = icell.axon[8]) \n",
    "\n",
    "\n",
    "        for section in icell.dend:\n",
    "            if section.diam >= 0 and section.diam <= 6:\n",
    "                icell.dendall.append(sec=section)\n",
    "            if section.diam >= 1.4 and section.diam <= 6:\n",
    "                icell.dendtrunk.append(sec=section)\n",
    "            if section.diam >= 2.05 and section.diam <= 6:  \n",
    "                icell.dendsodium.append(sec=section)\n",
    "            if section.diam >= 0 and section.diam <= 1.4:\n",
    "                section.cm = ((11.510294 * math.exp( - 1.376463 * section.diam) + 2.120503))*1.5\n",
    "       \n",
    "        sim.neuron.h.celsius = 37"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Purkinje cell are valid when, in the absence of sodium channel in the the AIS, the cell is unable to generate spontaneuos firing.\n",
    "\n",
    "The numbers printed in the output of this cell refers to the simulation of each results validated in the previous part of the usecase.\n",
    "\n",
    "***Warning***: **The next step can take up to 15m to be perfomed.**\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Purkinje cell optimization with BluePyOpt\n",
    "\n",
    "#import matplotlib.pyplot as plt\n",
    "import os\n",
    "import numpy as np\n",
    "\n",
    "#Morphology location\n",
    "morph = nrnPCmorphology('purkinje_mouse.swc', do_replace_axon = True)\n",
    "\n",
    "#Locations\n",
    "\n",
    "dend_all = ephys.locations.NrnSeclistLocation('dendall', seclist_name='dendall')\n",
    "dend_trunk = ephys.locations.NrnSeclistLocation('dendtrunk', seclist_name='dendtrunk')\n",
    "dend_sodium = ephys.locations.NrnSeclistLocation('dendsodium', seclist_name='dendsodium')\n",
    "\n",
    "somatic_loc = ephys.locations.NrnSeclistLocation('somatic', seclist_name='somatic')\n",
    "\n",
    "axon0 = ephys.locations.NrnSeclistLocation('AIS', seclist_name='AIS')\n",
    "axon1 = ephys.locations.NrnSeclistLocation('AISK', seclist_name='AISK')\n",
    "axon2 = ephys.locations.NrnSeclistLocation('myelin0', seclist_name='myelin0')\n",
    "axon3 = ephys.locations.NrnSeclistLocation('rn0', seclist_name='rn0')\n",
    "axon4 = ephys.locations.NrnSeclistLocation('myelin1', seclist_name='myelin1')\n",
    "axon5 = ephys.locations.NrnSeclistLocation('rn1', seclist_name='rn1')\n",
    "axon6 = ephys.locations.NrnSeclistLocation('myelin2', seclist_name='myelin2')\n",
    "axon7 = ephys.locations.NrnSeclistLocation('rn2', seclist_name='rn2')\n",
    "axon8 = ephys.locations.NrnSeclistLocation('myelin3', seclist_name='myelin3')\n",
    "\n",
    "#Channels list\n",
    "\n",
    "#DENDRITES.\n",
    "\n",
    "#Passive properties - Frozen for single parameters and not frozen for ranges of parameters ie channels.\n",
    "cm_dend_all = ephys.parameters.NrnSectionParameter(\n",
    "        name='cm_dend_all', #- must be unique\n",
    "        param_name = 'cm',\n",
    "        value = 1,\n",
    "        locations=[dend_trunk, dend_sodium],\n",
    "        frozen=True)\n",
    "\n",
    "cm_soma = ephys.parameters.NrnSectionParameter(\n",
    "        name='cm_soma', #- must be unique\n",
    "        param_name='cm',\n",
    "        value=0.77,\n",
    "        locations=[somatic_loc],\n",
    "        frozen=True)\n",
    "\n",
    "cm_rn_only = ephys.parameters.NrnSectionParameter(\n",
    "        name='cm_axon_rn', #- must be unique\n",
    "        param_name='cm',\n",
    "        value=0.77,\n",
    "        locations=[axon0, axon1, axon3, axon5, axon7],\n",
    "        frozen=True)\n",
    "\n",
    "cm_myelin_only = ephys.parameters.NrnSectionParameter(\n",
    "        name='cm_axon_myelin', #- must be unique\n",
    "        param_name='cm',\n",
    "        value=1.87e-11,\n",
    "        locations=[axon2, axon4, axon6, axon8],\n",
    "        frozen=True)\n",
    "\n",
    "\n",
    "ena = ephys.parameters.NrnSectionParameter(\n",
    "        name='ena_all', #- must be unique\n",
    "        param_name='ena',\n",
    "        value=80,\n",
    "        locations=[dend_sodium, somatic_loc, axon0, axon3, axon5, axon7],\n",
    "        frozen=True)\n",
    "\n",
    "ek = ephys.parameters.NrnSectionParameter(\n",
    "        name='ek_all', #- must be unique\n",
    "        param_name='ek',\n",
    "        value=-88,\n",
    "        locations=[dend_all, somatic_loc, axon0, axon1, axon3, axon5, axon7],\n",
    "        frozen=True)\n",
    "\n",
    "eleak = ephys.parameters.NrnSectionParameter(\n",
    "        name='eleak_all', #- must be unique\n",
    "        param_name='e_Leak',\n",
    "        value = -61,\n",
    "        locations=[dend_all, somatic_loc, axon0, axon1, axon3, axon5, axon7],\n",
    "        frozen=True)\n",
    "\n",
    "leak_soma = ephys.parameters.NrnSectionParameter(\n",
    "        name='leak_soma', #- must be unique\n",
    "        param_name='gmax_Leak',\n",
    "        value = 1E-3,\n",
    "        locations=[somatic_loc],\n",
    "        frozen=True)\n",
    "\n",
    "leak_rest = ephys.parameters.NrnSectionParameter(\n",
    "        name='leak_rest', #- must be unique\n",
    "        param_name='gmax_Leak',\n",
    "        value = 0.0003,\n",
    "        locations=[dend_all, axon0, axon1, axon3, axon5, axon7],\n",
    "        frozen=True)\n",
    "\n",
    "eca = ephys.parameters.NrnSectionParameter(\n",
    "        name='eca_all', #- must be unique\n",
    "        param_name='eca',\n",
    "        value = 137.5,\n",
    "        locations=[dend_all, somatic_loc, axon0, axon3, axon5, axon7],\n",
    "        frozen=True)\n",
    "\n",
    "Ra = ephys.parameters.NrnSectionParameter(\n",
    "        name='ra_all', #- must be unique\n",
    "        param_name='Ra',\n",
    "        value = 122,\n",
    "        locations=[dend_all, somatic_loc, axon0, axon1, axon2, axon3, axon4, axon5, axon6, axon7, axon8],\n",
    "        frozen=True)\n",
    "\n",
    "gpas_myelin_only = ephys.parameters.NrnSectionParameter(\n",
    "        name='g_pas_axon_myelin', #- must be unique\n",
    "        param_name='g_pas',\n",
    "        value=5e-8,\n",
    "        locations=[axon2, axon4, axon6, axon8],\n",
    "        frozen=True)\n",
    "\n",
    "epas_myelin_only = ephys.parameters.NrnSectionParameter(\n",
    "        name='e_pas_axon_myelin', #- must be unique\n",
    "        param_name='e_pas',\n",
    "        value=-63,\n",
    "        locations=[axon2, axon4, axon6, axon8],\n",
    "        frozen=True)\n",
    "\n",
    "ca_pump_dend = ephys.parameters.NrnSectionParameter(\n",
    "        name='TotalPump_dend', \n",
    "        param_name='TotalPump_cdp5',\n",
    "        value=5e-8,\n",
    "        locations=[dend_all],\n",
    "        frozen=True)\n",
    "\n",
    "ca_pump_soma = ephys.parameters.NrnSectionParameter(\n",
    "        name='TotalPump_soma', \n",
    "        param_name='TotalPump_cdp5',\n",
    "        value=2e-8,\n",
    "        locations=[somatic_loc, axon0],\n",
    "        frozen=True)\n",
    "\n",
    "ca_pump_axon = ephys.parameters.NrnSectionParameter(\n",
    "        name='TotalPump_axon', \n",
    "        param_name='TotalPump_cdp5',\n",
    "        value=5e-7,\n",
    "        locations=[axon3, axon5, axon7],\n",
    "        frozen=True)\n",
    "\n",
    "#CHANNELS\n",
    "#DEND\n",
    "\n",
    "#Mixed\n",
    "PC_Leak = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Leak_all',\n",
    "        prefix='Leak',\n",
    "        locations=[dend_all, somatic_loc, axon0, axon1, axon3, axon5, axon7])\n",
    "\n",
    "PC_Nav16 = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Nav1_6',\n",
    "        prefix='Nav1_6',\n",
    "        locations=[dend_sodium, somatic_loc, axon0, axon3, axon5, axon7])\n",
    "\n",
    "PC_Kv1_1 = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Kv1_1',\n",
    "        prefix='Kv1_1',\n",
    "        locations=[dend_all, somatic_loc, axon1])\n",
    "\n",
    "PC_Kv1_5 = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Kv1_5',\n",
    "        prefix='Kv1_5',\n",
    "        locations=[dend_all])\n",
    "\n",
    "PC_Kv3_3 = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Kv3_3',\n",
    "        prefix='Kv3_3',\n",
    "        locations=[dend_all])\n",
    "\n",
    "PC_Kv3_4 = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Kv3_4',\n",
    "        prefix='Kv3_4',\n",
    "        locations=[somatic_loc, axon0, axon3, axon5, axon7])\n",
    "\n",
    "PC_Kv4_3 = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Kv4_3',\n",
    "        prefix='Kv4_3',\n",
    "        locations=[dend_all])\n",
    "\n",
    "PC_Kir2_3 = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Kir2_3',\n",
    "        prefix='Kir2_3',\n",
    "        locations=[somatic_loc, dend_trunk])\n",
    "\n",
    "PC_Kca1_1 = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Kca1_1',\n",
    "        prefix='Kca1_1',\n",
    "        locations=[dend_all, somatic_loc])\n",
    "\n",
    "PC_Kca2_2 = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Kca2_2',\n",
    "        prefix='Kca2_2',\n",
    "        locations=[dend_trunk, somatic_loc])\n",
    "\n",
    "PC_Kca3_1 = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Kca3_1',\n",
    "        prefix='Kca3_1',\n",
    "        locations=[dend_trunk, somatic_loc])\n",
    "\n",
    "PC_Cav2_1 = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Cav2_1',\n",
    "        prefix='Cav2_1',\n",
    "        locations=[dend_all, somatic_loc, axon0, axon3, axon5, axon7])\n",
    "\n",
    "PC_Cav3_1 = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Cav3_1',\n",
    "        prefix='Cav3_1',\n",
    "        locations=[dend_trunk, somatic_loc, axon0, axon3, axon5, axon7])\n",
    "\n",
    "PC_Cav3_2 = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Cav3_2',\n",
    "        prefix='Cav3_2',\n",
    "        locations=[dend_trunk, somatic_loc])\n",
    "\n",
    "PC_Cav3_3 = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='Cav3_3',\n",
    "        prefix='Cav3_3',\n",
    "        locations=[dend_all, somatic_loc])\n",
    "\n",
    "PC_HCN1 = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='HCN1',\n",
    "        prefix='HCN1',\n",
    "        locations=[dend_all, somatic_loc])\n",
    "\n",
    "PC_cdp5 = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='cdp5',\n",
    "        prefix='cdp5',\n",
    "        locations=[dend_all, somatic_loc, axon0, axon3, axon5, axon7])\n",
    "\n",
    "PC_pass = ephys.mechanisms.NrnMODMechanism(\n",
    "        name='pas',\n",
    "        prefix='pas',\n",
    "        locations=[axon2, axon4, axon6, axon8])\n",
    "\n",
    "#CHANNELS PARAM\n",
    "\n",
    "#DEND_all\n",
    "PC_Kv1_1_param_dend_all = ephys.parameters.NrnSectionParameter(\n",
    "        name='gbar_Kv1_1_dend_all',\n",
    "        param_name='gbar_Kv1_1',\n",
    "        bounds=[0.001, 0.0015],\n",
    "        locations=[dend_all],\n",
    "        frozen=False)\n",
    "\n",
    "PC_Kv1_5_param_dend_all = ephys.parameters.NrnSectionParameter(\n",
    "        name='gKur_Kv1_5_dend_all',\n",
    "        param_name='gKur_Kv1_5',\n",
    "        bounds=[0.8e-4, 1.5e-4],\n",
    "        locations=[dend_all],\n",
    "        frozen=False)\n",
    "\n",
    "PC_Kv3_3_param_dend_all = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gbar_Kv3_3_dend_all',\n",
    "        param_name='gbar_Kv3_3',\n",
    "        locations=[dend_all],\n",
    "        bounds=[0.009, 0.015],\n",
    "        frozen=False) \n",
    "\n",
    "PC_Kv4_3_param_dend_all = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gkbar_Kv4_3_dend_all',\n",
    "        param_name='gkbar_Kv4_3',\n",
    "        locations=[dend_all],\n",
    "        bounds=[0.0008, 0.0015],\n",
    "        frozen=False) \n",
    "\n",
    "PC_Kca1_1_param_dend_all = ephys.parameters.NrnSectionParameter(\n",
    "        name='gbar_Kca1_1_dend_all',\n",
    "        param_name='gbar_Kca1_1',\n",
    "        bounds=[2e-2, 5e-2],\n",
    "        locations=[dend_all],\n",
    "        frozen=False)\n",
    "\n",
    "PC_Cav2_1_param_dend_all = ephys.parameters.NrnSectionParameter(\n",
    "        name='pcabar_Cav2_1_dend_all',\n",
    "        param_name='pcabar_Cav2_1',\n",
    "        bounds=[2.5e-3, 4.5e-3],\n",
    "        locations=[dend_all],\n",
    "        frozen=False)\n",
    "\n",
    "PC_Cav3_3_param_dend_all = ephys.parameters.NrnSectionParameter(\n",
    "        name='pcabar_Cav3_3_dend_all',\n",
    "        param_name='pcabar_Cav3_3',\n",
    "        bounds=[0.00008, 0.0002],\n",
    "        locations=[dend_all],\n",
    "        frozen=False)\n",
    "\n",
    "PC_HCN1_param_dend_all = ephys.parameters.NrnSectionParameter(\n",
    "        name='gbar_HCN1_dend_all',\n",
    "        param_name='gbar_HCN1',\n",
    "        bounds=[0.000002, 0.000005],\n",
    "        locations=[dend_all],\n",
    "        frozen=False)\n",
    "\n",
    "\n",
    "#dend_trunk\n",
    "PC_Kir2_3_param_dend_trunk = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gkbar_Kir2_3_dend_trunk',\n",
    "        param_name='gkbar_Kir2_3',\n",
    "        locations=[dend_trunk],\n",
    "        bounds=[0.000008, 0.00002],\n",
    "        frozen=False)\n",
    "\n",
    "PC_Kca2_2_param_dend_trunk = ephys.parameters.NrnSectionParameter(\n",
    "        name='gkbar_Kca2_2_dend_trunk',\n",
    "        param_name='gkbar_Kca2_2',\n",
    "        bounds=[7e-6, 2e-3],\n",
    "        locations=[dend_trunk],\n",
    "        frozen=False)\n",
    "\n",
    "PC_Kca3_1_param_dend_trunk = ephys.parameters.NrnSectionParameter(\n",
    "        name='gkbar_Kca3_1_dend_trunk',\n",
    "        param_name='gkbar_Kca3_1',\n",
    "        bounds=[0.0015, 0.004],\n",
    "        locations=[dend_trunk],\n",
    "        frozen=False)\n",
    "\n",
    "PC_Cav3_1_param_dend_all = ephys.parameters.NrnSectionParameter(\n",
    "        name='pcabar_Cav3_1_dend_all',\n",
    "        param_name='pcabar_Cav3_1',\n",
    "        bounds=[3e-6, 7e-6],\n",
    "        locations=[dend_trunk],\n",
    "        frozen=False)\n",
    "\n",
    "PC_Cav3_2_param_dend_all = ephys.parameters.NrnSectionParameter(\n",
    "        name='gcabar_Cav3_2_dend_all',\n",
    "        param_name='gcabar_Cav3_2',\n",
    "        bounds=[0.00095, 0.002],\n",
    "        locations=[dend_trunk],\n",
    "        frozen=False)\n",
    "\n",
    "#dend_sodium\n",
    "PC_Nav16_param_dend_sodium = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gbar_Nav1_6_dend_sodium',\n",
    "        param_name='gbar_Nav1_6',\n",
    "        locations=[dend_sodium],\n",
    "        bounds=[0.013, 0.017],\n",
    "        frozen=False)  \n",
    "\n",
    "#SOMA\n",
    "PC_Nav16_param_soma = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gbar_Nav1_6_soma',\n",
    "        param_name='gbar_Nav1_6',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[0.185, 0.21],\n",
    "        frozen=False)  \n",
    "\n",
    "PC_Kv1_1_param_soma = ephys.parameters.NrnSectionParameter(\n",
    "        name='gbar_Kv1_1_soma',\n",
    "        param_name='gbar_Kv1_1',\n",
    "        bounds=[0.001, 0.003],\n",
    "        locations=[somatic_loc],\n",
    "        frozen=False)\n",
    "\n",
    "PC_Kv3_4_param_soma = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gkbar_Kv3_4_soma',\n",
    "        param_name='gkbar_Kv3_4',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[0.04, 0.07],\n",
    "        frozen=False) \n",
    "\n",
    "PC_Kir2_3_param_soma = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gkbar_Kir2_3_soma',\n",
    "        param_name='gkbar_Kir2_3',\n",
    "        locations=[somatic_loc],\n",
    "        bounds=[0.00001, 0.00007],\n",
    "        frozen=False)\n",
    "\n",
    "PC_Kca1_1_param_soma = ephys.parameters.NrnSectionParameter(\n",
    "        name='gbar_Kca1_1_soma',\n",
    "        param_name='gbar_Kca1_1',\n",
    "        bounds=[0.008, 0.015],\n",
    "        locations=[somatic_loc],\n",
    "        frozen=False)\n",
    "\n",
    "PC_Kca2_2_param_soma = ephys.parameters.NrnSectionParameter(\n",
    "        name='gkbar_Kca2_2_soma',\n",
    "        param_name='gkbar_Kca2_2',\n",
    "        bounds=[8e-4, 1.5e-3],\n",
    "        locations=[somatic_loc],\n",
    "        frozen=False)\n",
    "\n",
    "PC_Kca3_1_param_soma = ephys.parameters.NrnSectionParameter(\n",
    "        name='gkbar_Kca3_1_soma',\n",
    "        param_name='gkbar_Kca3_1',\n",
    "        bounds=[0.008, 0.015],\n",
    "        locations=[somatic_loc],\n",
    "        frozen=False)\n",
    "\n",
    "PC_Cav2_1_param_soma = ephys.parameters.NrnSectionParameter(\n",
    "        name='pcabar_Cav2_1_soma',\n",
    "        param_name='pcabar_Cav2_1',\n",
    "        bounds=[1.5e-4, 3e-4],\n",
    "        locations=[somatic_loc],\n",
    "        frozen=False)\n",
    "\n",
    "PC_Cav3_1_param_soma = ephys.parameters.NrnSectionParameter(\n",
    "        name='pcabar_Cav3_1_soma',\n",
    "        param_name='pcabar_Cav3_1',\n",
    "        bounds=[4e-6, 9e-6],\n",
    "        locations=[somatic_loc],\n",
    "        frozen=False)\n",
    "\n",
    "PC_Cav3_2_param_soma = ephys.parameters.NrnSectionParameter(\n",
    "        name='gcabar_Cav3_2_soma',\n",
    "        param_name='gcabar_Cav3_2',\n",
    "        bounds=[0.0006, 0.001],\n",
    "        locations=[somatic_loc],\n",
    "        frozen=False)\n",
    "\n",
    "PC_Cav3_3_param_soma = ephys.parameters.NrnSectionParameter(\n",
    "        name='pcabar_Cav3_3_soma',\n",
    "        param_name='pcabar_Cav3_3',\n",
    "        bounds=[0.000095, 0.0004],\n",
    "        locations=[somatic_loc],\n",
    "        frozen=False)\n",
    "\n",
    "PC_HCN1_param_soma = ephys.parameters.NrnSectionParameter(\n",
    "        name='gbar_HCN1_soma',\n",
    "        param_name='gbar_HCN1',\n",
    "        bounds=[0.0008, 0.002],\n",
    "        locations=[somatic_loc],\n",
    "        frozen=False)\n",
    "\n",
    "#AXON\n",
    "#AIS\n",
    "PC_Nav16_param_axon0 = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gbar_Nav1_6_axon0',\n",
    "        param_name='gbar_Nav1_6',\n",
    "        locations=[axon0],\n",
    "        bounds=[0, 0.7],\n",
    "        frozen=False)  \n",
    "\n",
    "PC_Kv3_4_param_axon0 = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gkbar_Kv3_4_axon0',\n",
    "        param_name='gkbar_Kv3_4',\n",
    "        locations=[axon0],\n",
    "        bounds=[0.009, 0.02],\n",
    "        frozen=False) \n",
    "\n",
    "PC_Cav2_1_param_axon0 = ephys.parameters.NrnSectionParameter(\n",
    "        name='pcabar_Cav2_1_axon0',\n",
    "        param_name='pcabar_Cav2_1',\n",
    "        bounds=[1e-4, 4e-4],\n",
    "        locations=[axon0],\n",
    "        frozen=False)\n",
    "\n",
    "PC_Cav3_1_param_axon0 = ephys.parameters.NrnSectionParameter(\n",
    "        name='pcabar_Cav3_1_axon0',\n",
    "        param_name='pcabar_Cav3_1',\n",
    "        bounds=[5e-6, 1e-5],\n",
    "        locations=[axon0],\n",
    "        frozen=False)\n",
    "\n",
    "#AISK\n",
    "PC_Kv1_1_param_axon1 = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gbar_Kv1_1_axon1',\n",
    "        param_name='gbar_Kv1_1',\n",
    "        locations=[axon1],\n",
    "        bounds=[0.008, 0.02],\n",
    "        frozen=False)  \n",
    "\n",
    "#RN\n",
    "PC_Nav16_param_axon = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gbar_Nav1_6_axon',\n",
    "        param_name='gbar_Nav1_6',\n",
    "        locations=[axon3, axon5, axon7],\n",
    "        bounds=[0.025, 0.04],\n",
    "        frozen=False)  \n",
    "\n",
    "PC_Kv3_4_param_axon = ephys.parameters.NrnSectionParameter(                                    \n",
    "        name='gkbar_Kv3_4_axon',\n",
    "        param_name='gkbar_Kv3_4',\n",
    "        locations=[axon3, axon5, axon7],\n",
    "        bounds=[0.015, 0.03],\n",
    "        frozen=False) \n",
    "\n",
    "PC_Cav2_1_param_axon = ephys.parameters.NrnSectionParameter(\n",
    "        name='pcabar_Cav2_1_axon',\n",
    "        param_name='pcabar_Cav2_1',\n",
    "        bounds=[1e-4, 3e-4],\n",
    "        locations=[axon3, axon5, axon7],\n",
    "        frozen=False)\n",
    "\n",
    "PC_Cav3_1_param_axon = ephys.parameters.NrnSectionParameter(\n",
    "        name='pcabar_Cav3_1_axon',\n",
    "        param_name='pcabar_Cav3_1',\n",
    "        bounds=[1e-5, 2e-5],\n",
    "        locations=[axon3, axon5, axon7],\n",
    "        frozen=False)\n",
    "\n",
    "PC2018 = ephys.models.CellModel(\n",
    "        name='Purkinje2018StefanoMasoli',\n",
    "        morph=morph,\n",
    "        mechs=[PC_Leak,\n",
    "           PC_pass,\n",
    "           PC_Nav16, \n",
    "           PC_Kv1_1,\n",
    "           PC_Kv1_5,\n",
    "           PC_Kv3_3,        \n",
    "           PC_Kv3_4, \n",
    "           PC_Kv4_3,\n",
    "           PC_Kir2_3,       \n",
    "           PC_Kca1_1,\n",
    "           PC_Kca2_2,   \n",
    "           PC_Kca3_1,     \n",
    "           PC_Cav2_1,\n",
    "           PC_Cav3_1,\n",
    "           PC_Cav3_2,\n",
    "           PC_Cav3_3,\n",
    "           PC_HCN1,\n",
    "           PC_cdp5],\n",
    "\n",
    "        params=[cm_dend_all,\n",
    "        cm_soma,\n",
    "        cm_rn_only,\n",
    "        cm_myelin_only,\n",
    "        ena,\n",
    "        ek,\n",
    "        eleak,\n",
    "        leak_soma,\n",
    "        leak_rest,\n",
    "        eca,\n",
    "        Ra,\n",
    "        gpas_myelin_only,\n",
    "        epas_myelin_only,\n",
    "        ca_pump_dend,\n",
    "        ca_pump_soma,\n",
    "        ca_pump_axon,\n",
    "        PC_Nav16_param_dend_sodium,\n",
    "        PC_Kv1_1_param_dend_all,\n",
    "        PC_Kv1_5_param_dend_all,\n",
    "        PC_Kv3_3_param_dend_all,\n",
    "        PC_Kv4_3_param_dend_all,\n",
    "        PC_Kir2_3_param_dend_trunk,\n",
    "        PC_Cav2_1_param_dend_all,\n",
    "        PC_Cav3_1_param_dend_all,\n",
    "        PC_Cav3_2_param_dend_all,\n",
    "        PC_Cav3_3_param_dend_all,\n",
    "        PC_Kca1_1_param_dend_all,\n",
    "        PC_Kca2_2_param_dend_trunk,\n",
    "        PC_Kca3_1_param_dend_trunk,\n",
    "        PC_HCN1_param_dend_all,\n",
    "        PC_Nav16_param_soma,\n",
    "        PC_Kv1_1_param_soma,\n",
    "        PC_Kv3_4_param_soma,\n",
    "        PC_Kir2_3_param_soma,\n",
    "        PC_Cav2_1_param_soma,\n",
    "        PC_Cav3_1_param_soma,\n",
    "        PC_Cav3_2_param_soma,\n",
    "        PC_Cav3_3_param_soma,\n",
    "        PC_Kca1_1_param_soma,\n",
    "        PC_Kca2_2_param_soma,\n",
    "        PC_Kca3_1_param_soma,\n",
    "        PC_HCN1_param_soma,\n",
    "        PC_Nav16_param_axon0,\n",
    "        PC_Kv3_4_param_axon0,\n",
    "        PC_Cav2_1_param_axon0,\n",
    "        PC_Cav3_1_param_axon0,\n",
    "        PC_Kv1_1_param_axon1,\n",
    "        PC_Nav16_param_axon,\n",
    "        PC_Kv3_4_param_axon,\n",
    "        PC_Cav2_1_param_axon,\n",
    "        PC_Cav3_1_param_axon])  \n",
    "\n",
    "PC2018.seclist_names.append('dendall')\n",
    "PC2018.seclist_names.append('dendtrunk')\n",
    "PC2018.seclist_names.append('dendsodium')\n",
    "\n",
    "PC2018.seclist_names.append('AIS')\n",
    "PC2018.seclist_names.append('AISK')\n",
    "PC2018.seclist_names.append('myelin0')\n",
    "PC2018.seclist_names.append('rn0')\n",
    "PC2018.seclist_names.append('myelin1')\n",
    "PC2018.seclist_names.append('rn1')\n",
    "PC2018.seclist_names.append('myelin2')\n",
    "PC2018.seclist_names.append('rn2')\n",
    "PC2018.seclist_names.append('myelin3')\n",
    "\n",
    "print( PC2018)\n",
    "\n",
    "\n",
    "soma_loc = ephys.locations.NrnSeclistCompLocation(\n",
    "        name='soma',\n",
    "        seclist_name='somatic',\n",
    "        sec_index=0,\n",
    "        comp_x=0.5)\n",
    "\n",
    "\n",
    "sweep_protocols = []\n",
    "for protocol_name, amplitude in [('step1', 0)]:\n",
    "    stim = ephys.stimuli.NrnSquarePulse(\n",
    "                step_amplitude=amplitude,\n",
    "                step_delay=100,\n",
    "                step_duration=1000,\n",
    "                location=soma_loc,\n",
    "                total_duration=1000)\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)\n",
    "\n",
    "nrn = ephys.simulators.NrnSimulator()\n",
    "\n",
    "\n",
    "for i in (good_results):\n",
    "    print( 'Simulation of validated individual:',i)\n",
    "    default_params = {'gbar_Nav1_6_dend_sodium' : conductvalues[i,0],\n",
    "              'gbar_Kv1_1_dend_all': conductvalues[i,1], \n",
    "              'gKur_Kv1_5_dend_all' : conductvalues[i,2],\n",
    "              'gbar_Kv3_3_dend_all' : conductvalues[i,3],\n",
    "              'gkbar_Kv4_3_dend_all' : conductvalues[i,4],\n",
    "              'gkbar_Kir2_3_dend_trunk' : conductvalues[i,5],\n",
    "              'pcabar_Cav2_1_dend_all' : conductvalues[i,6],\n",
    "              'pcabar_Cav3_1_dend_all' : conductvalues[i,7],\n",
    "              'gcabar_Cav3_2_dend_all' : conductvalues[i,8],\n",
    "              'pcabar_Cav3_3_dend_all' : conductvalues[i,9],\n",
    "              'gbar_Kca1_1_dend_all' : conductvalues[i,10], \n",
    "              'gkbar_Kca2_2_dend_trunk' : conductvalues[i,11],\n",
    "              'gkbar_Kca3_1_dend_trunk' : conductvalues[i,12],\n",
    "              'gbar_HCN1_dend_all' : conductvalues[i,13],\n",
    "              'gbar_Nav1_6_soma' : conductvalues[i,14],\n",
    "              'gbar_Kv1_1_soma' : conductvalues[i,15],\n",
    "              'gkbar_Kv3_4_soma': conductvalues[i,16],\n",
    "              'gkbar_Kir2_3_soma': conductvalues[i,17],\n",
    "              'pcabar_Cav2_1_soma' : conductvalues[i,18],\n",
    "              'pcabar_Cav3_1_soma' : conductvalues[i,19],\n",
    "              'gcabar_Cav3_2_soma' : conductvalues[i,20],\n",
    "              'pcabar_Cav3_3_soma' : conductvalues[i,21],\n",
    "              'gbar_Kca1_1_soma' : conductvalues[i,22], \n",
    "              'gkbar_Kca2_2_soma' : conductvalues[i,23],\n",
    "              'gkbar_Kca3_1_soma' : conductvalues[i,24],\n",
    "              'gbar_HCN1_soma' : conductvalues[i,25],\n",
    "              'gbar_Nav1_6_axon0':0,\n",
    "              'gkbar_Kv3_4_axon0':conductvalues[i,27],\n",
    "              'pcabar_Cav2_1_axon0': conductvalues[i,28],\n",
    "              'pcabar_Cav3_1_axon0' : conductvalues[i,29], \n",
    "              'gbar_Kv1_1_axon1':conductvalues[i,30],\n",
    "              'gbar_Nav1_6_axon':conductvalues[i,31],\n",
    "              'gkbar_Kv3_4_axon':conductvalues[i,32],\n",
    "              'pcabar_Cav2_1_axon':conductvalues[i,33],\n",
    "              'pcabar_Cav3_1_axon' : conductvalues[i,34]} \n",
    "\n",
    "\n",
    "    responses = threestep_protocol.run(cell_model=PC2018, param_values=default_params, sim=nrn)\n",
    "    step1time = responses['step1.soma.v']['time']\n",
    "    step1vm = responses['step1.soma.v']['voltage']\n",
    "    sommastep1 = list(zip(step1time,step1vm))\n",
    "    saveval = str(i)\n",
    "\n",
    "\n",
    "    np.savetxt(str(saveval) + '_sim_step_noais' + '.txt', sommastep1, delimiter = ' ')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell plots the validated results based on the frequency absence of firing without sodium chanells in the Axon Initial Segment AIS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot single result.\n",
    "\n",
    "for i in good_results: #Number of the current step used in the optimization. 1, first step, 2 second step and so on.\n",
    "    filetxt = str(i) + '_sim_step_noais.txt'\n",
    "\n",
    "    trace = np.genfromtxt(filetxt)\n",
    "\n",
    "    time, voltage = zip(*trace)\n",
    "        \n",
    "    fig, ax = plt.subplots()\n",
    "    ax.plot(time, voltage, 'b', label='spikes indiv. %d' % (i))\n",
    "\n",
    "    legend = ax.legend(loc='upper right', shadow=True)\n",
    "\n",
    "    frame = legend.get_frame()\n",
    "    frame.set_facecolor('0.90')\n",
    "\n",
    "    plt.xlabel(\"Time (ms)\")\n",
    "    plt.ylabel(\"Membrane voltage (mv) \")                                          "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot single result.\n",
    "import os\n",
    "\n",
    "# Location of the recording. For the Purkinje cell there is only one.\n",
    "location = '' \n",
    "\n",
    "\n",
    "for z in good_results:\n",
    "    for current_inj_step in range (1,4): \n",
    "        filetxt = valid_dir + str(z) + '_sim_step' + str(current_inj_step) + location +'.txt'\n",
    "        if os.path.isfile(filetxt) == False:\n",
    "            print( 'file not found')\n",
    "            pass\n",
    "    \n",
    "        elif os.path.isfile(filetxt):\n",
    "            trace = np.genfromtxt(filetxt)\n",
    "\n",
    "            time, voltage = zip(*trace)\n",
    "        \n",
    "            fig, ax = plt.subplots()\n",
    "            ax.plot(time, voltage, 'b', label='spikes indiv. %d, %s' % (z, steps[current_inj_step]))\n",
    "\n",
    "            legend = ax.legend(loc='upper right', shadow=True)\n",
    "\n",
    "            frame = legend.get_frame()\n",
    "            frame.set_facecolor('0.90')\n",
    "\n",
    "            plt.xlabel(\"Time (ms)\")\n",
    "            plt.ylabel(\"Membrane voltage (mv) \")"
   ]
  },
  {
   "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
}