NMODL Kinetic Scheme

This notebook describes the reaction kinetics & mass action laws that apply within KINETIC blocks, and the implementation of the KineticBlockVisitor in NMODL which transforms KINETIC blocks into DERIVATIVE blocks containing an equivalent system of ODEs.

For a higher level overview of the approach to solving ODEs in NMODL, please see the nmodl-odes-overview notebook.

For a more general tutorial on using the NMODL python interface, please see the tutorial notebook.

Reaction Kinetics

We consider a set of reaction species $A_j$ , with corresponding Molar concentrations $Y_j$
They react according to a set of reaction equations:

$\sum_j \nu_{ij}^L A_j \overset{k_i}{\rightarrow} \sum_j \nu_{ij}^R A_j$

where - $k_i$ is the rate coefficient - $\nu_{ij}^L$ , $\nu_{ij}^R$ are stoichiometric coefficients - must be positive integers (including zero)

Law of Mass Action

This allows us to convert these reaction equations to a set of ODEs

$\frac{dY_j}{dt} = \sum_i \Delta \nu_{ij} r_i$

where $\Delta \nu_{ij} = \nu_{ij}^R - \nu_{ij}^L$ , and

$r_i = k_i \prod_j Y_j^{\nu_{ij}^R}$

KINETIC block format

A reaction equation is specifed in the mod file as

~ A0 + 3A1 + 2A2 + ... <-> 2A0 + A1 + ... (kf, kb)

where - A0 etc are the species $A_j$ - the integer preceeding a species (with or without a space) is the corresponding stochiometric coefficient $\nu_{ij}$ (implicitly 1 if not specified) - kf is the forwards reaction rate $k^{(f)}_j$ - kb is the backwards reaction rate $k^{(b)}_j$ , i.e. the reaction rate for the same reaction with the LHS and RHS exchanged *** We can convert these statements to a system of ODEs using the law of Mass Action:

$\frac{dY_j}{dt} = \sum_i \Delta \nu_{ij} (r^{(f)}_i - r^{(b)}_i)$

where $\Delta \nu_{ij} = \nu_{ij}^R - \nu_{ij}^L$ , and

$\begin{aligned} r^{(f)}_i &= k^{(f)}_i \prod_j Y_j^{\nu_{ij}^{L}} \\ r^{(b)}_i &= k^{(b)}_i \prod_j Y_j^{\nu_{ij}^{R}} \end{aligned}$

***

Other types of reaction statement

There is also have a reaction statement of the form

~ h << (a)

where the LHS must be a state variable, and the RHS is an expression inside parentheses.

The meaning of this statement is to add a to the differential equation for h, i.e.

$\frac{dh}{dt} += a$

*** Finally there is a statement of the form

~ x + 2y + ... -> (a)

which is a one-way reaction statement with no reaction products.

This is just syntactic sugar for a special case of the standard <-> reaction equation, where the backwards rate is set to zero and there are no states on the RHS of the reaction, so the above is equivalent to

~ x + 2y + ... <-> (a, 0)

f_flux / b_flux variables

Within the KINETIC block in the MOD file, the user can make use of the f_flux and b_flux variables, which refer to the forwards $r^{(f)}$ and backwards $r^{(b)}$ fluxes of the preceeding reaction statement.

The KineticBlockVisitor substitutes the current expression for these fluxes for these variables within the KINETIC block.

If these variables are referenced before a reaction statement then they are assumed to be zero.

[TODO] CONSERVE

Converse statement allows specification of a conservation law, e.g.

CONSERVE h + m + z = 1

In NEURON, the ODE for the last state variable on the rhs of this expression is replaced with this algebraic expression, so in this case instead of replacing $z' = ...$ with the forwards or backwards Euler equation, it would be replaced with $z = 1 - h - m$

In order to be consistent with NEURON, in particular the way STEADYSTATE is implemented, we should do this in the same way. *** #### Implementation Tests

The unit tests may be helpful to understand what these functions are doing
KineticBlockVisitor tests are located in test/visitor/visitor.cpp ***

