# 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}$$.

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.vprobeset_id = (0, 0)
self.iprobeset_id = (0, 1)
self.cprobeset_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:

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)


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
.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)
.paint("(all)", arbor.density("pas/e=-65", {"g": 0.0001}))
)

# 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, decor, labels)

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.vprobeset_id, schedule, arbor.sampling_policy.exact)
i_handle = sim.sample(recipe.iprobeset_id, schedule, arbor.sampling_policy.exact)
c_handle = sim.sample(recipe.cprobeset_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:

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

--------
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]])
CV_ind = np.r_[CV_ind, i]

super().__init__(x=x, y=y, z=z, d=d)


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)

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


## 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.

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