JDG Lab Notebook

Reproducing Deco et al. 2013

Entry Format: results first

The following are my notes on reproducing the bifurcation diagrams in Deco et al. JN 2013 'Resting State Functional Connectivity Emerges from Structurally and Dynamically Shaped Slow Linear Fluctuations', using TVB.

Contents

Overview and summary of results

Background

We are interested in replicating two figures from Deco et al. 2013.

The first is panel B in figure 2:

(0_0)
broken link Deco et al. 2013, Figure 2. The figure caption in the paper reads as follows: "A, B, Mean-field analyses of the attractor landscape of the cortical spiking network (A) and the reduced DMF model (B) as a function of the global coupling strength. Each point represents the maximum firing rate activity among all nodes. In both cases, depending on the coupling strength (G) and the initial conditions, the network converges to one of three possible dynamical regimens: for low values of G, the network converges to a single stable state of low firing activity (i.e., the spontaneous state); for increasing values of G, new stable states of high activity together with the spontaneous state coexist; for even larger values of G, the spontaneous state becomes unstable. Insets, For a particular set of parameters in the multistability regimen (with spontaneous state), the set of possible stable states (colored dots). For visualization, the states (firing rate vectors of dimension, N = 66) are projected in the two first principal components of the state space using the PCA method. C, D, These states are represented by the specific configuration of firing rates shown for the spiking model (C) and its dynamic mean field reduction (D). These are as follows: the spontaneous state of low activity, plus three states for which some nodes have a high firing rate. spont., Spontaneous."

Panel B shows the network bifurcations diagrams of the reduced wong-wang / dynamic mean field model. These consist of scatter plots of the maximum (over tie) of the projection of the network activity on to its first and second principal components.

The DMF bifurcation diagram is identical to that of the spiking network model; an important result for the purposes of the paper. Here we are not concerned at all with the spiking network model, just the DMF model.

The key feature that is shown in panels A and B is the existence of 3 separate regimes:

  1. 'Single state' - Globally low firing rate

  2. 'Multistability (with spont. state)' - Low firing rate coexists but a new high firing rate regime has also appeared. The PCA separates them as two clean, contigous lines

  3. 'Multistability (high firing rate)' - Same as 2., but now the low firing rate regime disppears

Global coupling strength (G) serves as a control parameter that separates the transitions from regime 1 to regime 2 and regime 2 to regime 3. These 2 bifurcations occur at G = 1.5 and G = 2.5, respectively.

The main result of the paper is then shown in panel A of figure 3:

(0_0)
broken link Deco et al. 2013 Fig 3. The figure caption in the paper reads as follows: "A, The fitting of empirical data as measured by the correlation between simulated and empirical FC, as a function of the global coupling parameter G, was calculated for both the spiking (red stars) and the reduced DMF models for two levels of noise (blue points: σ = 0.001, light blue trace: σ = 0.010). In the two models, the best fit is achieved at the edge of the second bifurcation (black vertical line) For σ = 0.010, the gray vertical indicates the critical value of G, for which the spontaneous state is destabilized. B, Comparison of pairwise correlation coefficients obtained using empirical BOLD activity and simulated BOLD activity from the spiking network. C, The subdiagonal elements represent the empirical FC, in warm color code, whereas upper-diagonal elements represent the FC generated by the spiking model (cold colors). D, Comparison of pairwise correlation coefficients obtained using empirical BOLD activity and simulated BOLD activity from the DMF model. E, The subdiagonal elements represent the empirical FC, in warm color code, whereas upper-diagonal elements represent the FC generated by the DMF model (cold colors). dyn., Dynamics."

The scatter plot in panel A shows fitting between simulated and empirical FC as a function of G. The key result is that this fitting is maximal at exactly the location of the second bifurcation.

I want to replicate the plots in panel B of figure 2 and panel A of figure 3 above, using the ReducedWongWang model in TVB.

It took quite a lot of experimentation and head scratching, but I am now confident I am now able to get a reasonably close replication of Deco figure 2, for a number of different connectivity matrices (TVB default, original Deco/Hagmann matrix, and some HCP subjects).

I can also indicate clearly the behaviour shown in the PCA plots in a few different complementary visual representations of the simulated time series.

Replication of Deco figure 3 panel A (which is the main result of the paper) is not quite so compelling at this point; although still work in progress. The main issues are a) how closely the bifurcation point and maximal fit point line up, and b) the magnitude of the fits, which are a lot lower than the maxmium of 0.5 as in the figure.

This is a 'results first' format notebook entry, so we look at some figures in the next few sections, then summarize and sign offf; all the analysis and simulation code used to generate the figures is then given in the second half of the notebook.

Right, let's get to it:

Connectivities

Shown below are structural connectivity and, when available, functional connectivity matrices for our three test cases.

The first two columns are before and after normalizing the connectivity weights - which is done using the TVB 'region-based' weights scaling option in TVB - which rescales all weights such that the maximum absolute value of the cumulative input to any region is unity. The third column is FC.

(0_0)
broken link Structural and functional connectivity matrices. Row 1 = TVB default connectivity (no FC available for this data). Row 2 - Hagmann/Deco connectivity. Row 3 - HCP WUMinn example subject (148941).

Bifurcation diagrams

For the bifurcation diagrams (unlike the BOLD FC correlations), we can do many, relatively short simulation runs, relatively quickly. Here we use 10 second runs.

Calculating firing rates

First, we calculate firing rates $H(x)$ from synaptic activation variable $S$ returned by TVB, using the equation

\begin{array} . H(x_i) &= \frac{ax_i - b}{1 - exp( -d(ax_i-b))} \\ x_i &= w J_N S_i + G J_N \sum_j C_{ij} S_j + I_0 \end{array}

We have to do this manually as TVB doesn't give these numbers by default, although it does calculate them on each integration step using the aquation abotvee an other.

Visualize time series

When we look at the time series (either synaptic activation or firing rates), we can see the transition through the three regimes described above very clearly in the WWD model time series, given the correct G parameter values:

(0_0)
broken link State variable time series. Shown are time series of the synaptic activation state variable S, over all network nodes, for various values of global coupling G, for the three connectivities used here.

Basically, from the low to high Gs, the majority of nodes move from being along the thick lower black lines, up to the thick upper black lines with increasing Gs.

Compute PCA

To make the bifurcation figures as per Deco figure 2, we compute PCA on the firing rate (and, for completeness, also synaptic activation) time series, and project on to the first two components, as follows:

\begin{array} _ U,S,V &= svd(Y) \\ \hat{U} &= U_{1:N} \\ P &= ((\hat{U} \hat{U}^T )Y^T)^T \\ \end{array}

Where $N=2$ is the number of components used for the projection.

These projections now give us N time series for each simulation run. Again, we can see the transition through the three regimes by looking at the projected time series

Bifurcation scatter plots

Now, we take the maximum over time and for each of the two component projections, and represent on a scatter plot. These correspond to the bifurcation plots in Deco figure 2

