Skip to content
Snippets Groups Projects
Commit 7233d1a3 authored by Jan Fousek's avatar Jan Fousek
Browse files

rc1 everything should work

parent f21b3f53
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
1_TVB_First_steps.ipynb
+ 73
− 70
View file @ 7233d1a3
%% Cell type:code id:ec97e74a tags: %% Cell type:code id:68ccd22b tags:
   
``` python ``` python
%load_ext autoreload %load_ext autoreload
%autoreload 2 %autoreload 2
``` ```
   
%% Cell type:markdown id:e8fd8f47 tags: %% Cell type:markdown id:e5a89d1a tags:
   
# First steps using TVB # First steps using TVB
   
**Summer School in Nonlinear Dynamics for the Life Sciences with Applications to Neuroscience and Psychology May 31-June 11, 2021** **Summer School in Nonlinear Dynamics for the Life Sciences with Applications to Neuroscience and Psychology May 31-June 11, 2021**
   
%% Cell type:markdown id:3f6bc59a tags: %% Cell type:markdown id:c669ec16 tags:
   
# Brain network modeling with TVB # Brain network modeling with TVB
   
*** ***
   
Components to be covered Components to be covered
* node dynamics * node dynamics
* connectivity * connectivity
* coupling functions * coupling functions
* monitors * monitors
* integrator * integrator
* stimulus * stimulus
   
%% Cell type:markdown id:dfe4517d tags: %% Cell type:markdown id:727cf9bf tags:
   
# Objectives # Objectives
   
Here we: Here we:
   
* Build a brain network model using subject-specific structural connectivity, * Build a brain network model using subject-specific structural connectivity,
* Simulate resting-state activity, * Simulate resting-state activity,
* Characterize the resting-state activity by calculating the functional connectivity (FC), * Characterize the resting-state activity by calculating the functional connectivity (FC).
* Extract the resting-state networks (RSNs).
   
%% Cell type:markdown id:329f9428 tags: %% Cell type:markdown id:48605a2c tags:
   
# How to do it with TVB? # How to do it with TVB?
   
*** ***
   
In the first part of this tutorial, we presents the basic anatomy of a region simulation using The Virtual Brain's scripting interface. In the first part of this tutorial, we presents the basic anatomy of a region simulation using The Virtual Brain's scripting interface.
   
The first thing we want to do is to import the modules we will need for a simulation. The first thing we want to do is to import the modules we will need for a simulation.
   
%% Cell type:code id:bac045da tags: %% Cell type:code id:48088d0e tags:
   
``` python ``` python
%%capture %%capture
import os import os
import time as tm import time as tm
   
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
   
from tvb.simulator.lab import * from tvb.simulator.lab import *
from utils import plot_connectivity from utils import plot_connectivity
from phase_plane import phase_plane_interactive from phase_plane import phase_plane_interactive
``` ```
   
%% Cell type:markdown id:8dd86a7b tags: %% Cell type:markdown id:1a410d7a tags:
   
A basic simulation of TVB consists of **5 main components**. Each of these components is an object within TVB: A basic simulation of TVB consists of **5 main components**. Each of these components is an object within TVB:
   
%% Cell type:markdown id:ec330e11 tags: %% Cell type:markdown id:ae302f73 tags:
   
### Connectivity ### Connectivity
   
   
We start by loading and visualizing the structural connectivity matrix that represents the set of all existing connections between brain areas. Having loaded the default dataset, we can then alter the speed of signal propagation through the network: We start by loading and visualizing the structural connectivity matrix that represents the set of all existing connections between brain areas. Having loaded the default dataset, we can then alter the speed of signal propagation through the network:
   
%% Cell type:code id:6c563ce4 tags: %% Cell type:code id:710d21e4 tags:
   
``` python ``` python
# Import the anatomical structural connectivity. # Import the anatomical structural connectivity.
conn = connectivity.Connectivity().from_file( conn = connectivity.Connectivity().from_file(
os.path.abspath('dataset/connectivity_76.zip') os.path.abspath('dataset/connectivity_76.zip')
) )
nregions = len(conn.region_labels) # Number of regions nregions = len(conn.region_labels) # Number of regions
conn.speed = np.array(np.inf) # Set the conduction speed conn.speed = np.array(np.inf) # Set the conduction speed
conn.configure() conn.configure()
``` ```
   
%% Cell type:markdown id:229ebbef tags: %% Cell type:markdown id:de6dafab tags:
   
*Take a look at some of the properties of the `conn` object: `weights`, `delays`, `region_labels`, etc.* *Take a look at some of the properties of the `conn` object: `weights`, `delays`, `region_labels`, etc.*
   
%% Cell type:code id:8e7e9be7 tags: %% Cell type:code id:a301f70b tags:
   
``` python ``` python
conn conn
``` ```
   
%% Output %% Output
   
<tvb.datatypes.connectivity.Connectivity at 0x7fb69414a6d8> <tvb.datatypes.connectivity.Connectivity at 0x7fb69414a6d8>
   
%% Cell type:code id:04438045 tags: %% Cell type:code id:c82cd475 tags:
   
``` python ``` python
plot_connectivity(conn) plot_connectivity(conn)
``` ```
   
%% Output %% Output
   
   
%% Cell type:markdown id:ae6af06e tags: %% Cell type:markdown id:e4d3660a tags:
   
### Model ### Model
   
%% Cell type:markdown id:96407494 tags: %% Cell type:markdown id:0adf2b72 tags:
   
A set of differential equations describing the local neural dynamics. There are a number of predefined models available in TVB, as an example here we will use the **Generic 2-dimensional Oscillator** model: A set of differential equations describing the local neural dynamics. There are a number of predefined models available in TVB, as an example here we will use the **Generic 2-dimensional Oscillator** model:
   
\begin{align} \begin{align}
\dot{V} &= d \, \tau (-f V^3 + e V^2 + g V + \alpha W + \gamma I), \\ \dot{V} &= d \, \tau (-f V^3 + e V^2 + g V + \alpha W + \gamma I), \\
\dot{W} &= \dfrac{d}{\tau}\,\,(c V^2 + b V - \beta W + a). \dot{W} &= \dfrac{d}{\tau}\,\,(c V^2 + b V - \beta W + a).
\end{align} \end{align}
   
%% Cell type:markdown id:431f293c tags: %% Cell type:markdown id:e3f828d7 tags:
   
We can explore it using the **phase-plane** tool in TVB. We can explore it using the **phase-plane** tool in TVB.
   
%% Cell type:code id:4ccbd26b tags: %% Cell type:code id:a7ae0fec tags:
   
``` python ``` python
# Create and launch the phase-plane tool # Create and launch the phase-plane tool
g2d = models.Generic2dOscillator() # The model to investigate g2d = models.Generic2dOscillator() # The model to investigate
heunint = integrators.HeunDeterministic(dt=1.0) # Integrator to use. heunint = integrators.HeunDeterministic(dt=1.0) # Integrator to use.
# dt is the integration step # dt is the integration step
   
phase_plane_interactive(model=g2d, integrator=heunint) phase_plane_interactive(model=g2d, integrator=heunint)
   
   
# Try changing the default parameters to # Try changing the default parameters to
# 1) a = 3 (limit cycle) # 1) a = 3 (limit cycle)
# 2) b = 1 (two stable and unstable fixed points) # 2) b = 1 (two stable and unstable fixed points)
``` ```
   
%% Output %% Output
   
   
<function phase_plane.phase_plane_interactive.<locals>.plot_phase_plane(**param_kwargs)> <function phase_plane.phase_plane_interactive.<locals>.plot_phase_plane(**param_kwargs)>
   
%% Cell type:markdown id:f806ac21 tags: %% Cell type:markdown id:9e915bee tags:
   
For our simulation we set it to the limit-cycle regime. For our simulation we set it to the limit-cycle regime.
   
%% Cell type:markdown id:2174002c tags: %% Cell type:markdown id:939696ba tags:
   
