This page was generated from docs/exa_lfp.ipynb. Interactive online version:

.

Example: Local field potential¶

In this example, we show how one can sample the Local Field Potential around a neuron, by extracting the necessary data from the simulation and applying the LFP calculations for the desired sample points.

If one adds a population of neurons, does the relevant changes and focuses the sample points closer, the same methods can also yield the multi-unit activity (mind the simulation timestep though).

First of all, let’s install the required software, if on certain platforms like Colab that run the bare notebooks:

import os; from pathlib import Path
if 'COLAB_BACKEND_VERSION' in os.environ:
  !TMP=$(mktemp -d); git clone https://eden-simulator.org/repo --depth 1 -b development "$TMP"; cp -r "$TMP/." .; rm -rf "$TMP"
  exec(Path('.binder/install_livenb.py').read_text())
if 'DEEPNOTE_PROJECT_ID' in os.environ: exec(Path('../.binder/install_livenb.py').read_text())

Get the model¶

Let’s see what the LFP looks like for the L5 pyramidal cell from Mainen et al. 1995:

# Download the zip file with the model
import urllib.request
urllib.request.urlretrieve('https://codeload.github.com/OpenSourceBrain/MainenEtAl_PyramidalCell/zip/refs/heads/master', 'Mainen95.zip')
# and unpack it
from zipfile import ZipFile
with ZipFile('Mainen95.zip', 'r') as zipp: zipp.extractall()
# Check which files are present now
# from os import listdir; listdir()

Open the folder (through the GUI or os.listdir) and locate the NeuroMLv2 folder.

nml_dir = 'MainenEtAl_PyramidalCell-master/neuroConstruct/generatedNeuroML2/'
# listdir(nml_dir)

Prepare and run the simulation¶

We’ll run the existing setup with a lone neuron and a DC pulse clamp on it, but we’ll make a whole new <Simulation> file to capture all points of the neuron.

First analyse the cell to extract the compartments this neuron contains:

net_filename = nml_dir+"MainenEtAl_PyramidalCell.net.nml"
celltype_name = 'MainenNeuroML'
import eden_simulator
cells_info = eden_simulator.experimental.explain_cell(net_filename)
cell_info = cells_info[celltype_name]

comp_mid_seg = cell_info['comp_midpoint_segment']
comp_mid_fra = cell_info['comp_midpoint_fractionAlong']
n_comps = len(comp_mid_seg) # Number of actual compartments in this cell, may be retrieved from most cell_info arrays
print(f'Compartments: {n_comps}')

Compartments: 739

And then generate the new simulation file with its traces:

network_name = 'network_MainenEtAl_PyramidalCell';
population_name = 'NeuroMLBased'; N_cells = 1
# Set a new routine for making the <Simulation> file with new recording options
def NmlSimParms(network_name, run_time, run_timestep=25e-6, sampling_period=1e-3,
            rec_lines='', href='file://results.gen.txt'):
    return f'''
        <Simulation id="sim1" length="{run_time}s" step="{run_timestep}s" target="{network_name}">
            <EdenOutputFile id="first" href="{href}" format="ascii_v0" sampling_interval="{sampling_period} s">
            '''+'\n\t        '.join(rec_lines)+'''
            </EdenOutputFile>
        </Simulation>
        <Target component="sim1"/>'''

# Make the traces to record each compartment, and save the file
traces = [f"{population_name}[{neu}]/{celltype_name}/{comp_mid_seg[i]}{('%.9f'%comp_mid_fra[i])[1:]}/v"
    for neu in range(N_cells) for i in range(len(comp_mid_seg))]
rec_lines = [f'<OutputColumn id="v_{i}" quantity="{x}"  output_units="mV"/>' for i,x in enumerate(traces) ]
print(f"recording {len(rec_lines)} waveforms this time !")
sim_filename = nml_dir+"/SimLFP.xml"
with open(sim_filename, "w") as f:
    f.write('<?xml version="1.0" encoding="UTF-8"?>\n<Lems>\n<include file="MainenEtAl_PyramidalCell.net.nml"/>\n')
    f.write(NmlSimParms(network_name, run_time=0.06, run_timestep=10e-6, sampling_period=0.1e-3, href='moresults.gen.txt', rec_lines=rec_lines))
    f.write('\n</Lems>\n')
# !cat $nml_dir/SimMore.xml

recording 739 waveforms this time !

And run the simulation:

%time moresults = eden_simulator.runEden(sim_filename)
timevec = moresults['t']

CPU times: user 54.4 ms, sys: 12.1 ms, total: 66.4 ms
Wall time: 838 ms

Since we’ll be doing physics with the simulation’s results so that we can calculate the LFP, let’s convert the $V_m$ data to a physical Quantity with the help of the pint Python package for unit math.

!pip install -q pint

import numpy as np
from pint.registry import Quantity as Q
neuron_membrane_voltage = np.array([moresults[x] for i,x in enumerate(traces)])
neuron_membrane_voltage = Q(neuron_membrane_voltage,'mV')
print("%d points in space, sampled over %d points in time." % neuron_membrane_voltage.shape)

739 points in space, sampled over 601 points in time.

Let’s see a raster of $V_m$ per compartment:

from matplotlib import pyplot as plt
im = plt.imshow(neuron_membrane_voltage.m_as('mV'),
    extent=[ timevec[0], timevec[-1], n_comps, 0 ],
    aspect='auto', interpolation='none', cmap ="viridis")
cbar = plt.colorbar(im); cbar.set_label('Voltage (mV)')
plt.clim((-50, +50)); plt.ylabel('Compartment #'); plt.xlabel('Time (sec)');
plt.xlim([0.03,0.035]);

And an animation of said $V_m$:

from eden_simulator.display.animation import subsample_trajectories
samples_vm, anim_axis_vm, sampled_time_vm, (sampled_voltage,) = subsample_trajectories(
    timevec, [neuron_membrane_voltage.m_as('mV').T], animation_speed=0.0048, animation_frames_per_second=24)

import k3d; from eden_simulator.display.spatial.k3d import Plot, plot_neuron
Vm_plot = plot = Plot()
plot += plot_neuron(cell_info, sampled_voltage, time_axis_sec=anim_axis_vm,
                             color_range=[-70, -20],color_map='rainbow');
plot.camera = plot.get_auto_camera(pitch=30, yaw=10)[:6]+[0,1,0] # set 'y' to up ! vecs are pos, tgt, up
k3d_label_dict = { str(real_time): f't = {sim_time*1000:.1f} ~ ms' for (real_time, sim_time) in zip(anim_axis_vm, sampled_time_vm)  }
plt_label = k3d.text2d(k3d_label_dict, (0.,0.)); plot += plt_label # add 2d elements AFTER setting auto camera
plot.show_html()

We have recorded $V_m$ all over the cell and we’ll use it (along with the cell’s properties) to obtain the resulting LFP.

Post-process simulation data to get $I_m$¶

In this tutorial we’ll be looking at the LFP, hence we need the trans-membrane current $I_m\,$ which flows to/from extracellular space. Each compartment’s electric dynamics are:

\[C\dot{V}_m = I_{axial} + I_{transport} + I_{injected}~;~I_m = C\dot{V}_m - I_{transport}\]

where $I_{transport}$ flows between the cytoplasm and extracellular space (e.g. ion channels and passive leaks), $I_{injected}$ flows between the cytoplasm and places other than extracellular space or adjacent compartments (e.g. patch clamps), and $V_m is measured from the inner to the outer side of the membrane.

To avoid supreme confusion, in the following we will consistently measure positive capacitive current $C\dot{V}_m$ as coming OUT of the compartment and into extracellular space, and all other currents coming INTO the compartment as positive. Hence capacitive current equals (is opposed to) everything else. We will also measure the total membrane current $Im$ coming OUT of the compartment; current locally flowing OUT will cause a positive LFP in the vicinity by Ohm’s law (even though voltage drop goes along the current flow, potential at infinity = ground = zero, remember that), and current locally flowing IN will make it negative.

Keep in mind the following important points:

It might look strange but $I_m$ includes the capacitive current of the membrane (dis)charging. If you’re not convinced, draw a circuit diagram of the compartments and circle the ion channels and capacitor in parallel. What current comes into the intracellular node (axial currents and patch clamps) must equal what goes out to extracellular space, following Kirchhoff’s law of currents.
Although direct injections cross the membrane (topologically), they don’t count towards the LFP because charge travels from inside a body to inside another; thus the extracellular environment is not affected. Consider: If we bring the neuron into a steady state using a current clamp, $I_{transport} + I_{\mathrm{injected}} = 0$ (modulo spatial details), $I_{m} = -I_{patch}$ and the LFP is very real: the neuron itself becomes the probe we pump current through.
Whether a stimulus from the experimental rig should count as trans-membrane current or direct injection is ultimately up to the modeller’s intention; due caution is advised! For example: is a “virtual synapse” stimulus imposed by a patch clamp, or a simulated synapse?

Having said all that, we can measure $I_m$ by either measuring $C\dot{V}_m - I_{transport}$, or $I_{\mathrm{axial}} + I_{\mathrm{injected}}$ . Which one is easier?

Well, presently, NeuroML doesn’t provide for directly recording any of these. But not all is lost: thanks to EDEN’s cell analysis feature we can deduce capacitive and axial current. And we should know clamp current, since we are the ones who put them (right?); even if it’s a dynamic clamp that we defined with LEMS we can record the resulting current. This leaves out the current from trans-membrane mechanisms (say, channels) which takes the bulk of computing. Let’s leave it then to EDEN and work with axial + injected current.

One last wrinkle in our plans is then, where in the world are injected currents located? One could fish the exact <segment> out of the NeuroML file - but remember, we discretise those into compartments so we want the actual compartment that the clamping sites fall on. Appropriate helper routines for that will be available soon, but for now let’s exploit the facts that a) most people inject in the soma, which is usually 2) a single compartment, the very first one.

Calculating axial current¶

Just multiply membrane potential by the conductance matrix joining neurons. This matrix is so generally useful (see also the imposed field example) that a helper for Python will soon provide it.

# First, get the axial resistance matrix: Ia = G*V. It follows the tree structure
comp_ga = cell_info['comp_conductance_to_parent'] # provided in nS
iii = []; jjj = []; vvv = [] # Construct the lists (i,j,v)
for i,p in enumerate(cell_info['comp_parent']):
    g = comp_ga[i]
    if p >= 0:
        iii += [ i, p, i, p] # add 4 sparse elms in one go
        jjj += [ p, i, i, p]
        vvv += [+g,+g,-g,-g] # same as:
        # axial_conductance_matrix[i,p] += g
        # axial_conductance_matrix[p,i] += g
        # axial_conductance_matrix[i,i] -= g
        # axial_conductance_matrix[p,p] -= g

import scipy
G = scipy.sparse.coo_matrix((vvv,(iii,jjj))).tocsr() # needs csr for some reason
G = Q(G,'nS')

axial_current = (G.m @ neuron_membrane_voltage.m)*G.units*neuron_membrane_voltage.units

Calculating injected current¶

As explained previously, we’ll have to figure that out from our model for the time being. Oh well…

comp_Iclamp = Q(np.zeros_like(neuron_membrane_voltage), 'pA')
# facts of life for github.com/OpenSourceBrain/MainenEtAl_PyramidalCell - change accordingly.
comp_Iclamp[0,(3e-3 < timevec)&(timevec <= 43e-3+11e-6)] = Q(160,'pA')
# NOTE: Check the inequalities.
# The effect of Im will appear on the timestep AFTER it nominally starts and ends,
#  as the capacitive and axial currents can't have reacted to it before it started!
# Might as well apply a highpass or median filter since the simulation timestep should be less than 25 μs.

And check that the total membrane current makes sense on, say, the soma.

Calculating membrane current¶

Add them up and check that the total membrane current makes sense on, say, the soma.

# Icap + Iion + Iax + Iclamp = 0 (outward); Im = Icap + Iion = -Iax -Iclamp (outward)
neuron_membrane_current = -(-comp_Iclamp + -axial_current) # going outward!

# Let's see the breakdown
comp_capa = Q(cell_info['comp_capacitance'], 'pF')
capa_current = comp_capa[:,None] * np.diff(neuron_membrane_voltage, axis=-1)/Q(np.diff(timevec),'s') # outward!
ion_current = (neuron_membrane_current[:,1:]+neuron_membrane_current[:,:-1])/2 - capa_current # outward!

comp = 0; xlim = [0.049,0.056]
fig, (ax1, ax2) = plt.subplots(2,1,figsize=(12,8), dpi=100)
ax1.plot(timevec,-comp_Iclamp[comp].m_as('uA'), '-',lw=1,label='Clamp current');
ax1.plot(timevec,-axial_current[comp].m_as('uA'), '-',lw=1,label='Axial current');
ax1.plot((timevec[1:]+timevec[:-1])/2,capa_current[comp].m_as('uA'), '-',lw=1, label='Capacitive current');
ax1.plot((timevec[1:]+timevec[:-1])/2,ion_current[comp].m_as('uA'), '--',lw=1, label='Ion+Leak (deduced)');
ax1.plot(timevec,neuron_membrane_current[comp].m_as('uA'),'k--',lw=1,label='Total Im (deduced)');
ax1.set_ylabel('uA'); ax1.legend()
# ax1.set_xlim(xlim)
ax2.plot(timevec,neuron_membrane_voltage[comp].m_as('mV')); ax2.set_ylabel('Vm (mV)'); ax2.set_xlabel('Time (sec)');
# ax2.set_xlim(xlim)

Let’s see the total current flow for each compartment… On second thought, that’s absolute numbers for unevenly sized compartments. The density over the membrane is perhaps more informative.

fig, axs = plt.subplots(1,2, figsize=(18,5))
im = axs[0].imshow(neuron_membrane_current.m_as('pA'),
    extent=[ moresults['t'][0], timevec[-1], n_comps, 0 ],
    aspect='auto', interpolation='none', cmap ="bwr");
cbar = plt.colorbar(im); cbar.set_label('pA')
im.set_clim((-50, +50)); axs[0].set_ylabel('Compartment #'); axs[0].set_xlabel('Time (sec)');
axs[0].set_title("Comp. current");

comp_area = Q(cell_info['comp_area'], 'um*um')
im = axs[1].imshow((neuron_membrane_current/comp_area[:,None]).m_as('pA/um**2'),
    extent=[ moresults['t'][0], timevec[-1], n_comps, 0 ],
    aspect='auto', interpolation='none', cmap ="bwr");
cbar = plt.colorbar(im); cbar.set_label('pA/μm²')
im.set_clim((-1, +1)); axs[1].set_ylabel('Compartment #'); axs[1].set_xlabel('Time (sec)');
axs[1].set_title("Current density");

Observe that some compartments are much more affected by axial currents (\(\approx\) membrane currents!) than others.
Generally, the less area they have, and the tighter they are coupled to adjacent compartments, the more affected they are by axial current. (Cytoplasmic resistivity Ra and membrane spec. capacitance Cm can also matter when non-uniform.)
Feel free to plot comp_ga and comp_capa (and the ratio thereof 1/RC) to see the sensitivity of each component.

# intentionally left blank

Visualisation¶

Let’s see then how the trans-membrane current evolves on the neuron. We’ll use the same Red-White-Blue colormap we’ll display the LFP with, in the following.

# resample the anim_axis because Im has 1 time point less than Vm, due to np.diff
samples_im, anim_axis_im, sampled_time_im, [sampled_current_density] = subsample_trajectories(
    timevec[:-1], [(neuron_membrane_current / comp_area[:,None]).m_as('pA/um**2').T], animation_speed=0.0048, animation_frames_per_second=24)

Im_plot = plot = Plot()
plot += plot_neuron(cell_info, sampled_current_density, time_axis_sec=anim_axis_im,
                             color_range=[-1, +1],color_map='bwr');
