.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_tutorials/tutorial_08_opm_data.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_tutorials_tutorial_08_opm_data.py: Laminar inference with OPM data =============================== This tutorial demonstrates how to perform laminar inference using model comparison based on free energy using simulated OPM data, as shown in [Helbling, 2025, Inferring laminar origins of MEG signals with optically pumped magnetometers (OPMs): A simulation study](https://direct.mit.edu/imag/article/doi/10.1162/imag_a_00410/125546). As in tutorial 02, a 20Hz oscillation is simulated at a particular cortical location in various layers. Source reconstruction is performed using the Empirical Bayesian Beamformer on the simulated sensor data using forward models based on different layer meshes. These models are then compared using free energy. .. GENERATED FROM PYTHON SOURCE LINES 15-19 Setting up the OPM simulations ------------------------------ Simulations are based on a synthetic OPM 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 19-48 .. code-block:: default import os import shutil import tempfile from IPython.display import Image import base64 import numpy as np import matplotlib.pyplot as plt from lameg.invert import coregister, invert_ebb from lameg.laminar import model_comparison from lameg.simulate import run_current_density_simulation, setup_opm_simulation from lameg.surf import LayerSurfaceSet from lameg.util import get_fiducial_coords from lameg.viz import color_map, show_surface import spm_standalone # Subject information for data to use subj_id = 'sub-104' ses_id = 'ses-01' # Fiducial coil coordinates fid_coords = get_fiducial_coords(subj_id, '../test_data/participants.tsv') # Where to put simulated data tmp_dir = tempfile.mkdtemp() spm = spm_standalone.initialize() .. GENERATED FROM PYTHON SOURCE LINES 49-50 The simulations will be based on a forward model using the multilayer mesh, and the model comparison will use each layer mesh .. GENERATED FROM PYTHON SOURCE LINES 50-54 .. code-block:: default surf_set_bilam = LayerSurfaceSet(subj_id, 2) surf_set = LayerSurfaceSet(subj_id, 11) .. GENERATED FROM PYTHON SOURCE LINES 55-66 Instead of using an existing MEG dataset as a template, we will now generate a new synthetic OPM dataset directly using :func:`lameg.simulate.setup_opm_simulation`. This function configures and runs SPM's OPM simulation utility (`spm_opm_sim`), which creates a valid SPM MEEG object containing empty data arrays with the desired sampling rate, number of trials, and sensor configuration. We will generate 200 trials of 1 second each (600 samples per trial), using triaxial sensors distributed over the superior half of the scalp. The simulated sensor array will be placed 3.5 mm above the scalp surface, with 35 mm spacing between adjacent sensors. The resulting `.mat` and `.dat` files will serve as the base dataset for subsequent simulations, laminar model inversion and, free energy comparison. .. GENERATED FROM PYTHON SOURCE LINES 66-82 .. code-block:: default s_rate = 600 base_fname = setup_opm_simulation( os.path.join(tmp_dir, 'test_opm'), surf_set, s_rate=s_rate, wholehead=False, sensor_spacing=35, sensor_offset=3.5, n_samples=s_rate, n_trials=200, axes=3, spm_instance=spm ) .. GENERATED FROM PYTHON SOURCE LINES 83-87 The rest of the tutorial is the same as tutorial 02 - we just now work with this base synthetic OPM dataset. 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 87-111 .. 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, n_spatial_modes=60, spm_instance=spm ) .. GENERATED FROM PYTHON SOURCE LINES 112-115 Simulating a signal on the pial surface --------------------------------------- We're going to simulate 1s of a 20Hz sine wave with a dipole moment of 10nAm .. GENERATED FROM PYTHON SOURCE LINES 115-129 .. code-block:: default # Frequency of simulated sinusoid (Hz) freq = 20 # Strength of simulated activity (nAm) dipole_moment = 10 # Generate 1s of a sine wave at a sampling rate of 600Hz (to match the data file) time = np.linspace(0,1,s_rate) sim_signal = np.sin(time*freq*2*np.pi).reshape(1,-1) plt.plot(time,dipole_moment*sim_signal[0,:]) plt.xlabel('Time (s)') plt.ylabel('Amplitude (nAm)') .. GENERATED FROM PYTHON SOURCE LINES 130-133 .. image:: ../_static/tutorial_08_sim_signal.png :width: 800 :alt: .. GENERATED FROM PYTHON SOURCE LINES 135-136 We need to pick a location (mesh vertex) to simulate at .. GENERATED FROM PYTHON SOURCE LINES 136-149 .. code-block:: default sim_vertex=10561 cam_view = [40, -240, 25, 60, 37, 17, 0, 0, 1] plot = show_surface( surf_set, marker_vertices=sim_vertex, marker_size=5, camera_view=cam_view ) .. GENERATED FROM PYTHON SOURCE LINES 150-153 .. image:: ../_static/tutorial_08_sim_location.png :width: 800 :alt: .. GENERATED FROM PYTHON SOURCE LINES 155-156 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 156-178 .. code-block:: default # Simulate at a vertex on the pial surface pial_vertex = surf_set.get_multilayer_vertex('pial', sim_vertex) prefix = f'sim_{sim_vertex}_pial_' # Size of simulated patch of activity (mm) sim_patch_size = 5 # SNR of simulated data (dB) SNR = -5 # Generate simulated data pial_sim_fname = run_current_density_simulation( base_fname, prefix, pial_vertex, sim_signal, dipole_moment, sim_patch_size, SNR, spm_instance=spm ) .. GENERATED FROM PYTHON SOURCE LINES 179-184 Model comparison (pial - white matter) -------------------------------------- Now we can run model comparison between source models based on the pial and white matter surfaces using free energy. Specifically, we'll look at the difference in free energy between the two models (pial - white matter). This should be positive (more model evidence for the pial surface model) because we simulated activity on the pial surface .. GENERATED FROM PYTHON SOURCE LINES 184-198 .. code-block:: default # Run model comparison between the first layer (pial) and the last layer (white matter) [F,_] = model_comparison( fid_coords, pial_sim_fname, surf_set_bilam, spm_instance=spm, invert_kwargs={ 'patch_size': patch_size, 'n_temp_modes': n_temp_modes, 'n_spatial_modes': 60 } ) .. GENERATED FROM PYTHON SOURCE LINES 199-201 The difference in free energy is an approximation of the Bayes factor between the two models This value should be positive (more model evidence for the pial layer model) .. GENERATED FROM PYTHON SOURCE LINES 201-203 .. code-block:: default F[0]-F[1] .. GENERATED FROM PYTHON SOURCE LINES 204-209 White matter surface simulation with pial - white matter model comparison ------------------------------------------------------------------------- Let's simulate the same pattern of activity, in the same location, but on the white matter surface. This time, model comparison should yield greater model evidence for the white matter surface, and therefore the difference in free energy (pial - white matter) should be negative. .. GENERATED FROM PYTHON SOURCE LINES 209-238 .. code-block:: default # Simulate at the corresponding vertex on the white matter surface white_vertex = surf_set.get_multilayer_vertex('white', sim_vertex) prefix = f'sim_{sim_vertex}_white_' # Generate simulated data white_sim_fname = run_current_density_simulation( base_fname, prefix, white_vertex, sim_signal, dipole_moment, sim_patch_size, SNR, spm_instance=spm ) # Run model comparison between the first layer (pial) and the last layer (white matter) [F,_] = model_comparison( fid_coords, white_sim_fname, surf_set_bilam, spm_instance=spm, invert_kwargs={ 'patch_size': patch_size, 'n_temp_modes': n_temp_modes, 'n_spatial_modes': 60 } ) .. GENERATED FROM PYTHON SOURCE LINES 239-241 The difference in free energy is an approximation of the Bayes factor between the two models This value should be negative (more model evidence for the white matter layer model) .. GENERATED FROM PYTHON SOURCE LINES 241-243 .. code-block:: default F[0]-F[1] .. GENERATED FROM PYTHON SOURCE LINES 244-249 Simulation in each layer with model comparison across layers ------------------------------------------------------------ That was model comparison with two candidate models: one based on the white matter surface, and one on the pial. Let's now simulate on each layer, and for each simulation, run model comparison across all layers. We'll turn off SPM visualization here. .. GENERATED FROM PYTHON SOURCE LINES 249-285 .. code-block:: default # Now simulate at the corresponding vertex on each layer, and for each simulation, run model comparison across # all layers all_layerF = [] for l in range(surf_set.n_layers): print(f'Simulating in layer {l}') l_vertex = surf_set.get_multilayer_vertex(l, sim_vertex) prefix = f'sim_{sim_vertex}_{l}_' l_sim_fname = run_current_density_simulation( base_fname, prefix, l_vertex, sim_signal, dipole_moment, sim_patch_size, SNR, spm_instance=spm ) [layerF,_] = model_comparison( fid_coords, l_sim_fname, surf_set, viz=False, spm_instance=spm, invert_kwargs={ 'patch_size': patch_size, 'n_temp_modes': n_temp_modes, 'n_spatial_modes': 60 } ) all_layerF.append(layerF) all_layerF = np.array(all_layerF) .. GENERATED FROM PYTHON SOURCE LINES 286-288 For each simulation, we can plot the free energy for all models relative to the worst model. The layer model with the highest free energy should correspond to the layer that the activity was simulated in. .. GENERATED FROM PYTHON SOURCE LINES 288-317 .. code-block:: default col_r = plt.cm.cool(np.linspace(0,1, num=surf_set.n_layers)) plt.figure(figsize=(10,4)) # For each simulation, plot the free energy of each layer model relative to that of the worst # model for that simulation plt.subplot(1,2,1) for l in range(surf_set.n_layers): layerF = all_layerF[l,:] plt.plot(layerF-np.min(layerF), label=f'{l}', color=col_r[l,:]) plt.legend() plt.xlabel('Eval layer') plt.ylabel(r'$\Delta$F') # For each simulation, find which layer model had the greatest free energy plt.subplot(1,2,2) peaks = [] for l in range(surf_set.n_layers): layerF = all_layerF[l,:] pk = np.argmax(layerF) peaks.append(pk) plt.plot(peaks) plt.xlim([-0.5,10.5]) plt.ylim([-0.5,10.5]) plt.plot([0,10],[0,10],'k--') plt.xlabel('Sim layer') plt.ylabel(r'Peak $\Delta$F') plt.tight_layout() .. GENERATED FROM PYTHON SOURCE LINES 318-321 .. image:: ../_static/tutorial_08_results.png :width: 800 :alt: .. GENERATED FROM PYTHON SOURCE LINES 321-342 .. code-block:: default # Normalization step norm_layerF = np.zeros(all_layerF.shape) for l in range(surf_set.n_layers): norm_layerF[l,:] = all_layerF[l,:] - np.min(all_layerF[l,:]) # Transpose for visualization im=plt.imshow(norm_layerF.T,cmap='Spectral_r') # Find the indices of the max value in each column max_indices = np.argmax(norm_layerF, axis=1) # Plot an 'X' at the center of the square for each column's maximum for idx, max_idx in enumerate(max_indices): plt.text(idx, max_idx, 'X', fontsize=12, ha='center', va='center', color='black', weight='bold') plt.xlabel('Simulated layer', fontsize=14) plt.ylabel('Evaluated layer', fontsize=14) cb=plt.colorbar(im) cb.set_label(r'$\Delta F$', fontsize=14) .. GENERATED FROM PYTHON SOURCE LINES 343-346 .. image:: ../_static/tutorial_08_results_matrix.png :width: 800 :alt: .. GENERATED FROM PYTHON SOURCE LINES 348-352 .. 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_08_opm_data.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_08_opm_data.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: tutorial_08_opm_data.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_