Parameter exploration with neurolib
In this example, we will draw a bifurcation diagram of a neural mass model that we load using the brain simulation framework neurolib. Please visit the Github repo to learn more about this library or read the gentle introduction to neurolib to learn more about the neuroscience background of neural mass models and whole-brain simulations.
What we will do
We will scan through a 2-dimensional parameter space and record all outputs of the model to a hdf file. We will then load the simulated results from the hdf file and condense the output to a single scalar so we can plot it. Let's get to it!
We'll first start with some unnecessary but useful imports for logging and auto-reloading code. You can skip this block if you don't know what it does.
# change into the root directory of the project
import os
if os.getcwd().split("/")[-1] == "examples":
os.chdir('..')
%load_ext autoreload
%autoreload 2
import logging
logger = logging.getLogger()
logger.setLevel(logging.INFO)
We have to install matplotlib and neurolib for this example, which we can do using the shell command syntax of jupyter using !pip install matplotlib.
# install matplotlib
!pip install matplotlib
import matplotlib.pyplot as plt
import numpy as np
# a nice color map
plt.rcParams['image.cmap'] = 'plasma'
# install neurolib
!pip install neurolib
from neurolib.models.aln import ALNModel
import mopet
We can now load the ALNModel model from neurolib.
model = ALNModel()
model.params.duration = 1*1000
Like in other examples, we need to define an evalFunction that mopet can call. The parameters params that are passed to that function are the parameters of the model to run. We load the parameters to the model using a simple dictionary update old_dict.update(new_dict). We then return the output of the model.
def evalFunction(params):
model.params.update(params)
model.run()
return model.outputs
Now we define the parameter ranges to explore (the parameters here are called mue_ext_mean and mui_ext_mean which represent the background input to two neural populations).
# NOTE: These values are low for testing!
explore_params = {"mue_ext_mean" : np.linspace(0, 3, 3),
"mui_ext_mean" : np.linspace(0, 3, 3)}
# For a real run, use these values:
# explore_params = {"mue_ext_mean" : np.linspace(0, 3, 31),
# "mui_ext_mean" : np.linspace(0, 3, 31)}
# we need this random filename to avoid testing clashes
hdf_filename = f"exploration-{np.random.randint(99999)}.h5"
ex = mopet.Exploration(evalFunction, explore_params, default_params=model.params, hdf_filename=hdf_filename)
Everything is ready and we can now simply run the exploration.
ex.run()
After the exploration is done, we can load the results into a pandas Dataframe. By adding the argument arrays=True, we also tell mopet to load all simulated output (including arrays!) of the exploration into the Dataframe.
ex.load_results(arrays=True, as_dict=True)
Because we have used as_dict=True, the outputs of each run are also stored in the results dictionary.
ex.results
Reducing the results
As you can see, the results above are time series. For each simulation, we've got multiple arrays. We would like to visualize the results somehow. However, it can be quite challenging to plot many time series in a single figure and still understand what's happening. Therefore, to reduce the dimensionality of the data to a single scalar number per simulation, we will loop through the entire results Dataframe, grab the time series called rates_exc and simply compute its maximum after a transient time of t>500 time steps. We then store this number in the Dataframe and simply call the new column result.
ex.df["result"] = None
for r in ex.df.index:
t = ex.results[r]['t']
rates_exc = ex.results[r]['rates_exc']
ex.df.loc[r, "result"] = np.max(rates_exc[:, t>500])
We can now inspect the updated dataframe using the attribute .df.
ex.df
As you can see, the column result only has a single number for each run. Perfect for plotting! In order to plot this long-format Dataframe, we need to pivot it to create a new Dataframe pivoted which has the correct x and y coordinates for plotting the results value.
pivoted = ex.df.pivot_table(values='result', index = 'mui_ext_mean', columns='mue_ext_mean', aggfunc='first')
Let's have a look at the new Dataframe that we're about to plot.
pivoted
Perfect, that's exactly what we need. We have two indices, which are both equal to the paramters that we have explored before. The only entry of the Dataframe is the results value that we've put in before. Now, let's plot this Dataframe using matplotlib.imshow().
plt.imshow(pivoted, \
extent = [min(ex.df.mue_ext_mean), max(ex.df.mue_ext_mean),
min(ex.df.mui_ext_mean), max(ex.df.mui_ext_mean)], origin='lower')
plt.colorbar(label='Maximum firing rate')
plt.xlabel("Input to E")
plt.ylabel("Input to I")
