Interacting with External Fields

In this tutorial we are going to discuss how to use Arbor and an external simulation of smooth fields to facilitate interaction at the microscopic level, e.g.

• stimulation by electric fields as induced by direct currents or magnetic fields [1]

• simulations of the extra-cellular medium [2]

As a concrete example, we will model the influence of a time-varying electric field on a single cell.

This is an experimental technique leaning into the implementation details of Arbor and modcc and as such might be changed sooner or later. For a valuable outcome, readers should be familiar with modeling cable cells in Arbor, NMODL and its use in Arbor, the basics of electro-dynamics, and some low-level (C++) programming.

All source code for all intermediate steps can be in the directory python/example/reading-external-fields of the Arbor source tree.

Single Cell Model

As the first step, we create a single cell model with a sufficiently complex geometry, following the standard method outlined in for example A simple single-cell recipe

import arbor as A
from arbor import units as U

import matplotlib.pyplot as plt
from pathlib import Path

here = Path(__file__).parent

ctr = "(location 0 0.5)"

mrf = A.load_swc_neuron(here / "Acker2008.swc")

dec = (
    A.decor()
    .set_property(Vm=-55 * U.mV)
    .paint('"soma"', A.density("hh"))
    .paint('"dend"', A.density("pas/e=-70"))
    .place(ctr, A.synapse("expsyn"), "syn")
)

cvp = A.cv_policy("(every-segment)")

prp = A.neuron_cable_properties()


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

    def num_cells(self):
        return 1

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

    def global_properties(self, kind):
        return prp

    def cell_description(self, _):
        return A.cable_cell(mrf.morphology, dec, mrf.labels, cvp)

    def event_generators(self, gid):
        return [
            A.event_generator(
                target="syn", weight=50, sched=A.regular_schedule(50 * U.ms)
            )
        ]

    def probes(self, gid):
        return [A.cable_probe_membrane_voltage(ctr, tag="Um")]


if __name__ == "__main__":
    T = 1000
    dt = 0.01
    rec = recipe()
    sim = A.simulation(rec)
    hdl = sim.sample((0, "Um"), A.regular_schedule(10 * dt * U.ms))
    sim.run(T * U.ms, dt * U.ms)
    fg, ax = plt.subplots()
    for data, meta in sim.samples(hdl):
        ax.plot(data[:, 0], data[:, 1], label=meta)
    ax.legend()
    ax.set_xlabel("Time $(ms)$")
    ax.set_ylabel("Membrane potential $(mV)$")
    ax.set_xlim(0, T)
    fg.savefig("single-cell.svg", bbox_inches="tight")

For convenience and the purposes of this tutorial, we store the actual cable cell data as global variables outside the recipe

ctr = "(location 0 0.5)"

mrf = A.load_swc_neuron(here / "Acker2008.swc")

dec = (
    A.decor()
    .set_property(Vm=-55 * U.mV)
    .paint('"soma"', A.density("hh"))
    .paint('"dend"', A.density("pas/e=-70"))
    .place(ctr, A.synapse("expsyn"), "syn")
)

cvp = A.cv_policy("(every-segment)")

prp = A.neuron_cable_properties()

We add a spiking input to the soma centre

ctr = "(location 0 0.5)"

with a regular rate of \(\nu = 1/50ms\)

    def event_generators(self, gid):
        return [
            A.event_generator(
                target="syn", weight=50, sched=A.regular_schedule(50 * U.ms)
            )
        ]

For reference, this is the voltage trace after simulation:

Theory: Injecting Currents by Electric Fields

Now, we want to apply the effects of an electric field \(\vec E(\vec x, t)\) to each segment of the morphology. We follow the quasi-potential approach [3]

\[\phi = - \int \mathrm{d}\vec l\cdot \vec E\]

and add an extra term to the cable equation

\[I_\mathrm{ext} = -\sigma \frac{\partial^2}{\partial x^2} \phi\]

This term can then be absorbed into the cable equation by adding an extra mechanism per segment, feeding it the segment’s endpoints \(\vec x\) and \(\vec x'\) in 3d space and a handle to the vector field \(\vec E(\vec x, t)\). Then, we compute

\[\begin{split}I_\mathrm{ext} \simeq -\frac{1}{r_\mathrm{L}} \frac{\vec E \cdot \vec a}{||a||}\\ \vec a = \vec x - \vec x'\end{split}\]

which converges to the proper result iff the segment lengths are sufficiently small [4], where \(r_\mathrm{L}\) is the familiar axial resistivity.

Let’s step back for a moment and take a look at the morphology

which has been coloured according to the angle to a uniform vector field \(\vec E = (1, 0, 0)\)

Stepping Stone: Implicit Electric Fields

We will now show how to manipulate a standard NMODL mechanism into an injected current for a homogeneous, varying electric field \(\vec E = (E_0 \sin(\omega t), 0, 0)\). This allows us to demonstrate the technique without delving into the implementation details of transferring data from and to C++.

First, we will create a point mechanism in NMODL, then compile it to C++ using modcc. In the next section, we use the compiled C++ code and change the implementation to produce the current we want. Without further ado

NEURON {
    POINT_PROCESS efield
    RANGE xp, yp, zp, xd, yd, zd
    NONSPECIFIC_CURRENT i
}
    
PARAMETER {
    xp yp zp : proximal coordinates
    xd yd zd : distal coordinates
    e0       : amplitude
    omega    : frequency
    r_axial  : axial resistivity
}

ASSIGNED { da }

STATE { t }

INITIAL {
    LOCAL dx, dy, dz
    : vector spanned by segment endopints
    dx = xd - xp
    dy = yd - yp
    dz = zd - zp
    : unit length of the segment
    da = 1.0/(r_axial * sqrt(dx*dx + dy*dy + dz*dz))
    : time
    t = 0
}

BREAKPOINT {
    SOLVE state METHOD cnexp
    LOCAL e
    : electrical field along x
    e = e0 * sin(omega*t)
    : induced current
    i = e*(xd - xp)*da
}

DERIVATIVE state { t' = 1 }

NET_RECEIVE(weight) {}

This code formalises what we discussed above, if you are unfamiliar with writing mechanisms, please refer to How to use NMODL to extend Arbor’s repertoire of Ion Channels. We pass the desired amplitude and frequency of the E-field, the axial resistivity, the distal and proximal coordinates as parameters

PARAMETER {
    xp yp zp : proximal coordinates
    xd yd zd : distal coordinates
    e0       : amplitude
    omega    : frequency
    r_axial  : axial resistivity
}

The time is modelled via an ODE

DERIVATIVE state { t' = 1 }

which is solved to compute the induced current as outlined above

BREAKPOINT {
    SOLVE state METHOD cnexp
    LOCAL e
    : electrical field along x
    e = e0 * sin(omega*t)
    : induced current
    i = e*(xd - xp)*da
}

Compile the catalogue

$ arbor-build-catalogue efields mod
Building catalogue 'efields' from mechanisms in [...]/mod
* NMODL
  * efield
Catalogue has been built and copied to [...]/efields-catalogue.so

and it is ready to use. We create a new simulation derived from the one above

import arbor as A
from arbor import units as U

import matplotlib.pyplot as plt
from pathlib import Path
import numpy as np

here = Path(__file__).parent

ctr = "(location 0 0.5)"

mrf = A.load_swc_neuron(here / "Acker2008.swc")

dec = (
    A.decor()
    .set_property(Vm=-55 * U.mV)
    .paint('"soma"', A.density("hh"))
    .paint('"dend"', A.density("pas/e=-70"))
)

cvp = A.cv_policy("(every-segment)")

prp = A.neuron_cable_properties()
cat = A.load_catalogue(here / "efields-catalogue.so")
prp.catalogue.extend(cat, "")

e0 = 30.0e-3
omega = 0.05  # 1/ms

for id, seg in enumerate(mrf.segment_tree.segments):
    loc = f"(support (on-components 0.5 (segment {id})))"
    dec.place(
        loc,
        A.synapse(
            f"efield/e0={e0},omega={omega},r_axial={1}",
            xp=seg.prox.x,
            xd=seg.dist.x,
            yp=seg.prox.y,
            yd=seg.dist.y,
            zp=seg.prox.z,
            zd=seg.dist.z,
        ),
        label=f"ef{id}",
    )


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

    def num_cells(self):
        return 1

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

    def global_properties(self, kind):
        return prp

    def cell_description(self, _):
        return A.cable_cell(mrf.morphology, dec, mrf.labels, cvp)

    def probes(self, gid):
        return [A.cable_probe_membrane_voltage(ctr, tag="Um")]


if __name__ == "__main__":
    T = 1000
    dt = 0.01
    rec = recipe()
    sim = A.simulation(rec)
    hdl = sim.sample((0, "Um"), A.regular_schedule(10 * dt * U.ms))
    sim.run(T * U.ms, dt * U.ms)

As promised, we add one point process per segment, passing in the endpoints and field parameters, after loading our catalogue into the model

prp = A.neuron_cable_properties()
cat = A.load_catalogue(here / "efields-catalogue.so")
prp.catalogue.extend(cat, "")

e0 = 30.0e-3
omega = 0.05  # 1/ms

for id, seg in enumerate(mrf.segment_tree.segments):
    loc = f"(support (on-components 0.5 (segment {id})))"
    dec.place(
        loc,
        A.synapse(
            f"efield/e0={e0},omega={omega},r_axial={1}",
            xp=seg.prox.x,
            xd=seg.dist.x,
            yp=seg.prox.y,
            yd=seg.dist.y,
            zp=seg.prox.z,
            zd=seg.dist.z,
        ),
        label=f"ef{id}",
    )

Note, that we set the axial resistivity to the unit value here to simplify the parameter setup. In production simulations, this needs to be adjusted to the electrophysiology, possibly varying in space along the dendrite.

We plot the membrane potential and induced current over time

and clearly observe the influence of the driving current and non-linearity of the HH-mechanism.

If you can – analytically or otherwise – determine the value of \(\vec E(\vec r, t)\) at the coordinates of the segment centres without reliance on external tools, this approach is sufficient and there is no need for the techniques demonstrated in the rest of this tutorial. The same holds for other quantities interacting with the simulation.

(Ab)using NMODL to Create Custom Mechanisms

Now, this gives us a current on every segment, but for interfacing with an external simulation, we need to enable the mechanism to read from that source. For this, we create a new point mechanism with a similar layout as above

NEURON {
    POINT_PROCESS template
    RANGE xp, yp, zp, xd, yd, zd
    NONSPECIFIC_CURRENT i
}
    
PARAMETER {
    field
    xp yp zp : proximal coordinates
    xd yd zd : distal coordinates
    r_axial  : axial resistivity
}

ASSIGNED { da }

STATE {}

INITIAL {
    LOCAL dx, dy, dz
    : vector spanned by segment endopints
    dx = xd - xp
    dy = yd - yp
    dz = zd - zp
    : unit length of the segment
    da = 1.0/(r_axial * sqrt(dx*dx + dy*dy + dz*dz))
}

BREAKPOINT {
    LOCAL e
    : electrical field along x
    e = 42
    : induced current
    i = e*(xd - xp)*da
}


NET_RECEIVE(weight) {}

This will not be used in a simulation, but rather we will compile it to C++

$ arbor-build-catalogue --debug tmp efields mod
Building catalogue 'efields' from mechanisms in [...]/mod
* NMODL
  * template
  * efield
Building debug code in '[...]/tmp'.
Catalogue has been built and copied to [...]/efields-catalogue.so

By using the debug mode, the intermediate C++ files will be preserved in the tmp directory, and we can edit and add them back to the catalogue. Open tmp/build/generated/efields/template_cpu.cpp and take a look at the compute_currents function

static void compute_currents(arb_mechanism_ppack* pp) {
    PPACK_IFACE_BLOCK;
    for (arb_size_type i_ = 0; i_ < _pp_var_width; ++i_) {
        auto node_indexi_ = _pp_var_node_index[i_];
        arb_value_type current_ = 0;
        arb_value_type i = 0;
        arb_value_type e;
        e =  42.0;
        i = e*(_pp_var_xd[i_]-_pp_var_xp[i_])*_pp_var_da[i_];
        current_ = i;
        _pp_var_vec_i[node_indexi_] = fma(_pp_var_weight[i_], current_, _pp_var_vec_i[node_indexi_]);
    }
}

Note that if you have SIMD enabled, this might look differently. We recommend disabling SIMD for this tutorial, since the generated C++ is simpler to follow.

We are going to continue, as before, under the assumption that the electric field is homogeneous, so we need only to concern ourselves with three double numbers. Generalizing to complete fields is straightforward, if tedious [5].

Copy the two files tmp/build/generated/efields/template_cpu.cpp and tmp/build/generated/efields/template.hpp into the mod directory and rename them to mod/reader.hpp and mod/reader_cpu.hpp.

cp tmp/build/generated/efields/template.hpp mod/reader.hpp
cp tmp/build/generated/efields/template_cpu.cpp mod/reader_cpu.cpp

then we adjust the interface names so we can use them in our catalogue:

sed -i '' 's/template/reader/g' mod/reader_cpu.cpp
sed -i '' 's/template/reader/g' mod/reader.hpp
sed -i '' 's/make__reader_/make_arb_efields_catalogue_reader/g' mod/reader_cpu.cpp
sed -i '' 's/make__reader_/make_arb_efields_catalogue_reader/g' mod/reader.hpp

on MacOS while Linux users drop the ''. Now, you should be able to compile the catalogue again

$ arbor-build-catalogue fields mod --raw reader
Building catalogue 'efields' from mechanisms in [...]/mod
* NMODL
  * template
  * efield
* Raw
  * reader
Building debug code in '[...]/tmp'.
Catalogue has been built and copied to [...]/efields-catalogue.so

This has built an extra point process reader from a raw C++ implementation. Let’s take it for a test drive using a new simulation

import arbor as A
from arbor import units as U

import matplotlib.pyplot as plt
from pathlib import Path
import numpy as np

here = Path(__file__).parent

E = np.zeros(shape=(3,), dtype=np.float64)

ctr = "(location 0 0.5)"

mrf = A.load_swc_neuron(here / "Acker2008.swc")

dec = (
    A.decor()
    .set_property(Vm=-55 * U.mV)
    .paint('"soma"', A.density("hh"))
    .paint('"dend"', A.density("pas/e=-70"))
)

cvp = A.cv_policy("(every-segment)")

prp = A.neuron_cable_properties()
cat = A.load_catalogue(here / "efields-catalogue.so")
prp.catalogue.extend(cat, "")

e0 = 30.0e-3
omega = 0.05  # 1/ms

for id, seg in enumerate(mrf.segment_tree.segments):
    loc = f"(support (on-components 0.5 (segment {id})))"
    dec.place(
        loc,
        A.synapse(
            f"reader/field={0xC0FFEE},r_axial={1}",
            xp=seg.prox.x,
            xd=seg.dist.x,
            yp=seg.prox.y,
            yd=seg.dist.y,
            zp=seg.prox.z,
            zd=seg.dist.z,
        ),
        label=f"ef{id}",
    )


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

    def num_cells(self):
        return 1

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

    def global_properties(self, kind):
        return prp

    def cell_description(self, _):
        return A.cable_cell(mrf.morphology, dec, mrf.labels, cvp)

    def probes(self, gid):
        return [A.cable_probe_membrane_voltage(ctr, tag="Um")]


if __name__ == "__main__":
    T = 1000
    dt = 0.01
    rec = recipe()
    sim = A.simulation(rec)
    hdl = sim.sample((0, "Um"), A.regular_schedule(10 * dt * U.ms))
    t = 0
    while t < T:
        E[0] = e0 * np.sin(omega * t)
        sim.run((t + T) * U.ms, dt * U.ms)
        t += 1

We added the electric field vector

E = np.zeros(shape=(3,), dtype=np.float64)

and swapped the point process for reader passing nono-sense data to field

for id, seg in enumerate(mrf.segment_tree.segments):
    loc = f"(support (on-components 0.5 (segment {id})))"
    dec.place(
        loc,
        A.synapse(
            f"reader/field={0xC0FFEE},r_axial={1}",
            xp=seg.prox.x,
            xd=seg.dist.x,
            yp=seg.prox.y,
            yd=seg.dist.y,
            zp=seg.prox.z,
            zd=seg.dist.z,
        ),
        label=f"ef{id}",
    )

and our simulation is now stepping in increments of 1ms after which we change the x-component of the electric field

    while t < T:
        E[0] = e0 * np.sin(omega * t)
        sim.run((t + T) * U.ms, dt * U.ms)
        t += 1

Final Step: Getting Data from NumPy into C++

Now we need to wire our ‘simulation’ of the electric field into the custom mechanism. To that end, we need to acquire the pointer to the backing data of the NumPy array.

import ctypes

E = np.zeros(shape=(3,), dtype=np.float64)
buf = ctypes.c_double * 3
addr = ctypes.addressof(buf.from_buffer(E))

We use ctypes to fetch the pointer from the array. However, this is a somewhat brittle process and the machinery relies on this pointer never changing again. This means, you can write to the array like this

E[0] = 31.0

which does not change the backing memory, just the values contained within. This, however, does

x = np.array(...)
y = np.array(...)
E = x + y

and thus breaks the programm, as do all operations that assign a new array to E.

Next, we pass the pointer to the mechanism

for id, seg in enumerate(mrf.segment_tree.segments):
    loc = f"(support (on-components 0.5 (segment {id})))"
    dec.place(
        loc,
        A.synapse(
            f"reader/field={addr},r_axial={1}",
            xp=seg.prox.x,
            xd=seg.dist.x,
            yp=seg.prox.y,
            yd=seg.dist.y,
            zp=seg.prox.z,
            zd=seg.dist.z,
        ),
        label=f"ef{id}",
    )

Note that we pass an integral value as a floating point number, which is questionable, but a fix requires tweaking the NMODL language and might happen in the future. For demonstration purposes, we shift the sine wave by 50ms versus the implicit field example

    t = 0
    while t < T:
        E[0] = e0 * np.sin(omega * (t + 50))
        sim.run((t + 1) * U.ms, dt * U.ms)
        t += 1

All that’s left to do is to read the field value from the reader mechanism

static void compute_currents(arb_mechanism_ppack* pp) {
    PPACK_IFACE_BLOCK;
    auto e_field_ptr = (double*)((uint64_t) _pp_var_field);
    for (arb_size_type i_ = 0; i_ < _pp_var_width; ++i_) {
        auto node_indexi_ = _pp_var_node_index[i_];
        arb_value_type e = e_field_ptr[0];
        arb_value_type i = e*(_pp_var_xd[i_]-_pp_var_xp[i_])*_pp_var_da[i_];
        _pp_var_vec_i[node_indexi_] = fma(_pp_var_weight[i_], i, _pp_var_vec_i[node_indexi_]);
    }
}

First, we cast the floating point value to a double-valued array

    auto e_field_ptr = (double*)((uint64_t) _pp_var_field);

then we read the first value \(E_x\), use it to compute the current i as \(I = \frac{E_x a_x}{r_\mathrm{L}|a|}\)

        arb_value_type e = e_field_ptr[0];
        arb_value_type i = e*(_pp_var_xd[i_]-_pp_var_xp[i_])*_pp_var_da[i_];

which is written back the data consumed by Arbor’s cable equation solver

        _pp_var_vec_i[node_indexi_] = fma(_pp_var_weight[i_], i, _pp_var_vec_i[node_indexi_]);

Conclusions

We have seen how custom C++ mechanisms can source data from external data and they can be created, relatively simply, from NMODL templates. There are many ways to use this, examples include driving two simulations exchanging data over a shared memory segment, using a network connection to obtain values from external sources, and many more. The largest hurdle here is the synchronisation between Arbor and the external data source.

It might be prudent to note that this approach will be slow, both in terms of setup and execution. We are explicitly adding one point process to each segment and setting its coordinates. Also, technically, setting both proximal and distal coordinates is redundant as segments (mostly) share their endpoints. The way the simulation is set up, the time is stepped at the granularity of the electric field evolution. We would much rather use asynchronous simulations and a lightweight signalling mechanism.

References and Remarks

[1]

Yonezawa et al Biophysical Simulation of the Effects of Transcranial Magnetic Stimulation (TMS) on a Cortical Column, ICNIP 2025

[2]

McDougal et al Reaction-diffusion in the NEURON simulator, FNINF 2013

[3]

Wang et al, Coupling Magnetically Induced Electric Fields to Neurons: Longitudinal and Transverse Activation, Biophys 2018

[4]

Pashut et al, Mechanisms of Magnetic Stimulation of Central Nervous System Neurons, Plos CompBio 2011

[5]

A general outline might be

1. Add parameters for the extents both in terms of colocation points and physical dimensions.

2. Shift the segment endpoints into the reference frame of the field layout above.

3. Compute the indices into the array and retrieve the value at index.

4. Continue as outlined in the tutorial.