Simulate Raw EEG using a Forward Model

Background

The source space contains thousands of possible source locations which makes selecting the right one challenging. This is where a brain atlas becomes useful. An atlas parcellates the source locations, or vertices, into a number of anatomically defined regions. In this section, you are going to use the Desikan-Killiany atlas to select anatomical regions and visualize them on the brain surface. When rendering anatomical regions, we can select different ways for displaying the surface. For example, we can view the brain as it looks from the outside or inflate the cortical surface to visualize regions that would otherwise be buried in a sulcus.

Exercises

In this section, you are going to select and visualize sources in the brain. You are going to select whole anatomical regions using the labels from the atlas parcellation as well as single vertices based on their coordinates. Here are the code snippets for solving the exercises:

The cell below loads the labels from the Desikan-Killiany parcellation atlas, which splits each hemisphere into 34 anatomical regions and prints the first 10 labels. Each label reflects an anatomical region with a varying number of vertices.

labels = mne.read_labels_from_annot('fsaverage', parc='aparc', subjects_dir=fs_dir.parent)
labels[:10]

Reading labels from parcellation...
   read 35 labels from /home/atleeri/mne_data/MNE-fsaverage-data/fsaverage/label/lh.aparc.annot
   read 34 labels from /home/atleeri/mne_data/MNE-fsaverage-data/fsaverage/label/rh.aparc.annot

[<Label | fsaverage, 'bankssts-lh', lh : 2137 vertices>,
 <Label | fsaverage, 'bankssts-rh', rh : 2196 vertices>,
 <Label | fsaverage, 'caudalanteriorcingulate-lh', lh : 1439 vertices>,
 <Label | fsaverage, 'caudalanteriorcingulate-rh', rh : 1608 vertices>,
 <Label | fsaverage, 'caudalmiddlefrontal-lh', lh : 3736 vertices>,
 <Label | fsaverage, 'caudalmiddlefrontal-rh', rh : 3494 vertices>,
 <Label | fsaverage, 'cuneus-lh', lh : 1630 vertices>,
 <Label | fsaverage, 'cuneus-rh', rh : 1638 vertices>,
 <Label | fsaverage, 'entorhinal-lh', lh : 1102 vertices>,
 <Label | fsaverage, 'entorhinal-rh', rh : 902 vertices>]

Example: Get the label for the left insula "insula-lh".

label_name = "insula-lh"
label = [l for l in labels if l.name == label_name][0]
label

<Label | fsaverage, 'insula-lh', lh : 5229 vertices>

Exercise: Get the label for the right insula "insula-rh".

label_name = "insula-rh"
label = [l for l in labels if l.name == label_name][0]
label

<Label | fsaverage, 'insula-rh', rh : 5090 vertices>

Example: Render the labeled region in the brain using a "pial" surface (shows the brain as it looks from the outside with visible gyri and sulci, good for visualizing regions that sit on the brain surface).

brain = mne.viz.Brain('fsaverage', hemi='both', subjects_dir=fs_dir.parent, surf='pial')
brain.add_label(label)

Exercise: Render the labeled region in the brain using a "white" surface (shows the white matter boundary, good for showing the exact source locations, especially for deep sources).

Solution

brain = mne.viz.Brain('fsaverage', hemi='both', subjects_dir=fs_dir.parent, surf='white')
brain.add_label(label)

Exercise: Get the label for the right transverse temporal cortex, "transversetemporal-rh" (also called Heschl’s gyrus) and render it on the brain using a "pial" surface.

Solution

label_name = "transversetemporal-rh"
label = [l for l in labels if l.name == label_name][0]
label

brain = mne.viz.Brain('fsaverage', hemi='both', subjects_dir=fs_dir.parent, surf='pial')
brain.add_label(label)

Exercise: Render the same region using an "inflated" surface (good for showing regions buried in sulci).

Solution

brain = mne.viz.Brain('fsaverage', hemi='both', subjects_dir=fs_dir.parent, surf='inflated')
brain.add_label(label)

Exercise: Choose a cortical region you are interested in from labels and render it on the brain using an appropriate surface.

Solution

label_name = "caudalmiddlefrontal-rh"
label = [l for l in labels if l.name == label_name][0]
label

brain = mne.viz.Brain('fsaverage', hemi='both', subjects_dir=fs_dir.parent, surf='pial')
brain.add_label(label)

Example: Highlight the vertex at coordinates [0, 0, 0] in the left hemisphere and set the opacity of the brain surface to alpha=0.5 to make deep sources visible.

brain = mne.viz.Brain('fsaverage', hemi='both', subjects_dir=fs_dir.parent, surf='pial', alpha=0.5)
brain.add_foci([0, 0, 0], hemi='lh', color="yellow")

Exercise: Highlight the vertex at coordinates [-44, -26, 50] in the left hemisphere (left primary motor cortex) and set the opacity of the brain surface to alpha=0.5 to make deep sources visible.

Solution

brain = mne.viz.Brain('fsaverage', hemi='both', subjects_dir=fs_dir.parent, surf='pial', alpha=0.5)
brain.add_foci([-44, -26, 50], hemi='lh', color="yellow")

Background

Now that we know how to select and visualize brain regions, we can simulate neural activity at specific source locations and project it onto the EEG sensors using the forward model. First, we need to compute the forward model from the source space and a boundary element method (BEM) solution. Then, the function simulate_sparse_stc generates source time courses at a given number of dipoles — either placed randomly or constrained to specific brain regions using the labels parameter. The result is a SourceEstimate object that stores the activation waveform at each dipole location over time. We can then use simulate_raw to project this source activity through the forward model, producing simulated EEG channel recordings.

Exercises

In this section, you are going to compute the forward model, simulate source activity with varying numbers of dipoles, visualize the source time courses on the brain and project them onto the EEG sensors. Here are the code examples required to solve the exercises:

To generate the forward solution, we need a model of how electric fields propagate from source to sensors. Such a model can be generated using the boundary element method (BEM). Instead of modeling every point in the volume, which would be computationally expensive, BEM focuses on the boundaries between tissues of different conductivity: brain, skull and skin. Here we use the BEM solution from FreeSurfer’s average template. The image generated below shows the boundaries between the different tissue layers.

bem = fs_dir / "bem" / "fsaverage-5120-5120-5120-bem-sol.fif"
mne.viz.plot_bem(subject="fsaverage", subjects_dir=fs_dir.parent, src=src);

Using surface: /home/atleeri/mne_data/MNE-fsaverage-data/fsaverage/bem/inner_skull.surf
Using surface: /home/atleeri/mne_data/MNE-fsaverage-data/fsaverage/bem/outer_skull.surf
Using surface: /home/atleeri/mne_data/MNE-fsaverage-data/fsaverage/bem/outer_skin.surf

With the source space and the BEM solution we can compute the forward model which essentially is a large matrix of shape (n_channels, n_sources*3) (since each source has a column for the x, y and z orientation) that captures the effect of every possible source on all of the sensors.

fwd = mne.make_forward_solution(
    raw.info, trans="fsaverage", src=src, bem=bem, eeg=True
)
raw.info["dev_head_t"] = fwd["info"]["dev_head_t"]

Source space          : <SourceSpaces: [<surface (lh), n_vertices=163842, n_used=10242>, <surface (rh), n_vertices=163842, n_used=10242>] MRI (surface RAS) coords, subject 'fsaverage', ~21.9 MiB>
MRI -> head transform : /home/atleeri/repositories/new-learning-platform/notebooks/mne/.pixi/envs/default/lib/python3.12/site-packages/mne/data/fsaverage/fsaverage-trans.fif
Measurement data      : instance of Info
Conductor model   : /home/atleeri/mne_data/MNE-fsaverage-data/fsaverage/bem/fsaverage-5120-5120-5120-bem-sol.fif
Accurate field computations
Do computations in head coordinates
Free source orientations

Read 2 source spaces a total of 20484 active source locations

Coordinate transformation: MRI (surface RAS) -> head
    0.999994 0.003552 0.000202      -1.76 mm
    -0.003558 0.998389 0.056626      31.09 mm
    -0.000001 -0.056626 0.998395      39.60 mm
    0.000000 0.000000 0.000000       1.00

Read  64 EEG channels from info
Head coordinate coil definitions created.
Source spaces are now in head coordinates.

Setting up the BEM model using /home/atleeri/mne_data/MNE-fsaverage-data/fsaverage/bem/fsaverage-5120-5120-5120-bem-sol.fif...

Loading surfaces...

Loading the solution matrix...

Three-layer model surfaces loaded.
Loaded linear collocation BEM solution from /home/atleeri/mne_data/MNE-fsaverage-data/fsaverage/bem/fsaverage-5120-5120-5120-bem-sol.fif
Employing the head->MRI coordinate transform with the BEM model.
BEM model fsaverage-5120-5120-5120-bem-sol.fif is now set up


Source spaces are in head coordinates.
Checking that the sources are inside the surface (will take a few...)
Checking surface interior status for 10242 points...
    Found  2433/10242 points inside  an interior sphere of radius   47.7 mm
    Found     0/10242 points outside an exterior sphere of radius   98.3 mm
    Found     0/ 7809 points outside using surface Qhull
    Found     0/ 7809 points outside using solid angles
    Total 10242/10242 points inside the surface
Interior check completed in 7953.0 ms
Checking surface interior status for 10242 points...
    Found  2241/10242 points inside  an interior sphere of radius   47.7 mm
    Found     0/10242 points outside an exterior sphere of radius   98.3 mm
    Found     0/ 8001 points outside using surface Qhull
    Found     0/ 8001 points outside using solid angles
    Total 10242/10242 points inside the surface
Interior check completed in 7982.4 ms
Setting up for EEG...
Computing EEG at 20484 source locations (free orientations)...

Finished.

Example: Create an array of time points for a 1 second long recording. The distance between neighboring time points is the inverse of the sampling rate 1/raw.info["sfreq"].

times = np.arange(0, 1, step=1 / raw.info["sfreq"])
times

array([0.     , 0.00625, 0.0125 , ..., 0.98125, 0.9875 , 0.99375],
      shape=(160,))

Exercise: Create an array of time points for a 5 second long recording.

times = np.arange(0, 5, step=1 / raw.info["sfreq"])
times

array([0.     , 0.00625, 0.0125 , ..., 4.98125, 4.9875 , 4.99375],
      shape=(800,))

Example: Simulate source activity at 1 randomly placed dipole for the array of times and plot it. This opens a GUI where you can visualize the source activity by moving the “Time(s)” slider.

Note: Depending on the source location, you may have to rotate the brain to view the source.

stc = simulate_sparse_stc(fwd["src"], n_dipoles=1, times=times)
stc.plot(hemi="both");

Using control points [1.e-07 1.e-07 1.e-07]

Exercise: Simulate source activity at 2 randomly placed dipoles for the array of times and plot it.

Note: By default, all sources have the same 10 Hz oscillation so their time courses will overlap.

Solution

stc = simulate_sparse_stc(fwd["src"], n_dipoles=2, times=times)
stc.plot(hemi="both");

Using control points [1.e-07 1.e-07 1.e-07]

Example: Instead of placing dipoles randomly, use the labels parameter to place 2 dipoles at the center of specific brain regions — left primary motor cortex ("precentral-lh") and right auditory cortex ("transversetemporal-rh") and plot the source activity.

regions = ["precentral-lh", "transversetemporal-rh"]
src_labels = [l for l in labels if l.name in regions]
stc = simulate_sparse_stc(fwd["src"], n_dipoles=2, times=times, labels=src_labels, location="center")
stc.plot(hemi="both")

Using control points [1.e-07 1.e-07 1.e-07]

Exercise: Place 2 dipoles at the center of the brain regions defined below — left visual association cortex ("lateraloccipital-lh") and right prefrontal cortex ("superiorfrontal-rh") and plot the source activity.

regions = ["lateraloccipital-lh", "superiorfrontal-rh"]

Solution

src_labels = [l for l in labels if l.name in regions]
stc = simulate_sparse_stc(fwd["src"], n_dipoles=2, times=times, labels=src_labels, location="center")
stc.plot(hemi="both")

Using control points [1.e-07 1.e-07 1.e-07]

Example: Simulate source activity at 2 randomly placed dipoles and use simulate_raw to project the source activity through the forward model onto the EEG sensors and plot the resulting simulated raw data.

stc = simulate_sparse_stc(fwd["src"], n_dipoles=2, times=times)
raw_sim = simulate_raw(raw.info, stc, forward=fwd, verbose=False)
raw_sim.plot();

Exercise: Simulate source activity at 15 randomly placed dipoles and use simulate_raw to project the source activity through the forward model onto the EEG sensors and plot the resulting simulated raw data. How does this compare to the simulation with 2 dipoles?

Solution

stc = simulate_sparse_stc(fwd["src"], n_dipoles=15, times=times)
raw_sim = simulate_raw(raw.info, stc, forward=fwd, verbose=False)
raw_sim.plot();

Background

The simulated data so far contains only the signal from the neural sources, which makes it look unrealistically clean. Real EEG data always contains noise from various sources. This noise isn’t purely random — it has a certain correlational structure because nearby channels will pick up similar noise patterns. In MNE, this can be modeled by creating a noise covariance matrix (make_ad_hoc_cov) where the standard deviation of the covariance matrix determines the strength of the noise. In addition to sensor noise, EEG recordings are contaminated by physiological artifacts such as eye blinks, which can be added using add_eog. Adding noise is important for testing analysis methods because it lets you study how well they perform at different signal-to-noise ratios.

Exercises

In this section, you are going to add different types and levels of noise to simulated data and observe the effect on the signal. Here are the code examples required to solve the exercises:

Run the cell below to create a fresh simulation with 2 dipoles that we will use as a baseline for adding noise. We store this clean data in raw_clean so we can create noisy copies of it in the exercises below.

times = np.arange(0, 30, step=1 / raw.info["sfreq"])
stc = simulate_sparse_stc(fwd["src"], n_dipoles=2, times=times, random_state=42)
raw_clean = simulate_raw(raw.info, stc, forward=fwd, verbose=False)
raw_clean.plot();

Example: Add a small amount of colored noise (1 µV, i.e. 1e-6) to a copy of the clean data and plot the noisy raw data. Compare the raw data and power spectrum to the clean version above.

raw_noisy = raw_clean.copy()
cov = make_ad_hoc_cov(raw_noisy.info, std=dict(eeg=1e-6))
add_noise(raw_noisy, cov)
raw_noisy.plot();

Adding noise to 64/64 channels (64 channels in cov)

Exercise: Add a large amount of noise (10 µV, i.e. 10e-6) to a copy of the clean data and plot the noisy raw data. Can you still identify the 10 Hz signal of the simulated sources?

Solution

raw_noisy = raw_clean.copy()
cov = make_ad_hoc_cov(raw_noisy.info, std=dict(eeg=10e-6))
add_noise(raw_noisy, cov)
raw_noisy.plot();

Adding noise to 64/64 channels (64 channels in cov)

Exercise: Add a medium amount of noise (3 µV, i.e. 3e-6) to a copy of the clean data and plot the noisy raw data. How does the signal-to-noise ratio compare to the previous example?

Solution

raw_noisy = raw_clean.copy()
cov = make_ad_hoc_cov(raw_noisy.info, std=dict(eeg=3e-6))
add_noise(raw_noisy, cov)
raw_noisy.plot();

Adding noise to 64/64 channels (64 channels in cov)

Exercise: Add simulated eye blink artifacts to the noisy raw data using add_eog and plot the result. Eye blinks produce large, slow deflections that are most prominent in frontal channels.

Note: since blinks are random events you may have to rerun the simulation if you can’t see any blink artifacts.

Solution

