Modeling resting-state MEG-Data
In this example we will learn how to use neurolib to simulate resting state functional connectivity of MEG recordings.
In the first part of the notebook, we will compute the frequency specific functional connectivity matrix of an examplary resting state MEG recording from the YouR-Study Uhlhaas, P.J., Gajwani, R., Gross, J. et al. The Youth Mental Health Risk and Resilience Study (YouR-Study). BMC Psychiatry 17, 43 (2017).
To this end we will:
- Band-Pass filter the signal
- Apply the
hilbert-transformation to extract the signal envelope - Orthogonalize the signal envelopes of two examplary regions
- Low-Pass filter the signal envelopes
- and compute the pairwise envelope correlations which yields the
functional connectivitymatrix.
We follow the approach presented in Hipp, J., Hawellek, D., Corbetta, M. et al., Large-scale cortical correlation structure of spontaneous oscillatory activity. Nat Neurosci 15, 884–890 (2012)
In the second part of this notebook, we will use a whole-brain model to simulate brain activity and compute functional connectivity matrix of the simulated signal envelope, as was done for the empirical MEG data. The parameters of this model have been previously optimized with neurolib's evolutionary algorithms (not shown here).
Finally, we will compute the fit (Pearson correlation) of the simulated functional connectivity to the empirical MEG data, which was used as a fitting objective in a previous optimization procedure.
# change to the root directory of the project
import os
if os.getcwd().split("/")[-1] == "examples":
os.chdir('..')
# This will reload all imports as soon as the code changes
%load_ext autoreload
%autoreload 2
import os
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import seaborn as sns
import ipywidgets as widgets
from IPython.utils import io
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
import time
import pandas as pd
Empirical Functional Connectivity
Load MEG-Data
First off, let's load the MEG data using the Signal class from neurolib. Our example data has already been preprocessed and projected into source space using the AAL2 atlas.
from neurolib.utils.signal import Signal
signal = Signal.from_file(os.path.join('examples', 'data','rs-meg.nc'))
region_labels = signal.data.regions.values
nr_regions = len(region_labels)
display(signal.data)
Band-Pass filter and Hilbert transform
We will now filter the signal into the desidered frequency band and apply the hilbert transform on the band-passed filtered signal. This will provide us with the analytic representation of the signal, which we can then use to extract the signal's envelope and its phase.
In the following, we plot each processing step for an example target region that you can chose using the widgets below (default: left Precentral Gyrus). Furthermore, we can also choose the frequency range that we'd like to filter the signal in (default: alpha (8-12Hz)).
print('Select a region from the AAL2 atlas and a frequency range')
# Select a Region
target = widgets.Select(options=region_labels, value='PreCG.L', description='Regions',
tooltips=['Description of slow', 'Description of regular', 'Description of fast'],
layout=widgets.Layout(width='50%', height='150px'))
display(target)
# Select Frequency Range
freq = widgets.IntRangeSlider(min=1, max=46, description='Frequency (Hz)', value=[8, 12], layout=widgets.Layout(width='80%'),
style={'description_width': 'initial'})
display(freq)
# Define how many timepoints you'd like to plot
plot_timepoints = 1000
# Plot unfiltered Signal
fig, ax = plt.subplots(2, 1, figsize=(12,8), sharex=True)
sns.lineplot(x=signal.data.time[:plot_timepoints].values, y=signal.data.sel(regions=target.value)[:plot_timepoints].values,
ax=ax[0], color='k', alpha=0.6)
ax[0].set_title(f'Unfiltered Signal ({target.value})');
# Band Pass Filter the Signal
signal.filter(freq.value[0], freq.value[1], inplace=True);
# Apply hilbert-transform to extract the signal envelope
complex_signal = signal.hilbert_transform('complex', inplace=False)
signal_env = np.abs(complex_signal.data)
# Plot filtered Signal and Signal Envelope
sns.lineplot(x=signal.data.time[:plot_timepoints].values, y=signal.data.sel(regions=target.value)[:plot_timepoints].values,
ax=ax[1], label='Bandpass-Filtered Signal')
sns.lineplot(x=signal_env.time[:plot_timepoints].values, y=signal_env.sel(regions=target.value)[:plot_timepoints].values,
ax=ax[1], label='Signal Envelope')
ax[1].set_title(f'Filtered Signal ({target.value})');
ax[1].legend(bbox_to_anchor=(1.2, 1),borderaxespad=0)
sns.despine(trim=True)
Orthogonalized signal envelope
Now we are going to address the main methodological issue of MEG when it comes to the analysis of the cortical functional connectivity structure, i.e. its low spatial resolution. The electric field generated by any given neural source spreads widely over the cortex so that the signal captured at the MEG sensors is a complex mixture of signals from multiple underlying neural sources.
To account for the effect of electric field spread on our MEG connectivity measures, we adapted the orthogonalization approach by Hipp, J., Hawellek, D., Corbetta, M. et al. Large-scale cortical correlation structure of spontaneous oscillatory activity. Nat Neurosci 15, 884–890 (2012) link.
The basic idea here is that a signal generated by one neural source and measured at two separate sensors must have exactly the same phase at both sensors. In contrast, signals from different neural sources have different phases. And thus it is possible to eliminate the effect of a reference signal on the target signal by removing the signal component that has the same phase as a reference region.
Formally, this can be expressed as: \(Y_{\perp X}(t,f) = imag\big(\ Y(t,f)\ \frac{X(t,f)^\star}{|X(t,f)|}\ \big)\ \label{eq:orth}\). Here, \(Y\) represents the analytic signal from our target regions that is being orthogonalized with respect to the signal from region \(X\).
Using the widgets below, you can choose the reference region \(X\) (default: right Precentral Gyrus)
print('Select a reference region for the orthogonalization')
# Select a Region
referenz = widgets.Select(options=region_labels, value='PreCG.R', description='Regions',
tooltips=['Description of slow', 'Description of regular', 'Description of fast'],
layout=widgets.Layout(width='50%', height='150px'))
display(referenz)
# Perform Orthogonalization
signal_conj = complex_signal.data.conj()
conj_div_env = signal_conj/signal_env
orth_signal = (complex_signal.data.sel(regions=target.value) * conj_div_env.sel(regions=referenz.value)).imag
orth_env = np.abs(orth_signal)
# Plot
fig, ax = plt.subplots(2,1,figsize=(12,8), sharex=True)
sns.lineplot(x=signal.data.time[:plot_timepoints].values, y=signal.data.sel(regions=referenz.value)[:plot_timepoints].values, ax=ax[0])
sns.lineplot(x=signal_env.time[:plot_timepoints].values, y=signal_env.sel(regions=referenz.value)[:plot_timepoints].values, ax=ax[0])
ax[0].set_title(f'Referenz Region X ({referenz.value})');
sns.lineplot(x=signal.data.time[:plot_timepoints].values, y=signal.data.sel(regions=target.value)[:plot_timepoints].values,
ax=ax[1], label='Bandpass-Filtered Signal')
sns.lineplot(x=signal_env.time[:plot_timepoints].values, y=signal_env.sel(regions=target.value)[:plot_timepoints].values,
ax=ax[1], label='Signal Envelope')
sns.lineplot(x = orth_env.time[:plot_timepoints].values, y=orth_env[:plot_timepoints].values, ax=ax[1], label='Orthogonalized Envelope')
ax[1].legend(bbox_to_anchor=(1.2, 1),borderaxespad=0)
ax[1].set_title(f'Target Region Y ({target.value})');
sns.despine(trim=True)
Low-Pass filtering of the envelopes
As a last step, before calculating the envelope correlations, we need to low-pass filter the signal envelopes since the connectivity measures of (ultra)-low frequency components of the MEG-signal correspond best to the functional connectivity as measured using fMRI.
Below, you can choose the low-pass frequency (default: 0.2 Hz).
low_pass = widgets.FloatSlider(value=0.2, min=0, max=2.0, step=0.1, description='Low-Pass Frequency (Hz)',
disabled=False, readout=True, readout_format='.1f', layout=widgets.Layout(width='80%'),
style={'description_width': 'initial'})
display(low_pass)
with io.capture_output() as captured:
low_orth_env = Signal(orth_env).filter(low_freq=None, high_freq=low_pass.value, inplace=False)
low_signal_env = Signal(signal_env.sel(regions=referenz.value)).filter(low_freq=None, high_freq=low_pass.value, inplace=False)
# Plot
fig, ax = plt.subplots(1,2,figsize=(15,4), sharey=True)
sns.lineplot(x=signal_env.time[:plot_timepoints].values, y=signal_env.sel(regions=referenz.value)[:plot_timepoints].values, ax=ax[0])
sns.lineplot(x=low_signal_env.data.time[:plot_timepoints].values, y=low_signal_env.data[:plot_timepoints].values, ax=ax[0])
ax[0].set_title(f'Referenz Region X ({referenz.value})');
sns.lineplot(x = orth_env.time[:plot_timepoints].values, y=orth_env[:plot_timepoints].values, ax=ax[1], label='Orthogonalized Envelope')
sns.lineplot(x = low_orth_env.data.time[:plot_timepoints].values, y=low_orth_env.data[:plot_timepoints].values, ax=ax[1], label='Low-Passed Orthogonalized Envelope')
ax[1].legend(bbox_to_anchor=(1, -0.18),borderaxespad=0)
ax[1].set_title(f'Target Region Y ({target.value})');
sns.despine(trim=True)
print(f'Orthogonalized envelope correlation between {referenz.value} and {target.value}: ', np.round(np.corrcoef(low_orth_env.data,low_signal_env.data)[0,1],2))
Computing the functional connectivity matrix
We will now define a function that iterates over each pair of brain regions and performs the previously presented processing steps, i.e. that extracts the envelopes, performs the orthogonalization, applies the low-pass filter, and returns the functional connectivity matrix that contains the pairwise envelope correlations.
This step may take a minute.
def orth_fc(signal, low_pass):
nr_regions = signal.data.shape[0]
progress = widgets.IntProgress(min=0, max=nr_regions, description=('Calculating FC Matrix'),
layout=widgets.Layout(width='80%'), style={'description_width': 'initial'})
display(progress)
complex_signal = signal.hilbert_transform('complex', inplace=False)
signal_env = signal.hilbert_transform('amplitude', inplace=False);
conj_div_env = complex_signal.data.conj()/signal_env.data
# Low-pass filter Signal envelope
with io.capture_output() as captured:
signal_env.filter(low_freq=None, high_freq=low_pass)
corr = []
for complex_region in complex_signal.data:
orth_signal = (complex_region * conj_div_env).imag
orth_env = np.abs(orth_signal).T
orth_env = Signal(orth_env)
with io.capture_output() as captured:
orth_env.filter(low_freq=None, high_freq=low_pass)
corr_mat = np.corrcoef(orth_env.data, signal_env.data)
corr.append(np.diag(corr_mat, k=nr_regions))
progress.value += 1
fc = np.array(corr)
# Since the orthogonalization process is not symmetric we take the mean of both directions.
fc = (fc.T + fc) / 2.
np.fill_diagonal(fc,0)
return fc
# Execute Function
fc = orth_fc(signal, low_pass.value)
Let's now plot the functional connectivity matrix. We label only every second row/column since right and left regions alternate in the AAL2 atlas.
fig, ax = plt.subplots(figsize=(10,8))
sns.heatmap(fc, square=True, ax=ax, cmap='YlGnBu', linewidth=0.005, cbar_kws={"shrink": .8})
ticks = [tick[:-2] for tick in region_labels[::2]]
ax.set_xticks(np.arange(0,94,2)); ax.set_yticks(np.arange(0,94,2))
ax.set_xticklabels(ticks, rotation=90, fontsize=8); ax.set_yticklabels(ticks, rotation=0, fontsize=8);
Exclude subcortical regions
For the following whole-brain simulation we are only interested in the cortical regions. So we'll now exclude all subcortical regions: * Hippocampus: 41 - 44 * Amygdala: 45-46 * Basal Ganglia: 75-80 * Thalamus: 81-82
Attention: AAL indices start with 1
exclude = list(range(40, 46)) + list(range(74, 82))
tmp = np.delete(fc, exclude, axis=0)
emp_fc = np.delete(tmp, exclude, axis=1)
# Exclude regions from the list of region labels
emp_labels = np.delete(region_labels, exclude)
Whole-brain model
In this part of the notebook, we will use neurolib to simulate the functional connectivity. We will therefore:
- Load structural connectivity matrices from the Human Connectome Project and initiate the whole-brain model using the Wilson-Cowan model to simulate each brain region
- Set the global coupling strength, exc. background input, and the noise strength parameters of the model
- Run the simulation
- Compute the functional connectivity using the signal envelopes
Please refer to the wc-minimal example for an introduction to the Wilson-Cowan model.
Initiate whole-brain model
# Let's import the neurolib
from neurolib.models.wc import WCModel
from neurolib.utils.loadData import Dataset
# First we load the structural data set from the Human Connectome Project
ds = Dataset("hcp")
# We initiate the Wilson-Cowan model
wc = WCModel(Cmat = ds.Cmat, Dmat = ds.Dmat, seed=0)
Parameter settings
You may now choose parameters settings for the global coupling, the excitatory background input, and the noise strength, which will be used when we run the model. The final fit between simulated and empirical connectivity matrices will depend on the parameters choosen here.
global_coupling = widgets.FloatSlider(value=6.55, min=0., max=20.0, step=0.01, description='Global Coupling',
disabled=False, readout=True, readout_format='.2f', layout=widgets.Layout(width='80%'),
style={'description_width': 'initial'})
exc_drive = widgets.FloatSlider(value=1.58, min=0.0, max=4.0, step=0.01, description='Exc. Background Drive',
disabled=False, readout=True, readout_format='.2f', layout=widgets.Layout(width='80%'),
style={'description_width': 'initial'})
inh_drive = widgets.FloatSlider(value=2.83, min=0.0, max=4.0, step=0.01, description='Inh. Background Drive',
disabled=False, readout=True, readout_format='.2f', layout=widgets.Layout(width='80%'),
style={'description_width': 'initial'})
noise_level = widgets.FloatSlider(value=0.02, min=0.001, max=0.05, step=0.001, description='Noise Level',
disabled=False, readout=True, readout_format='.3f', layout=widgets.Layout(width='80%'),
style={'description_width': 'initial'})
display(global_coupling)
display(exc_drive)
display(inh_drive)
display(noise_level)
Run the simulation
Let's now run the whole-brain model using the defined parameter settings. This may take some time since we're simulating a complete minute here.
# Let's set the previously defined parameters
# note: the duraiton here is short for testing:
wc.params['duration'] = 10*1000
# use longer simulation for real run:
#wc.params['duration'] = 1*60*1000
wc.params['K_gl'] = global_coupling.value
wc.params['exc_ext'] = exc_drive.value
wc.params['inh_ext'] = inh_drive.value
wc.params['sigma_ou'] = noise_level.value
# Run the model
wc.run()
Simulated functional connectivity
We'll now compute the functional connectivity matrix containing the pairwise envelope correlations between all cortical regions of the AAL2 atlas. We'll thus follow the processing steps as before, i.e. band-pass filter the signal, extract the signal envelopes using the hilbert transformation, low-pass filter the envelopes and compute the pairwise pearson correlations. Note that we don't apply the orthogonalization scheme here, since this was only done to account to the electric field spread in the empirical data.
# Create xr DataArray from the simulated excitatory timeseries (keeping the region labels)
sim_signal = xr.DataArray(wc.exc[:, int(1000/wc.params.dt):], dims=("regions", "time"), coords={"regions": emp_labels, "time": wc.t[int(1000/wc.params.dt):]/1000},
attrs={'atlas':'AAL2'})
# Initialize Figure
fig, ax = plt.subplots(figsize=(12,4))
# Filter signal
sim_signal = Signal(sim_signal)
sim_signal.resample(to_frequency=100)
with io.capture_output() as captured:
sim_signal.filter(freq.value[0], freq.value[1], inplace=True);
sns.lineplot(x=sim_signal.data.time[:plot_timepoints].values, y=sim_signal.data.sel(regions=target.value)[:plot_timepoints].values, ax=ax, label='Filtered Signal')
# Extract signal envelope
sim_signal.hilbert_transform('amplitude', inplace=True)
sns.lineplot(x=sim_signal.data.time[:plot_timepoints].values, y=sim_signal.data.sel(regions=target.value)[:plot_timepoints].values, ax=ax, label='Signal Envelope')
# Low-Pass Filter
with io.capture_output() as captured:
sim_signal.filter(low_freq=None, high_freq=low_pass.value, inplace=True)
sns.lineplot(x=sim_signal.data.time[:plot_timepoints].values, y=sim_signal.data.sel(regions=target.value)[:plot_timepoints].values, ax=ax, label='Low-Pass Signal Envelope')
ax.legend(bbox_to_anchor=(1.2, 1),borderaxespad=0)
ax.set_title(f'Simulated Signal of Target Region Y ({target.value})');
sns.despine(trim=True)
To compute the simulated functional connectivity matrix we use the fc functions from neurolib.
import neurolib.utils.functions as func
# Compute the functional connectivity matrix
sim_fc = func.fc(sim_signal.data)
# Set diagonal to zero
np.fill_diagonal(sim_fc, 0)
# Plot Empirical and simulated connectivity matrix
fig, ax = plt.subplots(1,2, figsize=(16,10))
sns.heatmap(emp_fc, square=True, ax=ax[0], cmap='YlGnBu', linewidth=0.005, cbar_kws={"shrink": .5})
ax[0].set_title('Empirical FC',pad=10);
sns.heatmap(sim_fc, square=True, ax=ax[1], cmap='YlGnBu', linewidth=0.005, cbar_kws={"shrink": .5})
ax[1].set_title('Simulated FC',pad=10);
ticks = [tick[:-2] for tick in emp_labels[::2]]
for ax in ax:
ax.set_xticks(np.arange(0,80,2)); ax.set_yticks(np.arange(0,80,2))
ax.set_xticklabels(ticks, rotation=90, fontsize=8); ax.set_yticklabels(ticks, rotation=0, fontsize=8);
Model fit
Lastly, we may evaluate the model fit by computing the pearson correlation between our simulated functional connectivity matrix and the empirical one. Additionally we'll also plot the correlation between structural and functional connectivity matrices to have a reference.
# Compare structural and simulated connectivity to the empirical functional connectivity
struct_emp = np.corrcoef(emp_fc.flatten(), ds.Cmat.flatten())[0,1]
sim_emp = np.corrcoef(emp_fc.flatten(), sim_fc.flatten())[0,1]
# Plot
fig, ax = plt.subplots(figsize=(6,6))
splot = sns.barplot(x=['Structural Connectivity', 'Simulated Connectivity'], y=[struct_emp, sim_emp], ax=ax)
ax.set_title('Correlation to Empiral Functional Connectivity', pad=10)
for p in splot.patches:
splot.annotate(format(p.get_height(), '.2f'),
(p.get_x() + p.get_width() / 2., p.get_height()),
ha = 'center', va = 'center',
size=20, color='white',
xytext = (0, -12),
textcoords = 'offset points')
sns.despine()
print(f"Parameters: \tGlobal Coupling: {wc.params['K_gl']}\n\t\tExc. Background Drive: {wc.params['exc_ext']}")
print(f"\t\tNoise Level: {wc.params['sigma_ou']}")
