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 {
    LOCAL kf0_, kb0_
    kf0_ = a
    kb0_ = b
    h' = (-1*(kf0_*h-kb0_*m))
    m' = (1*(kf0_*h-kb0_*m))
}

Ex. 2

Annihilation reaction statement

[4]:
ex2 = """
STATE {
    x
}
KINETIC kin {
    ~ x -> (a)
}
"""
print(run_kinetic_visitor_and_return_derivative(ex2))
DERIVATIVE kin {
    LOCAL kf0_
    kf0_ = a
    x' = (-1*(kf0_*x))
}

Ex. 3

<< reaction statement

[5]:
ex3 = """
STATE {
    x
}
KINETIC kin {
    ~ x << (a)
}
"""
print(run_kinetic_visitor_and_return_derivative(ex3))
DERIVATIVE kin {
    LOCAL source0_
    source0_ = a
    x' = (source0_)
}

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 {
    LOCAL source0_, kf1_
    source0_ = a
    kf1_ = b
    x' = (source0_)+(-1*(kf1_*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 {
    LOCAL kf0_, kb0_, kf1_
    kf0_ = a
    kb0_ = b
    f = kf0_*x-kb0_*y
    kf1_ = c
    g = kf1_*z
    h = 0
    x' = (-1*(kf0_*x-kb0_*y))
    y' = (1*(kf0_*x-kb0_*y))
    z' = (-1*(kf1_*z))
}