plot.camera = plot.get_auto_camera(pitch=30, yaw=10)[:6]+[0,1,0] # set 'y' to up ! vecs are pos, tgt, up
plot.show_html()

Post-process $I_m$ to get the LFP¶

The other half of calculating the LFP is how much each point outside the cell is being affected by the current flow on each part of the neuron (that is, mostly by the closest compartments, but also slightly from the more distant ones).

The usual assumptions apply for this example: The extracellular medium is uniform, highly conductive, isotropic and unobstructed until far away. If you need a more sophisticated model, apply the relevant adjustments to the following. (Perhaps a separate simulation process just for the LFP is called for, see the other articles TBD on “Pipelines” generally.)

Just for illustration, we’ll use one element per compartment. It would be even better if each NeuroML <segment> is applied separately, if it’s not too much for the computer.

Alternatively, the LFPykit could also be used to calculate the coupling matrix.

extracell_conductivity = Q(.3, 'S/m')

source_points = cell_info['comp_midpoint'][:]
x_space = np.linspace(-360, 000, 10); y_space = np.linspace(-180, 600, 30); z_space = [-30]
x, y, z = np.meshgrid(x_space, y_space, z_space, indexing='ij')
sampling_points = np.stack((x,y,z),axis=-1).reshape((-1,3))
# print("Sampling points:", sampling_points)

# Vectorize calculation of pairwise distance matrix https://jaykmody.com/blog/distance-matrices-with-numpy/
distmat =  source_points @ sampling_points.T
distmat = np.sum(sampling_points**2,axis=-1) - 2*distmat + np.sum(source_points**2,axis=-1)[:,None]
distmat = np.sqrt(distmat)
distmat = Q(distmat, 'um')

# LATER: Exclude current from points inside the cell, thereby avoiding singularities as a bonus.
coupling = (1 / (4*np.pi*extracell_conductivity*distmat)).to('uV/pA')
# NB: Coupling is positive *although* potential *drops* along the electric field's lines! Think of it this way:
#     Since potential is 0 at infinity, the potential any closer to the source has to be positive! Thus I/4πσr

print('Coupling is calculated for %d current sources times %d sampling points.' % coupling.shape)
def get_lfp(membrane_current): return ((coupling.m.T @ membrane_current.m)*coupling.units*membrane_current.units).to('uV')
# get_lfp(neuron_membrane_current[:,0])

Coupling is calculated for 739 current sources times 300 sampling points.

Visualisation¶

In 4-D¶

Now let’s show the membrane voltage and the observed LFP together in one animation.

# Animate the lfp, similar to label but with a vector instead of a string
lfp_uvolt = get_lfp(neuron_membrane_current).m_as('uV').T
k3d_lfp_dict = { str(real_time): x.astype(np.float32)
                for (real_time, x) in zip(anim_axis_vm, lfp_uvolt[samples_im])  }

LFP_plot = plot = Plot()

lfp_scale_range = [-10, 10] #other options to see transitions vs. absolute range: [-.2, .2], [-.5, +.5] ...
plt_points = k3d.points(positions=sampling_points.astype(np.float32),attribute=k3d_lfp_dict,
        color_range=lfp_scale_range,point_size=10,color_map=k3d.matplotlib_color_maps.Bwr); plot += plt_points

plot.camera = plot.get_auto_camera(pitch=30, yaw=10)[:6]+[0,1,0] # set 'y' to up ! vecs are pos, tgt, up
plt_mesh = plot_neuron(cell_info, sampled_voltage, time_axis_sec=anim_axis_vm,
                    color_range=[-70, -20],color_map='rainbow', compress_cells=False);
plot += plt_mesh
plt_label = k3d.text2d(k3d_label_dict, (0.,0.)); plot += plt_label # add 2d elements AFTER setting auto camera
plot.show_html()

Navigation

Related Topics

Example: Local field potential¶

Get the model¶

Prepare and run the simulation¶

Post-process simulation data to get \(I_m\)¶

Calculating axial current¶

Calculating injected current¶

Calculating membrane current¶

Visualisation¶

Post-process \(I_m\) to get the LFP¶

Visualisation¶

In 4-D¶