(0_0)
broken link Bifurcation diagram scatter plots. Shown are scatter plots of the maximal firing rate for the first two principal components of simulated network activity, as a function of global coupling G. The vertical lines indicate (qualitatively chosen) bifurcation points - first where a the high firing state emerges in some nodes, and second where the remaining nodes transition from low firing rate regime to the high firing rate regime.

So we can see here that both components start in the low firing rate regime, and then one component rises linearly with coupling while the other stays low, and then the second one leaves the low regimes and joins the first component in a high regime, climibing with coupling strength.

These figures differ a bit from the Deco figures in the extent of the low firing rate regime (i.e. the flat line at the bottom in the figures from the paper), the discontinuity in the jump to the high regime (whereas here we have a continuous climb, albeit a fast one); and also the values of G.

Removing the transient

Note that the projections above are calculated after a burn-in of 5 seconds to remove the initial conditions transient. If you don't do that, and compute the PCA over the full run, your PCA plots will look very different (basically you won't get low firing rate states for low G; the initial conditions will make there always be an initial high firing rate (on average, for some nodes), irrspective of the value of G).

Connectivity normalization

The above plots were donw using the connecitvity weights scaling described above. If that isn't done, then, unsurprisingly, the values of G at which the bifurcation occur can be wildly different.

Histogram plots

Another way of visualzing the transitions through the regimes is to compute a histogram of the node values for a single time slice (e.g. the final one), for the different G values, like so:

(0_0)
broken link KDE histograms of node states as a function of G. Shown are kernel density estimates (smooth histograms) over nodess for the final time slide in the short-run simulations. The colour scales move from low G (blues) to high G (reds), and indicate clearly the transition from the low to high firing rate regime as a function of global coupling.

Correlation with functional connectivity

As noted, the main result in Deco et al. is the observation that the fit of simulated FC to empirical FC is maximal at the bifurcation point.

Reproducing this is a little more involved than the bifurcation diagrams because a) we need an FC matrix to compare to, and b) we now need to do some very long simulation runs (in Deco et al. they used a rather excessive 20 minutes, which corresponded to the acquisition time of the empirical fMRI data).

So to do this we do exactly the same as in the previous section, but now for longer runs times and less G values. In the following I ran sims for 60 seconds. Longer simulation runs are, well, running.

(0_0)
broken link Simulated-empirical FC fits as a functio of G. Shown are correlations between simulated and empirical FC, and also simulated FC and empirical SC, for different values of global coupling. The correlations are computed for both the full matrices and for the lower triangle only ('tril').

Now if we put these FC correlation plots alongside the PCA bifurcation diagrams, we can see that the maximum FC is achieved around the vicinity of the second bifurcation.

(0_0)
broken link Simulated FC fits together with bifurcation diagrams. This figure combines the FC fit and bifurcation scatter plot figures above, so as to see clearly the correspondence between the maximal FC fit points and the bifurcation points. The vertical line in both sets of figures is the maximum of the BOLD FC fit, which occurs in the vicinity of the second bifurcation, as predicted.

For completeness, take a look at the simulated FC matrices also:

(0_0)
broken link SC,FC, and simulated FC matrices for Deco/Hagmann data. Shown are empirical SC, FC, and then simulated FC at the first, optimal, and last G value in the range used.
(0_0)
broken link SC,FC, and simulated FC matrices for HCP data. Shown are empirical SC, FC, and then simulated FC at the first, optimal, and last G value in the range used.

Summary

Some final quick comments on how these analyses were done and imporant points to remember:

Scaling

All the simulations done here used re-scaled connectivity matices, as described above. This has the handy property that it places the location of the bifurcation point in fairly similar locations for all matrices, even though the original values of the weights in each case were vastly different. One consequence of using scaling, however, is that the G values now aren't directly comparable with those from the paper.

I have done some explorations of sim runs with the unscaled connectivities, but they didn't make it into the final cut on this doc, so far at least. They could also be added to these notes for reference.

Run time

As noted, the long run FC fit simulations here were only run for 1 minute. Watch this space for full 20 minute runs. Hopefully the FC fits will be a lot higher for them.

Correlation fits

The simulated-empirical FC correlation fits are not as high here as in Deco et al. The most obvious reason for this is the run time, as above. Another is that I din't know for sure that the Hagmann/Deco SC/FC matrices used here are the same as the ones from the paper (I've kind of inherited them third hand). Actually the Hagmann/Deco FC matrix is a bit strange, as it has mostly negative values. Note also - when we just use the lower triangle of the matrix for the ocrrelations, they are substantially reduced. And, of course, if we use $R^2$ rather than $R$, then the fit values all drop by roughly half. So it would be nice if these numbers were a bit better, tbh.

Analysis documentation + code

Notebook Setup

Define some variables

# put some system-specific variables in the namespace ('le' dict)
%run ~/set_localenv_vars.py

nb_name = 'reproducing_deco2013'

outdir = '%s/notebooks/%s' %(le['data_dir'], nb_name)

# this is where stuff actually lives
outdir_loc = '%s/%s' %('/mnt/titania-hecuba/mcintosh_lab/john/Data/notebooks', nb_name)

# symlink the two
import os
if not os.path.isdir(outdir_loc): os.makedirs(outdir_loc)
if not os.path.isdir(outdir): os.system('ln -s %s %s' %(outdir_loc,outdir))


!mkdir -p $outdir/pics

# root path to tvb-library, tvb-data, etc. libraries. 
tvb_folder = le['code_dir'] + '/libraries_of_others/github/tvb-pack'

# stuff for workdocs-cloudfiles

aws_key = 'drowssaperucesyreva'
aws_secret = '?teytidettopsuoyevah'

Importage

# Generic imports

import  os,sys,glob,numpy as np,pandas as pd
from scipy.io import loadmat
from joblib import Parallel,delayed
from sklearn.decomposition import PCA
from scipy.linalg import svd
from datetime import datetime

# Visualization stuff

%matplotlib inline
from matplotlib import pyplot as plt
from IPython.display import Image,display as d
import seaborn as sns


# aws api and workdocs-cloudfiles stuff

sys.path.append(le['ipynb_workdocs_dir'])
from new_utils import nb_fig,cloudfiles_nb


# TVB stuff

nbso,nbse = sys.stdout,sys.stderr # hack part 1/2 to keep output printing properly
sys.path.append(tvb_folder +'/tvb-library')
sys.path.append(tvb_folder + '/tvb-data')
from tvb.simulator.lab import *
sys.stdout,sys.stderr = nbso,nbse  # ...hack part 2/2


from copy import copy # seems that this needs to be done last otherwise it gets 
                      # replaced with numpy copy function

Calico document extensions

%%javascript
IPython.load_extensions('calico-spell-check', 
                        'calico-document-tools', 'calico-cell-tools');

Initialize aws api and workdocs-cloudfiles folder

cnb = cloudfiles_nb('aws', [aws_key,aws_secret])
cnb.initialize_folder(nb_name)

Go to output folder

os.chdir(outdir)

Ok, let's get cracking

Define some functions

Simulation workhorse function