*Note that the parameters of the `Model` class has to be set as numpy arrays.* *Note that the parameters of the `Model` class has to be set as numpy arrays.*
   
%% Cell type:code id:2950caad tags: %% Cell type:code id:1b8ffc60 tags:
   
``` python ``` python
# Initialize the model # Initialize the model
g2d = models.Generic2dOscillator(a=np.array(1.7402)) g2d = models.Generic2dOscillator(a=np.array(1.7402))
g2d g2d
``` ```
   
%% Output %% Output
   
<tvb.simulator.models.oscillator.Generic2dOscillator at 0x7fb64c430b38> <tvb.simulator.models.oscillator.Generic2dOscillator at 0x7fb64c430b38>
   
%% Cell type:markdown id:c4d3b0d1 tags: %% Cell type:markdown id:5045a18b tags:
   
### Coupling function ### Coupling function
   
It is a function that is used to join the local `Model` dynamics at distinct spatial locations over the connections described in `Connectivity`. Proper setting of the parameters for this function requires some knowledge of the properties of both the model being used and the structure through which it is connected. For our present purposes, we happen to know that for this configuration of parameters of TVB's `Generic2dOscillator` connected through TVB's default connectivity matrix, a linear function with a slope of 0.0075 is a reasonable thing to use. It is a function that is used to join the local `Model` dynamics at distinct spatial locations over the connections described in `Connectivity`. Proper setting of the parameters for this function requires some knowledge of the properties of both the model being used and the structure through which it is connected. For our present purposes, we happen to know that for this configuration of parameters of TVB's `Generic2dOscillator` connected through TVB's default connectivity matrix, a linear function with a slope of 0.0075 is a reasonable thing to use.
   
%% Cell type:code id:966e422c tags: %% Cell type:code id:27292c27 tags:
   
``` python ``` python
# Initialise a Coupling function. # Initialise a Coupling function.
G = np.array(0.0075) G = np.array(0.0075)
con_coupling = coupling.Scaling(a=G) con_coupling = coupling.Scaling(a=G)
``` ```
   
%% Cell type:markdown id:8a0e4549 tags: %% Cell type:markdown id:d491e8c4 tags:
   
### Integrator ### Integrator
   
Now that we have defined our structure and dynamics, we need to select an integration scheme. While TVB supports a number of schemes, for most purposes you should use either `HeunDeterministic` or `HeunStochastic`. Now that we have defined our structure and dynamics, we need to select an integration scheme. While TVB supports a number of schemes, for most purposes you should use either `HeunDeterministic` or `HeunStochastic`.
   
Note that the most important thing here is to use a step size that is small enough for the integration to be numerically stable. Note that the most important thing here is to use a step size that is small enough for the integration to be numerically stable.
   
%% Cell type:code id:57e730bc tags: %% Cell type:code id:8b205225 tags:
   
``` python ``` python
# Initialise an Integrator scheme. # Initialise an Integrator scheme.
dt = 0.1 # Integration step [ms] dt = 0.1 # Integration step [ms]
   
# We can use the deteministic integrator: # We can use the deteministic integrator:
# heunint = integrators.HeunDeterministic(dt=dt) # heunint = integrators.HeunDeterministic(dt=dt)
   
# Or stochastic integrator: # Or stochastic integrator:
nsigma = 0.001 # Standard deviation of the noise nsigma = 0.001 # Standard deviation of the noise
hiss = noise.Additive(nsig=np.array([nsigma, 0])) hiss = noise.Additive(nsig=np.array([nsigma, 0]))
heunint = integrators.HeunStochastic(dt=dt, noise=hiss) heunint = integrators.HeunStochastic(dt=dt, noise=hiss)
``` ```
   
%% Cell type:markdown id:83b1f05a tags: %% Cell type:markdown id:23c2a35c tags:
   
### Monitors ### Monitors
   
The last component is `Monitor`: observer models. Although there are Monitors which apply a biophysical measurement process to the simulated neural activity, such as EEG, MEG, etc, here we will select two simple monitors to get the time series of the model state variables: The last component is `Monitor`: observer models. Although there are Monitors which apply a biophysical measurement process to the simulated neural activity, such as EEG, MEG, etc, here we will select two simple monitors to get the time series of the model state variables:
   
* the `Raw` monitor takes no arguments and simply returns all the simulated data at the time resolution of the integration step, * the `Raw` monitor takes no arguments and simply returns all the simulated data at the time resolution of the integration step,
* the `TemporalAverage` monitor averages over a time window of length period returning one time point every period (given in ms). * the `TemporalAverage` monitor averages over a time window of length period returning one time point every period (given in ms).
   
*Note: as a general rule, the `Raw` monitor shouldn't be used for anything but very short simulations as the amount of data returned can become prohibitively large* *Note: as a general rule, the `Raw` monitor shouldn't be used for anything but very short simulations as the amount of data returned can become prohibitively large*
   
%% Cell type:code id:a2512904 tags: %% Cell type:code id:66d578cd tags:
   
``` python ``` python
# Initialise some Monitors with period in physical time. # Initialise some Monitors with period in physical time.
mon_raw = monitors.Raw() mon_raw = monitors.Raw()
mon_tavg = monitors.TemporalAverage(period=1) # 1000 Hz mon_tavg = monitors.TemporalAverage(period=1) # 1000 Hz
   
#Bundle them #Bundle them
what_to_watch = (mon_raw, mon_tavg) what_to_watch = (mon_raw, mon_tavg)
``` ```
   
%% Cell type:markdown id:f15d3c85 tags: %% Cell type:markdown id:e3f3da9b tags:
   
*** ***
   
%% Cell type:markdown id:818487e4 tags: %% Cell type:markdown id:1a20a9a7 tags:
   
### Go! Simulate ### Go! Simulate
   
The last step is to bring all these components together into a `Simulator` object. We then need to run the configure method, which basically just acts to calculate information necessary for the simulation that draws on specific combinations of the components. The last step is to bring all these components together into a `Simulator` object. We then need to run the configure method, which basically just acts to calculate information necessary for the simulation that draws on specific combinations of the components.
   
%% Cell type:code id:0a98d0ab tags: %% Cell type:code id:95deeaa2 tags:
   
``` python ``` python
# Initialise the Simulator. # Initialise the Simulator.
sim = simulator.Simulator(model=g2d, sim = simulator.Simulator(model=g2d,
connectivity=conn, connectivity=conn,
conduction_speed=np.float(conn.speed), conduction_speed=np.float(conn.speed),
coupling=con_coupling, coupling=con_coupling,
integrator=heunint, integrator=heunint,
monitors=what_to_watch) monitors=what_to_watch)
sim.configure() sim.configure()
``` ```
   
%% Output %% Output
   
<tvb.simulator.simulator.Simulator at 0x7fb64c70bc18> <tvb.simulator.simulator.Simulator at 0x7fb64c70bc18>
   
%% Cell type:markdown id:7f418030 tags: %% Cell type:markdown id:4b408895 tags:
   
Now, we can run the simulation. All we need to do is iterate for some length, which we provide in *ms*, and collect the output. Now, we can run the simulation. All we need to do is iterate for some length, which we provide in *ms*, and collect the output.
   
The data returned by the simulator is in the form of a list of arrays. For most subsequent purposes it is much easier to deal with the data if it exists as a single contiguous array, hence the enumeration of the return values for the `sim.run`. The data returned by the simulator is in the form of a list of arrays. For most subsequent purposes it is much easier to deal with the data if it exists as a single contiguous array, hence the enumeration of the return values for the `sim.run`.
   
%% Cell type:code id:a02b73e4 tags: %% Cell type:code id:07822ba4 tags:
   
``` python ``` python
# Perform the simulation. # Perform the simulation.
tic = tm.time() tic = tm.time()
   
# this should take ~25s # this should take ~25s
(raw_time, raw_data), (tavg_time, tavg_data) = sim.run(simulation_length=10000.) (raw_time, raw_data), (tavg_time, tavg_data) = sim.run(simulation_length=10000.) # ~30s
   
'simulation required %0.3f seconds.' % (tm.time()-tic) 'simulation required %0.3f seconds.' % (tm.time()-tic)
``` ```
   
%% Output %% Output
   
'simulation required 33.790 seconds.' 'simulation required 33.790 seconds.'
   
%% Cell type:markdown id:b4c3d1d9 tags: %% Cell type:markdown id:d6e42579 tags:
   
For each monitor, the simulation returns the time points and the associated data array. For each monitor, the simulation returns the time points and the associated data array.
   
The data arrays have shapes $n_\text{timepoints} \times n_\text{variables of interest} \times n_\text{regions} \times n_\text{modes}$. The data arrays have shapes $n_\text{timepoints} \times n_\text{variables of interest} \times n_\text{regions} \times n_\text{modes}$.
   
%% Cell type:code id:e2aac7e0 tags: %% Cell type:code id:38e019da tags:
   
``` python ``` python
f"raw_time: {raw_time.shape} | raw_data: {raw_data.shape}" f"raw_time: {raw_time.shape} | raw_data: {raw_data.shape}"
``` ```
   
%% Output %% Output
   
'raw_time: (100000,) | raw_data: (100000, 1, 76, 1)' 'raw_time: (100000,) | raw_data: (100000, 1, 76, 1)'
   
%% Cell type:code id:1bfa54d6 tags: %% Cell type:markdown id:77422a89 tags:
***
%% Cell type:markdown id:a03c41a9 tags:
### Visualize our simulation
%% Cell type:markdown id:2eaa282c tags:
And finally, we can look at the results of our simulation in terms of time series.
%% Cell type:code id:73823cce tags:
   
``` python ``` python
# Remove the dimensions with one element for easier indexing # Remove the dimensions with one element for easier indexing
raw = np.squeeze(np.array(raw_data)) raw = np.squeeze(np.array(raw_data))
tavg = np.squeeze(np.array(tavg_data)) tavg = np.squeeze(np.array(tavg_data))
raw.shape, tavg.shape raw.shape, tavg.shape
``` ```
   
%% Output %% Output
   
((100000, 76), (10000, 76)) ((100000, 76), (10000, 76))
   
%% Cell type:markdown id:1c91b9ce tags: %% Cell type:code id:1cd9ceee tags:
***
%% Cell type:markdown id:8776e1c9 tags:
### Visualize our simulation
%% Cell type:markdown id:9e2c73af tags:
And finally, we can look at the results of our simulation in terms of time series.
%% Cell type:code id:ac65f763 tags:
   
``` python ``` python
# Normalize the time series for easier visualization # Normalize the time series for easier visualization
nraw = raw / (np.max(raw, 0) - np.min(raw, 0)) nraw = raw / (np.max(raw, 0) - np.min(raw, 0))
ntavg = tavg / (np.max(tavg, 0) - np.min(tavg, 0)) ntavg = tavg / (np.max(tavg, 0) - np.min(tavg, 0))
``` ```
   
%% Cell type:code id:ea95a7a0 tags: %% Cell type:code id:a04167a4 tags:
   
``` python ``` python
# Plot the raw time series # Plot the raw time series
fig1 = plt.figure(figsize=(10,8)) fig1 = plt.figure(figsize=(10,8))
plt.plot(raw_time[:], nraw[:, :10] + np.r_[:10]) plt.plot(raw_time[:], nraw[:, :10] + np.r_[:10])
plt.title('Raw Neuronal Activity', fontsize=20) plt.title('Raw Neuronal Activity', fontsize=20)
plt.xlabel('Time [ms]', fontsize=20) plt.xlabel('Time [ms]', fontsize=20)
plt.yticks(range(10), conn.region_labels[:10], fontsize=10) plt.yticks(range(10), conn.region_labels[:10], fontsize=10)
   