raw_eog = raw_noisy.copy()
add_eog(raw_eog)
raw_eog.plot();

Sphere                : origin at (0.0 0.0 0.0) mm
              radius  : 0.1 mm
Source location file  : dict()
Assuming input in millimeters
Assuming input in MRI coordinates

Positions (in meters) and orientations
2 sources
blink simulated and trace not stored
Setting up forward solutions

Background

So far, simulate_sparse_stc has been using its default waveform, a 10 Hz sinusoid. However, you can pass any custom function to the data_fun parameter. This function must take a time array as input and return an activation waveform of the same length. This allows you to simulate specific neural phenomena such as oscillations at different frequencies or evoked responses. When simulating multiple sources with different waveforms, each source must be simulated separately and then combined using utils.combine_stcs, which merges the source estimates before projecting them onto the EEG sensors.

Exercises

In this section, you are going to plot and simulate three predefined waveforms, combine them into a single source estimate, and optionally constrain the sources to specific brain regions. Here are the code examples required to solve the exercises:

Example: The function defined below returns an alpha_wave with a constant 10 Hz frequency. Generate an array of times from 0 to 5 seconds with a step size of 1/raw.info["sfreq"] and use it to plot the alpha_wave.

def alpha_wave(times):
    """Simulate a 10 Hz alpha oscillation."""
    return 1e-7 * np.sin(2 * np.pi * 10 * times)

times = np.arange(0, 5, 1/raw.info["sfreq"])
plt.plot(times, alpha_wave(times))

Exercise: The function defined below returns a theta_wave with oscillating amplitude. Generate an array of times from 0 to 10 seconds with a step size of 1/raw.info["sfreq"] and use it to plot the theta_wave.

def theta_wave(times):
    return 1e-7 * (1 + np.sin(2 * np.pi * 1 * times)) * np.sin(2 * np.pi * 6 * times)

Solution

times = np.arange(0, 10, 1/raw.info["sfreq"])
plt.plot(times, theta_wave(times))

Exercise: The function defined below returns an evoked_response with a negative and positive peak that repeats every second. Generate an array of times from 0 to 8 seconds with a step size of 1/raw.info["sfreq"] and use it to plot the evoked_response.

def evoked_response(times):
    """Simulate an evoked response with a single positive and negative peak repeating every 1 s."""
    t_local = (times % 1.0) - 0.1   # repeat every 1 s, peak at 100 ms
    return 0.5e-6 * np.sin(2 * np.pi * 5 * t_local) * np.exp(-t_local**2 / 0.005)

Solution

times = np.arange(0, 8, 1/raw.info["sfreq"])
plt.plot(times, evoked_response(times))

Example: Simulate source activity at 10 random dipoles using the alpha_wave function as data_fun. Then, simulate raw EEG from the source activity and plot it.

stc = simulate_sparse_stc(fwd["src"], n_dipoles=10, times=times, data_fun=alpha_wave)
raw = simulate_raw(raw.info, stc, forward=fwd, verbose=False)
raw.plot();

Exercise: Simulate source activity at 2 random dipoles using the theta_wave function as data_fun. Then, simulate raw EEG from the source activity and plot it.

Solution

stc = simulate_sparse_stc(fwd["src"], n_dipoles=2, times=times, data_fun=theta_wave)
raw = simulate_raw(raw.info, stc, forward=fwd, verbose=False)
raw.plot();

Exercise: Simulate source activity at 1 random dipole using the evoked_response function as data_fun. Then, simulate raw EEG from the source activity and plot it.

Solution

stc = simulate_sparse_stc(fwd["src"], n_dipoles=1, times=times, data_fun=evoked_response)
raw = simulate_raw(raw.info, stc, forward=fwd, verbose=False)
raw.plot();

Example: Generate two single-dipole source signals using the alpha_wave and theta_wave functions as data_fun and use utils.combine_stcs to combine them. Then simulate raw EEG from the combined source activity and plot it.

stc1 = simulate_sparse_stc(fwd["src"], n_dipoles=1, times=times, data_fun=alpha_wave)
stc2 = simulate_sparse_stc(fwd["src"], n_dipoles=1, times=times, data_fun=theta_wave)
stc = utils.combine_stcs(stc1, stc2)
raw = simulate_raw(raw.info, stc, forward=fwd, verbose=False)
raw.plot();

Exercise: Generate two single-dipole source signals using the alpha_wave and evoked_response functions as data_fun and use utils.combine_stcs to combine them. Then simulate raw EEG from the combined source activity and plot it.

Solution

stc1 = simulate_sparse_stc(fwd["src"], n_dipoles=1, times=times, data_fun=alpha_wave)
stc2 = simulate_sparse_stc(fwd["src"], n_dipoles=1, times=times, data_fun=evoked_response)
stc = utils.combine_stcs(stc1, stc2)
raw = simulate_raw(raw.info, stc, forward=fwd, verbose=False)
raw.plot();

Exercise: Generate three single-dipole source signals using the alpha_wave, theta_wave and evoked_response functions as data_fun and use utils.combine_stcs to combine them. Then simulate raw EEG from the combined source activity and plot it.

Solution

stc1 = simulate_sparse_stc(fwd["src"], n_dipoles=1, times=times, data_fun=alpha_wave)
stc2 = simulate_sparse_stc(fwd["src"], n_dipoles=1, times=times, data_fun=theta_wave)
stc3 = simulate_sparse_stc(fwd["src"], n_dipoles=1, times=times, data_fun=evoked_response)
stc = utils.combine_stcs(stc1, stc2, stc3)
raw = simulate_raw(raw.info, stc, forward=fwd, verbose=False)
raw.plot();

Bonus: Create your own custom data_fun (for example, combine an alpha and theta oscillation) and use it to generate source activity. Then simulate raw EEG from the source activity and plot the result.

Solution

def custom_wave(times):
    alpha = 1e-7 * np.sin(2 * np.pi * 10 * times)
    theta = 1e-7 * np.sin(2 * np.pi * 6 * times)
    return alpha + theta

stc = simulate_sparse_stc(fwd["src"], n_dipoles=1, times=times, data_fun=custom_wave)
raw = simulate_raw(raw.info, stc, forward=fwd, verbose=False)
raw.plot();

Bonus: Put everything you learned in this session together:

1. Simulate one source in the left auditory cortex ("transversetemporal-lh") using the evoked_response function as data_fun
2. Simulate one source in the right visual cortex ("pericalcarine-rh") using the alpha_wave function as data_fun
3. Use utils.combine_stcs to combine the source signals and simulate raw EEG from the combined source activity
4. Add sensor noise and EOG artifacts
5. Plot the simulated raw EEG signal

Solution

# auditory source
auditory_label = [l for l in labels if l.name == "transversetemporal-lh"][0]
stc_auditory = simulate_sparse_stc(
    fwd["src"], n_dipoles=1, times=times, data_fun=evoked_response,
    labels=[auditory_label], location="center"
)

# visual source
visual_label = [l for l in labels if l.name == "pericalcarine-rh"][0]
stc_alpha = simulate_sparse_stc(
    fwd["src"], n_dipoles=1, times=times, data_fun=alpha_wave,
    labels=[visual_label], location="center"
)

# combine sources and simulate raw EEG
stc = utils.combine_stcs(stc_auditory, stc_alpha)
raw_sim = simulate_raw(raw.info, stc, forward=fwd, verbose=False)

# add noise and EOG
cov = make_ad_hoc_cov(raw_sim.info, std=dict(eeg=3e-6))
add_noise(raw_sim, cov)
add_eog(raw_sim)

# plot result
raw_sim.plot();

Adding noise to 64/64 channels (64 channels in cov)

Sphere                : origin at (0.0 0.0 0.0) mm
              radius  : 0.1 mm
Source location file  : dict()
Assuming input in millimeters
Assuming input in MRI coordinates

Positions (in meters) and orientations
2 sources
blink simulated and trace not stored
Setting up forward solutions