# Run sim function
def run_wwd_sim(conn, cs=1., D=0.001, cv=3.0, dt=0.5, 
                simlen=1e3,int_type='heun_stoch',
                wwd_params={},initconds = None,return_B=True,return_SandH=True):
    
    if int_type=='heun_stoch':
        solver = integrators.HeunStochastic(dt=dt,noise=noise.Additive(nsig=array([D])))
    elif int_type=='euler_stoch':
        solver = integrators.EulerStochastic(dt=dt,noise=noise.Additive(nsig=array([D])))
    else: 
        solver = integrators.HeunDeterministic(dt=dt)
        
    sim = simulator.Simulator(
        model=models.ReducedWongWang(**wwd_params),
        connectivity=conn,
        coupling=coupling.Linear(a=cs),
        integrator=solver,
        monitors=(monitors.TemporalAverage(period=5.0),# 200 Hz
                  monitors.Bold())
    )
    sim.configure()
    
    if initconds:
        _ = sim.run(1)
        if initconds == 0:
            sim.history = np.zeros_like(test_sim.history)
        else: 
            sim.history = initconds
        
    res = sim.run(simulation_length=simlen)
    (S_ts,S_ys),(B_ts,B_ys) = res
    
    if return_B:
      df_B = pd.DataFrame(np.squeeze(B_ys), index=B_ts)        
    else: 
      df_B = None
        
        
    if return_SandH:

      df_S = pd.DataFrame(np.squeeze(S_ys), index=S_ts)
    
    
      _a = sim.model.a[0]
      _b = sim.model.b[0]
      _d = sim.model.d[0]
      _w = sim.model.w[0]
      _I_o = sim.model.I_o[0]
      _J_N = sim.model.J_N[0]
      _cs = sim.coupling.a[0]
      _SC = sim.connectivity.weights
    
      _x = rww_compute_x(df_S.values,w=_w,j=_J_N,io=_I_o,cs=_cs,SC=_SC)
      _h=1000*rww_compute_h(_x,a=_a,b=_b,d=_d)
                       
      df_H = pd.DataFrame(_h,index=df_S.index,columns=df_S.columns)
    
    else: 
   
      df_S,df_H = None,None
        
    return df_S,df_H,df_B

Functions to compute WWD X and H variables from state variable S

def rww_compute_h(x,a=None,b=None,d=None):
  h = (a*x - b) / (1 - np.exp(-d*(a*x - b)))
  return h
    
def rww_compute_x(S, w=None,j=None,io=None,cs=None,SC=None): 
  x = w*j*S + io + cs*j*np.dot(S,SC)
  return x

Wrapper to run parameter sweep over Gs. Parallellization over cores / threads with joblib. 

# Parameter sweep function
def run_wwd_Gsweep(Gs,fixed_params):
    
    paramslist = []
    for g in Gs: 
        newparams = copy(fixed_params)
        newparams['cs'] = np.array([g])
        paramslist.append(newparams)
    
    sim_res = Parallel(n_jobs=2)\
             (delayed(run_wwd_sim)(**p)\
             for p in paramslist)
    
    
    if fixed_params['return_B'] == True:
      Bs = {g: r[2] for g,r in zip(Gs,sim_res)}
      df_B_all = pd.concat(Bs)
      df_B_all.index.names = ['G', 't']
      df_B_all.columns.names = ['region']        
    else: 
      df_B_all = None
        
        
    if fixed_params['return_SandH'] == True:
        
      Ss = {g: r[0] for g,r in zip(Gs,sim_res)}
      df_S_all = pd.concat(Ss)
      df_S_all.index.names = ['G', 't']
      df_S_all.columns.names = ['region']    

      Hs = {g: r[1] for g,r in zip(Gs,sim_res)}
      df_H_all = pd.concat(Hs)
      df_H_all.index.names = ['G', 't']
      df_H_all.columns.names = ['region']    
        
    else: 
        
      df_S_all,df_H_all=None,None
    
    return df_S_all,df_H_all,df_B_all

Functions for computing correlations between empFC/simFC/empSC matrices

# Compute fits to FC
def corrmats(mat1,mat2):

    r_full = corrcoef(mat1.ravel(), mat2.ravel())[0, 1]

    tril_inds = np.tril_indices_from(mat1,-1)
    r_tril = corrcoef(mat1[tril_inds],mat2[tril_inds])[0,1]  
    
    return r_full,r_tril


def get_conn_corrs(dfdict,SC,FC):
    
    corr_res = {}
    
    for k,v in dfdict.items():

        fcmat = v.corr().values

        cmfc_full,cmfc_tril = corrmats(fcmat,FC)
        cmsc_full,cmsc_tril = corrmats(fcmat,SC)
    
        corr_res[k] = {'corrFC': cmfc_full,
                       'corrFCtril': cmfc_tril,
                       'corrSC': cmsc_full,
                       'corrSCtril': cmsc_tril}
             
    df_corr_res = pd.DataFrame(corr_res) 
        
    df_corr_res.index.names = ['corrvar']
    df_corr_res.columns.names = ['G']        
        
    return df_corr_res

PCA projection calculation, and wrapper over a sim dataframe with multiple Gs

def svdproj(df,mean_centre=False):
  
  df = df.copy()
  if mean_centre is True:   df = df - df.mean(axis=0)
  _u,_s,_v = svd(df.corr())
  proj= (np.dot(_u,_u.T)[:2,:].dot(df.T)).T
  return proj
    
    
def project_dfs_for_gs(df,gs,mean_centre=False,abs_val=False):

  pc1s,pc2s,pc1maxes,pc2maxes,pc1means,pc2means = {},{},{},{},{},{}

  for g in gs: 
        
    df_data = df.copy().ix[g]
    
    pc12= svdproj(df_data,mean_centre=mean_centre)
        
    pc1maxes[g] = np.max(pc12[:,0])
    pc2maxes[g] = np.max(pc12[:,1])

    pc1means[g] = np.mean(pc12[:,0])
    pc2means[g] = np.mean(pc12[:,1])  

  df_pc12means = pd.concat({'pc1': pd.DataFrame(pc1means.values(),index=pc1means.keys())[0],
                          'pc2': pd.DataFrame(pc2means.values(),index=pc2means.keys())[0]})
  df_pc12means = df_pc12means.unstack().T
  df_pc12means.index.names = ['G']
  df_pc12means.sort_index(inplace=True)

  df_pc12maxes = pd.concat({'pc1': pd.DataFrame(pc1maxes.values(),index=pc1maxes.keys())[0],
                          'pc2': pd.DataFrame(pc2maxes.values(),index=pc2maxes.keys())[0]})
  df_pc12maxes = df_pc12maxes.unstack().T
  df_pc12maxes.index.names = ['G']
  df_pc12maxes.sort_index(inplace=True)
  
  return df_pc12means,df_pc12maxes

Load and prepare data

Load in and prepare the SC and FC data from the various sources we are using:

TVB76 default data

# Original weights
TVB_conn = connectivity.Connectivity(load_default=True) 
TVB_conn.speed = np.inf
TVB_conn.tract_lengths = np.zeros_like(TVB_conn.weights)
TVB_conn.configure()
TVB_SC = TVB_conn.weights

