How to use NMODL to extend Arbor’s repertoire of Ion Channels

NMODL is Arbor’s way of expressing ion channel dynamics, which in turn can be added to cable cells via the Cable cell decoration interface. This tutorial will guide you through to create such an ion channel from scratch. Note that there are already large selections of channels built into Arbor in addition to online databases of ready-to-use models. Please check there first and only if you find nothing to use or adapt, try to build your own. If you want to learn how to use ion channels, consider other tutorials first, such as A simple single-cell recipe and the dynamic catalogue example.

Introducing NMODL

NMODL is a domain specific language (DSL) for design electro-chemical dynamics deriving from the earlier MODL language. It is used by Arbor, Neuron, and CoreNeuron. Although each software has its own dialect, the core ideas are identical. Unfortunately, documentation for both NMODL and MODL is scarce and outdated. This tutorial aims to help you become proficient in writing and adapting ion channels for Arbor. Note that while we focus on density mechanisms here, there is a host of behaviours customisable in Arbor through NMODL in addition to wholly new functionality, examples include gap junctions, synapses, and voltage clamps. NMODL can use stochastic differential equations and modify ion concentrations.

The running example: Hodgkin-Huxley

During this tutorial, we will re-implement the classic Hodgkin-Huxley ion channel model. You can use it directly in Arbor as hh from the default catalogue, but for this tutorial, we pretend that it is an entirely new model.

We begin by setting up a simple model using the default hh model

#!/usr/bin/env python3

import arbor as A
from arbor import units as U
import matplotlib.pyplot as plt

# Create a single segment morphology
tree = A.segment_tree()
tree.append(A.mnpos, (-3, 0, 0, 3), (3, 0, 0, 3), tag=1)

# Create (almost empty) decor
decor = (
    A.decor()
    .paint("(all)", A.density("hh"))
    .place("(location 0 0.5)", A.iclamp(10 * U.ms, 2 * U.ms, 0.8 * U.nA), "iclamp")
)

# Run the model, extracting the membrane voltage
model = A.single_cell_model(A.cable_cell(tree, decor))
model.probe("voltage", "(location 0 0.5)", tag="Um", frequency=10 * U.kHz)
model.run(tfinal=30 * U.ms)

# Create a basic plot
fg, ax = plt.subplots()
ax.plot(model.traces[0].time, model.traces[0].value)
ax.set_xlabel("t/ms")
ax.set_ylabel("U/mV")
plt.savefig("hh-01.svg")

Store it in step-01.py and — once run — is should produce a plot like this:

../_images/hh-01.svg

You can find all steps in the python/example/hh directory in Arbor’s source code.

Starting out: Loading our own catalogue

Next, we have to do multiple things in parallel

  1. create our own channel

  2. build a catalogue containing said channel

  3. update the example accordingly

We start by creating a new directory mod (the name is not important, but will be used throughout this example) and adding a file named hh02.mod to it. Its contents should be this:

NEURON {
    SUFFIX hh02
    NONSPECIFIC_CURRENT il
}

BREAKPOINT {
    il = 0
}

We will discuss this in detail below, but for now, we will just translate and use it in our model. Change to a shell, next to the mod directory and type:

arbor-build-catalogue cat mod

and an output like this should appear (again cat is an arbitrary name we will use from here on).

Building catalogue 'cat' from mechanisms in /Users/hater/src/arbor/python/example/hh/mod
 * NMODL
   * hh02
Catalogue has been built and copied to /Users/hater/src/arbor/python/example/hh/cat-catalogue.so

and the file cat-catalogue.so should appear in your current directory. Next, modify the Python file like this:

#!/usr/bin/env python3

import arbor as A
from arbor import units as U
import matplotlib.pyplot as plt

# Create a single segment morphology
tree = A.segment_tree()
tree.append(A.mnpos, (-3, 0, 0, 3), (3, 0, 0, 3), tag=1)

# Create (almost empty) decor
decor = (
    A.decor()
    .paint("(all)", A.density("hh03"))
    .place("(location 0 0.5)", A.iclamp(10 * U.ms, 2 * U.ms, 0.8 * U.nA), "iclamp")
)

# Run the model, extracting the membrane voltage
model = A.single_cell_model(A.cable_cell(tree, decor))
model.probe("voltage", "(location 0 0.5)", tag="Um", frequency=10 * U.kHz)

# add our catalogue
model.properties.catalogue = A.load_catalogue("cat-catalogue.so")

model.run(tfinal=30 * U.ms)

