.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_tutorials/tutorial_09_sliding_window_ebb_layer.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_tutorials_tutorial_09_sliding_window_ebb_layer.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_tutorials_tutorial_09_sliding_window_ebb_layer.py:


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.

.. GENERATED FROM PYTHON SOURCE LINES 23-29

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.

.. GENERATED FROM PYTHON SOURCE LINES 29-65

.. code-block:: default


    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)


.. GENERATED FROM PYTHON SOURCE LINES 66-67

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

.. GENERATED FROM PYTHON SOURCE LINES 67-72

.. code-block:: default


    surf_set = LayerSurfaceSet(subj_id, 11)

    verts_per_surf = surf_set.get_vertices_per_layer()


.. GENERATED FROM PYTHON SOURCE LINES 73-74

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

.. GENERATED FROM PYTHON SOURCE LINES 74-99

.. code-block:: default


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


.. GENERATED FROM PYTHON SOURCE LINES 100-101

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

.. GENERATED FROM PYTHON SOURCE LINES 101-124

.. code-block:: default


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


.. GENERATED FROM PYTHON SOURCE LINES 125-129

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

.. GENERATED FROM PYTHON SOURCE LINES 129-149

.. code-block:: default


    _, 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)')


.. GENERATED FROM PYTHON SOURCE LINES 150-151

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

.. GENERATED FROM PYTHON SOURCE LINES 151-159

.. code-block:: default


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


.. GENERATED FROM PYTHON SOURCE LINES 160-161

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

.. GENERATED FROM PYTHON SOURCE LINES 161-182

.. code-block:: default


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


.. GENERATED FROM PYTHON SOURCE LINES 183-186

.. image:: ../_static/tutorial_09_sim_consec.png
   :width: 800
   :alt:

.. GENERATED FROM PYTHON SOURCE LINES 188-195

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

.. GENERATED FROM PYTHON SOURCE LINES 195-217

.. code-block:: default


    [_,_, 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}')


.. GENERATED FROM PYTHON SOURCE LINES 218-219

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

.. GENERATED FROM PYTHON SOURCE LINES 219-249

.. code-block:: default


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


.. GENERATED FROM PYTHON SOURCE LINES 250-253

.. image:: ../_static/tutorial_09_localizer.png
   :width: 800
   :alt:

.. GENERATED FROM PYTHON SOURCE LINES 255-256

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

.. GENERATED FROM PYTHON SOURCE LINES 256-262

.. code-block:: default


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


.. GENERATED FROM PYTHON SOURCE LINES 263-264

Now we can compute and plot the laminar CSD

.. GENERATED FROM PYTHON SOURCE LINES 264-299

.. code-block:: default


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


.. GENERATED FROM PYTHON SOURCE LINES 300-303

.. image:: ../_static/tutorial_09_consec_ebb_csd.png
   :width: 800
   :alt:

.. GENERATED FROM PYTHON SOURCE LINES 305-310

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.

.. GENERATED FROM PYTHON SOURCE LINES 310-322

.. code-block:: default


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


.. GENERATED FROM PYTHON SOURCE LINES 323-324

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

.. GENERATED FROM PYTHON SOURCE LINES 324-327

.. code-block:: default


    load_and_plot_csd(sim_fname, surf_set, layer_verts, layer_dists)


.. GENERATED FROM PYTHON SOURCE LINES 328-331

.. image:: ../_static/tutorial_09_consec_sliding_ebb_csd.png
   :width: 800
   :alt:

.. GENERATED FROM PYTHON SOURCE LINES 333-343

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.

.. GENERATED FROM PYTHON SOURCE LINES 343-353

.. code-block:: default


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


.. GENERATED FROM PYTHON SOURCE LINES 354-357

.. image:: ../_static/tutorial_09_sim_simult.png
   :width: 800
   :alt:

.. GENERATED FROM PYTHON SOURCE LINES 357-373

.. code-block:: default


    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
    )


.. GENERATED FROM PYTHON SOURCE LINES 374-375

**Inversion with basic sliding window EBB**

.. GENERATED FROM PYTHON SOURCE LINES 375-388

.. code-block:: default


    [_,_] = 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)


.. GENERATED FROM PYTHON SOURCE LINES 389-392

.. image:: ../_static/tutorial_09_simult_sliding_ebb_csd.png
   :width: 800
   :alt:

.. GENERATED FROM PYTHON SOURCE LINES 394-404

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)

.. GENERATED FROM PYTHON SOURCE LINES 404-418

.. code-block:: default


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


.. GENERATED FROM PYTHON SOURCE LINES 419-422

.. image:: ../_static/tutorial_09_simult_sliding_ebb_layer_csd.png
   :width: 800
   :alt:

.. GENERATED FROM PYTHON SOURCE LINES 424-425

EBB layer can then reconstruct 2 simultaneous sources.

.. GENERATED FROM PYTHON SOURCE LINES 427-431

.. code-block:: default

    spm.terminate()

    # Delete simulation files
    shutil.rmtree(tmp_dir)


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.000 seconds)


.. _sphx_glr_download_auto_tutorials_tutorial_09_sliding_window_ebb_layer.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example




    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: tutorial_09_sliding_window_ebb_layer.py <tutorial_09_sliding_window_ebb_layer.py>`

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: tutorial_09_sliding_window_ebb_layer.ipynb <tutorial_09_sliding_window_ebb_layer.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_