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 deriving 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}\).

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 numpy as np
import arbor as A
from arbor import units as U
import lfpykit

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

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

        self.the_cell = cell

        self.vprobeset_id = (0, 0)
        self.iprobeset_id = (0, 1)
        self.cprobeset_id = (0, 2)

    def num_cells(self):
        return 1

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

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

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

    def probes(self, _):
        return [
            A.cable_probe_membrane_voltage_cell("Um-all"),
            A.cable_probe_total_current_cell("Itotal-all"),
            A.cable_probe_stimulus_current_cell("Istim-all"),
        ]

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.

Read the morphology from a file or use the default

# Read the SWC filename from input
if len(sys.argv) == 1:
    print("No SWC file passed to the program, using default.")
    filename = Path(__file__).parent / "single_cell_detailed.swc"
elif len(sys.argv) == 2:
    filename = Path(sys.argv[1])
else:
    print("Usage: single_cell_detailed.py [SWC file name]")
    sys.exit(1)

morphology = A.load_swc_arbor(filename).morphology

As a target for a current stimulus we define

clamp_location = A.location(4, 1 / 6)

Next, we define our cell model, based on the morphology and clamp location. We set various attributes (label_dict, decor) for the cell model, attaches a sinusoid current using iclamp, and sets discretization policy ats cv_policy_fixed_per_branch.

iclamp = A.iclamp(
    5 * U.ms,  # stimulation onset
    1e8 * U.ms,  # stimulation duration
    -0.001 * U.nA,  # stimulation amplitude
    frequency=100 * U.Hz,  # stimulation frequency
    phase=0 * U.rad,  # stimulation phase
)

decor = (
    A.decor()
    # set initial voltage, temperature, axial resistivity, membrane capacitance
    .set_property(
        Vm=-65 * U.mV,  # Initial membrane voltage (mV)
        tempK=300 * U.Kelvin,  # Temperature (Kelvin)
        rL=10 * U.kOhm * U.cm,  # Axial resistivity (Ω cm)
        cm=0.01 * U.F / U.m2,  # Membrane capacitance (F/m**2)
    )
    # set passive mech w. leak reversal potential (mV)
    .paint("(all)", A.density("pas/e=-65", g=0.0001))
    # attach the stimulus
    .place(str(clamp_location), iclamp, "iclamp")
    # use a fixed 3 CVs per branch
    .discretization(A.cv_policy_fixed_per_branch(3))
)

Next, we instantiate Recipe, configure the sampling, and execute the model for a few hundred ms:

cell = A.cable_cell(morphology, decor)

# instantiate recipe with cell
rec = Recipe(cell)

# instantiate simulation
sim = A.simulation(rec)

# set up sampling on probes with sampling every 1 ms
schedule = A.regular_schedule(1.0 * U.ms)
v_handle = sim.sample(0, "Um-all", schedule)
i_handle = sim.sample(0, "Itotal-all", schedule)
c_handle = sim.sample(0, "Istim-all", schedule)

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.

sim.run(tfinal=500 * U.ms)

# 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()

# prep recorded data for plotting and computation of extracellular potentials

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:

V_m = V_m_samples[:, 1:]

# Add stimulation current to transmembrane current to mimic sinusoid synapse

Compute extracellular potentials

Here we utilize the LFPykit library to map trans-membrane 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 with LFPykit, as this library is not solely written for Arbor. Starting with a class inherited from lfpykit.CellGeometry:

# ## Compute extracellular potentials
# First we define a couple of classes to interface the LFPykit
# library (https://LFPykit.readthedocs.io, https://github.com/LFPy/LFPykit):
class ArborCellGeometry(lfpykit.CellGeometry):
    """
    Class inherited from  ``lfpykit.CellGeometry`` for easier forward-model
    predictions in Arbor that keeps track of A.segment information
    for each CV.

    Parameters
    ----------
    p: ``A.place_pwlin`` object
        3-d locations and cables in a morphology (cf. ``A.place_pwlin``)
    cables: ``list``
        ``list`` of corresponding ``A.cable`` objects where transmembrane
        currents are recorded (cf. ``A.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 seg in 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]

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()
            )

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(ppwl, 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()

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