# Create a basic plot
fg, ax = plt.subplots()
ax.plot(model.traces[0].time, model.traces[0].value)
ax.set_xlabel("t/ms")
ax.set_ylabel("U/mV")
plt.savefig("hh-03.pdf")
plt.savefig("hh-03.svg")

This should — once run — produce a plot like this:

../_images/hh-02.svg

You can find all the steps in the python/example/hh directory in Arbor’s source code. Let’s return to what just happened, it’s quite a bit. First, we added our ion channel and used arbor-build-catalogue to translate it into a form Arbor can utilize. These collections of ion channels are — unsurprisingly — called catalogues, see Cable cell mechanisms. We pulled this into our model by loading and assigning it to the model.

Next, let’s look at the output graph. We observe a sudden jump in potential during the period the current clamp is active. As Arbor’s model for a single CV cable cell is \(\partial_t U_m = i_e - i_m\) (for multi-CV cells, we have additional terms that can be neglected here, see How the Cable Cell is made), this behaviour is expected. The current clamp provides a positive \(i_e\) and our ion channel model supplies the trans-membrane current \(i_m = 0\). To understand the latter part, consider the channel model file we just added

NEURON {
    SUFFIX hh02
    NONSPECIFIC_CURRENT il
}

This is the NEURON block declaring the channel’s name, here hh02, which is used when channels from a catalogue. Files that put SUFFIX in front of the name are converted to density channels, as opposed to synapses (POINT_PROCESS) and gap junctions (JUNCTION_PROCESS). In addition to naming the channel, we also need to set up all variables used to interface with Arbor, namely ion currents, ion concentrations, ion reversal potentials, and non-ion currents. While the ion variables follow a rigid naming scheme, which we will discuss later, non-ion currents can be freely named after NONSPECIFIC_CURRENT. We chose il here, alluding to ‘leak current’. Semantically, these currents are considered to be unassociated to any specific ion and thus can represent all ion currents we do not model explicitly as a lump sum. When computing i_m for the cable equation above, Arbor takes the sum over all non-specific and ion currents across all ion channels on the current CV. We will revisit the NEURON multiple times later on, but for now we turn to:

BREAKPOINT {
    il = 0
}

During the integration of the cable equation, Arbor will evaluate this block to update its internal picture of the currents, i.e. to calculate i_m. This occurs at an unspecified moment of the execution and might even be done multiple times, so we need to take care not to depend on the execution order. We are _expected_, yet not forced by the tooling, to update all such outputs, so, again, some care is needed.

Stepping Stone: Leak

As you might have anticipated, our next step is to produce a finite current to counteract any disturbance in the membrane potential. So, we start by adding a new mechanism to mod, called hh03, which is just a copy of hh02.mod. Next, adjust SUFFIX hh02 to SUFFIX hh03. Similarly copy step-02.py to step-03.py and change

# Create (almost empty) decor
decor = (
    A.decor()
    .paint("(all)", A.density("hh03"))
    .place("(location 0 0.5)", A.iclamp(10 * U.ms, 2 * U.ms, 0.8 * U.nA), "iclamp")
)

as well as plt.savefig('hh-03.pdf'). From on out, we’ll assume the following steps are completed at the beginning of each new section:

  1. Copy step-n.py to step-(n+1).py

    • update all references to hhn to hh(n+1)

    • update the output image to hh-(n+1).pdf

  2. Copy mod/hhn.mod to mod/hh(n+1).mod

    • change the name to SUFFIX hh(n+1)

  3. Start editing the new NMODL and Python files.

    • After each change to the NMODL file, you’ll need to call arbor-build-catalogue cat mod

Keep this in mind while we start altering the NMODL file to produce a more sensible current. Let’s start with the current itself:


