Extracellular signals (LFPykit)

This example takes elements from other tutorials to create a geometrically detailed single cell model from an SWC morphology file, and adds predictions of extracellular potentials using the LFPykit Python library. LFPykit provides a few different classes facilitating calculations of extracellular potentials and related electroencephalography (EEG) and magnetoencephalography (MEG) signals from geometrically detailed neuron models under various assumptions. These are signals that mainly stem from transmembrane currents.

Note

Concepts covered in this tutorial:

1. Recording of transmembrane currents using arbor.cable_probe_total_current_cell

2. Recording of stimulus currents using arbor.cable_probe_stimulus_current_cell

3. Using the arbor.place_pwlin API

4. Map recorded transmembrane currents to extracellular potentials by derivinging Arbor specific classes from LFPykit’s lfpykit.LineSourcePotential and lfpykit.CellGeometry

The line source approximation

First, let’s describe how one can compute extracellular potentials from transmembrane currents of a number of segments, assuming that each segment can be treated as an equivalent line current source using a formalism invented by Gary R. Holt and Christof Koch 1. The implementation used in this tutorial rely on lfpykit.LineSourcePotential. This class conveniently defines a 2D linear response matrix \(\mathbf{M}\) between transmembrane current array \(\mathbf{I}\) (nA) of a neuron model’s geometry (lfpykit.CellGeometry) and the corresponding extracellular electric potential in different extracellular locations \(\mathbf{V}_{e}\) (mV) so

\[\mathbf{V}_{e} = \mathbf{M} \mathbf{I}\]

