#!/usr/bin/env python
# coding: utf-8

# # 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[ ]:


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
get_ipython().run_line_magic('matplotlib', 'inline')
matplotlib.rcParams['figure.figsize'] = (20.0, 10.0)


# ## 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[ ]:


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


# In[ ]:


list_operation_results(import_op.id)


# In[ ]:


# Copy the id of the ConnectivityIndex obtained above
conn = load_dt(dt_gid='d9be05424e8011e69bf02c4138a1c4ef')


# In[ ]:


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)


# ## 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[ ]:


list_projects()


# In[ ]:


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


# In[ ]:


list_operation_results(launched_operation.id)


# The simulated bold signals can be visualized using matplotlib library.

# In[ ]:


time_series_region_index_id = get_operation_results(launched_operation.id)[1].id


# In[ ]:


#Load time series h5 file
ts = load_dt(time_series_region_index_id) # use the id of the TimeSeriesRegionIndex obtained above
bold_time = ts.time
bold_data = ts.data


# In[ ]:


# 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)


# ## 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[ ]:


# 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


# In[ ]:


list_operation_results(launched_operation.id)


# In[ ]:


fcd_index_id = get_operation_results(launched_operation.id)[0].id


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

# In[ ]:


# Plot the FCD matrix and the FCD matrix segmented in the epochs
FCD = load_dt(fcd_index_id).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)


# ### Functional hubs
# 
# The functional connectivity matrix of each epoch defines a functional network; for each functional network, TVB identifies 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 [7]). 
# 
# Here the functional hub regions of the mouse brain are defined as the regions with the largest eigenvector components, in absolute value, associated with the three largest eigenvalues of the FC matrix.
# 
# 
# The functional hubs are an output of the FCD function (that we have just run), so we can save the results and display them in the mouse brain sections.

# In[ ]:


# copy all the ids of the ConnectivityMeasureIndexes obtained before
connectivity_measure_ids = [i.id for i in get_operation_results(launched_operation.id)[1:]]


# In[ ]:


# Plot the functional hubs extracted in the first epoch of stable functional connectivity
# set visualization parameters
get_ipython().run_line_magic('matplotlib', 'inline')
matplotlib.rcParams['figure.figsize'] = (20.0, 10.0)

import h5py
from mpl_toolkits.axes_grid1 import make_axes_locatable
from tvb.basic.readers import try_get_absolute_path

fig, axes = plt.subplots(1,3)
slice_idy=73
j=0
for conn_measure_id in connectivity_measure_ids:
    f_path = try_get_absolute_path("tvb_data", "mouse/allen_2mm/RegionVolumeMapping.h5")
    f = h5py.File(f_path, 'r', libver='latest')
    Vol=f['array_data'][:,:,:]
    f_path = try_get_absolute_path("tvb_data", "mouse/allen_2mm/StructuralMRI.h5")
    f = h5py.File(f_path, 'r', libver='latest')
    template=f['array_data'][:,:,:]
    conn_measure = load_dt(conn_measure_id)
    eig=conn_measure.array_data
    for i in range(np.shape(eig)[0]):
        Vol = np.where(Vol==i, eig[i], Vol)
    Vol = np.ma.masked_where(Vol < (np.amax(eig)/2), Vol)
    im1 = axes[j].imshow((template[:,slice_idy,:].T)[::-1], cmap='gray', vmin=template.min(), vmax=template.max())
    cax = axes[j].imshow((Vol[:,slice_idy,:].T)[::-1], cmap='YlOrRd', alpha=1, vmin=np.amax(eig)/2., vmax=np.amax(eig))
    axes[j].axis('off')
    axes[j].set_title(conn_measure.title)
    divider = make_axes_locatable(axes[j])
    cax1 = divider.append_axes("right", size="5%", pad=0.05)
    axcb=plt.colorbar(cax,cax1,ticks=[np.amax(eig)/2.,np.amax(eig)],orientation='vertical')
    axcb.set_ticklabels(['Max/2', 'Max']) 
    axcb.set_label('Eigenvector components')
    j=j+1


# ## References
# 
# [1] Melozzi, Francesca, Marmaduke Woodman, Viktor Jirsa, and Christophe Bernard. "The Virtual Mouse Brain: A Computational Neuroinformatics Platform To Study Whole Mouse Brain Dynamics." bioRxiv (2017): 123406.
# 
# 
# [2] Oh, Seung Wook, Julie A. Harris, Lydia Ng, Brent Winslow, Nicholas Cain, Stefan Mihalas, Quanxin Wang et al. "A mesoscale connectome of the mouse brain." Nature 508, no. 7495 (2014): 207-214.
# 
# 
# [3] Deco Gustavo, Ponce Alvarez Adrian, Dante Mantini, Gian Luca Romani, Patric Hagmann and Maurizio Corbetta. Resting-State Functional Connectivity Emerges from Structurally and Dynamically Shaped Slow Linear Fluctuations. The Journal of Neuroscience 32(27), 11239-11252, 2013.
# 
# 
# [4] Hansen, Enrique CA, Demian Battaglia, Andreas Spiegler, Gustavo Deco, and Viktor K. Jirsa. "Functional connectivity dynamics: modeling the switching behavior of the resting state." Neuroimage 105 (2015): 525-535.
# 
# 
# [5] Friston, Karl J., Andrea Mechelli, Robert Turner, and Cathy J. Price. "Nonlinear responses in fMRI: the Balloon model, Volterra kernels, and other hemodynamics." NeuroImage 12, no. 4 (2000): 466-477.
# 
# 
# [6] Preti, Maria Giulia, Thomas AW Bolton, and Dimitri Van De Ville. "The dynamic functional connectome: State-of-the-art and perspectives." NeuroImage (2016).
# 
# 
# [7] Newman, Mark EJ. "The mathematics of networks." The new palgrave encyclopedia of economics 2, no. 2008 (2008): 1-12.

# In[ ]: