{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n\n# Simulate raw data using subject anatomy\n\nThis example illustrates how to generate source estimates and simulate raw data\nusing subject anatomy with the :class:`mne.simulation.SourceSimulator` class.\nOnce the raw data is simulated, generated source estimates are reconstructed\nusing dynamic statistical parametric mapping (dSPM) inverse operator.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Author: Ivana Kojcic \n# Eric Larson \n# Kostiantyn Maksymenko \n# Samuel Deslauriers-Gauthier \n\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\n\nimport mne\nfrom mne.datasets import sample\n\nprint(__doc__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this example, raw data will be simulated for the sample subject, so its\ninformation needs to be loaded. This step will download the data if it not\nalready on your machine. Subjects directory is also set so it doesn't need\nto be given to functions.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "data_path = sample.data_path()\nsubjects_dir = data_path / \"subjects\"\nsubject = \"sample\"\nmeg_path = data_path / \"MEG\" / subject" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we get an info structure from the sample subject.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fname_info = meg_path / \"sample_audvis_raw.fif\"\ninfo = mne.io.read_info(fname_info)\ntstep = 1 / info[\"sfreq\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To simulate sources, we also need a source space. It can be obtained from the\nforward solution of the sample subject.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fwd_fname = meg_path / \"sample_audvis-meg-eeg-oct-6-fwd.fif\"\nfwd = mne.read_forward_solution(fwd_fname)\nsrc = fwd[\"src\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To simulate raw data, we need to define when the activity occurs using events\nmatrix and specify the IDs of each event.\nNoise covariance matrix also needs to be defined.\nHere, both are loaded from the sample dataset, but they can also be specified\nby the user.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "fname_event = meg_path / \"sample_audvis_raw-eve.fif\"\nfname_cov = meg_path / \"sample_audvis-cov.fif\"\n\nevents = mne.read_events(fname_event)\nnoise_cov = mne.read_cov(fname_cov)\n\n# Standard sample event IDs. These values will correspond to the third column\n# in the events matrix.\nevent_id = {\n \"auditory/left\": 1,\n \"auditory/right\": 2,\n \"visual/left\": 3,\n \"visual/right\": 4,\n \"smiley\": 5,\n \"button\": 32,\n}\n\n\n# Take only a few events for speed\nevents = events[:80]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to simulate source time courses, labels of desired active regions\nneed to be specified for each of the 4 simulation conditions.\nMake a dictionary that maps conditions to activation strengths within\naparc.a2009s :footcite:`DestrieuxEtAl2010` labels.\nIn the aparc.a2009s parcellation:\n\n- 'G_temp_sup-G_T_transv' is the label for primary auditory area\n- 'S_calcarine' is the label for primary visual area\n\nIn each of the 4 conditions, only the primary area is activated. This means\nthat during the activations of auditory areas, there are no activations in\nvisual areas and vice versa.\nMoreover, for each condition, contralateral region is more active (here, 2\ntimes more) than the ipsilateral.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "activations = {\n \"auditory/left\": [\n (\"G_temp_sup-G_T_transv-lh\", 30), # label, activation (nAm)\n (\"G_temp_sup-G_T_transv-rh\", 60),\n ],\n \"auditory/right\": [\n (\"G_temp_sup-G_T_transv-lh\", 60),\n (\"G_temp_sup-G_T_transv-rh\", 30),\n ],\n \"visual/left\": [(\"S_calcarine-lh\", 30), (\"S_calcarine-rh\", 60)],\n \"visual/right\": [(\"S_calcarine-lh\", 60), (\"S_calcarine-rh\", 30)],\n}\n\nannot = \"aparc.a2009s\"\n\n# Load the 4 necessary label names.\nlabel_names = sorted(\n set(\n activation[0]\n for activation_list in activations.values()\n for activation in activation_list\n )\n)\nregion_names = list(activations.keys())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create simulated source activity\n\nGenerate source time courses for each region. In this example, we want to\nsimulate source activity for a single condition at a time. Therefore, each\nevoked response will be parametrized by latency and duration.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def data_fun(times, latency, duration):\n \"\"\"Generate source time courses for evoked responses.\"\"\"\n f = 15 # oscillating frequency, beta band [Hz]\n sigma = 0.375 * duration\n sinusoid = np.sin(2 * np.pi * f * (times - latency))\n gf = np.exp(\n -((times - latency - (sigma / 4.0) * rng.rand(1)) ** 2) / (2 * (sigma**2))\n )\n return 1e-9 * sinusoid * gf" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here, :class:`~mne.simulation.SourceSimulator` is used, which allows to\nspecify where (label), what (source_time_series), and when (events) event\ntype will occur.\n\nWe will add data for 4 areas, each of which contains 2 labels. Since add_data\nmethod accepts 1 label per call, it will be called 2 times per area.\n\nEvoked responses are generated such that the main component peaks at 100ms\nwith a duration of around 30ms, which first appears in the contralateral\ncortex. This is followed by a response in the ipsilateral cortex with a peak\nabout 15ms after. The amplitude of the activations will be 2 times higher in\nthe contralateral region, as explained before.\n\nWhen the activity occurs is defined using events. In this case, they are\ntaken from the original raw data. The first column is the sample of the\nevent, the second is not used. The third one is the event id, which is\ndifferent for each of the 4 areas.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "times = np.arange(150, dtype=np.float64) / info[\"sfreq\"]\nduration = 0.03\nrng = np.random.RandomState(7)\nsource_simulator = mne.simulation.SourceSimulator(src, tstep=tstep)\n\nfor region_id, region_name in enumerate(region_names, 1):\n events_tmp = events[np.where(events[:, 2] == region_id)[0], :]\n for i in range(2):\n label_name = activations[region_name][i][0]\n label_tmp = mne.read_labels_from_annot(\n subject, annot, subjects_dir=subjects_dir, regexp=label_name, verbose=False\n )\n label_tmp = label_tmp[0]\n amplitude_tmp = activations[region_name][i][1]\n if region_name.split(\"/\")[1][0] == label_tmp.hemi[0]:\n latency_tmp = 0.115\n else:\n latency_tmp = 0.1\n wf_tmp = data_fun(times, latency_tmp, duration)\n source_simulator.add_data(label_tmp, amplitude_tmp * wf_tmp, events_tmp)\n\n# To obtain a SourceEstimate object, we need to use `get_stc()` method of\n# SourceSimulator class.\nstc_data = source_simulator.get_stc()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simulate raw data\n\nProject the source time series to sensor space. Three types of noise will be\nadded to the simulated raw data:\n\n- multivariate Gaussian noise obtained from the noise covariance from the\n sample data\n- blink (EOG) noise\n- ECG noise\n\nThe :class:`~mne.simulation.SourceSimulator` can be given directly to the\n:func:`~mne.simulation.simulate_raw` function.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "raw_sim = mne.simulation.simulate_raw(info, source_simulator, forward=fwd)\nraw_sim.set_eeg_reference(projection=True)\n\nmne.simulation.add_noise(raw_sim, cov=noise_cov, random_state=0)\nmne.simulation.add_eog(raw_sim, random_state=0)\nmne.simulation.add_ecg(raw_sim, random_state=0)\n\n# Plot original and simulated raw data.\nraw_sim.plot(title=\"Simulated raw data\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Extract epochs and compute evoked responsses\n\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "epochs = mne.Epochs(raw_sim, events, event_id, tmin=-0.2, tmax=0.3, baseline=(None, 0))\nevoked_aud_left = epochs[\"auditory/left\"].average()\nevoked_vis_right = epochs[\"visual/right\"].average()\n\n# Visualize the evoked data\nevoked_aud_left.plot(spatial_colors=True)\nevoked_vis_right.plot(spatial_colors=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Reconstruct simulated source time courses using dSPM inverse operator\n\nHere, source time courses for auditory and visual areas are reconstructed\nseparately and their difference is shown. This was done merely for better\nvisual representation of source reconstruction.\nAs expected, when high activations appear in primary auditory areas, primary\nvisual areas will have low activations and vice versa.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "method, lambda2 = \"dSPM\", 1.0 / 9.0\ninv = mne.minimum_norm.make_inverse_operator(epochs.info, fwd, noise_cov)\nstc_aud = mne.minimum_norm.apply_inverse(evoked_aud_left, inv, lambda2, method)\nstc_vis = mne.minimum_norm.apply_inverse(evoked_vis_right, inv, lambda2, method)\nstc_diff = stc_aud - stc_vis\n\nbrain = stc_diff.plot(\n subjects_dir=subjects_dir, initial_time=0.1, hemi=\"split\", views=[\"lat\", \"med\"]\n)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n.. footbibliography::\n\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.2" } }, "nbformat": 4, "nbformat_minor": 0 }