The transmembrane current \(\mathbf{I}\) is an array of shape (# segments, # timesteps) with unit (nA), and each row indexed by \(j\) of \(\mathbf{V}_{e}\) represents the electric potential at each measurement site for every time step. The resulting shape of \(\mathbf{V}_{e}\) is (# locations, # timesteps). The elements of \(\mathbf{M}\) are computed as:

\[M_{ji} = \frac{1}{ 4 \pi \sigma L_i } \log \left| \frac{\sqrt{h_{ji}^2+r_{ji}^2}-h_{ji} }{ \sqrt{l_{ji}^2+r_{ji}^2}-l_{ji}} \right|~.\]

Here, segments are indexed by \(i\), segment length is denoted \(L_i\), perpendicular distance from the electrode point contact to the axis of the line segment is denoted \(r_{ji}\), longitudinal distance measured from the start of the segment is denoted \(h_{ji}\) and longitudinal distance from the other end of the segment is denoted \(l_{ji}= L_i + h_{ji}\).

Note

Assumptions:

1. The extracellular conductivity \(\sigma\) is infinite, homogeneous, frequency independent (linear) and isotropic

2. Each segment is treated as a straight line source with homogeneous current density between its start and end point coordinate. Although Arbor allows segments to be defined as conical frusta with varying radius, we shall assume that any variation in radius is small relative to overall segment length.

3. Each extracellular measurement site is treated as a point

4. The minimum distance to a line source is set equal to the average segment radius to avoid singularities.

The model

In this tutorial, the neuron model itself is kept deliberately simple with only passive (leaky) membrane dynamics, and it receives sinusoid synaptic current input in one arbitrary chosen control volume (CV).

First we import some required modules:

import sys
import numpy as np
import arbor
import lfpykit

We define a basic Recipe class, holding a cell and three probes (voltage, stimulus current and total current):

class Recipe(arbor.recipe):
    def __init__(self, cell):
        super().__init__()

        self.the_cell = cell

        self.vprobe_id = (0, 0)
        self.iprobe_id = (0, 1)
        self.cprobe_id = (0, 2)

        self.the_props = arbor.neuron_cable_properties()

    def num_cells(self):
        return 1

    def num_sources(self, gid):
        return 0

    def cell_kind(self, gid):
        return arbor.cell_kind.cable

    def cell_description(self, gid):
        return self.the_cell

    def global_properties(self, kind):
        return self.the_props

    def probes(self, gid):
        return [
            arbor.cable_probe_membrane_voltage_cell(),
            arbor.cable_probe_total_current_cell(),
            arbor.cable_probe_stimulus_current_cell(),
        ]

Then, load a morphology on SWC file format (interpreted according to Arbor’s specifications). Similar to the tutorial “A simple single cell model”, we use the morphology file in python/example/single_cell_detailed.swc:

A graphic depiction of single_cell_detailed.swc.

Pass the filename as an argument to the simulation script:

# Read the SWC filename from input
# Example from docs: single_cell_detailed.swc
if len(sys.argv) < 2:
    print("No SWC file passed to the program")
    sys.exit(0)

filename = sys.argv[1]

# define morphology (needed for ``arbor.place_pwlin`` and ``arbor.cable_cell`` below)
morphology = arbor.load_swc_arbor(filename)

As a target for a current stimuli we define an arbor.location :

# define a location on morphology for current clamp
clamp_location = arbor.location(4, 1 / 6)

Next, we let a basic function define our cell model, taking the morphology and clamp location as input. The function defines various attributes (label_dict, decor) for the cell model, sets sinusoid current clamp as stimuli using iclamp defines discretization policy (cv_policy_fixed_per_branch) and returns the corresponding place_pwlin and cable_cell objects for use later:

def make_cable_cell(morphology, clamp_location):
    # number of CVs per branch
    cvs_per_branch = 3

    # Label dictionary
    labels = arbor.label_dict()

    # decor
    decor = arbor.decor()

    # set initial voltage, temperature, axial resistivity, membrane capacitance
    decor.set_property(
        Vm=-65,  # Initial membrane voltage (mV)
        tempK=300,  # Temperature (Kelvin)
        rL=10000,  # Axial resistivity (Ω cm)
        cm=0.01,  # Membrane capacitance (F/m**2)
    )

    # set passive mechanism all over
    # passive mech w. leak reversal potential (mV)
    pas = arbor.mechanism("pas/e=-65")
    pas.set("g", 0.0001)  # leak conductivity (S/cm2)
    decor.paint("(all)", arbor.density(pas))

    # set number of CVs per branch
    policy = arbor.cv_policy_fixed_per_branch(cvs_per_branch)
    decor.discretization(policy)

    # place sinusoid input current
    iclamp = arbor.iclamp(
        5,  # stimulation onset (ms)
        1e8,  # stimulation duration (ms)
        -0.001,  # stimulation amplitude (nA)
        frequency=0.1,  # stimulation frequency (kHz)
        phase=0,
    )  # stimulation phase)
    decor.place(str(clamp_location), iclamp, '"iclamp"')

    # create ``arbor.place_pwlin`` object
    p = arbor.place_pwlin(morphology)

    # create cell and set properties
    cell = arbor.cable_cell(morphology, labels, decor)

    return p, cell

Store the function output:

# get place_pwlin and cable_cell objects
p, cell = make_cable_cell(morphology, clamp_location)

Next, we instantiate Recipe, and execute the model for a few hundred ms, sampling the different signals every 1 ms:

# instantiate recipe with cell
recipe = Recipe(cell)

# instantiate simulation
sim = arbor.simulation(recipe)

# set up sampling on probes with sampling every 1 ms
schedule = arbor.regular_schedule(1.0)
v_handle = sim.sample(recipe.vprobe_id, schedule, arbor.sampling_policy.exact)
i_handle = sim.sample(recipe.iprobe_id, schedule, arbor.sampling_policy.exact)
c_handle = sim.sample(recipe.cprobe_id, schedule, arbor.sampling_policy.exact)

# run simulation for 500 ms of simulated activity and collect results.
sim.run(tfinal=500)

Extract recorded membrane voltages, electrode and transmembrane currents. Note that membrane voltages at branch points and intersections between CVs are dropped as we only illustrate membrane voltages of segments with finite lengths.

# extract time, V_m, I_m and I_c for each CV
V_m_samples, V_m_meta = sim.samples(v_handle)[0]
I_m_samples, I_m_meta = sim.samples(i_handle)[0]
I_c_samples, I_c_meta = sim.samples(c_handle)[0]

# drop recorded V_m values and corresponding metadata of
# zero-sized CVs (branch-point potentials).
# Here this is done as the V_m values are only used for visualization purposes
inds = np.array([m.dist != m.prox for m in V_m_meta])
V_m_samples = V_m_samples[:, np.r_[True, inds]]
V_m_meta = np.array(V_m_meta)[inds].tolist()

# assert that the remaining cables comprising the metadata for each probe
# are identical, as well as the reported sample times.
assert V_m_meta == I_m_meta
assert (V_m_samples[:, 0] == I_m_samples[:, 0]).all()

Finally we sum the stimulation and transmembrane currents, allowing the stimuli to mimic a synapse current embedded in the membrane itself rather than an intracellular electrode current:

# prep recorded data for plotting and computation of extracellular potentials
time = V_m_samples[:, 0]
V_m = V_m_samples[:, 1:]

# Add stimulation current to transmembrane current to mimic sinusoid synapse
# current embedded in the membrane.
I_m = I_c_samples[:, 1:] + I_m_samples[:, 1:]  # (nA)

Compute extracellular potentials

Here we utilize the LFPykit library to map transmembrane currents recorded during the simulation to extracellular potentials in vicinity to the cell. We shall account for every segment in each CV using the so-called line-source approximation described above.

First we define a couple of inherited classes to interface LFPykit (as this library is not solely written for Arbor). Starting with a class inherited from lfpykit.CellGeometry:

class ArborCellGeometry(lfpykit.CellGeometry):
    """
    Class inherited from  ``lfpykit.CellGeometry`` for easier forward-model
    predictions in Arbor that keeps track of arbor.segment information
    for each CV.

    Parameters
    ----------
    p: ``arbor.place_pwlin`` object
        3-d locations and cables in a morphology (cf. ``arbor.place_pwlin``)
    cables: ``list``
        ``list`` of corresponding ``arbor.cable`` objects where transmembrane
        currents are recorded (cf. ``arbor.cable_probe_total_current_cell``)

    See also
    --------
    lfpykit.CellGeometry
    """

    def __init__(self, p, cables):
        x, y, z, d = [np.array([], dtype=float).reshape((0, 2))] * 4
        CV_ind = np.array([], dtype=int)  # tracks which CV owns segment
        for i, m in enumerate(cables):
            segs = p.segments([m])
            for j, seg in enumerate(segs):
                x = np.row_stack([x, [seg.prox.x, seg.dist.x]])
                y = np.row_stack([y, [seg.prox.y, seg.dist.y]])
                z = np.row_stack([z, [seg.prox.z, seg.dist.z]])
                d = np.row_stack([d, [seg.prox.radius * 2, seg.dist.radius * 2]])
                CV_ind = np.r_[CV_ind, i]

        super().__init__(x=x, y=y, z=z, d=d)
        self._CV_ind = CV_ind

Then, a class inherited from lfpykit.LineSourcePotential. Other use cases may inherit from any other parent class defined in lfpykit.models in a similar manner:

class ArborLineSourcePotential(lfpykit.LineSourcePotential):
    """subclass of ``lfpykit.LineSourcePotential`` modified for
    instances of ``ArborCellGeometry``.
    Each CV may consist of several segments , and this implementation
    accounts for their contributions normalized by surface area, that is,
    we assume constant transmembrane current density per area across each CV
    and constant current source density per unit length per segment
    (inherent in the line-source approximation).

    Parameters
    ----------
    cell: object
        ``ArborCellGeometry`` instance or similar.
    x: ndarray of floats
        x-position of measurement sites (µm)
    y: ndarray of floats
        y-position of measurement sites (µm)
    z: ndarray of floats
        z-position of measurement sites (µm)
    sigma: float > 0
        scalar extracellular conductivity (S/m)

    See also
    --------
    lfpykit.LineSourcePotential
    """

    def __init__(self, **kwargs):
        super().__init__(**kwargs)
        self._get_transformation_matrix = super().get_transformation_matrix

    def get_transformation_matrix(self):
        """Get linear response matrix

        Returns
        -------
        response_matrix: ndarray
            shape (n_coords, n_CVs) ndarray
        """
        M_tmp = self._get_transformation_matrix()
        n_CVs = np.unique(self.cell._CV_ind).size
        M = np.zeros((self.x.size, n_CVs))
        for i in range(n_CVs):
            inds = self.cell._CV_ind == i
            M[:, i] = M_tmp[:, inds] @ (
                self.cell.area[inds] / self.cell.area[inds].sum()
            )

        return M

With these two classes one may then compute extracellular potentials from transmembrane currents in space with a few lines of code:

# create ``ArborCellGeometry`` instance
cell_geometry = ArborCellGeometry(p, I_m_meta)

# define locations where extracellular potential is predicted in vicinity
# of cell.
# Axis limits [x-min, x-max, y-min, y-max] (µm)
axis = np.array([-110, 370, -80, 70])
dx = 2  # spatial resolution along x-axis (µm)
dz = 2  # spatial resolution along y-axis (µm)
X, Y = np.meshgrid(
    np.linspace(axis[0], axis[1], int(np.diff(axis[:2]) // dx) + 1),
    np.linspace(axis[2], axis[3], int(np.diff(axis[2:]) // dz) + 1),
)
Z = np.zeros_like(X)

# ``ArborLineSourcePotential`` instance, get mapping for all segments per CV
lsp = ArborLineSourcePotential(
    cell=cell_geometry, x=X.flatten(), y=Y.flatten(), z=Z.flatten()
)
M = lsp.get_transformation_matrix()

# Extracellular potential in x,y-plane (mV)
V_e = M @ I_m.T

The result

The visualization below of simulation results shows the cellular geometry and a contour plot of the extracellular potential (V_e) in a plane. Each part (CV) of the cell is shown with some color coding for the membrane potential (V_m). The stimulus location is denoted by the black marker.

The result.

Note

The spatial discretization is here deliberately coarse with only 3 CVs per branch. Hence the branch receiving input about 1/6 of the way from its root (from decor.place(str(arbor.location(4, 1/6)), iclamp, '"iclamp"')) is treated as 3 separate line sources with inhomogeneous current density per unit length. This inhomogeneity is due to the fact that the total transmembrane current per CV may be distributed across multiple segments with varying surface area. The transmembrane current is assumed to be constant per unit length per segment.

The full code

You can find the full code of the example at python/examples/probe_lfpykit.py.

Find the source of LFPykit here.

References

1

Holt, G.R., Koch, C. Electrical Interactions via the Extracellular Potential Near Cell Bodies. J Comput Neurosci 6, 169–184 (1999). https://doi.org/10.1023/A:1008832702585