Abstract
Connectome-based modeling of large-scale brain network dynamics enables causal in silico interrogation of the brain’s structure-function relationship, necessitating the close integration of diverse neuroinformatics fields. Here we extend the open-source simulation software The Virtual Brain (TVB) to whole mouse brain network modeling based on individual diffusion magnetic resonance imaging (dMRI)-based or tracer-based detailed mouse connectomes. We provide practical examples on how to use The Virtual Mouse Brain (TVMB) to simulate brain activity, such as seizure propagation and the switching behavior of the resting state dynamics in health and disease. TVMB enables theoretically driven experimental planning and ways to test predictions in the numerous strains of mice available to study brain function in normal and pathological conditions.
Significance Statement
The Virtual Mouse Brain (TVMB) represents a versatile and intuitive tool for virtualizing and simulating mouse whole-brain dynamics using a connectome-based model approach. TVMB enables the construction of individual mouse brains using diffusion magnetic resonance imaging (dMRI) experiments. TVMB also allows building detailed connectomes based on tracer experiments realized by the Allen Institute. Various modalities can be modeled [e.g., electroencephalography (EEG) and fMRI]. The platform can be used to generate predictions that can be tested experimentally.
Introduction
Dedicated software environments are available to simulate detailed neuronal dynamics such as Neuron, Genesis, and MOOSE, which model the complex dendrite geometry, reaction-diffusion processes and receptor distributions of individual neurons and smaller networks (Hines and Carnevale, 1997). To simulate larger networks, neuron models are reduced to point neurons (Izhikevich, 2003; Brette and Gerstner, 2005), which is the case of the simulator BRIAN (Goodman and Brette, 2009), NEST (Diesmann and Gewaltig, 2001), and NENGO (Eliasmith et al., 2012). However, scaling up for detailed models beyond an entire cortical column (Markram, 2012) and a few brain regions becomes quickly intractable even for networks of point neurons. Although neuromorphic computation offers interesting alternatives for the future (the SpiNNaker project http://apt.cs.manchester.ac.uk/projects/SpiNNaker/ and the BrainScaleS project http://brainscales.kip.uni-heidelberg.de/), macroscopic modeling using neural population approaches is the only viable whole-brain network modeling strategy nowadays.
The Virtual Brain (TVB) is an open-source simulation software designed to model whole-brain network dynamics, where the network’s connectivity is based on diffusion magnetic resonance imaging (dMRI)-based individual connectomes or adaptations of more precise primate connectomes (Sanz Leon et al., 2013). TVB comprises several generative neural population models, defined in physical 3D space and constrained by anatomy, allowing simulating neuroimaging signals [such as magneto- and electroencephalography (MEG, EEG), or functional MRI (fMRI)]. Whole-brain dynamics can also be manipulated in TVB, e.g., via stimulation. TVB provides a large set of tools for visualization and data analysis (Sanz-Leon et al., 2015). As such, TVB provides a conceptual framework to interpret neuroimaging data, offering promising diagnostic and therapeutic perspectives (Jirsa et al., 2017; Proix et al., 2017). Other groups demonstrate converging results using similar large-scale brain modeling approaches (Hutchison et al., 2013; Sinha et al., 2017). However, very few model predictions can be experimentally tested in humans for obvious ethical reasons. Thus, assessing causality and extracting general principles of brain dynamics in health and disease remains a challenge.
Rodent research enabled major advances in our understanding of brain function and dysfunction, but mostly at the microscopic scale. The advent of new generations of MRI machines now gives access to detailed anatomic, structural and functional information at the whole rodent brain scale (Stafford et al., 2014), thus providing a formidable opportunity to explore general principles of whole-brain dynamics. Indeed, hypotheses can be tested and causality can be assessed in the numerous transgenic mouse lines that have been generated to study neurologic disorders and to manipulate neuronal networks (e.g., with optogenetics and pharmacomogenetics). However, a conceptual framework is needed to interpret neuroimaging data and generate testable hypotheses. Such framework would considerably accelerate our understanding of the mechanisms controlling and affecting whole-brain dynamics.
Here, we present The Virtual Mouse Brain (TVMB), the first connectome-based simulation platform to study large-scale mouse brain dynamics.
TVMB is integrated into TVB and adapted to the mouse brain to enable its virtualization. It uses the validated TVB simulators to generate brain network activity, as well as analysis and visualization tools.
In what follows, we will show how the platform can be used to virtualize not only individual mouse brains (based on dMRI connectome) but also to construct very detailed connectome-based models using tracer data obtained by the Allen Institute for Brain Science (Oh et al., 2014).
As a worked example to show how the platform can be used to generate predictions or interpret data, we will simulate resting state dynamics in a control and “epileptic” mouse and seizure propagation. We will also show how to integrate TVMB in a research project in which theoretical and experimental approaches benefit from one another.
Materials and Methods
All the methods discussed in what follows are implemented in TVB and freely available to the community. The scripts to reproduce the same results presented in the paper are in Extended data (Figs. 1-1, 1-2, 3-1, 4-1); the tracer-based connectome used here are available in the TVB data folder (…/tvb_data/mouse/).
The Allen Connectivity Builder
The Allen Connectivity Builder is a pipeline that we have designed to build a complete mouse connectome based on tracer information.
Specifically we define the link between two brain regions according to the anterograde tracing information provided by the Allen Institute of Brain Science and presented in the work of Oh et al. (2014). In the latter, the axonal projections from a given region are mapped by injecting in adult male C57Bl/6J mice the recombinant adeno-associated virus, which expresses the EGFP anterograde tracer. The tracer migration signal is detected with a serial two-photon tomography system. This approach is repeated systematically to collect the information on the tracer migration from several injection sites in the right hemisphere to target regions in both ipsilateral and contralateral hemispheres; for each injection sites several experiments are run and distinct measures are accomplished. The Allen Institute provides its data through an Internet-accessible interface, namely the Allen Software Development kit (Allen SDK), from which TVB, through The Allen Connectivity Builder interface, is able to obtain a volumetric atlas as well as the raw experimental information necessary to build complete mouse brain connectomes. The platform allows to choose the main characteristics of the connectome; specifically the user can set:
The resolution of the grid volume in which the data are registered (25, 50, and 100 μm).
The definition of the connection strength between source region i and target region j. Specifically the connection strength can be defined as:
– The detected projection density (the number of detected pixels in the target region normalized on the total number of pixels belonging to that region).
– The detected projection energy (the intensity of detected pixels in the target region normalized on the total number of pixels belonging to that region).
– The ratio between the projection density, defined as explained above, and the injection density of the source region (the number of infected pixels in the source region normalized on the total number of pixels belonging to that region).
It is possible to choose the characteristics of the brain areas to be included in the parcellation using the two following criteria:
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).
Only brain areas that have a volume greater than a given threshold can be included.
The pipeline, once downloaded the raw data from the Allen dataset, cleans the data to obtain a set of experiments in which the injection structures are exactly the same as the target structures and vice versa; this step ensures that the connectome will be a square matrix. Then, the pipeline excludes from the experimental set the area that do not fulfill the criteria set by the user (minimum volume (3) and minimum number of voxels infected (4)).
The experiments of the Allen Institute consider source regions always located in the right hemisphere, thus we build a complete structural connectivity matrix, taking the mirror image of the right hemisphere to build the left one. Therefore, if we divide the SC matrix in four blocks R-R, R-L, L-R, and L-L (clockwise order starting from upper left), we will have the symmetries R-R = L-L and R-L = L-R. This assumption is justified by the fact that the mouse brain shows a high degree of lateral symmetry (Calabrese et al., 2015).
The connection strength between a given region and another one is averaged across all the experiments that use as source and target regions those particular brain areas.
The Allen Connectivity Builder approximates the length of the tracts as the Euclidean distance between the centers of the brain regions; the latter are calculated using the volume built from the Allen SDK.
Finally, The Allen Connectivity Builder creates a region volume mapping, i.e., a 3D matrix which represents the volume of the mouse brain, by modifying the annotation volume downloaded from the Allen SDK. In particular the volume is built so that the entries of the 3D volume matrix range from -1 (background) to N-1, where N is the total number of areas in the connectivity: entries in the volume equal to i – 1 label the brain region whose incoming and outgoing connections are organized in the i-th row and i-th column of the connectivity matrix.
The code named allen_creator.py is located in the folder TVB folder/tvb/adapters/creators/. The code is composed of a function “launch” that calls several functions (all written in the file allen_creator.py) that perform all the actions described above. To work outside TVB, users need to copy and paste all the functions that they need from the allen_creator.py.
The volumes and the connectivities used in the present work have a resolution of 100 μm and each connection strength is defined as the ratio between projection density and injection density. The areas included in the parcellation have a volume >2 mm3 and they have >50 voxels infected in at least one injection experiment. The connectome obtained is in the TVB data folder (/TVB_Distribution/tvb_data/lib/python2.7/site-packages/tvb_data/mouse/); the instructions to obtain it through the TVB Jupiter interface are in script Figure 1-1 (Extended data).
Resting state dynamics
Brain model
The mean activity of each brain region, composing the mouse brain network, is described by the reduced Wong Wang model (Wong and Wang, 2006). In this approach, the dynamics of a brain region is given by the whole dynamics of excitatory and inhibitory populations of leaky integrate-and-fire neurons interconnected via NMDA synapses. In this work, we take into account this model with a further reduction performed in Deco et al. (2013b): the dynamics of the output synaptic NMDA gating variable S of a local brain area i is strictly bound to the collective firing rate Hi. The resulting model is given by the following coupled equations:
(1)
(2)
(3)
where xi is the synaptic input to the i-th region. γ is a kinetic parameter fixed to 0.641, τs is the NMDA decay time constant and its value is 100 ms; a, b, and d are the parameters of the input and output function H and are, respectively, equal to 270 nC– 1, 108 Hz, 0.154 s. JN = 0.2609 nA is an intensity scale for the synaptic input current. w is the local excitatory recurrence and I0 is the external input current. G is the coupling strength, i.e., a scalar parameter which scales all the connection strengths Cij without altering the global topology of the network.
The value of the local excitatory recurrence, w, and the external input current, I0, are set, respectively, to 0.3 nA and 1 to enrich the nonlinearity of the dynamics of each brain region. Indeed in this case, studying the dynamics of isolated brain areas (G = 0 in Eq. 3), it is possible to notice that each brain area is in a bistable state and it oscillates between high and low activity fixed points (Hansen et al., 2015). When the brain areas are connected in the network (i.e., G ≠ 0 in Eq. 3), it is not possible to find the analytic solution of 
(Eq. 1), and thus the attractors of the system. The behavior of the connected brain network can be assessed by simulating the system in a deterministic fashion until it reaches the stationary state; repeating this approach several times it is possible to obtain a rich repertoire of stationary states (i.e., attractors of the system) that can be classified through clustering techniques as described previously (Golos et al, 2015). It has been noticed (Hansen et al., 2015) that enriching the nonlinearity of each brain areas has an impact on the global network by introducing attractors that are not in trivial relation with the anatomic connectivity. We refer to this model as the enhanced nonlinearity mean-field model (eMFM). The implementation of the eMFM in a brain network offers the chance to study the nonstationary features of the functional connectivity (FC) patterns.
The values of G, together with the value of the noise amplitude σ of the normally distributed stochastic variable ηi, are tuned, respectively, to 0.096 and 5.1 × 10– 3. In absence of experimental data, the optimal values of the model parameters cannot be assessed through data fitting techniques; we have choose these model parameters since they allow the system to satisfy the two following conditions. (1) The system is in unstable state so that it is able to explore several brain configurations and it has been shown that this is the optimal range for simulating resting state activity (Deco et al., 2013a). The range of coupling strengths for which the system displays multistability is identified by building its bifurcation diagram as described in Deco et al. (2013b). (2) The system is able to reproduce the checkboard pattern of the FC dynamics (FCD) observed in experimental data; this condition is satisfied for a smallest subset of coupling strength values identified with the first condition (Hansen et al., 2015; Deco et al., 2017). We have use the FCD metrics as clue of goodness of our simulation, since it has been identified in the literature as an important resting state feature (Hutchison et al., 2013; Preti et al., 2016).
Integration scheme and BOLD signals
Model equations are numerically solved using the Euler integration method with a fixed integration step of 0.1 ms.
Simulated BOLD signal is obtained by converting the simulated synaptic activity (Eq. 1) using the Balloon-Windkessel method (Friston et al., 2000) with the default value implemented in TVB (Sanz-Leon et al., 2015).
The BOLD time series are down-sampled to 2 s and 20 min total length.
Functional connections
Functional connections in the simulated time series are explored from both spatial and temporal point of views using, respectively, the FC and the FCD.
The ij-th element of the FC matrix is calculated as the Pearson correlation between the BOLD signal of the brain region i and of the brain region j.
To estimate the FCD, 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 ti were used to calculate FC(ti). The ij-th element of the FCD matrix is calculated as the Pearson correlation between the upper triangular part of the FC(ti) matrix arranged as a vector and the upper triangular part of the FC(tj) matrix arranged as a vector.
To observe signal correlations at frequency greater than the typical one of the BOLD signals, the sliding window length is fixed to 3 min, since, as demonstrated by Leonardi and Van De Ville (2015), the nonspurious correlations in the FCD are limited by high-pass filtering of the signals with a cutoff equal to the inverse of the window length.
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).
FCD segmentation: spectral embedding
To identify the epochs of stable FC configurations, we used the spectral embedding method, that permits to group together the nodes of the FCD, i.e., the different time windows, in clusters.
The spectral embedding is a general cluster technique founded on the possibility to map the nodes of the network in the Euclidean space such that the Euclidean distance between the nodes in the space corresponds with the distance between the nodes in the network.
To implement this idea, it is necessary to define the notion of distance between nodes in a network; this is made introducing the concept of the commute distance cij between the nodes i and j that is defined as the expected number of steps in a random walk starting to travel from node i to node j, and back (Von Luxburg, 2007).
To mathematically define cij it is necessary to introduce some quantity. Let us consider a graph that has an adjacency matrix W, i.e., a matrix whose element wij is the weight of the link between node i and j, that in our case is the FCD matrix; it is possible to define the laplacian of the graph as:
(4)
Let us denote with |ei > the eigenvector i of the Laplacian; if the matrix U is the matrix whose columns are the eigenvectors of L, and Λ the diagonal matrix with the eigenvalues λi on the diagonal, thus it is possible to decompose the Laplacian as: L = UΛUT.
The generalized inverse of the Laplacian is defined as 
, where 
is the diagonal matrix with on the diagonal 1/λi when λi is different from zero, otherwise zero. Thanks to 
it is possible to express the commute distance between node i and j as:
(5)
The variable |zi > maps the vertex vi in the Euclidean space (
) such that the Euclidean distance between node i and j is equal to the commute distance cij of the nodes in the graph if and only if:
(6)
from which it follows that < zi| corresponds to the i–th row of the matrix 
.
Functional hubs
The FC matrix of each epoch defines a functional network; for each functional network, we identify the hub regions with an approach analogous to the one used in graph theory for defining the eigenvector centrality of a network node (Newman, 2008).
Let us define the functional centrality 
of a brain region i as the sum of the functional centralities of the neighboring brain regions weighted on the functional connection strength fcij:
(7) x where λ is a constant. Defining the vector 
as the column vector whose components are the functional centrality of each network region, we can rewrite the previous equation in matrix form:
(8)
It is simple to notice that 
is the eigenvector of the FC matrix associated with the eigenvalue λ. Since the FC is a real symmetric matrix (thus diagonalizable), we can decompose it as:
(9)
It follows that the magnitude of the eigenvalue gives a measure of the role of the corresponding eigenvector in reproducing the original matrix.
Taking into account all these observations, we identify the functional hub regions of the mouse brain as the regions with the largest eigenvector components, in absolute value, associated with the three largest eigenvalues of the FC matrix.
The script to run resting state simulation with the tracer-based connectome is in Figure 1-2 (Extended data).
Modeling altered connectomes in pathologic contexts
TVB allows manipulating the structural connectivity by selectively changing the strength of the connections between brain areas to mimic structural lesions. Using this tool, we have simulated mouse brain dynamics mimicking some aspects of the anatomic reorganization found in mesial temporal lobe epilepsy: the neuronal connection lost in hippocampal regions, in particular fields CA1 and CA3 (Esclapez et al., 1999). To reproduce this feature in silico, we have removed all the in-coming and out-coming connections of fields CA1 and CA3 of the hippocampus and then scaled all the connection strengths by a constant factor, so that the total weights of the modified SC is equal to the one of the original matrix. We simulated the resting state BOLD activity and we calculated the FCD matrix as described in the previous sections.
The script to run resting state simulation in pathologic condition is in Figure 3-1 (Extended data).
Epileptic spread in silico
The epileptic network node model
The Epileptor (Jirsa et al., 2014) is a model describing the onset (through a saddle-node bifurcation), the time course and the offset (through a homoclinic bifurcation) of seizures with five state variables that operate at three different time scales. The variable that guides the neural population through the bifurcations is the slow permittivity variable, z, which operates at the slowest time scale. Ensemble 1, comprising the variables x1 and y1, describes the fast discharges registered during ictal states and stable state observed during interictal states; it operates at the fastest timescale. Finally, ensemble 2 (x2, y2) operates at the intermediate time scale and accounts for spike-and-wave events. The interaction between the variables of the system is the following: ensemble 1, through the function g(x1), excites ensemble 2, which in turn inhibits ensemble 1 through 
; both the ensembles are coupled to the slow variable, and the first ensemble acts directly on z. Proix et al. (2014) propose a permittivity coupling between brain areas via a linear difference coupling function that links the fast subsystems with the slow variable z with the weights given by the distance cij.
The full model equations read:
(10)
where:
(11)
(12)
(13)
with 
, τ = 10, γ = 0.01; the permittivity coupling term Ks is fixed to -60.
The degree of epileptogenicity x0 of a brain region i is a parameter that establishes if the region generates seizures autonomously.
Integration scheme and epileptogenicity zone
The epileptogenic zones in the model are nodes of the network that are implemented in the simulation with an epileptogenicity value, x0, so that, for those nodes, the transition between the preictal and the ictal state occurs spontaneously (Proix et al., 2014). An isolated brain area with x0 ≤ 2.06 is epileptogenic, otherwise the area is in its equilibrium state and it can generate seizures only if an external stimulus pushes it through the transition and makes it fall in the propagation zone (Proix et al., 2014).
In the simulated mouse brain, the classification of brain areas in epileptogenic and propagation zone follows the experimental results of the work of Toyoda et al. (2013), in which the authors, using recording electrodes, evaluate the seizures propagation in rats with spontaneous seizures. The authors observe that the earliest seizure activity is recorded most frequently within the hippocampal formation and then spreads, in chronological order, in the subiculum, the entorhinal cortex, the olfactory cortex, the neocortex and the striatum; in 7 over 10 rats analyzed in the paper the epileptogenic region is likely identified in either hemisphere. Accordingly we set as epileptogenic, x0 = –1.9 the left hippocampal regions (field CA1, field CA3, and dentate gyrus), and we set all other regions as propagation zones, x0 = –2.1.
The differential equations of the model are integrated with the Heun stochastic method with an integration step equal to 0.04 ms; we use additive white Gaussian noise in the fast variables (x2 and y2) with mean zero and variance 0.0025. The signals are down-sampled to 1 ms. We set the pre-expression monitor to keep track of the local field potential, defined in the Epileptor as –x1 + x2, as well as the slow permittivity value z.
We define the time at which seizure initiates in a large brain region, as for example the olfactory cortex, as the mean of the seizure onset time of all the network nodes composing that region (see Table 1); to evaluate the chronological order of areas recruitment we define the seizure onset latency of a region as the difference between the time at which the seizure initiates in that region and the time at which the seizure has started in the epileptogenic zone, i.e., the hippocampal regions.
The script to simulate epileptic activity in mouse brain is in Figure 4-1 (Extended Data).
Results
Virtualizing the mouse brain
Tracer-based connectome
To exploit present (and future) high-resolution structural information of the Allen Institute, we designed The Allen Connectivity Builder, a pipeline, which uploads their raw data and processes it to create a connectome and its brain volume representation. The user chooses four sets of parameters: the resolution of the tracing data (1); the way the connection strengths are calculated (2); and the criteria used to include or not a given injected region based on its volume (3); and its experimental significance (4). The pipeline then computes automatically the averaged connection strength between any two regions. Since injections were only performed in the right hemisphere and since the mouse brain shows a high degree of lateral symmetry (Calabrese et al., 2015), the pipeline uses the mirror image to build the left hemisphere. If time delays are considered as an important variable to simulate whole-brain activity, the length of each axonal tract becomes a key parameter. The Allen Connectivity Builder approximates the length of the tracts as the Euclidean distance between the centers of two regions. Finally, the pipeline automatically builds the brain volume using the same parcellation as used to build the connectome. TVMB includes a region volume mapping visualizer to display the brain volume as sections and the results of the computations in the brain sections.
An example of a structural connectivity matrix obtained through the Allen connectivity builder is shown in Figure 1A and the corresponding volume sections in Figure 1D.
A, We used The Allen Connectivity Builder to build the structural connectivity matrix. The color map represents the connections strengths with a base-ten logarithmic scale. The resolution of the grid is 100 μm; the weights of the matrices are defined as the ratio between the projection and the injection density; all the areas in the parcellation have at least one injection experiment that has infected >50 voxels in those areas; the matrix contains only regions with a volume >2 mm3. B, Simulated resting state BOLD time series using the connectome built in A, and the eMFM to model the dynamics of each brain area. C, FCD matrix obtained from the time series. The three black segments (I, II, and III) correspond to epochs of stability of the FCD identified with the spectral embedding technique. D, Functional hubs detected in silico mapped on brain sections using the brain region volume visualizer. Images in the same row represent the plotting of the eigenvectors components, in absolute value, of the FC belonging to the same epoch. Images organized in different columns refer to eigenvectors belonging to different eigenvalues of the matrices. The scale used allows highlighting only the brain area associated to component of the eigenvector greater than the half of the maximum component. Such scale permits to efficiently visualize the relative difference between eigenvectors. According to our definition (see Materials and Methods), the areas with warm colors are the hub regions of the brain network defined by the FC matrices calculated over the relative epoch; the importance of each hub region is proportional to the corresponding eigenvalue. E, Experimental resting state networks and the corresponding functional hubs detected in Mechling et al. (2014). The Ipython scripts to obtain the results in the figure are Figure 1-1 and Figure 1-2 (Extended data).
Figure 1-1
The script shows how to build a tracer-based connectome of mouse brain, using TVMB. Download Figure 1-1, ZIP file.
Figure 1-2
The script shows how to simulate and analyze resting state dynamics in a tracer-based mouse connectome, as described in the paper. Download Figure 1-2, ZIP file.
dMRI-based connectome
TVMB can also make use of user-based dMRI data, enabling the virtualization of individual mouse brains. to use the analysis tools and the visualizer above, the brain volume should be uploaded in nifti format with the same parcellation as the connectome.
As an example, we have used here the high-resolution open-source mouse connectome of Calabrese et al. (2015; Fig. 2A), which we have embedded in the Allen volume (Fig. 2D). In the general case, the user needs to upload the following files: (1) a weight matrix, i.e., a square matrix whose rows and columns label the areas in the parcellation and whose entry (i, j) represents the values of the connection strength between region i and region j; (2) a file containing the labels of the brain regions; and (3) the list of Cartesian triplets that specify the spatial location of each region (Sanz-Leon et al., 2015). As exhaustively explained in the TVB documentation (http://docs.thevirtualbrain.org/index.html), it is possible to provide additional information as the lengths of the tracts connecting the brain areas, or a file containing a vector providing a way of distinguishing cortical from subcortical regions, or the volumes where the connectome is embedded in nifti format, etc.
A, Connectivity matrix obtained from Calabrese et al. (2015). B, Simulated resting state BOLD time series using the connectome shown in A, and the eMFM to model the dynamics of each brain area. C, FCD matrix obtained from the time series. The black segment identifies the epoch of stability of the FCD identified with the spectral embedding technique. D, Functional hubs detected in silico mapped in mouse brain sections using the brain region volume visualizer, as in Figure 1C.
Simulated brain activity
Once a virtual brain is constructed, the TVB environment generates a large-scale brain network equation (Jirsa, 2009) offering multiple ways to produce electrophysiological and neuroimaging signals, and analyze their dynamics.
We present three illustrative examples on how mouse brain network simulations can be accomplished. The scripts and the data necessary to reproduce all the simulations and results presented here are in Extended data.
Resting state activity in the “healthy” brain
Since recent studies highlighted the importance of studying FCD (Allen et al., 2014, Hansen et al., 2015) and the functional hubs of rodents brain (Mechling et al., 2014; Liska et al., 2015) during resting state activity, we introduce an analyzer able to calculate the FCD and to extract the functional hubs (details of the algorithms in Materials and Methods).
We focus on the nonstationary nature of the fMRI FC in resting state observed both in humans (Chang and Glover, 2010; Allen et al., 2014) and in rodents (Keilholz et al., 2013, Liang et al., 2015). Using the simulator tool of TVB we simulate the resting state activity using the reduced Wong Wang model (Wong and Wang, 2006) in the dynamical regime studied by Hansen et al. (2015). The model differs from previous resting state models (Deco and Jirsa, 2012, Deco et al., 2013b) by having a richer dynamical repertoire for each brain region, which results in a greater number of attractors for the global system.
The BOLD signals and the corresponding FCD matrix are shown in Figure 1B,C, respectively. The blocks around the diagonal of the FCD matrix correspond to time intervals during which the FC(t)s are strongly correlated; following the work of Hansen et al. (2015), we call these periods epochs of stability. The FCD analyzer, using the spectral embedding algorithm, detects three epochs of stability (Fig. 1C, black lines) in the FCD matrix. As explained in Materials and Methods, it is possible to identify the central nodes of the i-th network (i = 1,2,3), i.e., the functional hub regions of the i-th epoch, as the nodes linked to the largest components associated with the largest eigenvalues of the FC matrix computed over the i-th epoch. The functional hubs identified using this argument by the FCD analyzer are plotted in the mouse brain sections in Figure 1D.
It is possible to notice several analogies between the simulated functional hubs and the ones previously reported in literature. In particular, the hypothalamus, the visual and somatosensory cortex have been identified as hubs when analyzing resting state networks in mice (Fig. 1E; Mechling et al., 2014). In addition, we find that the agranular insular area is associated to the largest component of the first eigenvector of all epochs (i.e., the most salient hub) in keeping with the experimental data of Liska et al. (2015). Finally, it is interesting to note that the agrunar insular area, the cingulate and temporal cortex are also hub regions in humans (van den Heuvel and Sporns, 2013). Since these hubs are altered in neurologic disorders (Buckner et al., 2009; Crossley et al., 2014), it is straightforward to predict functional consequences after modifying hubs in silico. This illustrates a potential use of TVMB.
Brain activity can be simulated also in individual virtual mouse brain built from fMRI diffusion data as explained before. As an example, we uploaded the detailed dMRI connectome from Calabrese et al. (2015) and simulated subsequent resting state activity of its BOLD signals (Fig. 2).
Resting state activity in epilepsy
TVMB can be used to assess the functional consequence of the anatomic reorganization that takes place in most, if not all, neurologic disorders. Temporal Lobe Epilepsy is a prototypical example of neurologic disorder with well-described anatomic alterations (Esclapez et al., 1999, Chen and Buckmaster, 2005) and functional reorganizations (Centeno and Carmichael, 2014).
Using the tracer-based connectome described above, we removed the connections from the hippocampal CA3 and CA1 regions known to be lost in some forms of medial temporal lobe epilepsy. The simulated BOLD and the corresponding FCD are shown in Figure 3.
A, B, BOLD signals and the corresponding FCD matrix, respectively, obtained by simulating the mouse brain in which some links are removed to mimic epilepsy conditions. C, Functional hubs detected in the epileptic mouse brain after removing links as seen in some forms of epilepsy. The hubs displayed here are extracted from the FC matrix calculated over all the simulated BOLD signals (20 min), i.e., the global FC, since the FCD simulated in the epileptic mouse brain does not present evident sign of nonstationarity and consequently the epoch of stability cannot be detected. The Ipython script to obtain the results in the figure is in Figure 3-1 (Extended data).
Figure 3-1
The script shows how to simulate and analyze resting state dynamics in an epileptic mouse brain, as described in the paper. Download Figure 3-1, ZIP file.
The comparison between the activity of the healthy and epileptic brain at the level of a single region (Figs. 1B, 3A, respectively), does not provide any particular insight. However the differences in brain activity between the two conditions are revealed at the network level when computing the FCD (Fig. 3B). The functional connections that emerge in the epileptic brain are not correlated in time resulting in a suppression of the switching behavior of the FCD, as compared with the control connectome (Fig. 1C). As a result the functional hubs are modified. Since there is no switching, only hubs of global FC can be identified (Fig. 3C).
Seizure propagation
TVB also contains numerous models to generate EEG-like activity, including the Epileptor to simulate seizure genesis and propagation (Jirsa et al., 2014, Proix et al., 2014).
As an experimental reference, we used the electrophysiological recordings performed by Toyoda et al. (2013) in a rat model of temporal lobe epilepsy. Based on the latter results, we used the left hippocampal regions as epileptogenic zones, and analyzed how and where seizures propagated in silico.
The results of the simulation are shown in Figure 4. Each region is characterized by a specific time of seizure onset. The chronological order of the different areas recruited during seizure propagation is shown in Figure 4B,C. The brain areas in abscissa in Figure 4B⇓are sorted according to the seizure onset latency rank found by Toyoda et al. (2013) in rats. Despite the difference in species (rat versus mouse), there is a remarkable analogy with experimental results, suggesting that the structural connectome (and the time delays it imposes) plays a key role in the spatiotemporal pattern of seizure propagation as already reported in humans (Jirsa et al., 2017, Proix et al., 2017).
Simulating epileptiform activity in the mouse brain. A, The time series show simulated seizure genesis and propagation (direct current recording) in silico. B, The graph shows the propagation pattern. Time 0 corresponds to seizure onset in the left hippocampus. On the x-axis, regions are ordered as they are progressively recruited in Toyoda et al. (2013). The y-axis shows the average time of recruitment in arbitrary units of these regions after triggering a seizure in the left hippocampus in silico. Note the good match between simulated and experimental data. Extensive names of the region composing each group are illustrated in the table in Table 1. C, The time distance from seizure onset in the left hippocampus is given by the color scale and plotted in the brain volume for each region. The Ipython script to obtain the results in the figure is Figure 4-1 (Extended data).
Figure 4-1
The script shows how to simulate epileptic seizure propagation in mouse brain, as described in the paper. Download Figure 4-1, ZIP file.
List of regions used for the simulation.
Interpreting and planning experiments with TVMB
Interpreting experimental data with TVMB
Physiologic (e.g., normal aging) and pathologic processes (e.g., neurologic disorders) are associated with both structural (connectome) and functional (resting state networks) alterations. A central issue in neuroscience research is to understand how much structural alterations can account for functional ones. At present, both observations remain at the correlation level. In the case of aging, DTI and rsfMRI can be obtained at different times in a given animal (Figure 5). A virtual brain can be constructed at each time step, to simulate whole-brain dynamics. Following data fitting, alterations found specifically at time t + 1 experimentally can be introduced in the connectome measured at time t. If the resulting in silico rsfMRI reproduces that experimentally measured at time t + 1, it is possible to propose that these structural alterations are sufficient to explain the changes in whole-brain dynamics.
The cartoon illustrates how it is possible to use TVMB to do predictions when studying aging. A mouse can be scanned at different times t extracting anatomic and functional brain information. The anatomic information can be processed to obtain a connectome that can be used in TVMB to create a virtual mouse at each time step. The functional experimental information can be compared with the predictions done in TVMB, investigating how, for example, anatomic modifications during aging affect whole-brain dynamics. Multiple other testable predictions can be done. For example, explore in silico which types of neurones can be stimulated (or silenced) to activate specific resting state networks. The predictions can then be tested in ad hoc transgenic mice with optogenetics.
Planning experiments with TVMB
We present two of the many possibilities offered by the platform. Brain surgery and stimulation are two common procedures used to treat patients, e.g., for epilepsy and Parkinson’s disease. After virtualizing a mouse model of these pathologies at a specific stage of their evolution, researchers can perform neurosurgery in silico and predict the efficacy of the procedure. Likewise, in silico stimulation of brain regions is straightforward in TVMB, which allows studying how resting state dynamics can be manipulated (Spiegler et al., 2016). The predictions thus generated can then be tested experimentally in vivo in the same mouse that was used to make them. Novel preclinical strategies may thus be tested in mice, before their possible clinical transfer.
Many brain functions require dynamical interactions and information transfer between numerous brain regions. The contribution of a given region is thus difficult to evaluate a priori. Using a parameteric study in TVMB, it is possible to predict which regions play a key role by successively activating and inactivating them. Then, one can plan the experiment, choosing the appropriate transgenic mouse line to control the identified region with optogenetics or pharmacogenetics. Such a priori knowledge provided by the in silico approach would considerably accelerate research.
Discussion
TVMB opens a new set of research possibilities: it allows researchers, from different fields, to easily build specific/individual mouse brains (using various resolutions, weighting definitions and parcellations), to simulate different dynamical behaviors (using diverse neural population models, numerical integration schemes, and simulated neuroimaging modalities) and finally to analyze the results.
However, while TVMB is a highly generic framework, its underlying mathematical framework and simulation techniques make standard assumptions, among which the two most essential are that (1) the average activity of large populations of neurons is a meaningful quantification of the phenomena to be modeled and (2) the statistics of white matter fibers sufficiently describes how regions interact. Both resting state dynamics and seizure propagation, as demonstrated above, satisfy these assumptions. On the other hand, for example, fine grained spike timing effects would not well described within TVMB’s mathematical framework.
The wide range of possibility offered by rodent experiments will easily accommodate the validation of the parameterization required by all the modeling approaches contained in the software. This validation sometimes can proceed at a qualitative level for phenomenological models, such as the Kuramoto model of synchronization, but many detailed biophysical models allow for quantitive comparison with empirical data, such as spike timing (Brette and Gerstner, 2005). At the whole-brain level, TVMB allows for direct comparison with common modalities such as EEG, MEG, and fMRI, or common statistics thereupon such as FC; these comparisons allow for the characterization of parameter values in terms of their fit with empirical data and thus biological validity. The experimentally observed functional characteristics of the mouse brain (e.g., a functional hub during resting state or the effects of specific connections removal) can be easily imposed in the output of the virtual system, and through data fitting algorithms, it will be possible to retrieve the parameters of the model that give rise to that particular functional behavior. In this way, closing the circle, the reliability of the new predictions accomplished with the fitted parameter set will be improved; additionally the knowledge of the key features responsible of the different functional behavior allows to control and manipulate the system in silico, and, going a step further, also in vivo.
TVMB thus offers not only a conceptual framework to interpret neuroimaging data but, combined with experimental approaches, it also offers an operative framework to investigate the causal links between structure and function in the brain.
It is important to note that, to build the Allen connectome, we took the mirror image of the right hemisphere to build the left one, since injections of tracers were performed in the right hemisphere only (Oh et al., 2014). However, contrary to humans (Toga and Thompson, 2003), there is a high degree of similarity in terms of connections between the right and the left hemispheres in rodents (Calabrese et al., 2015).
The detailed tracer connectome, that can be built with TVMB, can be used to build a precise brain model, and serve as a reference or template for dMRI-based modeling, since dMRI data suffer from three major limitations, that are: (1) our ignorance about the directionality of the connections, (2) the indirect nature of the measures of connectivity based on water diffusion in white matter, and (3) the fact that complex fiber pathways (as crossing fibers, sharp change in directionality and long brain wirings) cannot be properly detected. How such limitations affect the simulations is not yet established. On the other hand tracer-based connectome are built by averaging experiments over many mice, and the definition of the seed region is based on stereotaxic coordinate rather than on anatomy. A comparison of the effects in simulating brain network using one or the other kind of connectome is needed especially since the only data available in human are dMRI data.
TVMB is an actively developed software, with new versions released regularly with new features. Among those targeted specifically for the mouse, the module which builds connectivities from the public Allen data will continue to evolve as the available dataset becomes richer. For example, when cortical layer annotations become available, it will be possible to construct mouse connectivities in which the cortical layers are distinct, allowing for example manipulations of inter-layer interactions. As dMRI protocols and tractography techniques become more established for rodent datasets, TVMB has potential to include modules which automate, with visual inspection across each step, the generation of connectomes for individual rodent data, directly from the DICOM slices provided by the acquisition equipment.
Acknowledgments
Acknowledgements: This research was supported by ANR MOTION ANR-13-NEUC-0005-01, and the European Union's Horizon 2020 research and innovation programme under grant agreement No. 720270. We thank E. Calabrese, A. Badea, G. Cofer, Y. Cofer, G. A. Johnson, and the Allen Institute for sharing their data with the community.
Footnotes
The authors declare no competing financial interests.
This work was supported by Inserm.
↵* V.K.J. and C.B. are equally contributing last authors.
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.
References
In this issue