Examples

[1]:

%%capture
! pip install nmodl

[2]:

import nmodl.dsl as nmodl


def run_kinetic_visitor_and_return_derivative(mod_string):
    # parse NMDOL file (supplied as a string) into AST
    driver = nmodl.NmodlDriver()
    AST = driver.parse_string(mod_string)
    # run SymtabVisitor to generate Symbol Table
    nmodl.symtab.SymtabVisitor().visit_program(AST)
    # constant folding, inlining & local variable renaming passes
    nmodl.visitor.ConstantFolderVisitor().visit_program(AST)
    nmodl.visitor.InlineVisitor().visit_program(AST)
    nmodl.visitor.LocalVarRenameVisitor().visit_program(AST)
    # run KINETIC block visitor
    nmodl.visitor.KineticBlockVisitor().visit_program(AST)
    # return new DERIVATIVE block
    return nmodl.to_nmodl(
        nmodl.visitor.AstLookupVisitor().lookup(
            AST, nmodl.ast.AstNodeType.DERIVATIVE_BLOCK
        )[0]
    )

Ex. 1

Given the following KINETIC statement

~ h <-> m (a,b)

We have 2 state variables, and 1 reaction statement, so - $i = \{0\}$ - $j = \{0, 1\}$

i.e. - the state vector $Y$ has 2 elements

$Y = \left( \begin{aligned} h \\ m \end{aligned} \right)$

the stoichiometric coefficients are 1x2 matrices

$\nu^L = \left( \begin{aligned} 1 && 0 \end{aligned} \right)$

$\nu^R = \left( \begin{aligned} 0 && 1 \end{aligned} \right)$

the rate vectors contain 1 element:

$k^{(f)} = a$

$k^{(b)} = b$

Using these we can construct the forwards and blackwards fluxes:

$r^{(f)} = a h$

$r^{(b)} = b m$

and finally we find the ODEs in matrix form:

$\frac{dY}{dt} = \left( \begin{aligned} -1 && 1 \end{aligned} \right) (ah - bm)$

which in terms of the state variables can be written:

$\begin{aligned} \frac{dh}{dt} &= bm - ah \\ \frac{dm}{dt} &= ah - bm \end{aligned}$

[3]:

ex1 = """
STATE {
    h m
}
KINETIC kin {
    ~ h <-> m (a,b)
}
"""
print(run_kinetic_visitor_and_return_derivative(ex1))

DERIVATIVE kin {
    h' = (-1*(a*h-b*m))
    m' = (1*(a*h-b*m))
}

Ex. 2

Annihilation reaction statement

[4]:

ex2 = """
STATE {
    x
}
KINETIC kin {
    ~ x -> (a)
}
"""
print(run_kinetic_visitor_and_return_derivative(ex2))

DERIVATIVE kin {
    x' = (-1*(a*x))
}

Ex. 3

<< reaction statement

[5]:

ex3 = """
STATE {
    x
}
KINETIC kin {
    ~ x << (a)
}
"""
print(run_kinetic_visitor_and_return_derivative(ex3))

DERIVATIVE kin {
    x' = (a)
}

Ex. 4

Annihilation & << reaction statement for the same state variable

[6]:

ex4 = """
STATE {
    x
}
KINETIC kin {
    ~ x << (a)
    ~ x -> (b)
}
"""
print(run_kinetic_visitor_and_return_derivative(ex4))

DERIVATIVE kin {
    x' = (a)+(-1*(b*x))
}

Ex. 5

Reaction statements and use of f_flux, b_flux variables

[7]:

ex5 = """
STATE {
    x y z
}
KINETIC kin {
    ~ x <-> y (a,b)
    f = f_flux - b_flux
    ~ z -> (c)
    g = f_flux
    h = b_flux
}
"""
print(run_kinetic_visitor_and_return_derivative(ex5))

DERIVATIVE kin {
    f = a*x-b*y
    g = c*z
    h = 0
    x' = (-1*(a*x-b*y))
    y' = (1*(a*x-b*y))
    z' = (-1*(c*z))
}