# Scaled weights
TVB_conn_sw = connectivity.Connectivity(load_default=True) 
TVB_conn_sw.weights = TVB_conn_sw.scaled_weights('region')    
TVB_conn_sw.speed = np.inf
TVB_conn_sw.tract_lengths = np.zeros_like(TVB_conn_sw.weights)
TVB_conn_sw.configure()
TVB_SC_sw = TVB_conn_sw.weights

Hagmann/Deco data

# Original weights
Hag_SC = loadmat('Hagmann_SC.mat')['C']
Hag_conn = connectivity.Connectivity()
Hag_conn.weights = Hag_SC.copy()
Hag_conn.speed = np.inf
Hag_conn.tract_lengths = np.zeros_like(Hag_SC)
Hag_conn.configure()

# Scaled weights
Hag_conn_sw = connectivity.Connectivity()
Hag_conn_sw.weights = Hag_SC.copy()
Hag_conn_sw.weights = Hag_conn_sw.scaled_weights('region')    
Hag_conn_sw.speed = np.inf
Hag_conn_sw.tract_lengths = np.zeros_like(Hag_SC)
Hag_conn_sw.configure()
Hag_SC_sw = Hag_conn_sw.weights

# Functional connectivity
Hag_FC = loadmat('Hagmann_FC.mat')['FC_emp']

As an aside: the FC matrix is very odd. Most of the entries are negative:

HCP data

# Original weights 
HCP_conn = connectivity.Connectivity.from_file(outdir + '/HCP_148941_sc33.zip')
HCP_conn.speed = np.inf
HCP_conn.tract_lengths = np.zeros_like(HCP_conn.weights)
HCP_conn.configure()
HCP_SC = HCP_conn.weights


# Scaled weights
HCP_conn_sw = connectivity.Connectivity.from_file(outdir + '/HCP_148941_sc33.zip')
HCP_conn_sw.weights = np.log1p(HCP_conn_sw.weights)
HCP_conn_sw.weights = HCP_conn_sw.scaled_weights('region')    
HCP_conn_sw.speed = np.inf
HCP_conn_sw.tract_lengths = np.zeros_like(HCP_conn_sw.weights)
HCP_conn_sw.configure()
HCP_SC_sw = HCP_conn_sw.weights


# Functional connectivity
f = '148941_rfMRI_REST1_LR_hpc200_clean__l2k8_sc33_ts.pkl'
HCP_ts = pd.read_pickle(f)
HCP_FC = HCP_ts.corr().values

Plot all connectivities together

figf = outdir + '/pics/all_SC_FC_conmats.png' 



fig, ax  =plt.subplots(ncols=3,nrows=3,figsize=(16,12)) 


# TVB data

tit = 'TVB SC orig'
a = ax[0][0]
sns.heatmap(TVB_SC, ax=a,xticklabels='',yticklabels='')
a.set_title(tit)

tit = 'TVB SC sw'
a = ax[0][1]
sns.heatmap(TVB_SC_sw, ax=a,xticklabels='',yticklabels='')
a.set_title(tit)

tit = '(no TVB FC)'
a = ax[0][2]
a.set_title(tit)
a.axis('off')

# Hagmann/Deco data

tit = 'Hag SC orig'
a = ax[1][0]
sns.heatmap(Hag_SC, ax=a,xticklabels='',yticklabels='')
a.set_title(tit)

tit = 'Hag SC sw'
a = ax[1][1]
sns.heatmap(Hag_SC_sw, ax=a,xticklabels='',yticklabels='')
a.set_title(tit)

tit = 'Hag FC'
a = ax[1][2]
sns.heatmap(Hag_FC, ax=a,xticklabels='',yticklabels='')
a.set_title(tit)



# HCP data

tit = 'HCP SC orig'
a = ax[2][0]
sns.heatmap(HCP_SC, ax=a,xticklabels='',yticklabels='')
a.set_title(tit)

tit = 'HCP SC sw'
a = ax[2][1]
sns.heatmap(HCP_SC_sw, ax=a,xticklabels='',yticklabels='')
a.set_title(tit)

tit = 'HCP FC'
a = ax[2][2]
sns.heatmap(HCP_FC, ax=a,xticklabels='',yticklabels='')
a.set_title(tit)



plt.savefig(figf, bbox_inches='tight')
plt.close()

Connectivity matrix weight histograms

figf = outdir + '/pics/all_SC_FC_conmat_histograms.png' 


fig, ax  =plt.subplots(ncols=3,nrows=3,figsize=(12,12)) 



# TVB data

tit = 'TVB SC orig'
a = ax[0][0]
pd.DataFrame(TVB_SC.ravel())[0].hist(ax=a,bins=6);
a.set_yscale('log')
a.set_xlim([0,3]);
a.set_title(tit)

tit = 'TVB SC sw'
a = ax[0][1]
pd.DataFrame(TVB_SC_sw.ravel())[0].hist(ax=a,bins=6);
a.set_yscale('log')
a.set_xlim([0,1])
a.set_title(tit)

tit = '(no TVB FC)'
a = ax[0][2]
#a.set_title(tit)
a.axis('off')

# Hagmann/Deco data

tit = 'Hag SC orig'
a = ax[1][0]
pd.DataFrame(Hag_SC.ravel())[0].hist(ax=a,bins=50)
a.set_yscale('log')
a.set_xlim([0,0.1])
a.set_title(tit)

tit = 'Hag SC sw'
a = ax[1][1]
pd.DataFrame(Hag_SC_sw.ravel())[0].hist(ax=a,bins=50)
a.set_yscale('log')
a.set_xlim([0,0.2])
a.set_title(tit)

tit = 'Hag FC'
a = ax[1][2]
pd.DataFrame(Hag_FC.ravel())[0].hist(ax=a,bins=50);
a.set_xlim([-0.2,0.5])
a.set_title(tit)



# HCP data

tit = 'HCP SC orig'
a = ax[2][0]
pd.DataFrame(HCP_SC.ravel())[0].hist(ax=a,bins=50);
a.set_yscale('log')
a.set_xlim([0,10000]);
a.set_title(tit)

tit = 'HCP SC sw'
a = ax[2][1]
pd.DataFrame(HCP_SC_sw.ravel())[0].hist(ax=a,bins=50);
a.set_yscale('log')
a.set_xlim([0,0.1]);
a.set_title(tit)

tit = 'HCP FC'
a = ax[2][2]
pd.DataFrame(HCP_FC.ravel())[0].hist(ax=a,bins=50);
a.set_xlim([-1,1]);
a.set_title(tit)


plt.savefig(figf, bbox_inches='tight')
plt.close()

Short sim runs

(for bifurcation diagrams)

TVB data

ndps = {'w': 0.9, 'I_o': 0.33} # non-default params
fixed_params = {'D': 1E-6, 'wwd_params': ndps,'return_SandH': True, 'return_B': True}
fixed_params['conn'] = TVB_conn_sw

# Run G sweep with short runs
Gs = np.linspace(0.0001,0.1,50)
fixed_params['simlen'] = 10000
df_TVB_sw_S1,df_TVB_sw_H1,df_TVB_sw_B1 = run_wwd_Gsweep(Gs,fixed_params)    

