Simulate resting state dynamics in mouse brain

This demo shows how to simulate and analyze resting state dynamics in mouse brain using as connectome a tracer-based connectome built thanks to the Allen Connectivity Builder.

The results showed here are discussed in Melozzi et al., 2016 [1]

First, we import all the required dependencies

In [1]:
from tvb.interfaces.command.lab import *
from tvb.simulator.lab import *
LOG = get_logger('demo')
from tvb.simulator.plot.tools import *
import numpy as np
import pylab
import matplotlib.pyplot as plt
%matplotlib inline
matplotlib.rcParams['figure.figsize'] = (20.0, 10.0)
2022-05-12 11:49:17,404 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_adapters.DistanceDBIN'>
2022-05-12 11:49:17,407 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_adapters.DistanceDWEI'>
2022-05-12 11:49:17,408 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_adapters.DistanceNETW'>
2022-05-12 11:49:17,409 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_adapters.DistanceRDA'>
2022-05-12 11:49:17,410 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_adapters.DistanceRDM'>
2022-05-12 11:49:17,411 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_adapters.ModularityOCSM'>
2022-05-12 11:49:17,413 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_adapters.ModularityOpCSMU'>
2022-05-12 11:49:17,418 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_centrality_adapters.CentralityEdgeBinary'>
2022-05-12 11:49:17,419 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_centrality_adapters.CentralityEdgeWeighted'>
2022-05-12 11:49:17,421 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_centrality_adapters.CentralityEigenVector'>
2022-05-12 11:49:17,422 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_centrality_adapters.CentralityKCoreness'>
2022-05-12 11:49:17,423 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_centrality_adapters.CentralityKCorenessBD'>
2022-05-12 11:49:17,424 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_centrality_adapters.CentralityNodeBinary'>
2022-05-12 11:49:17,425 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_centrality_adapters.CentralityNodeWeighted'>
2022-05-12 11:49:17,426 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_centrality_adapters.CentralityShortcuts'>
2022-05-12 11:49:17,428 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_centrality_adapters.FlowCoefficients'>
2022-05-12 11:49:17,429 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_centrality_adapters.ParticipationCoefficient'>
2022-05-12 11:49:17,430 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_centrality_adapters.ParticipationCoefficientSign'>
2022-05-12 11:49:17,432 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_centrality_adapters.SubgraphCentrality'>
2022-05-12 11:49:17,437 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_clustering_adapters.ClusteringCoefficient'>
2022-05-12 11:49:17,438 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_clustering_adapters.ClusteringCoefficientBU'>
2022-05-12 11:49:17,439 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_clustering_adapters.ClusteringCoefficientWD'>
2022-05-12 11:49:17,441 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_clustering_adapters.ClusteringCoefficientWU'>
2022-05-12 11:49:17,443 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_clustering_adapters.TransitivityBinaryDirected'>
2022-05-12 11:49:17,444 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_clustering_adapters.TransitivityBinaryUnDirected'>
2022-05-12 11:49:17,445 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_clustering_adapters.TransitivityWeightedDirected'>
2022-05-12 11:49:17,447 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_clustering_adapters.TransitivityWeightedUnDirected'>
2022-05-12 11:49:17,455 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_degree_adapters.Degree'>
2022-05-12 11:49:17,457 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_degree_adapters.DegreeIOD'>
2022-05-12 11:49:17,458 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_degree_adapters.DensityDirected'>
2022-05-12 11:49:17,459 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_degree_adapters.DensityUndirected'>
2022-05-12 11:49:17,460 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_degree_adapters.JointDegree'>
2022-05-12 11:49:17,463 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_degree_adapters.MatchingIndex'>
2022-05-12 11:49:17,464 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_degree_adapters.Strength'>
2022-05-12 11:49:17,465 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_degree_adapters.StrengthISOS'>
2022-05-12 11:49:17,467 - WARNING - tvb.config.init.introspector_registry - Skipped Adapter(probably because MATLAB not found):<class 'tvb.adapters.analyzers.bct_degree_adapters.StrengthWeights'>
2022-05-12 11:49:19,898 - INFO - alembic.runtime.migration - Context impl SQLiteImpl.
2022-05-12 11:49:19,899 - INFO - alembic.runtime.migration - Will assume non-transactional DDL.
2022-05-12 11:49:19,924 - INFO - tvb.config.init.model_manager - Database already has some data, will not be re-created!
   INFO  log level set to INFO

The Connectome

In order to built the mouse brain network we used a tracer-based connectome.

In particular we used a structural connectivity matrix (stored in the data folder of TVB), which is built thanks to the Allen Connectivity Builder in TVB. The Allen Connectivity Builder is a tool that download and manipulate the open-source tracer experiments of the Allen Institute (Oh et al., 2014 [2]) in order to built a connectome and the corresponding parcelled volume according to the preferences of the user.

