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 ***

[ ]: