In [1]:
%pylab nbagg
import sys
from tvb.simulator.lab import *
LOG = get_logger('demo')
import scipy.stats
from sklearn.decomposition import FastICA
import time
import utils

Populating the interactive namespace from numpy and matplotlib

Introduction

Fluctuations in brain activity in non-task conditions are now a well-established phenomena in the literature. These fluctuations are not random but shown to exhibit spatial patterns, referred to as resting state networks. Despite being readily identifiable during rest, these networks are related to specific functions and on the other hand abnormalities in such RSNs have been associated with pathology.

In the following, we will demonstrate some starting points for modeling resting state networks in TVB, using the default data set.

Setting up the simulation

In the following, we'll use a basic region level simulation, with the generic oscillator set in an excitable regime, linear coupling with low strength, a stochastic integrator with low noise and a temporal average monitor at 200 Hz.

These settings are a good starting point for modeling resting state patterns because no particular factor dominates the dynamics and a balance between the structural connectivity, moderate intrinsic excitability and noise comes into play.

In [2]:
def run_sim(conn, cs, D, cv=3.0, dt=0.5, simlen=1e3):
    sim = simulator.Simulator(
        model=models.Generic2dOscillator(a=numpy.array([0.0])),
        connectivity=conn,
        coupling=coupling.Linear(a=numpy.array([cs])),
        integrator=integrators.HeunStochastic(dt=dt, noise=noise.Additive(nsig=array([D]))),
        monitors=(monitors.TemporalAverage(period=5.0),) # 200 Hz
    )
    sim.configure()
    (t, y), = sim.run(simulation_length=simlen)
    return t, y[:, 0, :, 0]

conn = connectivity.Connectivity.from_file()

WARNING  File 'hemispheres' not found in ZIP.

One of the common features of simulations is an initial transient, so we'll perform a one minute simulation, and as soon as the time series have been generated, we check that the transient has decayed:

In [3]:
tic = time.time()
t, y = run_sim(conn, 6e-2, 5e-4, simlen=10*60e3)
'simulation required %0.3f seconds.' % (time.time() - tic, )
Out[3]:
'simulation required 378.781 seconds.'

Functional Connectivity

Next, to quickly assess the presence of a network structure in the time series, we window the time series into 1 second non overlapping windows, obtain per-window correlation matrices

In [4]:
cs = []
for i in range(int(t[-1]/1e3)):
    cs.append(corrcoef(y[(t>(i*1e3))*(t<(1e3*(i+1)))].T))
cs = array(cs)
cs.shape
Out[4]:
(599, 76, 76)

The strength of correlation can be assessed statistically by Fisher Z transforming the coefficients and applying a t-test,

In [5]:
cs_z = arctanh(cs)
for i in range(cs.shape[1]):
    cs_z[:, i, i] = 0.0
_, p = scipy.stats.ttest_1samp(cs, 0.0)

/Library/Anaconda/anaconda3/envs/tvb-run/lib/python3.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: divide by zero encountered in arctanh
  """Entry point for launching an IPython kernel.

Which we then visualize the structural connectivity (left) and functional connectivity (right) as an adjacency matrices applying a threshold on significance:

In [12]:
figure(figsize=(10, 4))
subplot(121), imshow(conn.weights, cmap='binary', interpolation='none')
subplot(122), imshow(log10(p)*(p < 0.05), cmap='gray', interpolation='none');
show()

/Library/Anaconda/anaconda3/envs/tvb-run/lib/python3.7/site-packages/ipykernel_launcher.py:3: RuntimeWarning: divide by zero encountered in log10
  This is separate from the ipykernel package so we can avoid doing imports until

We can see there are significant deviations in the FC from the SC which are especially prominent in the interhemispheric FC, where interactions are found despite limited interhemispheric SC.

We can then ask what degree of similarity there is between the average functional connectivity and the structural connectivity, and how it varies over time:

In [13]:
figure()
plot(r_[1:len(cs)+1], [corrcoef(cs_i.ravel(), conn.weights.ravel())[0, 1] for cs_i in cs])
ylim([0, 0.5])
ylabel('FC-SC correlation')
xlabel('Time (s)');
show()

Seed-region correlation maps

A common visualization of FC specific to a given is to pull out its row of the FC matrix and plot a colormap on the anatomy. We can do this will the simulated functional connectivity to identify which regions are highly coordinated with the seed region.

In [14]:
ctx = cortex.Cortex.from_file()
ctx.region_mapping_data.connectivity = conn
rm = ctx.region_mapping

def plot_roi_corr_map(reg_name):
    roi = list(conn.ordered_labels).index(reg_name)
    cs_m = cs[2:].mean(axis=0)
    utils.multiview(cs_m[roi][rm], cortex=ctx, shaded=False, suptitle=reg_name, figsize=(10, 5))

As a few examples of such maps, seeding in the left motor cortex, right ventrolateral prefront cortex, and right superior parietal cortex:

In [16]:
for reg in 'lM1 rPFCVL rPCS'.split():
    plot_roi_corr_map(reg)

Seed-region maps are useful when one has a prior about which regions are implicated for a given network.

ICA

Another common exploratory tool in resting state data analysis, where the implicated regions or networks are not known a priori, is independent component analysis, which extracts components with unique or independent statistical properties.

For exapmle, we can perform an ICA keeping 5 components the above simulated data,

In [10]:
ica = FastICA(n_components=5, max_iter=250)
ica.fit(y[t>2e3])
Out[10]:
FastICA(max_iter=250, n_components=5)

And then visualize the different components

In [17]:
for i, comp in enumerate(ica.components_[:3]):
    utils.multiview(comp[rm], cortex=ctx, shaded=False, suptitle='ICA %d' % (i, ), figsize=(10, 5))

These components are not selected 'by hand', but represent independent subnetworks during the simulated resting state activity.

Finally, we point out, that commonly ICA analyses of fMRI are done at a group level to identify spatial patterns which are reproducible across subjects, whereas in this application to this simulation, spatial components may reflect as much non-robust, spontaneous fluctuations of the network passing through various state as the dominant rest state networks identified in human rest state.

Further exploration

These results are starting point, from which you can base you simulations, going in directions such as

  • Perform a parameter sweep to identify regions of improved FC-SC correlation
  • Introduce a subject specific structural connectivity and compare with subject specific functional connectivity
  • Threshold the calculated functional connectivities and apply graph theoretic measures

We hope this has been a useful tutorial and welcome any comments or questions on our mailing list (https://groups.google.com/forum/#!forum/tvb-users).

References

In [ ]: