Sliding time window EBBlayer

This tutorial demonstrates the need for a sliding window, layer-specific version of EBB reconstruction to accurately infer two distinct laminar sources across time and space. It builds on the laminar CSD tutorial and addresses limitations, as well as best practices for using EBB for laminar inference.

It is composed of two parts: one adressing the need for sliding windows and the other for EBBlayer, a version of EBB that accounts for correlated sources across layers.

Two sources are simulated, one in superficial and one in deep layers, either consecutively or simultaneously. Source reconstruction is performed on the simulated sensor data using a forward model based on the multilayer mesh, first on the whole time window using EBB, then using a sliding time window EBB, an finally using a sliding time window with a modified version of EBB with specific priors accounting for correlated sources across layers. Laminar CSD analysis is then run on the laminar source reconstructed signals at the location of the peak variance.

Part I: Temporal smoothing effect: the need for sliding time windows

Setting up the simulations

Simulations are based on an existing dataset, which is used to define the sampling rate, number of trials, duration of each trial, and the channel layout.

import os
import shutil
import numpy as np
import matplotlib.pyplot as plt
import tempfile
from IPython.display import Image
import base64

from lameg.invert import invert_ebb, invert_sliding_window_ebb, invert_sliding_window_ebb_layer, coregister, load_source_time_series
from lameg.laminar import compute_csd
from lameg.simulate import run_current_density_simulation
from lameg.surf import LayerSurfaceSet
from lameg.util import get_fiducial_coords, load_meg_sensor_data, init_spm_parallel_pool
from lameg.viz import show_surface, color_map, plot_csd
import spm_standalone

# Subject information for data to base the simulations on
subj_id = 'sub-104'
ses_id = 'ses-01'

# Fiducial coil coordinates
fid_coords = get_fiducial_coords(subj_id, '../test_data/participants.tsv')

# Data file to base simulations on
data_file = os.path.join(
    '../test_data',
    subj_id,
    'meg',
    ses_id,
    f'spm/pspm-converted_autoreject-{subj_id}-{ses_id}-001-btn_trial-epo.mat'
)

spm = spm_standalone.initialize()
init_spm_parallel_pool(spm, 16)

The simulations and source reconstructions will be based on a forward model using the multilayer mesh, with 11 layers

surf_set = LayerSurfaceSet(subj_id, 11)

verts_per_surf = surf_set.get_vertices_per_layer()

We’re going to copy the data file to a temporary directory and direct all output there.

# Extract base name and path of data file
data_path, data_file_name = os.path.split(data_file)
data_base = os.path.splitext(data_file_name)[0]

# Where to put simulated data
tmp_dir = tempfile.mkdtemp()

# Copy data files to tmp directory
if not os.path.exists(tmp_dir):
   os.mkdir(tmp_dir)

# Copy data files to tmp directory
shutil.copy(
    os.path.join(data_path, f'{data_base}.mat'),
    os.path.join(tmp_dir, f'{data_base}.mat')
)
shutil.copy(
    os.path.join(data_path, f'{data_base}.dat'),
    os.path.join(tmp_dir, f'{data_base}.dat')
)

# Construct base file name for simulations
base_fname = os.path.join(tmp_dir, f'{data_base}.mat')

Invert the subject’s data using the multilayer mesh. This step only has to be done once - this is just to compute the forward model that will be used in the simulations

# Patch size to use for inversion (in this case it matches the simulated patch size)
patch_size = 5
# Number of temporal modes to use for EBB inversion
n_temp_modes = 4

# Coregister data to multilayer mesh
coregister(
    fid_coords,
    base_fname,
    surf_set,
    spm_instance=spm
)

# Run inversion
[_,_] = invert_ebb(
    base_fname,
    surf_set,
    patch_size=patch_size,
    n_temp_modes=n_temp_modes,
    spm_instance=spm
)

Simulating two sources in different layers In other tutorials we explored the case of a single laminar source but in real data, we might expect several consecutive sources in different layers.

We’re going to simulate 200ms of 2 Gaussian with a dipole moment of 5nAm, a width of 25ms, the second one occuring 80ms after the first

_, time, _ = load_meg_sensor_data(base_fname)
print(f'Time from {time[0]} to {time[-1]}, with shape {time.shape}')

# Strength of simulated activity (nAm)
dipole_moment = 5
# Temporal width of the simulated Gaussian
signal_width= 10 # 25ms
# Sampling rate (must match the data file)
s_rate = 600

zero_time1=-40
zero_time2=40
sim_signal1=np.exp(-((time-zero_time1)**2)/(2*signal_width**2)).reshape(1,-1)
sim_signal2=np.exp(-((time-zero_time2)**2)/(2*signal_width**2)).reshape(1,-1)
plt.plot(time,dipole_moment*sim_signal1[0,:],label='sim1')
plt.plot(time,dipole_moment*sim_signal2[0,:],label='sim2')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude (nAm)')

