{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quickstart\n",
    "To run the code below:\n",
    "\n",
    "1. Click on the cell to select it.\n",
    "2. Press `SHIFT+ENTER` on your keyboard or press the play button\n",
    "   (<button class='fa fa-play icon-play btn btn-xs btn-default'></button>) in the toolbar above.\n",
    "\n",
    "Feel free to create new cells using the plus button\n",
    "(<button class='fa fa-plus icon-plus btn btn-xs btn-default'></button>), or pressing `SHIFT+ENTER` while this cell\n",
    "is selected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Appendix 4: parameter exploration\n",
    "\n",
    "In this example, we show how Brian can be used to do a parameter exploration. The most efficient way to implement this is to consider a group of neurons, where each of the neurons is described by the same model, but one or several model parameters systematically change across neurons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "from brian2 import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can speed up this \"embarassingly parallel\" problem by parallelizing it over CPUs, making use of the OpenMP interface for multithreading:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set_device('cpp_standalone')\n",
    "prefs.devices.cpp_standalone.openmp_threads = 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the [Brian2GeNN](https://brian2genn.readthedocs.io) interface, this simulation could also run on a GPU (needs a [CUDA](https://developer.nvidia.com/cuda-zone)-enabled NVIDIA GPU, and an installation of `brian2genn`): "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# These two lines would replace the above two lines:\n",
    "# import brian2genn\n",
    "# set_device('genn')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will first set some general constants for our model (a HH-type model with an injected current):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "area = 20000*umetre**2\n",
    "Cm = (1*ufarad*cm**-2) * area\n",
    "gl = (5e-5*siemens*cm**-2) * area\n",
    "\n",
    "El = -60*mV\n",
    "EK = -90*mV\n",
    "ENa = 50*mV\n",
    "g_kd = (30*msiemens*cm**-2) * area\n",
    "VT = -63*mV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will explore the effect of the $g_{Na}$ conductance, and the effect of the injected current $I$, varying each over 100 values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g_na_values = np.linspace(10, 100, num=100)*msiemens*cm**-2 * area\n",
    "I_values = np.linspace(0, 20, num=100)*pA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now define the model, with $g_{Na}$ and $I$ as parameters, set independently for each neuron, and create one neuron for each combination of parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eqs = Equations('''\n",
    "dv/dt = (gl*(El-v)-\n",
    "         g_na*(m*m*m)*h*(v-ENa)-\n",
    "         g_kd*(n*n*n*n)*(v-EK) + I)/Cm : volt\n",
    "dm/dt = alpha_m*(1-m)-beta_m*m : 1\n",
    "dn/dt = alpha_n*(1-n)-beta_n*n : 1\n",
    "dh/dt = alpha_h*(1-h)-beta_h*h : 1\n",
    "alpha_m = 0.32*(mV**-1)*(13*mV-v+VT)/\n",
    "         (exp((13*mV-v+VT)/(4*mV))-1.)/ms : Hz\n",
    "beta_m = 0.28*(mV**-1)*(v-VT-40*mV)/\n",
    "        (exp((v-VT-40*mV)/(5*mV))-1)/ms : Hz\n",
    "alpha_h = 0.128*exp((17*mV-v+VT)/(18*mV))/ms : Hz\n",
    "beta_h = 4./(1+exp((40*mV-v+VT)/(5*mV)))/ms : Hz\n",
    "alpha_n = 0.032*(mV**-1)*(15*mV-v+VT)/\n",
    "         (exp((15*mV-v+VT)/(5*mV))-1.)/ms : Hz\n",
    "beta_n = .5*exp((10*mV-v+VT)/(40*mV))/ms : Hz\n",
    "I : amp (constant)\n",
    "g_na : siemens (constant)\n",
    "''')\n",
    "neuron = NeuronGroup(len(g_na_values)*len(I_values), eqs,\n",
    "                     method='exponential_euler',\n",
    "                     threshold='v>-20*mV', refractory='v>-20*mV')\n",
    "neuron.v = El\n",
    "spike_mon = SpikeMonitor(neuron)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now create all combinations of the explored parameters, and assign them to the individual neurons. Running the \"network\" of neurons will then explore all the parameter combinations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_g_na_values, all_I_values = np.meshgrid(g_na_values, I_values)\n",
    "all_g_na_values = all_g_na_values.flat[:]\n",
    "all_I_values = all_I_values.flat[:]\n",
    "neuron.g_na = all_g_na_values\n",
    "neuron.I = all_I_values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start the run and determine the firing rate for each neuron:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run(10*second)\n",
    "rates = spike_mon.count/(10*second)/Hz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now plot the firing rate as a function of the two explored parameters (in this example here, vertical slices of the graph can also be interpreted as the f/I curve of the neuron for varying values of $g_{Na}$):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I_index = ((neuron.I - np.min(I_values)) / np.diff(I_values)[0]).round().astype('int')\n",
    "Na_index = ((neuron.g_na - np.min(g_na_values)) / np.diff(g_na_values)[0]).round().astype('int')\n",
    "matrix = np.full((len(I_values), len(g_na_values)), np.nan)\n",
    "matrix[I_index, Na_index] = rates\n",
    "norm = mpl.colors.BoundaryNorm(np.arange(0, 19), plt.cm.viridis.N)\n",
    "fig, ax = plt.subplots()\n",
    "m = ax.imshow(matrix, norm=norm)\n",
    "# We do manual ticks, easier than using extent and getting the scaling right\n",
    "ticks = [0, 49, 99]\n",
    "ax.set(xticks=ticks, xticklabels=['%.1f' % (g_na_values[i]/(mS*cm**-2 * area)) for i in ticks],\n",
    "       yticks=ticks, yticklabels=['%.1f' % (I_values[i]/pA) for i in ticks],\n",
    "       xlabel='$g_{Na}$ (mS/cm²)', ylabel='$I$ (nA)')\n",
    "ax.axes.invert_yaxis()\n",
    "cbar = fig.colorbar(m)\n",
    "cbar.set_label('number of spikes')"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}