{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### NMODL matexp solver\n", "\n", "This notebook describes the implementation of the `matexp` solver, which solves the systems of ODEs defined in `KINETIC` blocks when the ODEs are *linear* and *coupled*.\n", "\n", "For a more general tutorial on using the NMODL python interface, please see the [tutorial notebook](nmodl-python-tutorial.ipynb).\n", "\n", "***\n", "\n", "#### Mathematics\n", "ODEs are linear if they can be written in the form $\\frac{dx}{dt} = A \\cdot x$ where $x$ is the state vector of the model and $A$ is the Jacobian matrix. The analytic solution to this initial value problem is $x = e^{A \\Delta t} \\cdot x_0$. The matrix $e^{A \\Delta t}$ is known as the *propagator* matrix and it advances the state of the model by one time step. Since this is the exact solution, it is compatible with the Crank−Nicholson method which requires at least second-order accuracy.\n", "\n", "***\n", "\n", "#### Matrix Exponential Function\n", "The exponential function is defined over matrices as well as scalars by its Taylor series expansion $e^x = \\sum_{n=0}^{\\infty} \\frac{x^n}{n!}$. However, computing the exponential function by evaluating its Taylor series is slow and numerically unstable. There are many faster and more accurate methods of calculating the matrix exponential function. The NEURON simulator uses the [Eigen](https://eigen.tuxfamily.org) library for its advanced linear algebra functions, including for the matrix exponential function. For this function, the Eigen library implements the algorithm introduced in the publication:\n", "\n", "> \"The scaling and squaring method for the matrix exponential revisited,\" \\\n", "> By Nicholas J. Higham, \\\n", "> SIAM J. Matrix Anal. Applic., 26:1179–1193, 2005. \\\n", "> https://doi.org/10.1137/04061101X\n", "\n", "***\n", "\n", "#### Implementation\n", "The `MatexpVisitor` implements the solver method `matexp`.\n", "\n", "The visitor does the following:\n", "* make lists of all `SOLVE` statements, `STEADYSTATE` statements, and `KINETIC` blocks in the program\n", "* for each `SOLVE block STEADYSTATE matexp` statement\n", " * find the corresponding `KINETIC` block and convert it into a `MATEXP_BLOCK`\n", " * replace the `SOLVE ... STEADYSTATE` statement with the `MATEXP_BLOCK`\n", "* for each `SOLVE block METHOD matexp` statement\n", " * find the corresponding `KINETIC` block and convert it into a `MATEXP_BLOCK`\n", " * append the `MATEXP_BLOCK` to the program\n", "* remove the solved `KINETIC` blocks\n", "\n", "MATEXP_BLOCKs have the following fields:\n", "* steadystate (boolean flag)\n", "* kinetic block\n", "* list of conserve statements, which the visitor removed from the kinetic block\n", "\n", "During the code generation phase, all MATEXP_BLOCKs that were appended to the program are moved into the state update function. Then MATEXP_BLOCKs generate code in-place either in the initialization function or in the state update function. MATEXP_BLOCKs generate the following code:\n", "* define and zero-init the Jacobian matrix using the Eigen library\n", "* insert the body of the kinetic block\n", " * replace reaction statements with equivalent assignments to the Jacobian matrix\n", " * For example `~ X <-> Y (A, B)` would become:\n", " * `J[X, X] -= A`\n", " * `J[Y, X] += A`\n", " * `J[Y, Y] -= B`\n", " * `J[X, Y] += B`\n", "* initialize the state vector\n", "* define the solution vector as `solution = jacobian.exp() * state`\n", "* conserve the sum of the states, if applicable\n", " * multiply each state variable by the expected sum divided by the actual sum\n", "* assign the solution vector to the state variables\n", "\n", "***\n", "\n", "#### Tests Cases\n", "The unit tests may be helpful to understand what these functions are doing\n", " - `MatexpVisitor` tests are located in [test/nmodl/transpiler/unit/visitor/matexp.cpp](https://github.com/neuronsimulator/nrn/blob/master/test/nmodl/transpiler/unit/visitor/matexp.cpp) and have the tag `[matexp]`\n", "\n", "The usecase tests compare several forumations of the same Hodgkin-Huxley model\n", "- Located in [test/nmodl/transpiler/usecases/matexp](https://github.com/neuronsimulator/nrn/blob/master/test/nmodl/transpiler/usecases/matexp)\n", "- hh_cnexp.mod is the standard HH model\n", " - solved by `cnexp`\n", "- hh_matexp.mod is identical to hh_cnexp.mod except each gate (\"M\", \"H\", and \"N\") is a two state Markov model\n", " - solved by `matexp`\n", "- hhkin.mod is modified to enumerate each all 12 states of each ion channel\n", " - solved by `sparse`\n", "- kkkinmatexp.mod is identical to hhkin.mod except it is solved by `matexp`\n", "\n", "***\n", "\n", "#### Examples" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "%%capture\n", "! pip install neuron" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "import neuron.nmodl.dsl as nmodl\n", "\n", "\n", "def run_matexp_solver(mod_string):\n", " # Parse NMDOL file (supplied as a string) into AST\n", " driver = nmodl.NmodlDriver()\n", " AST = driver.parse_string(mod_string)\n", "\n", " # Run SymtabVisitor to generate Symbol Table\n", " nmodl.symtab.SymtabVisitor().visit_program(AST)\n", "\n", " # Run matexp solver\n", " nmodl.visitor.MatexpVisitor().visit_program(AST)\n", "\n", " # Return the solution as a string\n", " return nmodl.to_nmodl(AST)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Ex. 1\n", "Linear kinetic model with conserved sum of states" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "STATE {\n", " A\n", " B\n", "}\n", "\n", "BREAKPOINT {\n", "}\n", "\n", "MATEXP_SOLVE (0) {\n", " nmodl_eigen_j[0] = nmodl_eigen_j[0]-(0.123)*dt\n", " nmodl_eigen_j[1] = nmodl_eigen_j[1]+(0.123)*dt\n", " nmodl_eigen_j[3] = nmodl_eigen_j[3]-(0.456)*dt\n", " nmodl_eigen_j[2] = nmodl_eigen_j[2]+(0.456)*dt\n", "} CONSERVE = 0.789\n", "\n" ] } ], "source": [ "ex1 = \"\"\"\n", "STATE {\n", " A\n", " B\n", "}\n", "BREAKPOINT {\n", " SOLVE states METHOD matexp\n", "}\n", "KINETIC states {\n", " ~ A <-> B (0.123, 0.456)\n", " CONSERVE A + B = 0.789\n", "}\n", "\"\"\"\n", "print(run_matexp_solver(ex1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##### Ex. 2\n", "Non-linear ODEs\n", "\n", "The matexp solver checks for linear equations and raises exception if it finds any non-linearities.\n", "It asserts that each reaction has exactly one reactant and one product.\n", "This means that the following statements will cause errors:\n", "\n", "```\n", "KINETIC nonlinear_errors {\n", " ~ A <-> B + C (kf, kr) : Multiple products!\n", " ~ 2 A <-> C (kf, kr) : Multiple reactants!\n", " ~ A << (kf) : Zero reactants!\n", "}\n", "```\n", "\n", "However the solver does not check that the kinetic rates are independent from state variables, and so the following non-linear equation will be accepted by the program and solved inaccurately.\n", "\n", "```\n", "KINETIC nonlinear_silent_bug {\n", " ~ A <-> B (func(A), kr)\n", "}\n", "```\n" ] } ], "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.10.12" } }, "nbformat": 4, "nbformat_minor": 4 }