# Plot the temporally averaged time series # Plot the temporally averaged time series
fig2 = plt.figure(figsize=(10,8)) fig2 = plt.figure(figsize=(10,8))
plt.plot(tavg_time[:], ntavg[:, :10] + np.r_[:10]) plt.plot(tavg_time[:], ntavg[:, :10] + np.r_[:10])
plt.title('Temporally Averaged Neuronal Activity', fontsize=20) plt.title('Temporally Averaged Neuronal Activity', fontsize=20)
plt.xlabel('Time [ms]', fontsize=20) plt.xlabel('Time [ms]', fontsize=20)
plt.yticks(range(10), conn.region_labels[:10], fontsize=10) plt.yticks(range(10), conn.region_labels[:10], fontsize=10)
   
plt.show() plt.show()
``` ```
   
%% Output %% Output
   
   
   
%% Cell type:markdown id:e8da10ba tags: %% Cell type:markdown id:53dde102 tags:
   
*** ***
   
%% Cell type:markdown id:36973616 tags: %% Cell type:markdown id:ed77a3e0 tags:
   
# Analysis example: Functional Connectivity (FC) # Analysis example: Functional Connectivity (FC)
   
%% Cell type:markdown id:80d70041 tags: %% Cell type:markdown id:6c2533e4 tags:
   
**Functional Connectivity (FC)** describes the connectedness of two brain regions by means of the covariance between their time series. The classic and most widely used method to infer the strength of network interactions or functional connections consists in estimating the linear (Pearson) correlation coefficient between temporal signals. If two regions activate and deactivate at the same time, there is likely a functional connection. **Functional Connectivity (FC)** describes the connectedness of two brain regions by means of the covariance between their time series. The classic and most widely used method to infer the strength of network interactions or functional connections consists in estimating the linear (Pearson) correlation coefficient between temporal signals. If two regions activate and deactivate at the same time, there is likely a functional connection.
   
To calculate the FC, we ignore the first second of the simulation which might be affected by the initial conditions. To calculate the FC, we ignore the first second of the simulation which might be affected by the initial conditions.
   
%% Cell type:code id:24706c44 tags: %% Cell type:code id:3e7ecaf9 tags:
   
``` python ``` python
fc = np.corrcoef(tavg[1000:].T) fc = np.corrcoef(tavg[1000:].T)
``` ```
   
%% Cell type:markdown id:2412634a tags: %% Cell type:markdown id:76165909 tags:
   
And now we display the resulting FC matrix: And now we display the resulting FC matrix:
   
%% Cell type:code id:1c4d17b3 tags: %% Cell type:code id:6220213a tags:
   
``` python ``` python
# Visualize FC matrix # Visualize FC matrix
plt.figure(figsize=(8,8)) plt.figure(figsize=(8,8))
plt.imshow(fc, interpolation='nearest', cmap='jet') plt.imshow(fc, interpolation='nearest', cmap='jet')
plt.title('SIM FC', fontsize=20) plt.title('SIM FC', fontsize=20)
plt.xlabel('ROIs', fontsize=20); plt.ylabel('ROIs', fontsize=20) plt.xlabel('ROIs', fontsize=20); plt.ylabel('ROIs', fontsize=20)
plt.xticks([], fontsize=20); plt.yticks([], fontsize=20) plt.xticks([], fontsize=20); plt.yticks([], fontsize=20)
cb = plt.colorbar(shrink=0.7) cb = plt.colorbar(shrink=0.7)
cb.set_label('PCC', fontsize=20) cb.set_label('PCC', fontsize=20)
   
plt.show() plt.show()
``` ```
   
%% Output %% Output
   
   
%% Cell type:markdown id:b3811ffc tags: %% Cell type:markdown id:2189b99b tags:
   
# What is next? # What is next?
* take a look at the BOLD monitor
* introduce external stimulus
* epilepsy use case
......
2_TVB_BOLD_and_Stimulus.ipynb
+ 639
− 40
View file @ 7233d1a3
source diff could not be displayed: it is too large. Options to address this: view the blob.
3_TVB_Modelling_Epilepsy.ipynb
+ 40
− 40
View file @ 7233d1a3
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment