A Diffusion Process with Full Feedback

Originally we intended for the axial diffusion process to be fully integrated into the cable model (linked to the internal concentration, see here). However, for various technical reasons, this has proven infeasible and was replaced by introduction of an extra variable Xd which is subject to diffusion as opposed to the internal concentration.

This tutorial will show you how to build a simulation that fully integrates diffusive particle dynamics with a cable cell neuron. While we focus on a single neuron, this will naturally work as well in a network of cells.

Theory

In order to build a complete model, we will need to review the foundational equations as implemented in Arbor. We assume a single ion (in general: particle) species \(s\) with valence \(q_s\) and no leak currents such that \(I_m = I_s\). More ion species can be added, but do not change the fundamental mechanisms. First, the cable equation

\[\frac{1}{C}\dot U_m = \sigma \Delta U_m + I_m\]

where

\[I_s = g_s (U_m - E_s)\]

with the reversal potential \(E_s\). The current \(I_s\) – or actually the conductance \(g_s\) – is usually computed by a density mechanism.

Note that we use a slight simplifications here; \(\sigma\) and the cable radius \(r\) are taken as constant along the dendrite. Second, the diffusion equation under similar assumptions

\[\dot X_s = \beta \Delta X_s\]

where \(X_s\) is the diffusive concentration of species \(s\). Finally, the Nernst equation

\[E_s = \frac{R\cdot T}{q_s\cdot F}\log\left(\frac{X_o}{X_i}\right)\]

where \(X_{i,o}\) are the internal/external concentrations of species \(s\) If your model does not use the Nernst equation to compute \(E_s\), you can ignore this part, but be aware that this approach will not have a closed feedback loop.

Incorporating Diffusion

We will leave the cable equation and the current computation unaltered. However, from the intuition that a current is caused by movement of ions across the membrane, one would expect \(I_s\) to exact a change in \(X_s\). Thus, we will add an extra term to the diffusion equation

\[\dot X_s = \beta \Delta X_s + \frac{I_s}{q_s\cdot F\cdot V}\]

where Faraday’s constant \(F\) scales from current to change in molar amount, in conjunction with the charge per ion \(q_s\), which we scale again by the volume \(V\) to arrive at a change in concentration.

We implement this as an NMODL file intended to be added to the full cell that should look like this (cf. here):

NEURON {
    SUFFIX current_to_delta_x
    USEION x READ ix WRITE xd, xi VALENCE zx
    GLOBAL F
    RANGE xi0
}

PARAMETER {
    F = 96485.3321233100184 (coulomb/mole) : Faraday's constant
    diam                                   : compartment diameter
    xi0                                    : initial concentration
}

STATE { xi }

INITIAL {
    xi = xi0
    xd = xi0
}

BREAKPOINT  {
    : area   = pi * length * radius
    : volume = pi * length * radius^2
    : thus area/volume = radius
    xd = xd + 2*ix/(q*F*diam)
    xi = xd
}

where we exploit our knowledge about the cylindrical geometry of the CVs to compute \(\frac{A}{V} = \frac{1}{r}\). This mechanism also turns the internal concentration xi into a STATE variable, which is the standard way of handling concentration mechanisms, and couples it to xd directly. Note that this requires explicit intialisation of both xi and xd.

Setting up a Simulation

Having laid the foundation, adding this to a simulation is pretty simple. Save the NMODL file, add it to your local catalogue, and compile & load that via the usual method (cf. here).

First, declare the ion – we’ll use a new species X here, but any name will do – by calling set_ion("X", valence=1) on your global properties object. This object is part of the recipe or single_cell_model, the use of both is explained in other tutorials. When using an existing ion, this upfront declaration is not needed and you’ll get the default values for this ion.

dec = (A.decor()
    # Add our new ion to the cell; the `int_con` value has no effect.
    .set_ion("X", int_con=0.0, ext_con=42.0, diff=0.005, method="nernst/X") )
    # Place a synapse that _directly_ adds to the diffusive concentration
    .place("(location 0 0.5)", A.synapse("inject/x=X", {"alpha": 200.0}), "Zap")
    # also add an exponential decay to Xd
    .paint("(all)", A.density("decay/x=X"))
    # turn iX into a change in Xd and bind Xi to Xd
    .paint("(all)", A.density("current_to_delta_x/x=X", {"xi0": 10.0}))
    # ...
)

While simple, note some subtleties around our custom concentration mechanism: - The mechanism current_to_delta_x uses xi as a STATE and is thus

solely responsible for managing its value. This makes adding an explicit initialisation via xi0 necessary. Only one mechanism with this property should exist. See above for an alternative.

  • The change in xd due to events arriving at the synapse Zap will be synchronised with xi in our custom mechanism. If no concentration mechanism is used, the synapse needs to be modified to write to xi as well.

  • By using xi=xd, the Nernst mechanism will pick up the correct value for xi. If that is not your intention, you will have to provide a modified version of nernst in which xi is replaced with xd.

Final notes

Although the related theory is somewhat intricate, as this tutorial has shown, adding an ion (or another particle) with diffusion and full transmembrane current dynamics to a simulation is quite straightforward. You might also want to consider changing the external concentration Xo according to the ion current iX. This was not shown here for two reasons. First, Arbor does not handle extra-cellular dynamics and thus has no extra-cellular diffusion. Second, the method for handling this is identical to what we have done for xi, so including it doesn’t add any insight.