In [1]:
%matplotlib inline
4 - Dynamic Connectivity (Group Analysis)¶
In this short tutorial, we will build and expand on the previous tutorials by computing the dynamic connectivity, using Time-Varying Functional Connectivity Graphs.
In the near future, the standard method of "sliding window" will be supported.
Load data¶
In [2]:
import numpy as np
In [3]:
raw_eeg_eyes_open = np.load("data/eyes_opened.npy")
raw_eeg_eyes_closed = np.load("data/eyes_closed.npy")
num_trials, num_channels, num_samples = np.shape(raw_eeg_eyes_open)
read_trials = 10
In [4]:
eeg_eyes_open = raw_eeg_eyes_open[0:read_trials, ...]
eeg_eyes_closed = raw_eeg_eyes_closed[0:read_trials, ...]
Dynamic connectivity¶
Prepare and configure the estimator object¶
In [5]:
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
from dyconnmap import tvfcg
from dyconnmap.fc import IPLV
In [6]:
fb = [7.0, 13.0]
cc = 4.0
fs = 160.0
step = 80
In [7]:
estimator = IPLV(fb, fs)
Process condition "eyes open"¶
In [8]:
print('Working .', end='')
X = np.squeeze(eeg_eyes_open[0])
fcgs = tvfcg(X, estimator, fb, fs, cc, step)
fcgs_eyes_open = np.array(np.real(fcgs))
for i in range(1, read_trials):
print('.', end='')
X = np.squeeze(eeg_eyes_open[i])
fcgs = tvfcg(X, estimator, fb, fs, cc, step)
fcgs_eyes_open = np.vstack([fcgs_eyes_open, np.real(fcgs)])
print('')
Working ..........
Process condition "eyes closed"¶
In [9]:
print('Working .', end='')
X = np.squeeze(eeg_eyes_closed[0])
fcgs = tvfcg(X, estimator, fb, fs, cc, step)
fcgs_eyes_closed = np.array(np.real(fcgs))
for i in range(1, read_trials):
print('.', end='')
X = np.squeeze(eeg_eyes_closed[i])
fcgs = tvfcg(X, estimator, fb, fs, cc, step)
fcgs_eyes_closed = np.vstack([fcgs_eyes_closed, np.real(fcgs)])
print('')
Working ..........
FCμstates / Clustering¶
In [10]:
from dyconnmap.cluster import NeuralGas
In [11]:
num_fcgs_eo, _, _ = np.shape(fcgs_eyes_open)
num_fcgs_ec, _, _ = np.shape(fcgs_eyes_closed)
In [12]:
fcgs = np.vstack([fcgs_eyes_open, fcgs_eyes_closed])
num_fcgs, num_channels, num_channels = np.shape(fcgs)
triu_ind = np.triu_indices_from(np.squeeze(fcgs[0, ...]), k=1)
fcgs = fcgs[:, triu_ind[0], triu_ind[1]]
In [13]:
rng = np.random.RandomState(0)
mdl = NeuralGas(n_protos=5, rng=rng).fit(fcgs)
encoding, symbols = mdl.encode(fcgs)
Separate the encoded symbols based on their original groupings¶
In [14]:
grp_dist_eo = symbols[0:num_fcgs_eo]
grp_dist_ec = symbols[num_fcgs_eo:]
Plot¶
In [15]:
h_grp_dist_eo = np.histogram(grp_dist_eo, bins=mdl.n_protos, normed=True)
h_grp_dist_ec = np.histogram(grp_dist_ec, bins=mdl.n_protos, normed=True)
In [16]:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(12, 6))
ind = np.arange(mdl.n_protos)
p1 = ax.bar(ind - 0.125, h_grp_dist_ec[0], 0.25, label='Eyes Closed')
p2 = ax.bar(ind + 0.125, h_grp_dist_eo[0], 0.25, label='Eyes Open')
ax.legend()
ax.set_xlabel('Symbol Index')
ax.set_ylabel('Hits %')
ax.set_xticks(np.arange(mdl.n_protos))
plt.show()
Convert state prototypes to symmetric matrices and plot them
In [17]:
protos_mtx = np.zeros((mdl.n_protos, 64, 64))
for i in range(mdl.n_protos):
symbol_state = np.zeros((64, 64))
symbol_state[triu_ind] = mdl.protos[i, :]
symbol_state = symbol_state + symbol_state.T
np.fill_diagonal(symbol_state, 1.0)
protos_mtx[i, :, :] = symbol_state
In [18]:
mtx_min = np.min(protos_mtx)
mtx_max = np.max(protos_mtx)
f, ax = plt.subplots(ncols=mdl.n_protos, figsize=(12, 12))
for i in range(mdl.n_protos):
cax = ax[i].imshow(np.squeeze(protos_mtx[i,...]), vmin=mtx_min, vmax=mtx_max, cmap=plt.cm.Spectral)
ax[i].set_title('#{0}'.format(i))
# move the colorbar to the side ;)
f.subplots_adjust(right=0.8)
cbar_ax = f.add_axes([0.82, 0.445, 0.0125, 0.115])
cb = f.colorbar(cax, cax=cbar_ax)
cb.set_label('Imaginary PLV')
Separate symbols per subject¶
Now we would like to analyze the symbols per subject, per group.
In [19]:
grp_sym_eo = np.array_split(grp_dist_eo, 10, axis=0)
grp_sym_ec = np.array_split(grp_dist_ec, 10, axis=0)
Examine the first subject¶
In [20]:
subj1_eyes_open = grp_sym_eo[0]
subj1_eyes_closed = grp_sym_ec[0]
In [21]:
from dyconnmap.ts import markov_matrix
In [22]:
markov_matrix_eo = markov_matrix(subj1_eyes_open)
markov_matrix_ec = markov_matrix(subj1_eyes_closed)
In [23]:
from mpl_toolkits.axes_grid1 import ImageGrid
f = plt.figure(figsize=(8, 6))
grid = ImageGrid(f, 111,
nrows_ncols=(1,2),
axes_pad=0.15,
share_all=True,
cbar_location="right",
cbar_mode="single",
cbar_size="7%",
cbar_pad=0.15,
)
im = grid[0].imshow(markov_matrix_eo, vmin=0.0, vmax=1.0, cmap=plt.cm.Spectral)
grid[0].set_xlabel('Prototype')
grid[0].set_ylabel('Prototype')
grid[0].set_title('Eyes Open')
im = grid[1].imshow(markov_matrix_ec, vmin=0.0, vmax=1.0, cmap=plt.cm.Spectral)
grid[1].set_xlabel('Prototype')
grid[1].set_ylabel('Prototype')
grid[1].set_title('Eyes Close')
cb = grid[1].cax.colorbar(im)
cax = grid.cbar_axes[0]
axis = cax.axis[cax.orientation]
axis.label.set_text("Transition Probability")
plt.show()
In [24]:
from dyconnmap.ts import transition_rate, occupancy_time
In [55]:
tr_eo = transition_rate(subj1_eyes_open)
tr_ec = transition_rate(subj1_eyes_closed)
print("""
Transition rate
===============
Eyes open: {0:.3f}
Eyes closed: {1:.3f}
""".format(tr_eo, tr_ec))
Transition rate =============== Eyes open: 0.580 Eyes closed: 0.765
In [61]:
occ_eo = occupancy_time(subj1_eyes_open)[0]
occ_ec = occupancy_time(subj1_eyes_closed)[0]
print("""
Occupancy time
==============
State \t 0 \t 1 \t 2 \t 3 \t 4
-----
Eyes open \t {0:.3f} \t {1:.3f} \t {2:.3f} \t {3:.3f} \t {4:.3f}
Eyes closed \t {5:.3f} \t {6:.3f} \t {7:.3f} \t {8:.3f} \t {9:.3f}
""".format(*occ_eo, *occ_ec))
Occupancy time
==============
State 0 1 2 3 4
-----
Eyes open 0.000 0.336 0.017 0.059 0.000
Eyes closed 0.017 0.092 0.025 0.050 0.042