df_TVB_sw_S1_pc12means,df_TVB_sw_S1_pc12maxes = project_dfs_for_gs(df_TVB_sw_S1,Gs)    
df_TVB_sw_S1_nt = pd.concat({g: df_TVB_sw_S1.ix[g].ix[5000:] for g in Gs})
df_TVB_sw_S1_nt.index.names = ['G', 't']
df_TVB_sw_S1_nt_pc12means,df_TVB_sw_S1_nt_pc12maxes = project_dfs_for_gs(df_TVB_sw_S1_nt,Gs)    
df_TVB_sw_H1_pc12means,df_TVB_sw_H1_pc12maxes = project_dfs_for_gs(df_TVB_sw_H1,Gs)    
df_TVB_sw_H1_nt = pd.concat({g: df_TVB_sw_H1.ix[g].ix[5000:] for g in Gs})
df_TVB_sw_H1_nt.index.names = ['G', 't']
df_TVB_sw_H1_nt_pc12means,df_TVB_sw_H1_nt_pc12maxes = project_dfs_for_gs(df_TVB_sw_H1_nt,Gs)  
df_TVB_sw_S1_pc12means_nmc,df_TVB_sw_S1_pc12maxes_nmc = project_dfs_for_gs(df_TVB_sw_S1,Gs,mean_centre=False)
df_TVB_sw_S1_nt_pc12means_nmc,df_TVB_sw_S1_nt_pc12maxes_nmc = project_dfs_for_gs(df_TVB_sw_S1_nt,Gs,mean_centre=False)    
df_TVB_sw_H1_pc12means_nmc,df_TVB_sw_H1_pc12maxes_nmc = project_dfs_for_gs(df_TVB_sw_H1,Gs,mean_centre=False)
df_TVB_sw_H1_nt_pc12means_nmc,df_TVB_sw_H1_nt_pc12maxes_nmc = project_dfs_for_gs(df_TVB_sw_H1_nt,Gs,mean_centre=False)

Hagmann / Deco data

ndps = {'w': 0.9, 'I_o': 0.33} # non-default params
fixed_params = {'D': 1E-6, 'wwd_params': ndps,'return_SandH': True, 'return_B': True}
fixed_params['conn'] = Hag_conn_sw

# Run G sweep with short runs
Gs = np.linspace(0.0001,0.1,50)
fixed_params['simlen'] = 10000
df_Hag_sw_S1,df_Hag_sw_H1,df_Hag_sw_B1 = run_wwd_Gsweep(Gs,fixed_params)    

df_Hag_sw_S1_pc12means,df_Hag_sw_S1_pc12maxes = project_dfs_for_gs(df_Hag_sw_S1,Gs)    
df_Hag_sw_S1_nt = pd.concat({g: df_Hag_sw_S1.ix[g].ix[5000:] for g in Gs})
df_Hag_sw_S1_nt.index.names = ['G', 't']
df_Hag_sw_S1_nt_pc12means,df_Hag_sw_S1_nt_pc12maxes = project_dfs_for_gs(df_Hag_sw_S1_nt,Gs)    
df_Hag_sw_H1_pc12means,df_Hag_sw_H1_pc12maxes = project_dfs_for_gs(df_Hag_sw_H1,Gs)    
df_Hag_sw_H1_nt = pd.concat({g: df_Hag_sw_H1.ix[g].ix[5000:] for g in Gs})
df_Hag_sw_H1_nt.index.names = ['G', 't']
df_Hag_sw_H1_nt_pc12means,df_Hag_sw_H1_nt_pc12maxes = project_dfs_for_gs(df_Hag_sw_H1_nt,Gs)  
df_Hag_sw_S1_pc12means_nmc,df_Hag_sw_S1_pc12maxes_nmc = project_dfs_for_gs(df_Hag_sw_S1,Gs,mean_centre=False)
df_Hag_sw_S1_nt_pc12means_nmc,df_Hag_sw_S1_nt_pc12maxes_nmc = project_dfs_for_gs(df_Hag_sw_S1_nt,Gs,mean_centre=False)    
df_Hag_sw_H1_pc12means_nmc,df_Hag_sw_H1_pc12maxes_nmc = project_dfs_for_gs(df_Hag_sw_H1,Gs,mean_centre=False)
df_Hag_sw_H1_nt_pc12means_nmc,df_Hag_sw_H1_nt_pc12maxes_nmc = project_dfs_for_gs(df_Hag_sw_H1_nt,Gs,mean_centre=False)

HCP data

ndps = {'w': 0.9, 'I_o': 0.33} # non-default params
fixed_params = {'D': 1E-6, 'wwd_params': ndps,'return_SandH': True, 'return_B': True}
fixed_params['conn'] = HCP_conn_sw

# Run G sweep with short runs
Gs = np.linspace(0.0001,0.1,50)
fixed_params['simlen'] = 10000
df_HCP_sw_S1,df_HCP_sw_H1,df_HCP_sw_B1 = run_wwd_Gsweep(Gs,fixed_params)    

df_HCP_sw_S1_pc12means,df_HCP_sw_S1_pc12maxes = project_dfs_for_gs(df_HCP_sw_S1,Gs)    
df_HCP_sw_S1_nt = pd.concat({g: df_HCP_sw_S1.ix[g].ix[5000:] for g in Gs})
df_HCP_sw_S1_nt.index.names = ['G', 't']
df_HCP_sw_S1_nt_pc12means,df_HCP_sw_S1_nt_pc12maxes = project_dfs_for_gs(df_HCP_sw_S1_nt,Gs)    
df_HCP_sw_H1_pc12means,df_HCP_sw_H1_pc12maxes = project_dfs_for_gs(df_HCP_sw_H1,Gs)    
df_HCP_sw_H1_nt = pd.concat({g: df_HCP_sw_H1.ix[g].ix[5000:] for g in Gs})
df_HCP_sw_H1_nt.index.names = ['G', 't']
df_HCP_sw_H1_nt_pc12means,df_HCP_sw_H1_nt_pc12maxes = project_dfs_for_gs(df_HCP_sw_H1_nt,Gs)  
df_HCP_sw_S1_pc12means_nmc,df_HCP_sw_S1_pc12maxes_nmc = project_dfs_for_gs(df_HCP_sw_S1,Gs,mean_centre=False)
df_HCP_sw_S1_nt_pc12means_nmc,df_HCP_sw_S1_nt_pc12maxes_nmc = project_dfs_for_gs(df_HCP_sw_S1_nt,Gs,mean_centre=False)    
df_HCP_sw_H1_pc12means_nmc,df_HCP_sw_H1_pc12maxes_nmc = project_dfs_for_gs(df_HCP_sw_H1,Gs,mean_centre=False)
df_HCP_sw_H1_nt_pc12means_nmc,df_HCP_sw_H1_nt_pc12maxes_nmc = project_dfs_for_gs(df_HCP_sw_H1_nt,Gs,mean_centre=False)

Save to file

tvb_sw_ss_file = outdir + '/TVB_sw_short_sims.hdf'

df_TVB_sw_S1.to_hdf(tvb_sw_ss_file, 'df_TVB_sw_S1')
df_TVB_sw_H1.to_hdf(tvb_sw_ss_file, 'df_TVB_sw_H1')
df_TVB_sw_B1.to_hdf(tvb_sw_ss_file, 'df_TVB_sw_B1')

df_TVB_sw_H1_pc12maxes.to_hdf(tvb_sw_ss_file, 'df_TVB_sw_H1_pc12maxes')
df_TVB_sw_S1_pc12maxes.to_hdf(tvb_sw_ss_file, 'df_TVB_sw_S1_pc12maxes')

df_TVB_sw_H1_nt_pc12maxes.to_hdf(tvb_sw_ss_file, 'df_TVB_sw_H1_nt_pc12maxes')
df_TVB_sw_S1_nt_pc12maxes.to_hdf(tvb_sw_ss_file, 'df_TVB_sw_S1_nt_pc12maxes')

hag_sw_ss_file = outdir + '/Hag_sw_short_sims.hdf'

df_Hag_sw_S1.to_hdf(hag_sw_ss_file, 'df_Hag_sw_S1')
df_Hag_sw_H1.to_hdf(hag_sw_ss_file, 'df_Hag_sw_H1')
df_Hag_sw_B1.to_hdf(hag_sw_ss_file, 'df_Hag_sw_B1')

df_Hag_sw_H1_pc12maxes.to_hdf(hag_sw_ss_file, 'df_Hag_sw_H1_pc12maxes')
df_Hag_sw_S1_pc12maxes.to_hdf(hag_sw_ss_file, 'df_Hag_sw_S1_pc12maxes')

df_Hag_sw_H1_nt_pc12maxes.to_hdf(hag_sw_ss_file, 'df_Hag_sw_H1_nt_pc12maxes')
df_Hag_sw_S1_nt_pc12maxes.to_hdf(hag_sw_ss_file, 'df_Hag_sw_S1_nt_pc12maxes')

hcp_sw_ss_file = outdir + '/HCP_sw_short_sims.hdf'

df_HCP_sw_S1.to_hdf(hcp_sw_ss_file, 'df_HCP_sw_S1')
df_HCP_sw_H1.to_hdf(hcp_sw_ss_file, 'df_HCP_sw_H1')
df_HCP_sw_B1.to_hdf(hcp_sw_ss_file, 'df_HCP_sw_B1')

df_HCP_sw_H1_pc12maxes.to_hdf(hcp_sw_ss_file, 'df_HCP_sw_H1_pc12maxes')
df_HCP_sw_S1_pc12maxes.to_hdf(hcp_sw_ss_file, 'df_HCP_sw_S1_pc12maxes')

df_HCP_sw_H1_nt_pc12maxes.to_hdf(hcp_sw_ss_file, 'df_HCP_sw_H1_nt_pc12maxes')
df_HCP_sw_S1_nt_pc12maxes.to_hdf(hcp_sw_ss_file, 'df_HCP_sw_S1_nt_pc12maxes')

Make combined figures

time series plots

figf = outdir + '/pics/all_S_timeseries_vs_G.png'


fig, ax = plt.subplots(ncols=3,nrows=6, figsize=(16,20))


df = df_TVB_sw_S1
Gs = np.unique(df.index.get_level_values('G'))
Gstoplot = [Gs[0], Gs[5], Gs[15],
            Gs[25], Gs[35], Gs[45]]
for g_it,g in enumerate(Gstoplot):
  a = ax[g_it][0]
  df.ix[g].plot(legend=False,c='k', linewidth=0.5,alpha=0.5,ax=a,
                   title='TVB sw S1\ng=%1.3f' %g)
    

    
df = df_Hag_sw_S1
Gs = np.unique(df.index.get_level_values('G'))
Gstoplot = [Gs[0], Gs[5], Gs[15],
            Gs[25], Gs[35], Gs[45]]
for g_it,g in enumerate(Gstoplot):
  a = ax[g_it][1]
  df.ix[g].plot(legend=False,c='k', linewidth=0.5,alpha=0.5,ax=a,
                   title='Hag sw S1\ng=%1.3f' %g)
    


df = df_HCP_sw_S1
Gs = np.unique(df.index.get_level_values('G'))
Gstoplot = [Gs[0], Gs[5], Gs[15],
            Gs[25], Gs[35], Gs[45]]
for g_it,g in enumerate(Gstoplot):
  a = ax[g_it][2]
  df.ix[g].plot(legend=False,c='k', linewidth=0.5,alpha=0.5,ax=a,
                   title='HCP sw S1\ng=%1.3f' %g)
    

    
nrows,ncols = ax.shape
for row in range(nrows):
  for col in range(ncols):
    a = ax[row][col]
    if row != (nrows-1):
      a.set_xticklabels('')
      a.set_xlabel('')
    if col != 0:
      a.set_yticklabels('')
      a.set_ylabel('')
        
    
plt.tight_layout()


plt.savefig(figf, bbox_inches='tight')
plt.close()

bifurcation scatter plots

figf = outdir + '/pics/all_bif_scatter_plots.png'


fig, ax  = plt.subplots(ncols=3,figsize=(12,6))

bifpoint1,bifpoint2 = 0,0.025
tit = 'TVB sw nt \n(bif1: %s, bif2: %s' %(bifpoint1,bifpoint2)
a = ax[0]
df = df_TVB_sw_H1_nt_pc12maxes
df.reset_index().plot(kind='scatter', x='G', y='pc1',c='k',ax=a)
df.reset_index().plot(kind='scatter', x='G', y='pc2',c='k',ax=a)
a.axvline(x=bifpoint1,linewidth=2, color='k',alpha=0.5)
a.axvline(x=bifpoint2,linewidth=2, color='k',alpha=0.5)
a.set_title(tit)


bifpoint1,bifpoint2 = 0.01,0.02
tit = 'Hag sw nt \n(bif1: %s, bif2: %s' %(bifpoint1,bifpoint2)
a = ax[1]
df = df_Hag_sw_H1_nt_pc12maxes
df.reset_index().plot(kind='scatter', x='G', y='pc1',c='k',ax=a)
df.reset_index().plot(kind='scatter', x='G', y='pc2',c='k',ax=a)
a.axvline(x=bifpoint1,linewidth=2, color='k',alpha=0.5)
a.axvline(x=bifpoint2,linewidth=2, color='k',alpha=0.5)
a.set_title(tit)

bifpoint1,bifpoint2 = 0,0.02
tit = 'HCP sw nt \n(bif1: %s, bif2: %s' %(bifpoint1,bifpoint2)
a = ax[2]
df = df_HCP_sw_H1_nt_pc12maxes
df.reset_index().plot(kind='scatter', x='G', y='pc1',c='k',ax=a)
df.reset_index().plot(kind='scatter', x='G', y='pc2',c='k',ax=a)
a.axvline(x=bifpoint1,linewidth=2, color='k',alpha=0.5)
a.axvline(x=bifpoint2,linewidth=2, color='k',alpha=0.5)
a.set_title(tit)


plt.tight_layout()


plt.savefig(figf, bbox_inches='tight')
plt.close()

Note - can get the middle figure looking more similar to Deco et al. figure by changing zlims...

xlims = [0,0.03]


fig, ax  = plt.subplots(ncols=3,figsize=(12,6))

bifpoint1,bifpoint2 = 0,0.025
tit = 'TVB sw nt \n(bif1: %s, bif2: %s' %(bifpoint1,bifpoint2)
a = ax[0]
df = df_TVB_sw_H1_nt_pc12maxes
df.reset_index().plot(kind='scatter', x='G', y='pc1',c='k',ax=a)
df.reset_index().plot(kind='scatter', x='G', y='pc2',c='k',ax=a)
a.axvline(x=bifpoint1,linewidth=2, color='k',alpha=0.5)
a.axvline(x=bifpoint2,linewidth=2, color='k',alpha=0.5)
a.set_title(tit)
a.set_xlim(xlims)


bifpoint1,bifpoint2 = 0.01,0.02
tit = 'Hag sw nt \n(bif1: %s, bif2: %s' %(bifpoint1,bifpoint2)
a = ax[1]
df = df_Hag_sw_H1_nt_pc12maxes
df.reset_index().plot(kind='scatter', x='G', y='pc1',c='k',ax=a)
df.reset_index().plot(kind='scatter', x='G', y='pc2',c='k',ax=a)
a.axvline(x=bifpoint1,linewidth=2, color='k',alpha=0.5)
a.axvline(x=bifpoint2,linewidth=2, color='k',alpha=0.5)
a.set_title(tit)
a.set_xlim(xlims)

bifpoint1,bifpoint2 = 0,0.02
tit = 'HCP sw nt \n(bif1: %s, bif2: %s' %(bifpoint1,bifpoint2)
a = ax[2]
df = df_HCP_sw_H1_nt_pc12maxes
df.reset_index().plot(kind='scatter', x='G', y='pc1',c='k',ax=a)
df.reset_index().plot(kind='scatter', x='G', y='pc2',c='k',ax=a)
a.axvline(x=bifpoint1,linewidth=2, color='k',alpha=0.5)
a.axvline(x=bifpoint2,linewidth=2, color='k',alpha=0.5)
a.set_title(tit)
a.set_xlim(xlims)

plt.tight_layout()


#plt.savefig(figf, bbox_inches='tight')
plt.close()

KDE plots

figf= outdir+ '/pics/all_bif_kde_plots.png'



fig, ax = plt.subplots(ncols=3,figsize=(12,3))

tit = 'TVB sw S1'
a=ax[0]
df = df_TVB_sw_S1
Gs = np.unique(df.index.get_level_values('G'))
Gstoplot = [Gs[0],Gs[5],Gs[15],
            Gs[25], Gs[35], Gs[45]]
dfs = {}
for g_it,g in enumerate(Gstoplot):
  dfs['%1.4f' %g] = df.ix[g].iloc[1000:,:].max()    
dfs_cat = pd.concat(dfs)
dfs_cat.unstack().T.plot(kind='kde',cmap='coolwarm',ax=a,
                         title=tit)


tit = 'Hag sw S1'
a=ax[1]
df = df_Hag_sw_S1
Gs = np.unique(df.index.get_level_values('G'))
Gstoplot = [Gs[0],Gs[5],Gs[15],
            Gs[25], Gs[35], Gs[45]]
dfs = {}
for g_it,g in enumerate(Gstoplot):
  dfs['%1.4f' %g] = df.ix[g].iloc[1000:,:].max()    
dfs_cat = pd.concat(dfs)
dfs_cat.unstack().T.plot(kind='kde',cmap='coolwarm',ax=a,
                         title=tit)


tit = 'HCP sw S1'
a = ax[2]
df = df_HCP_sw_S1
Gs = np.unique(df.index.get_level_values('G'))
Gstoplot = [Gs[0],Gs[5],Gs[15],
            Gs[25], Gs[35], Gs[45]]
dfs = {}
for g_it,g in enumerate(Gstoplot):
  dfs['%1.4f' %g] = df.ix[g].iloc[1000:,:].max()    
dfs_cat = pd.concat(dfs)
dfs_cat.unstack().T.plot(kind='kde',cmap='coolwarm',ax=a,
                         title=tit)


plt.tight_layout()

plt.savefig(figf, bbox_inches='tight')

plt.close()

Long sims runs

(for FC fits)

Run sims

ndps = {'w': 0.9, 'I_o': 0.33} # non-default params
fixed_params = {'D': 1E-6, 'wwd_params': ndps,'return_SandH': False, 'return_B': True}

Gs = np.linspace(0.0001,0.1,25)
fixed_params['simlen'] = 1000*60*5

fixed_params['conn'] = TVB_conn_sw
df_TVB_sw_S2,df_TVB_sw_H2,df_TVB_sw_B2 = run_wwd_Gsweep(Gs,fixed_params)    

fixed_params['conn'] = Hag_conn_sw
df_Hag_sw_S2,df_Hag_sw_H2,df_Hag_sw_B2 = run_wwd_Gsweep(Gs,fixed_params)  

fixed_params['conn'] = HCP_conn_sw
df_HCP_sw_S2,df_HCP_sw_H2,df_HCP_sw_B2 = run_wwd_Gsweep(Gs,fixed_params)

Compute correlations with SC and FC

zs = np.zeros_like(TVB_SC_sw)
df_TVB_sw_B2_corrs = get_conn_corrs({g: df_TVB_sw_B2.ix[g] for g in Gs}, TVB_SC_sw,zs)
df_Hag_sw_B2_corrs = get_conn_corrs({g: df_Hag_sw_B2.ix[g] for g in Gs}, Hag_SC_sw,Hag_FC)
df_HCP_sw_B2_corrs = get_conn_corrs({g: df_HCP_sw_B2.ix[g] for g in Gs}, HCP_SC_sw,HCP_FC)

Make combined figures

correlation plots

figf = outdir + '/pics/all_FC_fit_plots.png'

fig, ax = plt.subplots(ncols=3,figsize=(12,3))


df_TVB_sw_B2_corrs.T.plot(marker='.',ax=ax[0],ylim=[-0.05,0.5],title='corr BOLD')
df_Hag_sw_B2_corrs.T.plot(marker='.',ax=ax[1],ylim=[-0.05,0.5],title='corr BOLD')
df_HCP_sw_B2_corrs.T.plot(marker='.',ax=ax[2],ylim=[-0.05,0.5],title='corr BOLD')

plt.savefig(figf, bbox_inches='tight')
plt.close()

correlation (R2) plots

figf = outdir + '/pics/all_FC_R2_fit_plots.png'


fig, ax = plt.subplots(ncols=3,figsize=(12,3))

(df_TVB_sw_B2_corrs**2).T.plot(marker='.',ax=ax[0],ylim=[-0.05,0.5],title='corr BOLD')
(df_Hag_sw_B2_corrs**2).T.plot(marker='.',ax=ax[1],ylim=[-0.05,0.5],title='corr BOLD')
(df_HCP_sw_B2_corrs**2).T.plot(marker='.',ax=ax[2],ylim=[-0.05,0.5],title='corr BOLD')



