NMODL integration of ODEs

This is an overview of - the different ways a user can specify the equations that define the system they want to simulate in the MOD file - how these equations can be related to each other - how these equations are solved in NMODL

The user can specify information about the system in a variety of ways: - A high level way to describe a system is to specify a Mass Action Kinetic scheme of reaction equations in a KINETIC block - Alternatively a system of ODEs can be specified in a DERIVATIVE block (any kinetic scheme can also be written as a system of ODEs) - Finally a system of linear/non-linear algebraic equations can be specified in a LINEAR/NONLINEAR block (a numerical integration scheme, such as backwards Euler, transforms a system of ODEs into a system of linear or non-linear equations)

To reduce duplication of functionality for dealing with these related systems, we implement a hierarchy of transformations: - KINETIC blocks of reaction statements are translated to DERIVATIVE blocks of equivalent ODE systems using the law of Mass Action - DERIVATIVE blocks of ODEs are translated to (NON)LINEAR blocks of algebraic equations using a numerical integration scheme

After these transformations we are left with only LINEAR/NONLINEAR blocks that are then solved numerically (by Gaussian Elimination, LU factorization or Newton iteration)

KINETIC block

  • Mass Action kinetics: a set of reaction equations with associated reaction rates

  • converted to a DERIVATIVE blocking containing an equivalent system of ODEs using the law of Mass Action

  • see the nmodl-kinetic-schemes notebook for more details

DERIVATIVE block

  • system of ODEs & associated solve method: cnexp, sparse, derivimplicit or euler

  • cnexp

    • applicable if ODEs are linear & independent

    • exact analytic integration

    • see the nmodl-sympy-solver-cnexp notebook for more details

  • sparse

    • applicable if ODEs are linear & coupled

    • backwards Euler numerical integration scheme: f(t+\Delta t) = f(t) + dt f'(t+ \Delta t)

    • results in a linear algebraic system to solve

    • numerically stable

    • \mathcal{O}(\Delta t) integration error

    • see the nmodl-sympy-solver-sparse notebook for more details

  • derivimplcit

    • always applicable

    • backwards Euler numerical integration scheme: f(t+\Delta t) = f(t) + dt f'(t+ \Delta t)

    • results in a non-linear algebraic system to solve

    • numerically stable

    • \mathcal{O}(\Delta t) integration error

    • see the nmodl-sympy-solver-sparse notebook for more details

  • euler

    • always applicable

    • forwards Euler numerical integration scheme: f(t+\Delta t) = f(t) + \Delta t f'(t)

    • numerically unstable

    • \mathcal{O}(\Delta t) integration error

    • not recommended due to instability of scheme

LINEAR block

  • system of linear algebraic equations

  • for small systems (N<=3)

    • solve at compile time by Gaussian elimination

  • for larger systems

    • solve at run-time by LU factorization with partial pivoting

  • see the nmodl-linear-solver notebook for more details

NONLINEAR block

  • system of non-linear algebraic equations

  • solve by Newton iteration

    • construct F(X), with Jacobian J(X)=\frac{\partial F_i}{\partial X_j}

    • such that desired solution X^* satisfies condition F(X^*) = 0

    • iterative solution given by X_{n+1} = X_n + J(X_n)^{-1} F(X_n)

  • see the nmodl-nonlinear-solver notebook for more details ***

[ ]: