Brain network exploration with neurolib
In this example, we will run a parameter exploration of a whole-brain 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.
# 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)
!pip install matplotlib
import matplotlib.pyplot as plt
import numpy as np
# a nice color map
plt.rcParams['image.cmap'] = 'plasma'
!pip install neurolib
from neurolib.models.aln import ALNModel
from neurolib.utils.loadData import Dataset
import neurolib.utils.functions as func
ds = Dataset("hcp")
import mopet
We load a model with parameters that generate interesting dynamics.
model = ALNModel(Cmat = ds.Cmat, Dmat = ds.Dmat)
model.params['duration'] = 0.2*60*1000
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 run it to see what kind of output it produces!
model.run(bold=True, chunkwise=True)
plt.plot(model.output.T);
We simualted the model with BOLD output, so let's compute the functional connectivity (fc) matrix:
plt.imshow(func.fc(model.BOLD.BOLD[:, model.BOLD.t_BOLD > 5000]))
This is our multi-stage evaluation function.
def evaluateSimulation(params):
model.params.update(params)
defaultDuration = model.params['duration']
invalid_result = {"fc" : [0]* len(ds.BOLDs)}
logging.info("Running stage 1")
# -------- stage wise simulation --------
# Stage 1 : simulate for a few seconds to see if there is any activity
# ---------------------------------------
model.params['duration'] = 3*1000.
model.run()
# check if stage 1 was successful
amplitude = np.max(model.output[:, model.t > 500]) - np.min(model.output[:, model.t > 500])
if amplitude < 0.05:
invalid_result = {"fc" : 0}
return invalid_result
logging.info("Running stage 2")
# Stage 2: simulate BOLD for a few seconds to see if it moves
# ---------------------------------------
model.params['duration'] = 20*1000.
model.run(bold = True, chunkwise=True)
if np.std(model.BOLD.BOLD[:, 5:10]) < 0.0001:
invalid_result = {"fc" : -1}
return invalid_result
logging.info("Running stage 3")
# Stage 3: full and final simulation
# ---------------------------------------
model.params['duration'] = defaultDuration
model.run(bold = True, chunkwise=True)
# -------- evaluation here --------
scores = []
for i, fc in enumerate(ds.FCs):#range(len(ds.FCs)):
fc_score = func.matrix_correlation(func.fc(model.BOLD.BOLD[:, 5:]), fc)
scores.append(fc_score)
meanScore = np.mean(scores)
result_dict = {"fc" : meanScore}
return result_dict
We test run the evaluation function.
model.params['duration'] = 20*1000.
evaluateSimulation(model.params)
# NOTE: These values are low for testing
model.params['duration'] = 10*1000.
explore_params = {"a": np.linspace(0, 40.0, 2)
,"K_gl": np.linspace(100, 400, 2)
,"sigma_ou" : np.linspace(0.1, 0.5, 2)
}
# we need this random filename to avoid testing clashes
hdf_filename = f"exploration-{np.random.randint(99999)}.h5"
ex = mopet.Exploration(evaluateSimulation, explore_params, default_params=model.params, hdf_filename=hdf_filename)
ex.run()
ex.load_results(as_dict=True)
ex.results
ex.params
ex.df
sigma_selectors = np.unique(ex.df.sigma_ou)
for s in sigma_selectors:
df = ex.df[(ex.df.sigma_ou == s)]
pivotdf = df.pivot_table(values='fc', index = 'K_gl', columns='a')
plt.imshow(pivotdf, \
extent = [min(df.a), max(df.a),
min(df.K_gl), max(df.K_gl)], origin='lower', aspect='auto')
plt.colorbar(label='Mean correlation to empirical rs-FC')
plt.xlabel("a")
plt.ylabel("K_gl")
plt.title("$\sigma_{ou}$" + "={}".format(s))
plt.show()
