The neural mass model
In this example, we will learn about the basic of neurolib. We will create a two-population mean-field model of exponential integrate-and-fire neurons called the aln model. We will learn how to create a Model, set some parameters and run a simulation. We will also see how we can easily access the output of each simulation.
aln - the adaptive linear-nonlinear cascade model
The adaptive linear-nonlinear (aln) cascade model is a low-dimensional population model of spiking neural networks. Mathematically, it is a dynamical system of non-linear ODEs. The dynamical variables of the system simulated in the aln model describe the average firing rate and other macroscopic variables of a randomly connected, delay-coupled network of excitatory and inhibitory adaptive exponential integrate-and-fire neurons (AdEx) with non-linear synaptic currents.
Ultimately, the model is a result of various steps of model reduction starting from the Fokker-Planck equation of the AdEx neuron subject to white noise input at many steps of input means \(\mu\) and variances \(\sigma\). The resulting mean firing rates and mean membrane potentials are then stored in a lookup table and serve as the nonlinear firing rate transfer function, \(r = \Phi(\mu, \sigma)\).
Basic use
# 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
try:
import matplotlib.pyplot as plt
except ImportError:
import sys
!{sys.executable} -m pip install matplotlib
import matplotlib.pyplot as plt
import numpy as np
import scipy
# Let's import the aln model
from neurolib.models.aln import ALNModel
# Some useful functions are provided here
import neurolib.utils.functions as func
# a nice color map
plt.rcParams['image.cmap'] = 'plasma'
Simulating a single aln node
To create a single node, we simply instantiate the model without any arguments.
# Create the model
model = ALNModel()
# Each model comes with a set of default parameters which are a dictionary.
# Let's change the parameter that controls the duration of a simulation to 10s.
model.params['duration'] = 10.0 * 1000
# For convenience, we could also use:
model.params.duration = 10.0 * 1000
# In the aln model an Ornstein-Uhlenbeck process is simulated in parallel
# as the source of input noise fluctuations. Here we can set the variance
# of the process.
# For more info: https://en.wikipedia.org/wiki/Ornstein%E2%80%93Uhlenbeck_process
# Let's add some noise.
model.params['sigma_ou'] = 0.1
# Finally, we run the model
model.run()
Accessing the outputs
Accessing the outputs is straight-forward. Every model's outputs
are stored in the model.outputs attribute. According to the specific
name of each of the model's outputs, they can also be accessed as a
key of the Model object, i.e. model['rates_exc'].
plt.plot(model['t'], model['rates_exc'].T, lw=2, c='k')
plt.xlabel("t [ms]")
plt.ylabel("Rate [Hz]")
plt.xlim(1000, 2000);
# Outputs are also available as an xr DataArray
xr = model.xr()
print(xr.dims)
# outputs can also be accessed via attributes in dot.notation
print("rates_exc", model.rates_exc)
Bifurcation diagram
Bifurcation diagrams can give us an overview of how different parameters of the model affect its dynamics. The simplest method for drawing a bifurcation diagram is to simply change relevant parameters step by step and record the model's behavior in response to these changes. In this example, we want to see how the model's dynamics changes with respect to the external input currents to the excitatory population. These input currents could be due to couplings with other nodes in a brain network or we could model other factors like external electrical stimulation.
Below, you can see a schematic of the aln model. As you can see, a single node consists of one excitatory (red) and one inhibitory population (blue). The parameter that controls the mean input to the excitatory population is \(\mu_{E}\) or model.params["mue_ext_mean"] .
Let's first decrease the duration of a single run so we can scan the parameter space a bit faster and let's also disable the noisy input.
model.params['duration'] = 2.0 * 1000
model.params['sigma_ou'] = 0.0
Let's fix the input to the inhibitory population:
model.params['mui_ext_mean'] = 0.5
We draw a one-dimensional bifurcation diagram, so it is enough to loop through different values of mue_ext_mean and record the minimum and maximum of the rate for each parameter.
max_rate_e = []
min_rate_e = []
# these are the different input values that we want to scan
mue_inputs = np.linspace(0, 2, 50)
for mue in mue_inputs:
# Note: this has to be a vector since it is input for all nodes
# (but we have only one node in this example)
model.params['mue_ext_mean'] = mue
model.run()
# we add the maximum and the minimum of the last second of the
# simulation to a list
max_rate_e.append(np.max(model.output[0, -int(1000/model.params['dt']):]))
min_rate_e.append(np.min(model.output[0, -int(1000/model.params['dt']):]))
Let's plot the results!
plt.plot(mue_inputs, max_rate_e, c='k', lw = 2)
plt.plot(mue_inputs, min_rate_e, c='k', lw = 2)
plt.title("Bifurcation diagram of the aln model")
plt.xlabel("Input to excitatory population")
plt.ylabel("Min / max firing rate")
Whole-brain model
neurolib comes with some example datasets for exploring its functionality. Please be aware that these datasets are not tested and should not be used for your research, only for experimentation with the software.
A dataset for whole-brain modeling can consists of the following parts:
- A structural connectivity matrix capturing the synaptic connection strengths between brain areas, often derived from DTI tractography of the whole brain. The connectome is then typically parcellated in a preferred atlas (for example the AAL2 atlas) and the number of axonal fibers connecting each brain area with every other area is counted. This number serves as a indication of the synaptic coupling strengths between the areas of the brain.
- A delay matrix which can be calculated from the average length of the axonal fibers connecting each brain area with another.
- A set of functional data that can act as a target for model optimization. Resting-state fMRI offers an easy and fairly unbiased way for calibrating whole-brain models. EEG data could be used as well.
We can load a Dataset by passing the name of it in the constructor.
from neurolib.utils.loadData import Dataset
ds = Dataset("gw")
We now create the aln model with a structural connectivity matrix and a delay matrix. In order to achieve a good fit of the BOLD activity to the empirical data, the model has to run for quite a while. A a rule of thumb, a simulation of resting-state BOLD activity should not be shorter than 3 minutes and preferably longer than 5 minutes real time. If the empirical recordings are for example 10 minutes long, ideally, a simulation of 10 minutes would be used to compare the output of the model to the resting state recording.
model = ALNModel(Cmat = ds.Cmat, Dmat = ds.Dmat)
model.params['duration'] = 0.4*60*1000
# Info: value 0.4*60*1000 is low for testing
# use 5*60*1000 for real simulation
After some optimization to the resting-state fMRI data of the dataset, we found a set of parameters that creates interesting whole-brain dynamics. We set the mean input of the excitatory and the inhibitory population to be close to the E-I limit cycle.
model.params['mue_ext_mean'] = 1.57
model.params['mui_ext_mean'] = 1.6
# We set an appropriate level of noise
model.params['sigma_ou'] = 0.09
# And turn on adaptation with a low value of spike-triggered adaptation currents.
model.params['b'] = 5.0
Let's have a look what the data looks like. We can access the
data of each model by calling its internal attributes.
Here, we plot the structural connectivity matrix by calling
model.params['Cmat'] and fiber length matrix by calling
model.params['lengthMat'].
Of course, we can also access the dataset using the Dataset
object itself. For example the functional connectivity matrices
of the BOLD timeseries in the datasets are given as list with
ds.FCs.
from matplotlib.colors import LogNorm
fig, axs = plt.subplots(1, 3, figsize=(12,8), dpi=75)
fig.subplots_adjust(wspace=0.28)
im = axs[0].imshow(model.params['Cmat'], norm=LogNorm(vmin=10e-5, vmax=np.max(model.params['Cmat'])))
axs[0].set_title("Cmat")
fig.colorbar(im, ax=axs[0],fraction=0.046, pad=0.04)
im = axs[1].imshow(model.params['lengthMat'], cmap='inferno')
axs[1].set_title("Dmat")
fig.colorbar(im, ax=axs[1],fraction=0.046, pad=0.04)
im = axs[2].imshow(ds.FCs[0], cmap='inferno')
axs[2].set_title("Empirical FC")
fig.colorbar(im, ax=axs[2],fraction=0.046, pad=0.04)
Run model
We run the model with bold simulation by using bold=True.
This simulates the Balloon-Windkessel BOLD model in parallel to
the neural population model in order to estimate the blood oxygen
levels of the underlying neural activity. The output of the bold
model can be used to compare the simulated data to empirical fMRI
data (resting-state fMRI for example).
To save (a lot of) RAM, we can run the simulation in chunkwise mode.
In this mode, the model will be simulated for a length of chunksize
steps (not time in ms, but actual integration steps!), and the output
of that chunk will be used to automatically reinitiate the model with
the appropriate initial conditions. This allows for a serial continuation
of the model without having to store all the data in memory and is
particularly useful for very long and many parallel simulations.
model.run(chunkwise=True, chunksize=60000, bold=True)
Results
The outputs of the model can be accessed using the attribute model.outputs
model.outputs
For convenience, they can also be accessed directly using attributes of the model with the outputs name, like model.rates_exc. The outputs are also available as xr DataArrays as model.xr().
Since we used bold=True to simulate BOLD, we can also access model.BOLD.BOLD for the actual BOLD activity, and model.BOLD.t for the time steps of the BOLD simulation (which are downsampled to 0.5 Hz by default).
Plot simulated activity
# Plot functional connectivity and BOLD timeseries (z-scored)
fig, axs = plt.subplots(1, 2, figsize=(7, 3), dpi=75, gridspec_kw={'width_ratios' : [1, 2]})
axs[0].imshow(func.fc(model.BOLD.BOLD[:, 5:]))
axs[1].imshow(scipy.stats.mstats.zscore(model.BOLD.BOLD[:, model.BOLD.t_BOLD>10000], axis=1), aspect='auto', extent=[model.BOLD.t_BOLD[model.BOLD.t_BOLD>10000][0], model.BOLD.t_BOLD[-1], 0, model.params['N']]);
axs[0].set_title("FC")
axs[0].set_xlabel("Node")
axs[0].set_ylabel("Node")
axs[1].set_xlabel("t [ms]")
# the results of the model are also accessible through an xarray DataArray
fig, axs = plt.subplots(1, 1, figsize=(6, 2), dpi=75)
plt.plot(model.xr().time[-10000:], model.xr().loc['rates_exc'].T[-10000:])
plt.xlabel("t [ms]");
Correlation of simulated BOLD to empirical data
We can compute the element-wise Pearson correlation of the functional connectivity matrices of the simulated data to the empirical data to estimate how well the model captures the inter-areal BOLD correlations found in empirical resting-state recordings.
scores = [func.matrix_correlation(func.fc(model.BOLD.BOLD[:, 5:]), fcemp) for fcemp in ds.FCs]
print("Correlation per subject:", [f"{s:.2}" for s in scores])
print(f"Mean FC/FC correlation: {np.mean(scores):.2}")