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

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

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

:

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.

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.