BREAKPOINT {
    il = gl*(v - el)

This will pull the membrane potential v towards a resting potential el since our reduced cable equation is now \(\partial_t U_m = i_e - g_l\cdot(U_m - E_l)\). The membrane potential is available in NMODL as a read-only built-in symbol v and can be used in any ion channel. However, we need a way to set the resting potential el and the conductivity gl. This is accomplished by adding a new block to the NMODL file:


PARAMETER {
    gl =   0.0003 (S/cm2)
    el = -54.3    (mV)
}

These parameters have an optional default value and a likewise optional unit. Both are helpful to have, though. The units chosen internally by Arbor come together such that the conductivity must have units S/cm2. Note that there is neither a check nor a conversion of units; the annotation serves purely as a reminder. Now, running the example step-03.py gives us the expected result of the membrane potential returning to the resting value:

../_images/hh-03.svg

We have now recreated the leak current from the HH neuron model, which is one of three currents needed. Before we turn to the other two, though, we’ll apply some polish. Variables declared in PARAMETER blocks can be set in the call to paint, like so:

decor = (
    A.decor()
    .paint('(all)', A.density('hh03', g=0.0005, el=-70))
    .place('(location 0 0.5)', A.iclamp(10 * U.ms, 2 * U.ms, 0.8 * U.nA), "iclamp")
)

To enable this, we need to tell NMODL that each CV will have its own value of gl and el, via

NEURON {
    SUFFIX hh03
    NONSPECIFIC_CURRENT il
    RANGE el, gl
}

Without this addition, there would be one global copy for each, which could be set by writing

decor = (
    A.decor()
    .paint('(all)', A.density('hh03/el=-70,gl=0.0005'))
    .place('(location 0 0.5)', A.iclamp(10 * U.ms, 2 * U.ms, 0.8 * U.nA), "iclamp")
)

instead. Parameters are either GLOBAL or RANGE, never both. The difference is subtle and non-existent for our single CV. The rule of thumb is that if you expect a parameter to vary smoothly across the neuron, make it RANGE and if you expect discrete, clearly delineated regions with discontinuous values, go for GLOBAL. Comparing the performance impact of parameters, GLOBAL is more efficient than RANGE. While RANGE parameters consume one memory location (and access) per CV, those tagged GLOBAL require only one location (and access) regardless of CV count. Thus, all things being equal, prefer GLOBAL.

Differential Equations in NMODL

After observing the ceremony of making copies of both Python and NMODL once more, we turn to the final task. There are currents left to handle in the HH model for potassium and sodium ions. Their formulations are quite similar, so we will discuss the potassium current here and leave the sodium current as an exercise

\[\begin{split}i_{na} = \bar g_{na} m^3 h (v - E_{na})\\ i_{k} = \bar g_{k} n^4(v - E_{k})\end{split}\]

In these equations, three new variables appear: \(m, h, n\), which are defined via differential equations

\[\begin{split}n' = \alpha_{n}(v) (1 - n) - \beta_{n}(v)n\\ m' = \alpha_{m}(v) (1 - m) - \beta_{m}(v)m\\ h' = \alpha_{h}(v) (1 - h) - \beta_{h}(v)h\end{split}\]

The coefficients \(\alpha_{m,h,n}\) and \(\alpha_{m,h,n}\) are in turn

\[\begin{split}\alpha_{x}(v) = \frac{x_{\infty}(v)}{\tau_{x}}\\ \beta_{x}(v) = \frac{1 - x_{\infty}(v)}{\tau_{x}}\end{split}\]

where the steady state activations \(m,h,n_\infty\) can be determined by fitting. We will simply use them in NMODL without further justification. Add this to (the end of) hh04.mod:

FUNCTION m_alpha(v) { m_alpha = exprelr(-0.1*v - 4.0) }
FUNCTION h_alpha(v) { h_alpha = 0.07*exp(-0.05*v - 3.25) }
FUNCTION n_alpha(v) { n_alpha = 0.1*exprelr(-0.1*v - 5.5) }

FUNCTION m_beta(v)  { m_beta  = 4.0*exp(-(v + 65.0)/18.0) }
FUNCTION h_beta(v)  { h_beta  = 1.0/(exp(-0.1*v - 3.5) + 1.0) }
FUNCTION n_beta(v)  { n_beta  = 0.125*exp(-0.0125*v - 0.8125) }

The FUNCTION construct introduces a function which can only access its parameters and can have no side-effects like writing to global variables. Its return value is set by formally assigning a value to the function’s name. Arbor provides some builtin functions like exprelr, which is used here, and returns

\[\mathrm{exprelr}(x) = \frac{x}{\exp(x) - 1}\]

and its smooth continuation over \(x=0\). exprelr is useful in many models similar to HH.

Now we turn to formulating the differential equations (ODE). Naively, we could be tempted to add something like this to the BREAKPOINT block (notice the LOCAL keyword to declare block-local variables)

BREAKPOINT {
    LOCAL alpha, beta

    alpha = n_alpha(v)
    beta  = n_beta(v)
    n     = n + dt*(alpha - n*(alpha + beta))

    ik = gkbar*n*n*n*n*(v - ek)
    il = gl*(v - el)
}

and attempt to solve the ODE manually via Euler’s method. Alas, this is inconvenient and cumbersome as we needed to adapt the ion channel’s paint call every time we change the time step dt (assuming we pass it as a parameter). It’s also less accurate than desirable. What’s more, as we don’t know the order and count of evaluation of these blocks, it’s also likely to be incorrect.

There is, though, a better way by using a new variable kind, the STATE variable. Add this in your NMODL file

STATE { n }

DERIVATIVE states {
    LOCAL alpha, beta

    : potassium activation system
    alpha = n_alpha(v)
    beta  = n_beta(v)
    n'    = (alpha - n*(alpha + beta))
}

Now we have told NMODL how to compute the derivative of n. For an initial value problem, we also need to add its initial value and actually solve the ODE. This is achieved by

INITIAL {
    LOCAL alpha, beta

    : potassium activation system
    alpha = n_alpha(v)
    beta  = n_beta(v)
    n     = alpha/(alpha + beta)
}

BREAKPOINT {
    SOLVE states METHOD cnexp
    LOCAL gk, n2

    gk = gkbar*n*n*n*n

    ik  = gk*(v - ek)
    il  = gl*(v - el)
}

For closing the loop, we need to adjust the NEURON block once more:

NEURON {
    SUFFIX hh04
    : a variable can be READ *or* WRITE, the latter granting read and write access
    USEION k READ ek WRITE ik
    NONSPECIFIC_CURRENT il
    : Note, no RANGE for STATE (these are implicitly unique to each CV)
    : gkbar is a new PARAMETER
    RANGE gl, el, gkbar
}

This is how we add ionic currents in NMODL. There’s a list of predefined ions (sodium na, potassium k, and calcium ca) and new ones can be added via set_ion. The variable names in use here follow a strict naming scheme: For a hypothetical ion X, the following variables are defined:

  • iX the current

  • eX the reversal potential

  • Xi the internal concentration

  • Xo the external concentration

Their default values can be set in the simulation using set_ion. We can run the simulation again to obtain:

../_images/hh-04.svg

Final Polish: Temperature

A common problem is that measurements are often taken at temperatures different from the natural environment of a neuron. This is fixed by adjusting for the difference. For us, it presents an opportunity to introduce the ASSIGNED construct. Like STATE variables, ASSIGNED variables persist across blocks. Unlike STATE variables, we cannot take their derivative. We reproduce hh05.mod here as a reference.

NEURON {
    SUFFIX hh05
    USEION k READ ek WRITE ik
    NONSPECIFIC_CURRENT il
    : qt is ASSIGNED, but needs to be RANGE
    RANGE gl, el, gkbar, q10
}

ASSIGNED { q10 }

STATE { n }

PARAMETER {
    gkbar  =   0.036  (S/cm2)
    gl     =   0.0003 (S/cm2)
    el     = -54.3    (mV)
    celsius           (degC)
    v                 (mV)
}


INITIAL {
    LOCAL alpha, beta

    q10 = 3^(0.1*celsius - 0.63)

    : potassium activation system
    alpha = n_alpha(v)
    beta  = n_beta(v)
    n     = alpha/(alpha + beta)
}

DERIVATIVE states {
    LOCAL alpha, beta

    : potassium activation system
    alpha = n_alpha(v)
    beta  = n_beta(v)
    n'    = (alpha - n*(alpha + beta))*q10
}

BREAKPOINT {
    SOLVE states METHOD cnexp
    LOCAL gk, gna, n2

    gk = gkbar*n*n*n*n

    ik  = gk*(v - ek)
    il  = gl*(v - el)
}

FUNCTION m_alpha(v) { m_alpha = exprelr(-0.1*v - 4.0) }
FUNCTION h_alpha(v) { h_alpha = 0.07*exp(-0.05*v - 3.25) }
FUNCTION n_alpha(v) { n_alpha = 0.1*exprelr(-0.1*v - 5.5) }

FUNCTION m_beta(v)  { m_beta  = 4.0*exp(-(v + 65.0)/18.0) }
FUNCTION h_beta(v)  { h_beta  = 1.0/(exp(-0.1*v - 3.5) + 1.0) }
FUNCTION n_beta(v)  { n_beta  = 0.125*exp(-0.0125*v - 0.8125) }

Things to take note of here is celsius, which contains the temperature in degrees Celsius. While it is listed as PARAMETER here, it is not a real parameter but rather a built-in variable. Adding it to PARAMETER makes it available in the NMODL file. Adding the potential here, too, is not required but is considered good form. Using ASSIGNED over LOCAL here is again a performance consideration. While costing memory capacity and operations, the exponentiation is expensive enough to warrant the expense. This simulation results in:

../_images/hh-05.svg

Your own Journey

Now, you are free to either explore implementing ina on your own or to look up how we did it in hh06.mod and step-06.py.

If you want to try it but need a hint, here’s a rough outline:

  1. Add a declaration of the sodium ion and its variables to the NEURON block

  2. Add a parameter for the conductivity

  3. Compute the sodium current in BREAKPOINT

  4. Add two new STATE variables h and m

  5. Formulate their ODEs in INITIAL and DERIVATIVE, one block can compute all derivatives simultaneously.

Conclusion

You have probably picked up some of the quirks in syntax and semantics of NMODL. Let us be blunt: NMODL isn’t anyone’s idea of a favourite language. So, why should you be investing time in learning it? The answer is simply that it serves as the basis for all complex networks as simulated by Arbor or Neuron. Even modern attempts at addressing the same problems, like NeuroML2, are translated into NMODL. Almost all existing models are rooted in NMODL.

We did introduce the most relevant constructs here, with some notable exceptions like KINETIC and NET_RECEIVE. We will possibly return to them in the future. Note that NMODL is significantly more complex than we have shown here, but these constructs can be avoided almost entirely and/or do not apply to Arbor.

Synapses in NMODL

If you have followed the tutorial above, there isn’t much more to know before you can write your own synapse models. Thus, we will just show and discuss the exponential synapse coming with Arbor.

NEURON {
    POINT_PROCESS expsyn
    RANGE tau, e
    NONSPECIFIC_CURRENT i
}

UNITS {
    (mV) = (millivolt)
    (uS) = (microsiemens)
}

PARAMETER {
    tau = 2.0 (ms) : the default for Neuron is 0.1
    e = 0   (mV)
}

ASSIGNED {}

STATE {
    g (uS)
}

INITIAL {
    g = 0
}

BREAKPOINT {
    SOLVE state METHOD cnexp
    i = g*(v - e)
}

DERIVATIVE state {
    g' = -g/tau
}

NET_RECEIVE(weight) {
    g = g + weight
}

The differences in density mechanisms like hh are fairly minimal on the surface; instead of SUFFIX, we write POINT_PROCESS in front of the name. The rest can be summarised as follows: We have a non-specific current i driving the model towards the resting potential e, which is proportional to the conductivity g. The conductivity g itself decays exponentially towards zero with time constant tau. One subtle difference is the unit for currents; where density mechanisms output a current density of mA/cm2 (thus g in hh above has unit S/cm2) point mechanisms directly produce a current in units of nA, resulting in units uS for g.

One new block NET_RECEIVE, which is exclusive to POINT_PROCESS, makes its appearance. It is evaluated every time an event is dispatched to the synapse, which might be never, once, or even multiple times during a time step, depending on the surrounding network state. This evaluation is done in unspecified order w.r.t. to BREAKPOINT, i.e. before, after, or both. However, as this is specifically to handle events, Arbor will evaluate this block exactly once per event, so the advice about BREAKPOINT and multiple evaluation does not apply to NET_RECEIVE.

The formal parameter weight contains the connection weight according to the connection along which the action potential being processed was transmitted. To put this into the larger context, we sketch the infrastructure used to tie synapses and connections into a model:

class recipe(A.recipe):
    def __init__(self, ncells):
        A.recipe.__init__(self)

    def num_cells(self):
        return 2

    def cell_description(self, gid):
        tree = A.segment_tree()
        tree.append(A.mnpos, (-3, 0, 0, 3), (3, 0, 0, 3), tag=1)
        # syn and det form the targets and sources for connections.
        decor = (
            A.decor()
            .paint('(all)', A.density("hh"))
            .place('(location 0 0.5)', A.synapse("expsyn"), "syn")
            .place('(location 0 0.5)', A.threshold_detector(-10 * U.mV), "det")
        )
        return A.cable_cell(tree, decor)

    def cell_kind(self, _):
        return A.cell_kind.cable

    def connections_on(self, gid):
        # gid 0 is purely a source
        if gid == 0:
            return []
        # gid 1 has an incoming connection from gid 0
        # any spike dispatched to syn will trigger NET_RECEIVE on expsyn with weight=0.01
        return [A.connection((0, "det"), "syn", 0.01, 5 * U.ms)]

    def global_properties(self, _):
        return A.neuron_cable_properties()

This concludes our look at synapses or POINT_PROCESS es. NMODL has more to offer, like the ability to formulate KINETIC reaction systems, and Arbor adds some extensions, e.g. voltage processes, longitudinal diffusion, support for STDP, and mechanisms for gap junctions.