plt.savefig(figf, bbox_inches='tight')

plt.close()

Correlation plots + bifurcation figs

figf = outdir + '/pics/all_bif_scatter_and_FC_plots.png'

fig, ax = plt.subplots(nrows=2,ncols=3, figsize=(12,8))

xlims = [0,0.1]

a = ax[0][0]
df = df_TVB_sw_B2_corrs.T
Gwithmaxcorr = df.idxmax()['corrSC']
tit = 'TVB sw nt \n(maxGcorr: %1.4f' %(Gwithmaxcorr)
df.plot(marker='.',ax=a,ylim=[-0.05,0.5],title=tit)
a.axvline(x=Gwithmaxcorr,linewidth=2, color='k',alpha=0.5)
a.set_xlim(xlims)
a.set_xticklabels(''); a.set_xlabel('')


a = ax[1][0]
df = df_TVB_sw_H1_nt_pc12maxes
df.reset_index().plot(kind='scatter', x='G', y='pc1',c='k',ax=a)
df.reset_index().plot(kind='scatter', x='G', y='pc2',c='k',ax=a)
a.axvline(x=Gwithmaxcorr,linewidth=2, color='k',alpha=0.5)
a.set_xlim(xlims)


a = ax[0][1]
df = df_Hag_sw_B2_corrs.T
Gwithmaxcorr = df.idxmax()['corrFC']
tit = 'Hag sw nt \n(maxGcorr: %1.4f' %(Gwithmaxcorr)
df.plot(marker='.',ax=a,ylim=[-0.05,0.5],title=tit)
a.axvline(x=Gwithmaxcorr,linewidth=2, color='k',alpha=0.5)
a.set_xlim(xlims)
a.set_xticklabels(''); a.set_xlabel('')


a = ax[1][1]
df = df_Hag_sw_H1_nt_pc12maxes
df.reset_index().plot(kind='scatter', x='G', y='pc1',c='k',ax=a)
df.reset_index().plot(kind='scatter', x='G', y='pc2',c='k',ax=a)
a.axvline(x=Gwithmaxcorr,linewidth=2, color='k',alpha=0.5)
a.set_xlim(xlims)



a = ax[0][2]
df = df_HCP_sw_B2_corrs.T
Gwithmaxcorr = df.idxmax()['corrFC']
tit = 'HCP sw nt \n(maxGcorr: %1.4f' %(Gwithmaxcorr)
df.plot(marker='.',ax=a,ylim=[-0.05,0.5],title=tit)
a.axvline(x=Gwithmaxcorr,linewidth=2, color='k',alpha=0.5)
a.set_xlim(xlims)
a.set_xticklabels(''); a.set_xlabel('')


a = ax[1][2]
df = df_HCP_sw_H1_nt_pc12maxes
df.reset_index().plot(kind='scatter', x='G', y='pc1',c='k',ax=a)
df.reset_index().plot(kind='scatter', x='G', y='pc2',c='k',ax=a)
a.axvline(x=Gwithmaxcorr,linewidth=2, color='k',alpha=0.5)
a.set_xlim(xlims);



plt.savefig(figf, bbox_inches='tight')
plt.close()

figf = outdir + '/pics/Hag_SC_FC_conmats.png'


fig, ax = plt.subplots(ncols=5,figsize=(16,3))


df = df_Hag_sw_B2_corrs.T
Gwithmaxcorr = df.idxmax()['corrFC']

tit = 'Hag SC'
a = ax[0]
sns.heatmap(Hag_SC_sw,ax=a, xticklabels='',yticklabels='',
             mask=np.eye(Hag_SC.shape[0]))
a.set_title(tit)


tit = 'Hag FC'
a = ax[1]
sns.heatmap(Hag_FC,ax=a,xticklabels='',yticklabels='',
             mask=np.eye(Hag_SC.shape[0]))
a.set_title(tit)


g = Gs[0]
tit = 'Hag sim FC G=%1.2f' %g
a = ax[2]
sns.heatmap(df_Hag_sw_B2.ix[g].corr(),
            mask=np.eye(Hag_SC.shape[0]),
            vmin=-0.3,vmax=0.3,ax=a,xticklabels='',yticklabels='')
a.set_title(tit)


g = Gwithmaxcorr
tit = 'Hag sim FC (best fit) G=%1.2f' %g
a = ax[3]
df = df_Hag_sw_B2_corrs.T
Gwithmaxcorr = df.idxmax()['corrFC']
sns.heatmap(df_Hag_sw_B2.ix[g].corr(),
            mask=np.eye(Hag_SC.shape[0]),
            vmin=-0.3,vmax=0.3,ax=a,xticklabels='',yticklabels='')
a.set_title(tit)

g = Gs[-1]
tit = 'Hag sim FC G=%1.2f' %g
a = ax[4]
sns.heatmap(df_Hag_sw_B2.ix[g].corr(),
            mask=np.eye(Hag_SC.shape[0]),
            vmin=-0.3,vmax=0.3,ax=a,xticklabels='',yticklabels='')
a.set_title(tit)


plt.tight_layout()


plt.savefig(figf, bbox_inches='tight')
plt.close()

figf = outdir + '/pics/HCP_SC_FC_conmats.png'


fig, ax = plt.subplots(ncols=5,figsize=(16,3))


df = df_HCP_sw_B2_corrs.T
Gwithmaxcorr = df.idxmax()['corrFC']

tit = 'HCP SC'
a = ax[0]
sns.heatmap(HCP_SC_sw,ax=a, xticklabels='',yticklabels='',
             mask=np.eye(HCP_SC.shape[0]))
a.set_title(tit)


tit = 'HCP FC'
a = ax[1]
sns.heatmap(HCP_FC,ax=a,xticklabels='',yticklabels='',
             mask=np.eye(HCP_SC.shape[0]))
a.set_title(tit)


g = Gs[0]
tit = 'HCP sim FC G=%1.2f' %g
a = ax[2]
sns.heatmap(df_HCP_sw_B2.ix[g].corr(),
            mask=np.eye(HCP_SC.shape[0]),
            vmin=-0.3,vmax=0.3,ax=a,xticklabels='',yticklabels='')
a.set_title(tit)


g = Gwithmaxcorr
tit = 'HCP sim FC (best fit) G=%1.2f' %g
a = ax[3]
df = df_HCP_sw_B2_corrs.T
Gwithmaxcorr = df.idxmax()['corrFC']
sns.heatmap(df_HCP_sw_B2.ix[g].corr(),
            mask=np.eye(HCP_SC.shape[0]),
            vmin=-0.3,vmax=0.3,ax=a,xticklabels='',yticklabels='')
a.set_title(tit)

g = Gs[-1]
tit = 'HCP sim FC G=%1.2f' %g
a = ax[4]
sns.heatmap(df_HCP_sw_B2.ix[g].corr(),
            mask=np.eye(HCP_SC.shape[0]),
            vmin=-0.3,vmax=0.3,ax=a,xticklabels='',yticklabels='')
a.set_title(tit)


plt.tight_layout()


plt.savefig(figf, bbox_inches='tight')
plt.close()