{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NMODL CONDUCTANCE\n",
    "\n",
    "This notebook described the `CONDUCTANCE` keyword in NEURON, how it is implemented in NMODL, and shows some examples of the output generated by NMODL in different situations.\n",
    "\n",
    "For a more general tutorial on using the NMODL python interface, please see the [tutorial notebook](nmodl-python-tutorial.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Introduction\n",
    "\n",
    "Motivation:\n",
    "  - during a NEURON simulation, a number of currents $I$ may be generated\n",
    "  - at each time step, for each current $I$ the corresponding conductance $dI/dV$ needs to be calculated\n",
    "  - by default in NEURON this is approximated by a forwards difference: $dI/dV \\simeq (I(v+\\Delta v)-I(v))/\\Delta v$, with $\\Delta v = 0.001$\n",
    "      - this introduces an $\\mathcal{O}(\\Delta v)$ numerical error\n",
    "      - it also requires two current calculations, which may be computationally inefficient\n",
    "\n",
    "Solution:\n",
    "  - the `CONDUCTANCE` keyword was added to the NMODL language\n",
    "  - this allows the user to manually specify the analytic expression for the conductance in the MOD file\n",
    "  - during the simulation, instead of the numerical differentiation, the user supplied expression is used\n",
    "    - this solves the problem, but requires additional effort from the user\n",
    "    - it also opens up room for user error: an incorrect expression will still run but the results will not be correct\n",
    "\n",
    "SymPy improvement:\n",
    "  - the currents in the mod file are differentiated analytically using SymPy\n",
    "  - the corresponding `CONDUCTANCE` statements are generated automatically\n",
    "    - no additional input required from the user\n",
    "    - avoids the possibility of user error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Implementation\n",
    "\n",
    "The `SympyConductanceVisitor` is defined in [src/visitors/sympy_conductance_visitor.hpp](https://github.com/BlueBrain/nmodl/blob/master/src/visitors/sympy_conductance_visitor.hpp), and it makes use of the python function [differentiate2c](https://github.com/BlueBrain/nmodl/blob/master/nmodl/ode.py#L377) to perform the analytic differentation using the [SymPy](https://www.sympy.org/en/index.html) symbolic math Python library.\n",
    "\n",
    "* For each ion write statement $i = \\dots$ in the BREAKPOINT block\n",
    "    * Differentiate to find the conductance $g_i=di/dv$\n",
    "    * If this $g_i$ coincides with an existing variable, e.g. $g$, add to BREAKPOINT the statement:\n",
    "        * CONDUCTANCE g USEION ion_name\n",
    "    * If not, also need to declare and asign a variable for the calculated conductance:\n",
    "        * LOCAL g_i_0\n",
    "        * CONDUCTANCE g_i_0 USEION ion_name\n",
    "        * g_i_0 = ...\n",
    "    * But if there is an existing CONDUCTANCE statement, then do not modify it\n",
    "\n",
    "\n",
    "* It may be the case that a variable in the write statement $i = \\dots$ itself depends on $v$, so to take this into account:\n",
    "    * an inlining visitor is first ran, after which all variable assignments occur within the BREAKPOINT block\n",
    "    * each preceeding expression is analysed in reverse order for $v$ dependence \n",
    "    * if it depends on $v$, the rhs of the expression is substituted for the lhs in all following statements\n",
    "    * the end result is a (complicated) expression $i = ...$ where all v dependence is explicit\n",
    "    * this is then differentiated w.r.t $v$ to give the conductance\n",
    "    * it then checks if this expression is equivalent to an existing variable\n",
    "    * for this step it is necessary to also substitute all non-$v$-dependent expressions on both sides & simplify\n",
    "\n",
    "#### Implementation Tests\n",
    "\n",
    "  - The unit tests may be helpful to understand what these functions are doing\n",
    "  - `SympyConductanceVisitor` tests are located in [test/visitor/visitor.cpp](https://github.com/BlueBrain/nmodl/blob/master/test/visitor/visitor.cpp#L3261)\n",
    "  - `differentiate2c` tests are located in [test/ode/test_ode.py](https://github.com/BlueBrain/nmodl/blob/master/test/ode/test_ode.py#L56)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Examples\n",
    "Here are some examples of generated CONDUCTANCE statements for a variety of sample mod files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "! pip install nmodl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nmodl.dsl as nmodl\n",
    "\n",
    "\n",
    "def run_conductance_visitor_and_return_breakpoint(mod_string):\n",
    "    # parse NMDOL file (supplied as a string) into AST\n",
    "    driver = nmodl.NmodlDriver()\n",
    "    AST = driver.parse_string(mod_string)\n",
    "    # run SymtabVisitor to generate Symbol Table\n",
    "    nmodl.symtab.SymtabVisitor().visit_program(AST)\n",
    "    # constant folding, inlining & local variable renaming passes\n",
    "    nmodl.visitor.ConstantFolderVisitor().visit_program(AST)\n",
    "    nmodl.visitor.InlineVisitor().visit_program(AST)\n",
    "    nmodl.visitor.LocalVarRenameVisitor().visit_program(AST)\n",
    "\n",
    "    # run CONDUCTANCE visitor\n",
    "    nmodl.visitor.SympyConductanceVisitor().visit_program(AST)\n",
    "    # return new BREAKPOINT block\n",
    "    return nmodl.to_nmodl(\n",
    "        nmodl.visitor.AstLookupVisitor().lookup(\n",
    "            AST, nmodl.ast.AstNodeType.BREAKPOINT_BLOCK\n",
    "        )[0]\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Ex. 1\n",
    "  - simple USEION statement, conductance equal to existing variable\n",
    "    - add CONDUCTANCE statement using existing variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BREAKPOINT {\n",
      "    CONDUCTANCE gna USEION na\n",
      "    ina = gna*(v-ena)\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "ex1 = \"\"\"\n",
    "NEURON {\n",
    "    USEION na READ ena WRITE ina\n",
    "    RANGE gna\n",
    "}\n",
    "BREAKPOINT {\n",
    "    ina = gna*(v - ena)\n",
    "}\n",
    "\"\"\"\n",
    "print(run_conductance_visitor_and_return_breakpoint(ex1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Ex. 2\n",
    "  - simple USEION statement, conductance not equal to existing variable\n",
    "    - declare new local variable\n",
    "    - assign conductance to it\n",
    "    - add CONDUCTANCE statement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BREAKPOINT {\n",
      "    LOCAL g_na_0\n",
      "    CONDUCTANCE g_na_0 USEION na\n",
      "    g_na_0 = 0.1*gna\n",
      "    ina = 0.1*gna*(v-ena)\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "ex2 = \"\"\"\n",
    "NEURON {\n",
    "    USEION na READ ena WRITE ina\n",
    "    RANGE gna\n",
    "}\n",
    "BREAKPOINT {\n",
    "    ina = 0.1*gna*(v - ena)\n",
    "}\n",
    "\"\"\"\n",
    "print(run_conductance_visitor_and_return_breakpoint(ex2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Ex. 3\n",
    "  - simple NONSPECIFIC_CURRENT statement, conductance equal to existing variable\n",
    "    - add CONDUCTANCE statement using existing variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BREAKPOINT {\n",
      "    CONDUCTANCE g\n",
      "    i = g*v\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "ex3 = \"\"\"\n",
    "NEURON {\n",
    "    NONSPECIFIC_CURRENT i\n",
    "    RANGE g\n",
    "}\n",
    "BREAKPOINT {\n",
    "    i = g*v\n",
    "}\n",
    "\"\"\"\n",
    "print(run_conductance_visitor_and_return_breakpoint(ex3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Ex. 4\n",
    "  - non-linear NONSPECIFIC_CURRENT statement, conductance not equal to existing variable\n",
    "    - declare new local variable\n",
    "    - assign conductance to it\n",
    "    - add CONDUCTANCE statement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BREAKPOINT {\n",
      "    LOCAL g__0\n",
      "    CONDUCTANCE g__0\n",
      "    g__0 = g+2*v\n",
      "    i = g*v+v*v\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "ex4 = \"\"\"\n",
    "NEURON {\n",
    "    NONSPECIFIC_CURRENT i\n",
    "    RANGE g\n",
    "}\n",
    "BREAKPOINT {\n",
    "    i = g*v + v*v\n",
    "}\n",
    "\"\"\"\n",
    "print(run_conductance_visitor_and_return_breakpoint(ex4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Ex. 5\n",
    "  - several current statements, conductance equal to existing variables\n",
    "    - add CONDUCTANCE statements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BREAKPOINT {\n",
      "    CONDUCTANCE gl\n",
      "    CONDUCTANCE gk USEION k\n",
      "    CONDUCTANCE gna USEION na\n",
      "    gna = gnabar*m*m*m*h\n",
      "    ina = gna*(v-ena)\n",
      "    gk = gkbar*n*n*n*n\n",
      "    ik = gk*(v-ek)\n",
      "    il = gl*(v-el)\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "ex5 = \"\"\"\n",
    "NEURON {\n",
    "    USEION na READ ena WRITE ina\n",
    "    USEION k READ ek WRITE ik\n",
    "    NONSPECIFIC_CURRENT il\n",
    "    RANGE gnabar, gkbar, gl, el, gna, gk\n",
    "}\n",
    "STATE {\n",
    "    m n h\n",
    "}\n",
    "BREAKPOINT {\n",
    "    gna = gnabar*m*m*m*h\n",
    "    ina = gna*(v - ena)\n",
    "    gk = gkbar*n*n*n*n\n",
    "    ik = gk*(v - ek)\n",
    "    il = gl*(v - el)\n",
    "}\n",
    "\"\"\"\n",
    "print(run_conductance_visitor_and_return_breakpoint(ex5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Ex. 6\n",
    "  - current contains variables that depend on $v$, conductance equal to existing variable `x3`\n",
    "    - substitute all variables with $v$-dependence, differentiate to find conductance\n",
    "    - compare result to each existing variable (after substituting all preceeding declarations on both sides)\n",
    "    - identify that expression for conductance is equivalent to `x3`\n",
    "    - add CONDUCTANCE statement using this existing variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BREAKPOINT {\n",
      "    CONDUCTANCE x3 USEION na\n",
      "    x1 = 0.2+3*v\n",
      "    x2 = v*v\n",
      "    x3 = 3*v*v+5*v-1.3\n",
      "    gna = x1+x2\n",
      "    ina = gna*(v-0.5)\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "ex6 = \"\"\"\n",
    "NEURON {\n",
    "    USEION na READ ena WRITE ina\n",
    "    RANGE gna, x1, x2, x3\n",
    "}\n",
    "BREAKPOINT {\n",
    "    x1 = 0.2+3*v\n",
    "    x2 = v*v\n",
    "    x3 = 3*v*v+5*v-1.3\n",
    "    gna = x1 + x2\n",
    "    ina = gna*(v-0.5)\n",
    "}\n",
    "\"\"\"\n",
    "print(run_conductance_visitor_and_return_breakpoint(ex6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Ex. 7\n",
    "  - current contains variables that depend on $v$, conductance not equal to existing variable\n",
    "    - substitute all variables with $v$-dependence, differentiate to find conductance\n",
    "    - compare result to each existing variable, no equivalent expression found\n",
    "    - declare new local variable, assign conductance to it\n",
    "    - add CONDUCTANCE statement using this new variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BREAKPOINT {\n",
      "    LOCAL g_na_0\n",
      "    CONDUCTANCE g_na_0 USEION na\n",
      "    g_na_0 = 3*pow(v, 2)+5*v-1.3\n",
      "    x1 = 0.2+3*v\n",
      "    x2 = v*v\n",
      "    gna = x1+x2\n",
      "    ina = gna*(v-0.5)\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "ex7 = \"\"\"\n",
    "NEURON {\n",
    "    USEION na READ ena WRITE ina\n",
    "    RANGE gna, x1, x2\n",
    "}\n",
    "BREAKPOINT {\n",
    "    x1 = 0.2+3*v\n",
    "    x2 = v*v\n",
    "    gna = x1 + x2\n",
    "    ina = gna*(v-0.5)\n",
    "}\n",
    "\"\"\"\n",
    "print(run_conductance_visitor_and_return_breakpoint(ex7))"
   ]
  },
  {
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}