Select the best vertex to simulate at (based on anatomical measures), and get the corresponding layer vertices. See other tutorials.

# Vertex to simulate activity at
sim_vertex=10561
sim_vertices = [
    surf_set.get_multilayer_vertex(1, sim_vertex), # simulating in superficial
    surf_set.get_multilayer_vertex(9, sim_vertex) # simulating in deep
]

We’ll simulate a 5mm patch of activity with -5 dB SNR at the sensor level. The desired level of SNR is achieved by adding white noise to the projected sensor signals

# Size of simulated patch of activity (mm)
sim_patch_size = 5
# SNR of simulated data (dB)
SNR = -5

prefix = f'sim_{sim_vertex}_'

# Generate simulated data
sim_fname = run_current_density_simulation(
    base_fname,
    prefix,
    sim_vertices,
    np.vstack([sim_signal1, sim_signal2]),
    [dipole_moment, dipole_moment],
    [sim_patch_size, sim_patch_size],
    SNR,
    average_trials=True,
    spm_instance=spm
)

Inversion

Now we’ll run a source reconstruction using the multilayer mesh, select the vertex to examine, extract the source signals at each layer in that location, and compute a laminar CSD

Inversion on the full time window

We will first look at the behavior of inverting the full window with the basic EBB algorithm

[_,_, MU] = invert_ebb(
    sim_fname,
    surf_set,
    patch_size=patch_size,
    n_temp_modes=n_temp_modes,
    return_mu_matrix= True,
    spm_instance=spm
)

pial_layer_vertices = np.arange(verts_per_surf)
pial_layer_ts, time, _ = load_source_time_series(
    sim_fname,
    vertices=pial_layer_vertices
)

# Peak
mean_pial_layer_ts = np.mean(pial_layer_ts,axis=-1)
peak = np.argmax(mean_pial_layer_ts)

print(f'Simulated vertex={sim_vertex}, Prior vertex={peak}')

We can see that the peak is at the same location we simulated at

# Plot colors and camera view
max_abs = np.max(np.abs(mean_pial_layer_ts))
c_range = [-max_abs, max_abs]

# Plot peak
colors,_ = color_map(
    mean_pial_layer_ts,
    "RdYlBu_r",
    c_range[0],
    c_range[1]
)
thresh_colors=np.ones((colors.shape[0],4))*255
thresh_colors[:,:3]=colors
thresh_colors[mean_pial_layer_ts<np.percentile(mean_pial_layer_ts,99.9),3]=0

cam_view = [40, -240, 25,
            60, 37, 17,
            0, 0, 1]

plot = show_surface(
    surf_set,
    vertex_colors=thresh_colors,
    info=True,
    camera_view=cam_view,
    marker_vertices=peak,
    marker_size=3,
    marker_color=[0,0,255]
)

We need the indices of the vertex at each layer for this location, and the distances between them

# first retreiving the distance between layers of the maximum peak
layer_verts = surf_set.get_layer_vertices(peak)
layer_dists = surf_set.get_interlayer_distance(peak)
print(layer_dists)

Now we can compute and plot the laminar CSD

# Set a function for loading the time series, computing CSD and plotting
def load_and_plot_csd(sim_fname, surf_set, layer_verts, layer_dists, MU=None):
    # Get source time series for each layer (already averaged)
    if MU is not None:
        mean_layer_ts, time, _ = load_source_time_series(sim_fname, mu_matrix=MU, vertices=layer_verts)
    else:
        mean_layer_ts, time, _ = load_source_time_series(sim_fname, vertices=layer_verts)

    csd = compute_csd(mean_layer_ts, np.sum(layer_dists), s_rate, method='KCSD1D')

    col_r = plt.cm.cool(np.linspace(0, 1, num=surf_set.n_layers))
    plt.figure(figsize=(15, 4))
    plt.subplot(1, 3, 1)
    for l in range(surf_set.n_layers):
        plt.plot(time, mean_layer_ts[l, :], label=f'{l}', color=col_r[l, :])
    plt.legend(loc='upper left')
    plt.xlabel('Time (ms)')
    plt.ylabel('Source (nAm)')

    plt.subplot(1, 3, 2)
    col_r = plt.cm.cool(np.linspace(0, 1, num=csd.shape[0]))
    for l in range(csd.shape[0]):
        plt.plot(time, csd[l, :], color=col_r[l, :])
    plt.xlabel('Time (ms)')
    plt.ylabel('CSD')

    ax = plt.subplot(1, 3, 3)
    plot_csd(csd, time, ax, n_layers=surf_set.n_layers)
    plt.xlabel('Time (ms)')
    plt.ylabel('Layer')
    plt.tight_layout()

