NeuroDyn Python model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NeuroDyn Python model\n",
    "\n",
    "## Overview\n",
    "\n",
    "This guide introduces the Python software environment for setting up and simulating a model of the NeuroDyn hardware equations. By the end of the guide you should get an understanding of the dynamics of a NeuroDyn neuron, be able to construct a NeuroDyn object, as well as understand how to fit the parameters of the model in order to emulate the biophysical conductance-based models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up the environment\n",
    "\n",
    "We first need to import the modules with the appropriate model definitions, as well as the Python packages that we will need:\n",
    "\n",
    "- `matplotlib.pyplot`: Python plotting module\n",
    "- `numpy`: Python library for numerical operations\n",
    "- `cb_models`: model definitions for both standard conductance-based models as well as NeuroDyn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from cb_models import HHModel, NeuroDynModel\n",
    "\n",
    "# **Ignore overflow warnings**\n",
    "import numpy as np\n",
    "np.seterr(all=\"ignore\")\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hodgkin-Huxley model\n",
    "\n",
    "We start with a short overview of the standard Hodgkin-Huxley model of squid giant axon and how we can define it using the environment.\n",
    "\n",
    "The Hodgkin-Huxley model of the neural membrane dynamics is given by the following set of equations:\n",
    "\n",
    "\\begin{align}\n",
    "C \\dot{V} &= - \\bar{g}_l (V - E_l) - \\bar{g}_{Na} m^3 h (V - E_{Na}) - \\bar{g}_K n^4 (V - E_K) + I_{app} \\\\\n",
    "\\dot{m} &= \\alpha_m (V) (1 - m) - \\beta_m (V) m \\\\\n",
    "\\dot{h} &= \\alpha_h (V) (1 - h) - \\beta_h (V) h \\\\\n",
    "\\dot{n} &= \\alpha_n (V) (1 - n) - \\beta_n (V) n\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "Let's start by defining a Hodgkin-Huxley model of a neuron that implements these equations - a NeuroDyn model will be done in very much the same way in the following sections."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hh_neuron = HHModel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By passing no arguments to the `HHModel` class, an instance of the neuron is created with default parameters taken from the original Hodgkin-Huxley equations. Note that the values of the constants are shifted so that the resting potential is at $0 mV$. If we instead wanted to initialize a model with, for example, different sodium current parameters, we could do that as:\n",
    "\n",
    "`hh_neuron = HHModel(gna = 100, Ena = 50)`\n",
    "\n",
    "which would set $g_{Na} = 100 mS$ and $E_{Na} = 50 mV$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's try simulating the model for which we only need to define the $I_{app}(t)$ function for the applied current. We can try different values for the magnitude of the constant current $I_0$ to showcase both resting and spiking behaviors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Iapp(t):\n",
    "    I0 = 10 # applied current in uA\n",
    "    return I0\n",
    "\n",
    "T = 100 # length of the simulation in ms\n",
    "trange = (0, T)\n",
    "\n",
    "x0 = [0,0,0,0] # initial conditions [V0, m0, h0, n0]\n",
    "\n",
    "sol = hh_neuron.simulate(trange, x0, Iapp)\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"Time [$ms$]\")\n",
    "plt.ylabel(\"Voltage [$mV$]\")\n",
    "plt.title(\"Membrane voltage\")\n",
    "plt.plot(sol.t, sol.y[0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Additional options\n",
    "\n",
    "##### Converting units to SI\n",
    "\n",
    "There is a few things to note here. For convenience, all of the parameters and variables are in non-SI units, so that the voltages are expressed in $mV$, currents in $\\mu A$, time in $m s $, etc. If we instead wanted to convert the model to SI units, we can initialize the neuron as:\n",
    "\n",
    "`hh_neuron = HHModel(SI_units = True)`\n",
    "\n",
    "which will convert all parameters of the model to SI units."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Voltage and time scaling\n",
    "\n",
    "A few other options that will be useful once we start defining NeuroDyn models are the parameters for scaling the time and the voltage scales of the model, `scl_t` and `scl_v`. Let's see how this works by instead defining a neuron as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hh_neuron_scaled = HHModel(scl_v = 3, scl_t = 2)\n",
    "\n",
    "def Iapp(t):\n",
    "    I0 = 10 # applied current in uA\n",
    "    return I0\n",
    "\n",
    "T = 100 # length of the simulation in ms\n",
    "trange = (0, T)\n",
    "\n",
    "x0 = [0,0,0,0] # initial conditions [V0, m0, h0, n0]\n",
    "\n",
    "sol = hh_neuron_scaled.simulate(trange, x0, Iapp)\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"Time [$ms$]\")\n",
    "plt.ylabel(\"Voltage [$mV$]\")\n",
    "plt.title(\"Membrane voltage\")\n",
    "plt.plot(sol.t, sol.y[0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can try changing the values for `scl_t` and `scl_v` and see how it changes the behavior compared to the non-scaled simulation we have previously produced."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we try fitting the conductance-based models to the NeuroDyn hardware we will see that it will be necessary to *scale* both voltage and time of the original model in order to fit them to the hardware constraints. This will change the timescale and the voltage range in which the NeuroDyn neuron operates, but importantly, will not change the underlying dynamical properties of the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NeuroDyn model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initializing a NeuroDyn model\n",
    "\n",
    "Creating a NeuroDyn object is done in exactly the same way as creating a Hodgkin-Huxley object by using the class `NeuroDynModel`. We can therefore create a NeuroDyn neuron using default parameters by calling:\n",
    "\n",
    "`nd_neuron = NeuroDynModel()`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that NeuroDyn uses SI units for all of its parameters and variables, so that the voltage output is in $V$, time is in $s$, etc.\n",
    "\n",
    "Before we create a NeuroDyn neuron and use it to simulate spiking behavior, let's first look at the NeuroDyn equations and how they connect to the Hodgkin-Huxley equations. The details of the section may be skipped and readers might go directly to the **Summary** that lists all of the controllable parameters and the mappings to the corresponding Hodgkin-Huxley parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NeuroDyn equations\n",
    "\n",
    "Let's now take a look at the full list of parameters that are required to define a NeuroDyn model. NeuroDyn circuit replicates the equations of Hodgkin-Huxley in hardware, and we will now see how the different parameters map to the biophysical constants of the conductance-based models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Membrane equation\n",
    "\n",
    "Let's start with the membrane dynamics equation:\n",
    "\n",
    "\\begin{align}\n",
    "C_{ND} \\dot{V} &= -I_l - I_{Na} - I_K + I_{app}\\\\\n",
    " &= - 2 I_{g_l} \\tanh \\bigg( \\frac{\\kappa_{lin} (V - E_l)}{2 V_T} \\bigg) - 2 I_{g_{Na}} \\tanh \\bigg( \\frac{\\kappa_{lin} (V - E_{Na})}{2 V_T} \\bigg) - 2 I_{g_{K}} \\tanh \\bigg( \\frac{\\kappa_{lin} (V - E_K)}{2 V_T} \\bigg)\n",
    "\\end{align}\n",
    "\n",
    "Note that capacitance of NeuroDyn $C_{ND} = 4 pF$ is significantly smaller than the capacitance of Hodgkin-Huxley $C_{HH} = 1 \\mu F$ which we will need to take into consideration when setting the other parameters of the circuit. \n",
    "\n",
    "The nonlinear $\\tanh$ relationship comes from the transconductance amplifier characteristic that implements the individual conductance relatinships. These amplifiers have additional circuitry to linearize this relationship which makes the constant $\\kappa_{lin}$ small. In the simulation, we set $\\kappa_{lin} = 0.2$. The constant $V_T$ is the temperature-dependent thermal voltage and is assumed to be $V_T \\approx 26 mV$.\n",
    "\n",
    "##### Maximal conductances\n",
    "\n",
    "Assuming the conductance circuit is in the linear range so that $\\tanh{x} \\approx x$, we can define the relationship between the circuit parameters and the original Hodgkin-Huxley equation by comparing the two directly:\n",
    "\n",
    "\\begin{equation}\n",
    "g_l = \\frac{\\kappa_{lin}}{V_T} I_{g_l}, \\quad\n",
    "g_{Na} = \\frac{\\kappa_{lin}}{V_T} I_{g_{Na}}, \\quad\n",
    "g_{K} = \\frac{\\kappa_{lin}}{V_T} I_{g_{K}}\n",
    "\\end{equation}\n",
    "\n",
    "It is clear from these equations that the conductances $g_l$, $g_{Na}$ and $g_K$ are effectively represented by three currents in the circuit $I_{g_l}$, $I_{g_{Na}}$ and $I_{g_K}$. The current $I_{g_l}$ is constant, reflecting the static nature of the leak current, while currents $I_{g_{Na}}$ and $I_{g_K}$ have the dynamics replicating the ion channel opening and closing processes.\n",
    "\n",
    "##### Reversal potentials\n",
    "\n",
    "The reversal potentials $E_l$, $E_{Na}$ and $E_K$ are defined using digital values using the following equation:\n",
    "\n",
    "\\begin{equation}\n",
    "E_{rev} = I_{voltage} \\cdot \\frac{dE_{rev}}{1024} \\cdot sign \\cdot R_{rev} + V_{ref}\n",
    "\\end{equation}\n",
    "\n",
    "The current $I_{voltage}$ is controllable and is shared by all the variables in the circuit. Similarly, $V_{ref}$ is controllable and is shared by all the circuits - it effectively shifts all of the parameters in the voltage range. The resistor $R_{rev} = 1.63 M\\Omega$ is fixed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ion currents\n",
    "\n",
    "Let's now look at the dynamics for the currents $I_{g_{Na}}$ and $I_{g_K}$:\n",
    "\n",
    "\\begin{align}\n",
    "I_{g_{Na}} &= I_{\\bar{g}_{Na}} \\bigg(\\frac{I_m}{I_{ref}} \\bigg)^3 \\bigg(\\frac{I_h}{I_{ref}} \\bigg)\\\\\n",
    "I_{g_K} &= I_{\\bar{g}_{K}} \\bigg(\\frac{I_n}{I_{ref}} \\bigg)^4\n",
    "\\end{align}\n",
    "\n",
    "Comparing to the equations of Hodgkin-Huxley, it is clear that the relationship between gating variables and currents is simply:\n",
    "\n",
    "\\begin{equation}\n",
    "m = \\frac{I_m}{I_{ref}}, \\quad\n",
    "h = \\frac{I_h}{I_{ref}}, \\quad\n",
    "n = \\frac{I_n}{I_{ref}}\n",
    "\\end{equation}\n",
    "\n",
    "The currents determining the maximal conductance parameters are set by digital values $dg$ that set the currents between $0$ and the current $I_{master}$:\n",
    "\n",
    "\\begin{equation}\n",
    "I_{\\bar{g}} = I_{master} \\frac{dg}{1024}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Gating variables\n",
    "\n",
    "Let's now take a look at the gating variable equations. The dynamical equations for the gating currents are given by:\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{d}{dt} \\bigg( \\frac{I_m}{I_{ref}} \\bigg) &= \\frac{I_{\\alpha_m}}{C_g V_T} \\bigg(1 - \\frac{I_m}{I_{ref}} \\bigg) - \\frac{I_{\\beta_m}}{C_g V_T} \\bigg( \\frac{I_m}{I_{ref}} \\bigg) \\\\\n",
    "\\frac{d}{dt} \\bigg( \\frac{I_h}{I_{ref}} \\bigg) &= \\frac{I_{\\alpha_h}}{C_g V_T} \\bigg(1 - \\frac{I_h}{I_{ref}} \\bigg) - \\frac{I_{\\beta_h}}{C_g V_T} \\bigg( \\frac{I_h}{I_{ref}} \\bigg) \\\\\n",
    "\\frac{d}{dt} \\bigg( \\frac{I_n}{I_{ref}} \\bigg) &= \\frac{I_{\\alpha_n}}{C_g V_T} \\bigg(1 - \\frac{I_n}{I_{ref}} \\bigg) - \\frac{I_{\\beta_n}}{C_g V_T} \\bigg( \\frac{I_n}{I_{ref}} \\bigg)\n",
    "\\end{align}\n",
    "\n",
    "Comparing again to the original Hodgkin-Huxley equations, we see that a simple substitution relates $I_{\\alpha}$ and $I_{\\beta}$ with $\\alpha (V)$ and $\\beta (V)$ for each gating variable $m$, $h$ and $n$:\n",
    "\n",
    "\\begin{equation}\n",
    "\\alpha (V) = \\frac{I_{\\alpha} (V)}{C_g V_T}, \\quad\n",
    "\\beta(V) = \\frac{I_{\\beta} (V)}{C_g V_T}\n",
    "\\end{equation}\n",
    "\n",
    "The constant $V_T$ is the same thermal voltage defined previously, while the capacitance used in the gating variable equations $C_g = 5 p F$ is different from the membrane capacitance $C_{ND}$.\n",
    "\n",
    "\n",
    "So far, we have seen that there is a one-to-one correspondence between the physical parameters of NeuroDyn and the Hodgkin-Huxley parameters. Let's now see how we can shape the $\\alpha(V)$ and $\\beta(V)$ functions in NeuroDyn by designing the corresponding currents $I_{\\alpha}$ and $I_{\\beta}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Designing $\\alpha (V)$ and $\\beta (V)$ functions\n",
    "\n",
    "Each current $I_{\\alpha}$ and $I_{\\beta}$ is built as a summation of **seven** sigmoidal IV currents, so that by changing the weights of these seven sigmoids we are able to fit the parameters to different $\\alpha (V)$ and $\\beta (V)$ functions. Therefore, the currents $I_{\\alpha}$ and $I_{\\beta}$ for each gating variable are represented by:\n",
    "\n",
    "\\begin{align}\n",
    "I_{\\alpha} &= \\sum_{j=1}^7 \\frac{I_{\\alpha, j}}{1 + \\exp{\\big( sign \\cdot \\kappa \\cdot (V_{b,j} - V)} \\big)} \\\\\n",
    "I_{\\beta} &= \\sum_{j=1}^7 \\frac{I_{\\beta, j}}{1 + \\exp{\\big( sign \\cdot \\kappa \\cdot (V_{b,j} - V)} \\big)}\n",
    "\\end{align}\n",
    "\n",
    "The process-dependent constant $\\kappa \\approx 0.7$ is constant throughout the circuit. Note that this is different from $\\kappa_{lin}$ that was encountered in the membrane equation, where the constant is decreased in order to increase the linear range of the $\\tanh$ nonlinearity. \n",
    "\n",
    "Importantly, note that all sigmoids in the summation have a $sign$ parameter that is either $+1$ or $-1$ . This is because $\\alpha(V)$ and $\\beta(V)$ are either monotonically increasing or decreasing functions - a monotonically increasing function will be built as a summation of monotonically increasing sigmoids, and vice-versa, a monotonically decreasing function will be built as a summation of monotonically decreasing sigmoids. So for example, the gating variable $m$ is an activating variable, so its $\\alpha(V)$ function is monotonically increasing, while the gating variable $h$ is an inactivating variable, and so its $\\alpha(V)$ function is instead monotonically decreasing.\n",
    "\n",
    "The sigmoid placements $V_{b,j}$ are equally spaced and shared by all currents in the NeuroDyn chip, and can be controlled by the current $I_{voltage}$.\n",
    "\n",
    "Each of the sigmoid weights $I_{\\alpha, j}$ and $I_{\\beta, j}$ are set by digital values in the same way as we have defined previously for maximal conductance currents:\n",
    "\n",
    "\\begin{align}\n",
    "I_{\\alpha, j} &= I_{master} \\frac{dI_{\\alpha,j}}{1024} \\\\ \n",
    "I_{\\beta, j} &= I_{master} \\frac{dI_{\\beta,j}}{1024}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary of the parameters\n",
    "\n",
    "The parameters that can be passed to the NeuroDyn class are the following:\n",
    "\n",
    "- `dg = [dgna, dgk, dgl]`: digital values for maximal conductances $\\bar{g}_{Na}$, $\\bar{g}_{K}$ and $\\bar{g}_l$\n",
    "- `dErev = [dEna, dEk, dEl]`: digital values for reversal potentials $E_{Na}$, $E_K$ and $E_l$\n",
    "- `dIb = [dIb_m, dIb_h, dIb_n]`: digital values for the weights of the spline functions\n",
    "    - each `dIb_m` `dIb_h` and `dIb_n` is an array consisting of `[dIb_alpha, dIb_beta]`\n",
    "    - each `dIb_alpha` and `dIb_beta` is an array of 7 digital values for the spline coefficients\n",
    "- `V_ref`: analog reference voltage (shifts all voltage values to be around $V_{ref}$)\n",
    "- `I_voltage`: analog current determining the range of reversal potentials, as well as spline positions\n",
    "- `I_master`: analog current determining the range of maximal conductances, as well as spline weights\n",
    "- `I_ref`: reference current for gating variables\n",
    "    - **Note**: the current has no effect on the model equations, but is a controllable parameter in the chip\n",
    "- `capacitance_scaling`: additional scaling factor that multiplies the membrane capacitance of the model\n",
    "    - **Note**: by default, the scaling factor is $1.0$, and the membrane capacitance is set to $4 pF$\n",
    "    - Setting `capacitance_scaling = 0.1` for instance would set the membrane capacitance to $0.4 pF$\n",
    "    \n",
    "Note that the currents `I_voltage`, `I_master` and `I_ref` are directly passed in the simulation, while in the chip we control these currents with corresponding voltage values. All digital values are integers between $0$ and $1023$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fitting the equations\n",
    "\n",
    "Let's now see how we can fit the parameters of the NeuroDyn model so that it replicates the behavior of a conductance-based model of a neuron. For this we will use the `fitting_utilities` module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from fitting_utilities import FitND"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The class `FitND` takes a single `HHModel` object at the initialization, defining the conductance-based neuron we would like to fit to.\n",
    "\n",
    "We can also pass additional arguments that are otherwise set to their default values:\n",
    "- `I_master`: the $I_{master}$ current that we have set for the NeuroDyn model. Default value is $150 nA$.\n",
    "- `I_voltage`: the $I_{voltage}$ current that we have set for the NeuroDyn model. Default value is $200 nA$.\n",
    "- `capacitance_scaling`: capacitance scaling for the NeuroDyn model we will use. Default values is $1$.\n",
    "- `vrange`: numpy array of voltage points at which we evaluate the fitting. Default is `np.arange(Ek, Ena/2, 5e-4).T`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first try to pass a regular Hodgkin-Huxley neuron to the fitting class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hh_neuron = HHModel(SI_units = True)\n",
    "\n",
    "fit_obj = FitND(hh_neuron)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we have set `SI_units = True` - this is important as all of the NeuroDyn values are defined in SI units, so we would like to first convert our original model to SI units before we continue fitting.\n",
    "\n",
    "If we now call the method `fitHH` of the class, the code will conduct a non-negative least squares fit in order to find the right sigmoid weights $I_{b,j}$ to fit the gating variables of the Hodgkin-Huxley model. The fitting first tries to find the optimal *analog* values for these currents, and then converts them to digital values based on the quantization discussed previously.\n",
    "\n",
    "We can pass optinal keyword arguments to the `fitHH` function that are by default set to `False`:\n",
    "\n",
    "- `plot_alpha_beta = True`: Plot the comparison of the original $\\alpha(V)$ and $\\beta(V)$ functions to the fitting results\n",
    "- `plot_tau_inf = True`: Plot the comparison of the original $\\tau(V)$ and $m_{\\infty}(V)$ functions to the fitting results\n",
    "\n",
    "Let's see how this works:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = fit_obj.fitHH(plot_alpha_beta=True) # weights stores the sigmoid coefficients"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the fits do not look very good initially.\n",
    "\n",
    "Firstly, it seems that the slopes of the sigmoids are too gentle to capture the steeper slopes of the $\\alpha(V)$ and $\\beta(V)$ functions, which we observe by comparing the blue and the orange curves. As the slopes of the sigmoid functions are fixed by the physical implementations, we can instead change the slopes of the Hodgkin-Huxley model instead! By changing the `scl_v` parameter of the Hodgkin-Huxley model, we can stretch the voltage range so that our sigmoids fit better to the stretched $\\alpha(V)$ and $\\beta(V)$ functions. This will effectively render the spikes of the NeuroDyn model larger.\n",
    "\n",
    "Secondly, we can see that there is a significant discrepancy between the orange line (continuous fit) and the green line (fit after quantization). This happens when the current coefficients are to small to be captured by our quantization range - remember that the smallest current that NeuroDyn can output here is\n",
    "\n",
    "\\begin{equation}\n",
    "\\Delta I_b = \\frac{I_{master}}{1024}\n",
    "\\end{equation}\n",
    "\n",
    "so that if the currents returned by the fitting procedure are smaller than this, they will be set to $0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Changing the voltage scaling\n",
    "\n",
    "Let's now try instead to fit the NeuroDyn model to a *scaled* version of the Hodgkin-Huxley object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scl_v = 3 # voltage scaling\n",
    "\n",
    "hh_neuron = HHModel(scl_v = scl_v, SI_units = True)\n",
    "\n",
    "fit_obj = FitND(hh_neuron)\n",
    "\n",
    "weights = fit_obj.fitHH(plot_alpha_beta = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results now look better.\n",
    "\n",
    "Still, we can see in the plot for $\\beta_n(V)$ that the currents are too small to be captured by the quantization. In order to improve the quantized results, we can try manually increasing the coefficients corresponding to the function $\\beta_n(V)$.\n",
    "\n",
    "You can also try changing the fitting parameters to see how they affect the procedure. The parameters that have the most effect are `I_voltage` as it controls the spread of the sigmoids in the voltage range, as well as `scl_v` as it effectively changes the relative slope of the sigmoids."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Scaling the time\n",
    "\n",
    "We need to determine what is the appropriate `scl_t` for the Hodgkin-Huxley model that will allow us to emulate it in hardware.\n",
    "\n",
    "First, let's address the fact that the capacitance of $C_{ND}$ is different from the capacitance of Hodgkin-Huxley $C_{HH}$. Looking at the membrane voltage equation of Hodgkin-Huxley, we can see that we can keep the dynamics intact if we multiple both sides of the equation by the appropriate constant so that on the left we have the capacitance of the NeuroDyn circuit. This gives the corresponding NeuroDyn maximal conductances and the applied current:\n",
    "\n",
    "\\begin{align}\n",
    "\\bar{g}^{(ND)} &= \\bar{g}^{(HH)} \\cdot \\frac{C_{ND}}{C_{HH}} \\\\\n",
    "I_{app}^{(ND)} &= I_{app}^{(HH)} \\cdot \\frac{C_{ND}}{C_{HH}}\n",
    "\\end{align}\n",
    "\n",
    "This effectively means that we need to decrease the conductances and the applied current of the NeuroDyn circuit to adjust for the difference in membrane capacitance.\n",
    "\n",
    "We now need to ensure that the required parameter values are realizable within the hardware constraints. Remember that the $I_{master}$ current in the NeuroDyn chip effectively sets the range of possible values for the maximal conductance values $\\bar{g}$ as well as the sigmoid weights $I_b$.\n",
    "\n",
    "After we have scaled the maximal conductances appropriately, and have received the analog values for the bias currents $I_b$ from the fitting procedure, we need to make sure that all current values are within the range $(0, I_{master})$. If not, we need to scale all currents so that the maximal value fits into this range, effectively determining the constant `scl_t`.\n",
    "\n",
    "This procedure is done automatically when we call the `fitHH` method!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Getting the parameters\n",
    "\n",
    "The previous procedure of finding the right time scaling is done automatically by the quantization method in the fitting module. The quantization takes into account the maximum value of the maximal conductance parameter we would like to fit to, as well as the value of the maximal sigmoid weight, and scales the time scale so that the parameters fit within the NeuroDyn physical range.\n",
    "\n",
    "Let's now obtain the digital parameters that we can then pass to our NeuroDyn class. For this, we can use the methods `get_digital_Ib(weights)`, `get_digital_g(g)` and `get_digital_E(E)` that will give us the digital values for the sigmoid weights, maximal conductances, and reversal potentials, respectively. The values will be automatically scaled, taking into account the `scl_v` parameter of the Hodgkin-Huxley model passed at the initialization, as well as the time scaling `scl_t` necessary for all the parameters to fit within the physical range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g0 = [120e-3,36e-3,0.3e-3] # maximal conductances\n",
    "E0 = [120e-3,-12e-3,10.6e-3] # reversal potentials\n",
    "dIb = fit_obj.get_digital_Ib(weights)\n",
    "dg = fit_obj.get_digital_g(g0)\n",
    "dE = fit_obj.get_digital_E(E0)\n",
    "\n",
    "I_master = fit_obj.I_master\n",
    "I_voltage = fit_obj.I_voltage\n",
    "\n",
    "# Manually increase coefficients if too low after fitting\n",
    "dIb[2][1] = dIb[2][1]*15\n",
    "\n",
    "scl_t = fit_obj.get_scl_t()\n",
    "print(\"The time scaling scl_t is\", scl_t) \n",
    "\n",
    "#%% Print the parameter values\n",
    "print('\\nImaster = ', I_master)\n",
    "print('Ivoltage = ', I_voltage)\n",
    "\n",
    "print('\\nDigital values for maximal conductances:')\n",
    "print('[gna, gk, gl] = ', dg, '\\n')\n",
    "\n",
    "print('Digital values for reversal potentials:')\n",
    "print('[Ena, Ek, El] = ', dE, '\\n')\n",
    "\n",
    "print('Digital values for gating variable kinetics:')\n",
    "print('alpha_m = ', dIb[0][0])\n",
    "print('beta_m = ', dIb[0][1], '\\n')\n",
    "print('alpha_h = ', dIb[1][0])\n",
    "print('beta_h = ', dIb[1][1], '\\n')\n",
    "print('alpha_n = ', dIb[2][0])\n",
    "print('beta_n = ', dIb[2][1], '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulating the fitted model\n",
    "\n",
    "Let's now try to use the parameters we have just calculated in order to define a new NeuroDyn object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I0 = 0e-6 # current in the original Hodgkin-Huxley model\n",
    "Iapp = lambda t : fit_obj.convert_I(I0) # scale the current based on scl_t, scl_v and C_HH/C_ND\n",
    "\n",
    "V_ref = 0.9\n",
    "\n",
    "ND = NeuroDynModel(dg, dE, dIb, V_ref, I_voltage, I_master)\n",
    "\n",
    "T = 0.02\n",
    "trange = (0, T)\n",
    "\n",
    "sol = ND.simulate(trange,[0.7,0,0,0],Iapp)\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel('t [s]')\n",
    "plt.ylabel('V [V]')\n",
    "plt.title('NeuroDyn simulation')\n",
    "plt.plot(sol.t, sol.y[0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compare this to the scaled Hodgkin-Huxley model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hh_neuron = HHModel(scl_v = 3, scl_t = scl_t, SI_units = True)\n",
    "\n",
    "I0 = 8e-6\n",
    "Iapp = lambda t: I0\n",
    "\n",
    "T = 0.02\n",
    "trange = (0, T)\n",
    "\n",
    "sol = hh_neuron.simulate(trange,[0,0,0,0],Iapp)\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('V')\n",
    "plt.title('Hodgkin-Huxley simulation')\n",
    "plt.plot(sol.t, sol.y[0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Passing analog values to NeuroDyn\n",
    "\n",
    "Instead of passing the digital parameter values `dg`, `dErev` and `dIb` we can also see how the NeuroDyn circuit would behave with the exact analog values obtained through the fitting process. This can allow us to get a better understanding of how the quantization error affects the circuit behavior.\n",
    "\n",
    "This can be easily done by passing a keyword argument `digital_values = False` when creating a NeuroDyn object. The object in this case expects analog values to be passed to `dg`, `dErev` and `dIb` variables instead of the standard integers between $0$ and $1023$.\n",
    "\n",
    "In order to get the analog values, we can call the corresponding `get_analog_Ib(weights)`, `get_analog_g(g)` and `get_analog_E(E)` instead of the digital methods that we have used previously.\n",
    "\n",
    "Let's see how this works using the results from our previous fitting procedure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ib = fit_obj.get_analog_Ib(weights)\n",
    "g = fit_obj.get_analog_g(g0)\n",
    "E = fit_obj.get_analog_E(E0)\n",
    "\n",
    "ND_analog = NeuroDynModel(g, E, Ib, V_ref, I_voltage, I_master, digital_values = False)\n",
    "\n",
    "I0 = 0e-6 # current in the original Hodgkin-Huxley model\n",
    "Iapp = lambda t : fit_obj.convert_I(I0) # scale the current based on scl_t and C_HH/C_ND\n",
    "\n",
    "T = 0.02\n",
    "trange = (0, T)\n",
    "\n",
    "sol = ND_analog.simulate(trange,[0.7,0,0,0],Iapp)\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('V')\n",
    "plt.title('NeuroDyn simulation')\n",
    "plt.plot(sol.t, sol.y[0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may notice in this case that even by increasing the applied current to the neuron, the circuit does not start regularly spiking. The main source of discrepancy that is left here is the non-linearity of the conductance relationship involving the $\\tanh$ input-output relationship. This non-linearity means that the same conductance values will lead to a somewhat different behavior for the same parameters, as observed by the lack of regular spiking here. Still, slight changes to the values may fix this - in this case, we can decrease the potassium conductance to *increase* the excitability of the neuron."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "g[1] = 0.6*g[1] # decrease potassium conductance\n",
    "\n",
    "ND_analog = NeuroDynModel(g, E, Ib, V_ref, I_voltage, I_master, digital_values = False)\n",
    "\n",
    "I0 = 7e-6 # current in the original Hodgkin-Huxley model\n",
    "Iapp = lambda t : fit_obj.convert_I(I0) # scale the current based on scl_t and C_HH/C_ND\n",
    "\n",
    "T = 0.02\n",
    "trange = (0, T)\n",
    "\n",
    "sol = ND_analog.simulate(trange,[0.7,0,0,0],Iapp)\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('V')\n",
    "plt.title('NeuroDyn simulation')\n",
    "plt.plot(sol.t, sol.y[0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modelling variability\n",
    "\n",
    "The model we have introduced in the previous sections of course presents an idealized picture compared to the physical realization in hardware. The difference stems from the inherent noise and variability present in electrical circuits, which is especially true for subthreshold analog techniques that we are utilizing within NeuroDyn. Due to this effect, setting the same set of external parameters to different NeuroDyn chips will yield different behavior and appropriate compensation techniques will need to be applied to tackle the variability.\n",
    "\n",
    "The basic source of variability is the *transistor mismatch*: applying the same gate voltage to different transistors within the circuit will generate different drain currents. We can model this variability in a simple way, by considering that these currents will obey a Gaussian distribution:\n",
    "\n",
    "\\begin{equation}\n",
    "I_{out} = I_0 \\cdot (1 + \\mathcal{N}(0, \\sigma^2))\n",
    "\\end{equation}\n",
    "\n",
    "where $I_0$ is the nominal output current and $\\mathcal{N}$ is a Gaussian distribution with the standard deviation $\\sigma$. A rule of thumb that we use here is to set $\\sigma \\approx 0.15$, which is the default value in the model.\n",
    "\n",
    "### Perturb function\n",
    "\n",
    "We can apply the transistor current variability to our NeuroDyn model by calling the function `perturb(sigma)`. This function will apply the Gaussian distribution variability to the main parameters affected: maximal conductances, reversal potentials, sigmoid weights of the $\\alpha(V)$ and $\\beta(V)$ functions and the gating variable exponents (e.g. $4$ in $m^4$ for the sodium current).\n",
    "\n",
    "Let's see how our previously defined model behaves after we add the variability (make sure you ran the previous code where we defined the initial digital NeuroDyn model):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I0 = 0e-6\n",
    "Iapp = lambda t : fit_obj.convert_I(I0)\n",
    "\n",
    "ND.perturb() # leaving out the keyword argument sets default sigma = 0.15\n",
    "\n",
    "sol = ND.simulate(trange,[0.7,0,0,0],Iapp)\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('V')\n",
    "plt.title('Perturbed NeuroDyn simulation')\n",
    "plt.plot(sol.t, sol.y[0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can rerun this code multiple times and see how each instance behaves. You should notice that each time we generate an instance of the neuron the behavior changes - the spike frequency and height may change, and sometimes we do not get spiking at all!\n",
    "\n",
    "You can try compensating for this by changing some parameters of the model to make the neuron excitable again. For example, maybe due to the perturbation the threshold voltage has increased and a simple increase in the applied current can get it back to spiking. Alternatively we may need to increase the maximal conductance of the sodium current relative to the maximal conductance of the potassium current in order to increase the excitability of the model. Sometimes though the perturbation might just be too large for these simple techniques to work!\n",
    "\n",
    "Let's try this out. First, let's check the current values of these parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"The current value of applied current is \", I0)\n",
    "print(\"The current values for the maximal conductances are\", dg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's see if our simple compensation technique can work:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First try changing the applied current\n",
    "I0 = 0e-6\n",
    "Iapp = lambda t : fit_obj.convert_I(I0)\n",
    "\n",
    "# If changing applied current didn't help, try adjusting the maximal conductances\n",
    "# dg = [1023, 307, 3]\n",
    "dg_comp = np.array([1023, 307, 3])\n",
    "ND.update_dg(dg_comp)\n",
    "\n",
    "sol = ND.simulate(trange,[0.7,0,0,0],Iapp)\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('V')\n",
    "plt.title('Compensated NeuroDyn simulation')\n",
    "plt.plot(sol.t, sol.y[0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can try changing other parameters as well and get a feel for how hard it is to compensate for the internal mismatch -  this inherent variability is common to both biological neurons and neuromorphic circuits and is the main challenge to designing and controlling bioinspired circuits."
   ]
  },
  {
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}