The user can choose:

  • the resolution of the grid volume in which the experimental data have been registered (here 100 $\mu m$).
  • The definition of the connection strength between source region $i$ and target region $j$. (here $w_{ij}=\frac{PD_j}{ID_i}$, where PD=projection density, ID=injection density)

It is possible to choose the characteristics of the brain areas to be included in the parcellation using the two following criteria:

  • Only brain areas where at least one injection has infected more than a given threshold of voxels. This kind of selection ensures that only the data with a certain level of experimental relevance is included in the connectome (Oh et al., 2014[2]), (here 50 voxels).
  • Only brain areas that have a volume greater than a given threshold can be included (here 2$mm^3$).

In the following the connectome is loaded and plotted.

In [3]:
from tvb.basic.readers import try_get_absolute_path
connectivity_path = try_get_absolute_path("tvb_data","mouse/allen_2mm/Connectivity.h5")
import_op = import_conn_h5(1, connectivity_path)
import_op = wait_to_finish(import_op)
import_op
2022-05-12 11:50:12,890 - INFO - tvb.core.services.operation_service - Starting operation TVBImporter
2022-05-12 11:50:13,116 - INFO - tvb.core.services.backend_clients.standalone_client - Start processing operation id:43
2022-05-12 11:50:32,330 - INFO - tvb.core.services.backend_clients.standalone_client - Finished with launch of operation 43
2022-05-12 11:50:32,331 - INFO - tvb.core.services.backend_clients.standalone_client - Return code: 0. Stopped: False
2022-05-12 11:50:32,332 - INFO - tvb.core.services.backend_clients.standalone_client - Thread: <OperationExecutor(Thread-117, initial)>
2022-05-12 11:50:32,345 - INFO - tvb.core.services.operation_service - Finished operation launch:TVBImporter
2022-05-12 11:50:33,367 - INFO - tvb.interfaces.command.lab - Operation finished successfully
Out[3]:
<Operation('80fbfebe462a4fb3990a581da01cac88', 89554d5c-d1d0-11ec-ac16-782b46fbd208, 2,'1','29','2022-05-12 11:50:12.911446','2022-05-12 11:50:20.336023', '2022-05-12 11:50:31.740440','5-FINISHED',True, 'None', '2022-05-12,11-50-31', '', 0)>
In [4]:
list_operation_results(import_op.id)
              id                     type                              gid         date
              39        ConnectivityIndex d9be05424e8011e69bf02c4138a1c4ef 2016-07-20 15:49:59.831038
In [5]:
# Copy the id of the ConnectivityIndex obtained above
conn = load_dt(39)
In [6]:
from tvb.basic.readers import try_get_absolute_path

