# 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:**

Recording of transmembrane currents using

`arbor.cable_probe_total_current_cell`

Recording of stimulus currents using

`arbor.cable_probe_stimulus_current_cell`

Using the

`arbor.place_pwlin`

APIMap 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

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:

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:**

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

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.

Each extracellular measurement site is treated as a point

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)
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
.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, labels, decor)
return p, cell
```

Store the function output:

```
# instantiate recipe with cell
```

Next, we instantiate `Recipe`

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

```
# 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 time, V_m, I_m and I_c for each CV
```

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.

```
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
# 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 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:

```
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:

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

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