load_and_plot_csd(sim_fname, surf_set, layer_verts, layer_dists, MU=MU)

We can observe that the two sources are not separated spatially, we do not retreive the layers we simulated at. Note that the amplitude difference in the second source comes from the source being deep and further away from the surface.

Inversion with a sliding window

We then reconstruct layer source time series using a sliding time window of 25ms.

# a modified version of invert_ebb using sliding time windows
[_,_] = invert_sliding_window_ebb(
    sim_fname,
    surf_set,
    patch_size=patch_size,
    n_temp_modes=n_temp_modes,
    win_size=25,
    win_overlap=True,
    spm_instance=spm
)

Then, we extract the time series and compute the laminar CSD

load_and_plot_csd(sim_fname, surf_set, layer_verts, layer_dists)

Here the spatiotemporal dynamics of layer activity are accurately reconstructed.

Explanation: The empirical Bayesian beamformer uses a sensor covariance matrix to compute its inversion. When inverting on the whole window, the covariance contains both deep and superficial layer activity, that have extremely similar leadfields, and it collapses into a compromise between the two. It then sees a dominant spatial pattern, for the whole time window. If we invert with a sliding time window, we get only one laminar source at a time, spatially decorrelated. But keep in mind that temporal and spatial discrimination of laminar sources strongly depends on the width of your reconstructed window. If you use too small a time window, you don’t get enough time samples to build an accurate covariance matrix.

Part II: Correlated sources across layers: the need for layer-specific spatial source priors

Simulating two simultaneous sources

Here we simulate two sources - one in a superficial layer and one in a deep layer (as before) but occuring at the same time. We keep the same SNR = -5bB and dipole moment 5nAm as previously.

zero_time_1new=0
zero_time_2new=0
sim_signal1_new=np.exp(-((time-zero_time_1new)**2)/(2*signal_width**2)).reshape(1,-1)
sim_signal2_new=np.exp(-((time-zero_time_2new)**2)/(2*signal_width**2)).reshape(1,-1)
plt.plot(time,dipole_moment*sim_signal1_new[0,:],label='sim1')
plt.plot(time,dipole_moment*sim_signal2_new[0,:],label='sim2')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude (nAm)')

prefix_2 = f'sim2_{sim_vertex}_'

# Generate simulated data
sim_fname_2 = run_current_density_simulation(
    base_fname,
    prefix_2,
    sim_vertices,
    np.vstack([sim_signal1_new, sim_signal2_new]),
    [dipole_moment, dipole_moment],
    [sim_patch_size, sim_patch_size],
    SNR,
    average_trials=True,
    spm_instance=spm
)

Inversion with basic sliding window EBB

[_,_] = invert_sliding_window_ebb(
    sim_fname_2,
    surf_set,
    patch_size=patch_size,
    n_temp_modes=n_temp_modes,
    win_size=25,
    win_overlap=True,
    spm_instance=spm
)

load_and_plot_csd(sim_fname_2, surf_set, layer_verts, layer_dists)

We can observe that with simultaneous (correlated) sources, we reconstruct a source in the middle layer. This issue also occurs with the free energy model comparison approach, with the surface in the middle of the two simulated layers being the model with the highest likelihood.

Inversion with EBBlayer

To counteract this behavior which is a common issue with beamformers, we developed a new source reconstruction algorithm based on the implementation of EBB for correlated sources O’Neill, G.C., Barry, D.N., Tierney, T.M. et al. Testing covariance models for MEG source reconstruction of hippocampal activity. Sci Rep 11, 17615 (2021).. Our mesh has a specific shape of Nlayers x Nvertices, and leadfields within cortical column are highly correlated because dipole orientations are the same across layers (link vectors). We therefore need to take this structure into account:

For this 3 spatial source priors will be used: - an independent prior: a non-zero covariance prior for vertices with the same surface index (but same “cortical column”) - a correlated sum prior: for each layer pair, we create a combined lead field (q+=la+lb) - a correlated depth contrast prior: for each pair of potential layer sources, we create the difference lead field (q-=la-lb)

# Inverting with this layer ebb is (about 2 times longer than standard EBB)
[_,_] = invert_sliding_window_ebb_layer(
    sim_fname_2,
    surf_set,
    patch_size=patch_size,
    n_temp_modes=n_temp_modes,
    win_size=25,
    win_overlap=True,
    spm_instance=spm
)

load_and_plot_csd(sim_fname_2, surf_set, layer_verts, layer_dists)

EBB layer can then reconstruct 2 simultaneous sources.

spm.terminate()

# Delete simulation files
shutil.rmtree(tmp_dir)