# Visualize the structural connectivity matrix
plt.subplots()
cs=plt.imshow(np.log10(conn.weights), cmap='jet', aspect='equal', interpolation='none')
plt.title('Structural connectivity matrix', fontsize=20)
axcb=plt.colorbar(cs)
axcb.set_label('Log10(weights)', fontsize=20)
C:\Users\romina.baila\.conda\envs\tvb-run\lib\site-packages\ipykernel_launcher.py:5: RuntimeWarning: divide by zero encountered in log10
  """

The simulation

Once the brain network is defined is possible to simulate its activity. Here we simulate resting state dynamics using the reduced Wong Wang model (Deco et al. 2013 [3], Hansen et al., 2015 [4]).

In order to convert the synaptic activity in BOLD signals we used the Balloon-Windkessel method (Friston et al., 200 [5]) using the default value implemented in The Virtual Brain.

In [7]:
list_projects()
                    name    id
         Default_Project     1
In [9]:
from tvb.core.services.algorithm_service import AlgorithmService
from tvb.core.services.simulator_service import SimulatorService
from tvb.core.entities.model.model_burst import BurstConfiguration
from tvb.config.init.introspector_registry import IntrospectionRegistry
from tvb.core.entities.file.simulator.view_model import SimulatorAdapterModel, EulerStochasticViewModel, BoldViewModel, AdditiveNoiseViewModel
from time import sleep

# define the neural mass model (here: reduced wong wang)
RWW = models.ReducedWongWang(w=numpy.array([1.0]), I_o=numpy.array([0.3]))

#define variables to monitor during simulation (here: BOLD activity)
monitor = BoldViewModel()
monitor.period=2e3

#define long range coupling parameter
longcoupling = coupling.Linear(a=numpy.array([0.096]))

#define duration of simulation in ms
duration=1200e3

#define integrator
integrator = EulerStochasticViewModel()
integrator.dt = 0.1
integrator.noise = AdditiveNoiseViewModel(nsig=np.array([0.00013]))

# Instantiate a SimulatorAdapterModel and configure it
simulator_model = SimulatorAdapterModel()
simulator_model.model=RWW
# Copy ConnectivityIndex gid from the result of the list_operation_results function
simulator_model.connectivity = "d9be05424e8011e69bf02c4138a1c4ef"
simulator_model.simulation_length = duration
simulator_model.coupling = longcoupling
simulator_model.integrator = integrator
simulator_model.monitors = [monitor]

# use id of your current project as first argument
launched_operation = fire_simulation(1, simulator_model)


launched_operation = wait_to_finish(launched_operation)

launched_operation
2022-05-12 11:51:16,821 - INFO - tvb.core.services.backend_clients.standalone_client - Start processing operation id:44
2022-05-12 11:51:16,851 - INFO - tvb.interfaces.command.lab - Operation launched ....
2022-05-12 13:12:29,810 - INFO - tvb.interfaces.command.lab - Operation finished successfully
Out[9]:
<Operation('14fe6b50ae954ae9937787c837f11793', ae723512-d1d0-11ec-851a-782b46fbd208, 2,'1','13','2022-05-12 11:51:15.176577','2022-05-12 11:51:31.351324', '2022-05-12 13:12:29.689230','5-FINISHED',True, 'None', '2022-05-12,13-12-29', '', 516.796875)>

The simulated bold signals can be visualized using matplotlib library.

In [10]:
list_operation_results(launched_operation.id)
              id                     type                              gid         date
              41   SimulationHistoryIndex 40587ce41c4442da8b11c7863bec7377 2022-05-12 13:12:29.240490
              42    TimeSeriesRegionIndex be01c9ee1b4444c2b4b2b1c955af3ca1 2022-05-12 13:12:29.529542
2022-05-12 13:12:30,201 - INFO - tvb.core.services.backend_clients.standalone_client - Finished with launch of operation 44
2022-05-12 13:12:30,203 - INFO - tvb.core.services.backend_clients.standalone_client - Return code: 0. Stopped: False
2022-05-12 13:12:30,204 - INFO - tvb.core.services.backend_clients.standalone_client - Thread: <OperationExecutor(Thread-151, started 15132)>
In [12]:
#Load time series h5 file
ts = load_dt(42) # use the id of the TimeSeriesRegionIndex obtained above
bold_time = ts.time
bold_data = ts.data
In [13]:
# Display the simulated bold timeseries
plt.subplots()
plt.plot(bold_time,bold_data[:,0,:,0])
plt.xlabel('Time (ms)', fontsize=20)
plt.ylabel('Amplitude (au)', fontsize=20)
plt.title('Simulated BOLD timeseries', fontsize=20)
Out[13]:
Text(0.5, 1.0, 'Simulated BOLD timeseries')

Analysis

The simulated BOLD signals can be analyzed in different way.

Functional Connectivity Dynamics

In particular here we focus on the Functional Connectivity Dynamics (FCD) a metric which is able to quantify the evolution of the functional states in time. There are several ways to estimate FCD (for a review Preti et al., 2016 [6]), TVB uses the sliding windows technique.

In order to estimate the FCD using the sliding window technique, the entire BOLD time-series is divided in time windows of a fixed length (3 min) and with an overlap of 176 s; the data points within each window centered at the time $t_i$ were used to calculate FC($t_i$). The \emph{ij}-th element of the FCD matrix is calculated as the Pearson correlation between the upper triangular part of the $FC(t_i)$ matrix arranged as a vector and the upper triangular part of the $FC(t_j)$ matrix arranged as a vector.

The FCD matrix allows identifying the epochs of stable FC configurations as blocks of elevated inter-$FC(t)$ correlation; these blocks are organized around the diagonal of the FCD matrix (Hansen et al., 2015 [4]).

In order to identify the epochs of stable FC configurations, TVB uses the spectral embedding method, that permits to group together the nodes of the FCD, i.e. the different time windows, in clusters.

In [15]:
# Run FCD Adapter in order to compute the FCD Matrix
from tvb.adapters.analyzers.fcd_adapter import FCDAdapterModel,FunctionalConnectivityDynamicsAdapter
from tvb.adapters.datatypes.db.time_series import TimeSeriesRegionIndex
from tvb.adapters.datatypes.h5.time_series_h5 import TimeSeriesRegionH5
from tvb.core.neocom import h5
from tvb.core.entities.storage import dao
# from tvb.core.entities.file.files_helper import FilesHelper
from tvb.core.adapters.abcadapter import ABCAdapterForm, ABCAdapter
from tvb.core.services.operation_service import OperationService
from time import sleep

adapter_instance = ABCAdapter.build_adapter_from_class(FunctionalConnectivityDynamicsAdapter)

# Create and evaluate the analysis
# build FCDAdapterModel
fcd_model = adapter_instance.get_view_model_class()()
fcd_model.time_series= ts.gid
fcd_model.sw=180e3 # windows length (ms)
fcd_model.sp=4e3 # spanning between sliding windows (ms)

# launch an operation and have the results stored both in DB and on disk
launched_operation = fire_operation(1, adapter_instance, fcd_model)

launched_operation = wait_to_finish(launched_operation)

launched_operation
2022-05-12 13:39:59,619 - INFO - tvb.core.services.operation_service - Starting operation FunctionalConnectivityDynamicsAdapter
2022-05-12 13:39:59,752 - INFO - tvb.core.services.backend_clients.standalone_client - Start processing operation id:45
2022-05-12 13:39:59,777 - INFO - tvb.core.services.operation_service - Finished operation launch:FunctionalConnectivityDynamicsAdapter
2022-05-12 13:39:59,778 - INFO - tvb.interfaces.command.lab - Operation launched....
2022-05-12 13:40:09,508 - INFO - tvb.core.services.backend_clients.standalone_client - Finished with launch of operation 45
2022-05-12 13:40:09,511 - INFO - tvb.core.services.backend_clients.standalone_client - Return code: 0. Stopped: False
2022-05-12 13:40:09,512 - INFO - tvb.core.services.backend_clients.standalone_client - Thread: <OperationExecutor(Thread-168, started 19452)>
2022-05-12 13:40:09,937 - INFO - tvb.interfaces.command.lab - Operation finished successfully
Out[15]:
<Operation('4e5aec8fe7014c1f8d8fb0e07711c2b2', df52b902-d1df-11ec-b253-782b46fbd208, 2,'1','3','2022-05-12 13:39:59.629542','2022-05-12 13:40:04.063964', '2022-05-12 13:40:09.018226','5-FINISHED',True, 'None', '2022-05-12,13-40-08', '', 0)>
In [16]:
list_operation_results(launched_operation.id)
              id                     type                              gid         date
              43                 FcdIndex 9fbc1f3d3ca14547b6d6773ad12815f3 2022-05-12 13:40:08.574846
              44 ConnectivityMeasureIndex a24c649466db4374b9cef900b5af3016 2022-05-12 13:40:08.684809
              45 ConnectivityMeasureIndex 239cdcb273aa41588f86dd9e78673c27 2022-05-12 13:40:08.801459
              46 ConnectivityMeasureIndex 55e2dc384174482cbb406663ecf7f686 2022-05-12 13:40:08.901247

The original and segmented FCD matrices can be visualized using the matplotlib library.

In [17]:
# Plot the FCD matrix and the FCD matrix segmented in the epochs
FCD = load_dt(43).array_data[:,:,0,0] # use the id of the FcdIndex obtained above

# If we have just one FCDIndex as a result of the FCD Adapter it means that the FCD Segmented is the same as FCD
FCD_SEGMENTED = FCD
    
plt.subplot(121)
cs=plt.imshow(FCD, cmap='jet', aspect='equal')
axcb =plt.colorbar(ticks=[0, 0.5, 1])
axcb.set_label(r'CC [FC($t_i$), FC($t_j$)]', fontsize=20)
cs.set_clim(0, 1.0)
for t in axcb.ax.get_yticklabels():
     t.set_fontsize(18)
plt.xticks([0,len(FCD)/2-1, len(FCD)-1],['0','10', '20'], fontsize=18)
plt.yticks([0,len(FCD)/2-1, len(FCD)-1],['0','10', '20'], fontsize=18)
plt.xlabel(r'Time $t_j$ (min)', fontsize=20)
plt.ylabel(r'Time $t_i$ (min)', fontsize=20)
plt.title('FCD', fontsize=20)

plt.subplot(122)
cs=plt.imshow(FCD_SEGMENTED, cmap='jet', aspect='equal')
axcb =plt.colorbar(ticks=[0, 0.5, 1])
axcb.set_label(r'CC [FC($t_i$), FC($t_j$)]', fontsize=20)
cs.set_clim(0, 1.0)
for t in axcb.ax.get_yticklabels():
     t.set_fontsize(18)
plt.xticks([0,len(FCD)/2-1, len(FCD)-1],['0','10', '20'], fontsize=18)
plt.yticks([0,len(FCD)/2-1, len(FCD)-1],['0','10', '20'], fontsize=18)
plt.xlabel(r'Time $t_j$ (min)', fontsize=20)
plt.ylabel(r'Time $t_i$ (min)', fontsize=20)
plt.title('FCD segmented', fontsize=20)
Out[17]:
Text(0.5, 1.0, 'FCD segmented')