diff --git a/examples/eg__griffiths2022.py b/examples/eg__griffiths2022.py new file mode 100644 index 00000000..7be22d0a --- /dev/null +++ b/examples/eg__griffiths2022.py @@ -0,0 +1,1598 @@ + +""" +.. _ex-griffiths2022: + +=================================================================== +DL based parameter estimation for neurophysiological models +=================================================================== + +This example contains the replication results of the paper: Griffiths, J. D., Wang, Z., Ather, S. H., Momi, D., Rich, S., Diaconescu, A., McIntosh, A. R., & Shen, K. Deep Learning-Based Parameter Estimation for Neurophysiological Models of Neuroimaging Data. bioRxiv (2022). DOI: 10.1101/2022.05.19.492664 + sphinx-build -b html -nT --keep-going -D plot_gallery=0 . _build/html +The code in this example allows for the replication of the findings in the paper. If you utilize this code, please cite the original publication. +""" + + +# %% +# 1. Setup +# -------------------------------------------------- +# +# Importage: +# General and whobpyt importage +import os +import sys +import numpy as np +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D +from scipy.optimize import fsolve +import sys +import json +import pandas as pd # for data manipulation +import seaborn as sns # for plotting +import time # for timer +import os +from numpy import exp,sin,cos,sqrt,pi, r_,floor,zeros +import scipy.io +from matplotlib.tri import Triangulation +#from nilearn import plotting, datasets +from matplotlib import cm +from matplotlib.pyplot import subplot +import warnings # for suppressing warnings and output +import nibabel as nib +import matplotlib.pyplot as plt # for plotting +warnings.filterwarnings('ignore') +from whobpyt.datasets.fetchers import fetch_eggriffiths2022 +import os +import numpy as np + +from whobpyt.depr.griffiths2022.plotting import ( + get_rotation_matrix, + get_combined_rotation_matrix, + plot_surface_mpl_mv, + plot_surface_mpl, + plot_sim_states_outputs, + plot_fit_parameters, + R_stat +) +from whobpyt.depr.griffiths2022.bifurcation_analysis import ( + h_tf, + dh_tf, + smooth_normalize, + derivative_orig, + derivative, + get_eig_sys, + regime_search_I0, + regime_search_gEE +) +from whobpyt.depr.griffiths2022.rww_pytorch_model_updated import ( + RNNWWD, + Model_fitting +) + +from whobpyt.depr.griffiths2022.wwd_test import WWD_test +# Load data +eg_dir = fetch_eggriffiths2022() +# Load a file, e.g., weights +f = os.path.join(eg_dir, 'hagmann', 'weights.txt') +W = np.loadtxt(f) + +# %% +# 2. bifurcation_analysis_singlenode +# -------------------------------------------------- +# + +param = { + +# Parameters for the integration +"dt" : 0.001, # integration step size +"ROI_num" : 1, # number of neural nodes + + +# Parameters for the ODEs +# Excitatory population +"W_E" : 1., # scale of the external input +"tau_E" : 100., # decay time +"gamma_E" : 0.641/1000., # other dynamic parameter (?) + +# Inhibitory population +"W_I" : 0.7, # scale of the external input +"tau_I" : 10., # decay time +"gamma_I" : 1./1000., # other dynamic parameter (?) + +# External input +"I_0" : 0.45, # external input +"I_external" : 0., # external stimulation + +# Coupling parameters +"g" : 20., # global coupling (from all nodes E_j to single node E_i) +"g_EE" : .5, # local self excitatory feedback (from E_i to E_i) +"g_IE" : .14, # local inhibitory coupling (from I_i to E_i) +"g_EI" : 0.15, # local excitatory coupling (from E_i to I_i) +"lamb" : 0., # scale of global coupling on I_i (compared to E_i) + +"aE":310, +"bE" :125, +"dE":0.16, +"aI":615, +"bI" :177, +"dI" :0.087, + + + +} + +# initial values of all model prameters + +json_path = os.path.join(eg_dir, 'regime_search_I0_result.json') + +# Load the JSON content +with open(json_path, 'r') as f: + data = json.load(f) + +# 2.1 Now access I0 and c_I0 +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# +I0 = data['I0'] +c_I0 = data['c_I0'] +print('I0 <', I0, 'bifurcation') + +# 2.2 fixed I_0 = 0.2 coastly search g_EE +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# +param_study = ['g_EE', 'g_IE', 'g_EI', 'I_0'] +num_param = 4 +num_trials = 20 +I0_rng = np.linspace(0.2,.25, 1) +gEE_rng = np.linspace(0.4,1, 10) +gIE_rng = np.linspace(0.1,0.5, num_trials) +gEI_rng = np.linspace(0.1,0.5, num_trials) + +c_gEE = regime_search_gEE(I0_rng, gEE_rng, gIE_rng, gEI_rng, param) #CS +gEE = gEE_rng[np.argmax(c_gEE)] +print('gEE >', gEE, 'bistable') + +# 2.3 Figures for the paper +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# +fig = plt.figure(figsize=(6,4)) +ax = fig.add_subplot(111) +im = ax.scatter(np.array(gIE_good), np.array(gEI_good), cmap ='bwr', c = np.array(c), vmax =2, vmin=1, alpha=0.3) +ax.set_xlabel('$g_{IE}$', fontsize=12) +ax.set_ylabel('$g_{EI}$', fontsize=12) +ax.text(1,1.25, 'I', fontsize=15, weight='bold') +ax.text(.1,.15, 'II', fontsize=15, weight='bold') +#fig.colorbar(im, ax = ax) +ax.set_title("Bifurcation analysis on $g_{IE}$ and $g_{EI}$ when $g_{EE} = 0.5$ and $I_0=0.2$") +#fig.savefig('/home/jwang/bifurcationanalysis_gEI_gIE.png') +plt.show() + +fig, ax = plt.subplots(1,1, figsize=(12,8)) +I0= 0.4167 +x1=np.linspace(0,I0- 0.001,50) +x2 =np.linspace(I0,2) +legend_text =['stable fixed points', 'unstable fixed points'] +y1 = np.sqrt(-x1+I0) +y2= - np.sqrt(-x1+I0) +ax.plot(x1, y1, 'r', linewidth=10) +ax.plot(x1, y2, 'r', linewidth=10) +h1, = ax.plot(x1, np.zeros_like(x1), 'r', linewidth=10) +h2, = ax.plot(x2,np.zeros_like(x2), 'g', linewidth=10) +ax.text(I0+0.05, 0.15, '$I_0=$'+str(I0), fontsize = 30) +ax.tick_params(labelsize = 30) +ax.legend([h2,h1], legend_text, fontsize = 20) +plt.show() + +fig, ax = plt.subplots(1,1, figsize=(12,8)) +gEE= 0.4526 +x1=np.linspace(0,gEE,50) +x2 =np.linspace(gEE+0.001,2) +legend_text =['stable fixed points', 'unstable fixed points'] +y1 = np.sqrt(x2-gEE) +y2= - np.sqrt(x2-gEE) +ax.plot(x2, y1, 'g', linewidth=10) +ax.plot(x2, y2, 'g', linewidth=10) +h1, = ax.plot(x1, np.zeros_like(x1), 'g', linewidth=10) +h2, = ax.plot(x2,np.zeros_like(x2), 'r', linewidth=10) +ax.text(gEE+0.05, 0.15, '$g_{EE}=$'+str(gEE), fontsize = 30) +ax.tick_params(labelsize = 30) +ax.legend([h1,h2], legend_text, fontsize = 20) +plt.show() + +# %% +# **Result Description:** +# The bifurcation analysis results show the dynamics of the RWW model that change as different parameters are tested. +# The dynamics of the single node network are shown using the eigenvalues that can be used to represent stability. + + + +# %% +# 3. Synthetic data model-fitting analysis +# -------------------------------------------------- +# + +#** Description:** +#The initial sections of the code handle downloading and importing the required files. Please ensure you manually download the necessary files and update the file paths accordingly. +#Additionally, for compatibility with newer versions of NumPy, if you encounter errors related to numpy.int, please update it to numpy.int64 in rww_pytorch_model.py + +# %% +# 3.1 Load and Sort Simulation Result Files +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# +mask = np.tril_indices(66,-1) +node_size = 66 +output_size = 66 +num_epoches = 1 +batch_size = 20 +step_size = .05 +input_size = 2 +tr = .75 + +G = 100 +gEE = 3.5 +gIE = 0.42 +gEI = 0.42 + +eg_dir = fetch_eggriffiths2022() +sc_file = os.path.join(eg_dir, 'hagmann', 'weights.txt') #CS +sc = np.loadtxt(sc_file) +sc = sc -np.diag(np.diag(sc)) +sc = 0.5*(sc.T+sc) +sc = np.log1p(sc)/np.linalg.norm(np.log1p(sc)) + +# make a fake ts +len_sim = 2400 +ts = np.ones((len_sim, node_size)) + + + +model = RNNWWD(input_size, node_size, batch_size, step_size, tr, sc, False, g_mean_ini=100, g_std_ini = .1, gEE_mean_ini=2.5, gEE_std_ini = .1) + + +# call model fit method +F = Model_fitting(model, ts, num_epoches) + + +output_test = F.test(240) + +fc_sim = np.corrcoef(output_test['simBOLD'][:,2400:]) +## connectivity 66 +fig, ax=plt.subplots(1,2, figsize=(12,8)) +ax[0].imshow(fc_sim - np.diag(np.diag(fc_sim)), cmap='bwr') +ax[1].imshow(sc, cmap='bwr') +plt.show() + +# %% +# 3.2 +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# +out_dir = 'results/' +if not os.path.isdir(out_dir): + os.makedirs(out_dir) + +base_dir = os.path.join(eg_dir, 'hcp') # CS + +subs =sorted([sc_file[-10:-4] for sc_file in os.listdir(base_dir) if sc_file[:8] == 'weights_']) + +start_time = time.time() + +#subs_s = ['562446', '257542', '154936' ] +eg_dir = fetch_eggriffiths2022() + +base_dir = os.path.join(eg_dir, 'hcp') +out_dir = os.path.join(eg_dir, 'hagmann') +for i in range(0,40): + + + node_size = 83 + mask = np.tril_indices(node_size, -1) + num_epoches = 20 + batch_size = 20 + step_size = 0.05 + input_size = 2 + tr = 0.75 + sub=subs[i] + print(i, sub) + sc_file = base_dir +'weights_'+sub+'.txt' + ts_file = base_dir +sub+'_rfMRI_REST1_LR_hpc200_clean__l2k8_sc33_ts.pkl'#out_dir+'sub_'+sub+'simBOLD_idt.txt'# + + if os.path.isfile(sc_file) and os.path.isfile(ts_file): + sc = np.loadtxt(sc_file) + SC =(sc+sc.T)*0.5 + + sc = np.log1p(SC)/np.linalg.norm(np.log1p(SC)) + + ts_pd = pd.read_pickle(ts_file) + ts = ts_pd.values + #ts = np.loadtxt(ts_file) + ts =ts/np.max(ts) + fc_emp = np.corrcoef(ts.T) + # Get the WWD model module for forward in a batch. + + + model = RNNWWD(input_size, node_size, batch_size, step_size, tr, sc, True, g_mean_ini=80, g_std_ini = .1, gEE_mean_ini=2.5, gEE_std_ini = .1) + + + + # call model fit method + F = Model_fitting(model, ts, num_epoches) + + # fit data(train) + + output_train = F.train() + + + output_test = F.test(120) + plot_fit_parameters(output_train) + plot_sim_states_outputs(ts, output_test) + + sc_mod = np.zeros_like(sc) + sc_mod[mask] = output_train['gains'][-100:].mean(0) + sc_mod = sc_mod+sc_mod.T + w = (1 + np.tanh(sc_mod))*sc + w_n = 0.5*(w + w.T)/np.linalg.norm(0.5*(w + w.T)) + np.savetxt(out_dir + 'bold_test_'+ sub +'.txt', output_test['simBOLD']) + np.savetxt(out_dir + 'bold_train_'+ sub +'.txt', output_train['simBOLD']) + np.savetxt(out_dir + 'sc_mod_'+ sub +'.txt', sc_mod) + np.savetxt(out_dir + 'sc_'+ sub +'.txt', sc) + g= output_train['g'][-100:].mean() + gEE = output_train['gEE'][-100:].mean() + gEI = output_train['gEI'][-100:].mean() + gIE = output_train['gIE'][-100:].mean() + + + np.savetxt(out_dir + 'parameters_'+ sub +'.txt', np.array([g,gEE, gIE, gEI]).T) +end_time = time.time() +print('running time is {0} \'s'.format(end_time - start_time )) + +Syn_ts_sim ={} +Syn_ts_test ={} +Syn_paras = {} +Syn_sc_mod = {} +Syn_sc = {} +for i in range(39): + sub = 'sub_'+str(i) + + + Syn_ts_test[sub] = np.loadtxt(os.path.join(hagmann_dir, f'bold_test_{sub}.txt')) + Syn_ts_sim[sub] = np.loadtxt(os.path.join(hagmann_dir, f'bold_train_{sub}.txt')) + Syn_sc_mod[sub] = np.loadtxt(os.path.join(hagmann_dir, f'sc_mod_{sub}.txt')) + Syn_sc[sub] = np.loadtxt(os.path.join(hagmann_dir, f'Hagmann', f'sc_{sub}.txt')) + Syn_paras[sub] = np.loadtxt(os.path.join(hagmann_dir, f'parameters_{sub}.txt')) + +# Save results +np.save(os.path.join(hagmann_dir, 'Syn_ts_sim.npy'), Syn_ts_sim) +np.save(os.path.join(hagmann_dir, 'Syn_ts_test.npy'), Syn_ts_test) +np.save(os.path.join(hagmann_dir, 'Syn_sc_mod.npy'), Syn_sc_mod) +np.save(os.path.join(hagmann_dir, 'Syn_sc.npy'), Syn_sc) +np.save(os.path.join(hagmann_dir, 'Syn_fitparas.npy'), Syn_paras) + +# %% +# 3.3 Plot Post Distributions of the fitted model parameters +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# +mask = np.tril_indices(66,-1) +node_size = 83 +step_size = 0.05 +Tr = 2.5 + +g = 82 +gEE = 2.5 +gIE = 0.4473 +gEI = 0.501 + +#CS here changed the direction of her own, it might be different than others +eg_dir = fetch_eggriffiths2022() +sc_file = os.path.join(eg_dir, 'hcp', 'weights.txt') +sc = np.loadtxt(sc_file) +sc = (sc + sc.T)*0.5 +SC = sc - np.diag(np.diag(sc)) +#SC[SC 0].mean() +1.1*SC[SC>0].std()] =0 # for connectivity 96 only +wo = np.log1p(SC)/np.linalg.norm(np.log1p(SC)) + +test_sim = WWD_test(g, gEE, gIE, gEI, wo, step_size, node_size, Tr) + +bold_test = test_sim.generate_bold(4800) #could be any number you want + +fc_sim = np.corrcoef(bold_test.T) + +## connectivity 83 +fig, ax=plt.subplots(1,2, figsize=(12,8)) +ax[0].imshow(fc_sim - np.diag(np.diag(fc_sim)), cmap='bwr') +ax[1].imshow(wo, cmap='bwr') +plt.show() + +ts_file = "SYN_bold_test_output.txt" # Change the name if needed + +# Save the array to a text file +np.savetxt(ts_file, bold_test) + +print(f"Saved bold_test results to {ts_file}") + +eg_dir = fetch_eggriffiths2022() +base_dir = os.path.join(eg_dir, 'hcp') +ts_file = os.path.join(eg_dir, 'SYN_bold_test_output.txt') # make sure this file exists +ts = np.loadtxt(ts_file) +print(ts.shape) + +# Set the structural connectivity file +sc_file = os.path.join(base_dir, 'weights.txt') + +sc = np.loadtxt(sc_file) +sc = (sc + sc.T)*0.5 +SC = sc - np.diag(np.diag(sc)) +#SC[SC 0].mean() + 1.*SC[SC>0].std()] = 0 +sc_true = np.log1p(SC)/np.linalg.norm(np.log1p(SC)) #SC/np.linalg.norm(SC) +fc_true = np.corrcoef(ts[2400:3600,:].T) +fc_cross = np.corrcoef(ts[0:2400,:].T) + +fig, ax = plt.subplots(1,1, figsize=(12,8)) +im = ax.imshow(sc_true, cmap='bwr') +mn = round(sc_true.min(),2) +mx = round(sc_true.max(), 2) +md = round((mx + mn)/2, 2) +mnd = round((md + mn)/2, 2) +mxd = round((mx + md)/2, 2) +cbar = plt.colorbar(im, ax=ax) +cbar.ax.tick_params(labelsize=30) +cbar.set_ticks([mn,mnd, md,mxd, mx]) +cbar.set_ticklabels([mn,mnd, md,mxd, mx]) +ax.set_xticklabels([]) +ax.set_yticklabels([]) +plt.show() +# plt.savefig(base_dir+'Syn_sc_true.png') + +fig, ax = plt.subplots(1,1, figsize=(12,8)) +fc_nd = fc_true - np.diag(np.diag(fc_true)) +im = ax.imshow(fc_nd, cmap='bwr') +mn = round(fc_nd.min(),2) +mx = round(fc_nd.max(), 2) +md = round((mx + mn)/2, 2) +mnd = round((md + mn)/2, 2) +mxd = round((mx + md)/2, 2) +cbar = plt.colorbar(im, ax=ax) +cbar.ax.tick_params(labelsize=30) +cbar.set_ticks([mn,mnd, md,mxd, mx]) +cbar.set_ticklabels([mn,mnd, md,mxd, mx]) +ax.set_xticklabels([]) +ax.set_yticklabels([]) +plt.show() +# plt.savefig(base_dir+'Syn_fc_true.png') + +# %% +# Read the model-fitting data with gains +eg_dir = fetch_eggriffiths2022() +base_dir = os.path.join(eg_dir, 'hagmann') +ts_sim_file = base_dir + 'Syn_ts_sim.npy' +ts_test_file = base_dir + 'Syn_ts_test.npy' +sc_mod_file = base_dir + 'Syn_sc_mod.npy' +sc_file = base_dir + 'Syn_sc.npy' +paras_file = base_dir + 'Syn_fitparas.npy' + +Syn_ts_sim_data = np.load(ts_sim_file, allow_pickle=True) +Syn_ts_test_data = np.load(ts_test_file, allow_pickle=True) +Syn_sc_mod_data = np.load(sc_mod_file, allow_pickle=True) +Syn_sc_data = np.load(sc_file, allow_pickle=True) +Syn_para_data = np.load(paras_file, allow_pickle=True) + +Syn_ts_test = Syn_ts_test_data.item() +Syn_ts_sim = Syn_ts_sim_data.item() +Syn_sc_mod = Syn_sc_mod_data.item() +Syn_sc = Syn_sc_data.item() +Syn_paras = Syn_para_data.item() + +start_time = time.time() + +plt.rcParams["axes.labelsize"] = 30 +mask = np.tril_indices(66,-1) +data_dict = {} +data_max = {} +data_dict['g'] = [] +data_dict['gEE'] = [] +data_dict['gIE'] = [] +data_dict['gEI'] = [] + +data_max['g'] = [] +data_max['gEE'] = [] +data_max['gIE'] = [] +data_max['gEI'] = [] + + +# Across the test subjects, take the maximum of the absolute value of the mean of the parameter +# values identified for the model and the same for the parameter values. Append the two values and +# add them to the `data_max` dictionary for the coupling strength parameters. + + +para_name = ['g', + 'gEE', + 'gIE', + 'gEI'] +for sub in Syn_paras: + + data_max['g'].append(max(np.abs(Syn_paras[sub][0].mean()), 100)) + data_max['gEE'].append(max(np.abs(Syn_paras[sub][1].mean()), 3.5)) + data_max['gIE'].append(max(np.abs(Syn_paras[sub][2].mean()), 0.42)) + data_max['gEI'].append(max(np.abs(Syn_paras[sub][3].mean()), 0.42)) + + +# Take these parameter values and use these maximum values to normalize the range of parameter values. + + +for sub in Syn_paras: + + #data_dict['g'].append((Syn_paras[sub][0].mean() - 100)/max(data_max['g'])) + data_dict['gEE'].append((Syn_paras[sub][1].mean() - 23)/max(data_max['gEE'])) #CS + data_dict['gIE'].append((Syn_paras[sub][2].mean() - 18)/max(data_max['gIE']))#CS + data_dict['gEI'].append((Syn_paras[sub][3].mean() - 17)/max(data_max['gEI']))#CS + data_dict['g'].append((Syn_paras[sub][0].mean() - 20)/max(data_max['g'])) + +data_dict['g'] = np.array(data_dict['g']) +data_dict['gEE'] = np.array(data_dict['gEE']) +data_dict['gIE'] = np.array(data_dict['gIE']) +data_dict['gEI'] = np.array(data_dict['gEI']) + +# %% +# Plot the results. + +import seaborn as sns +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +# Ensure ax is properly indexed +fig, ax = plt.subplots(2, 2, figsize=(12, 8), sharex=True) +ax = np.array(ax) # Ensure proper indexing +plt.rcParams["axes.labelsize"] = 30 + +# Verify data range and adjust if needed +for i, key in enumerate(para_name): + data = np.array(data_dict[key]) + + if len(data) > 1: # KDE requires at least two points + print(f"Plotting {key}: min={data.min()}, max={data.max()}, std={data.std()}") # Debug info + + # Use a larger bandwidth for smoother KDE + sns.kdeplot(ax=ax[i // 2, i % 2], data=data, color='r', linewidth=3, bw_adjust=0.5) + + # Fallback: if KDE fails, show histogram + if not np.any(data): # Check if data is effectively empty + print(f"Warning: {key} data has no variance, switching to histogram.") + ax[i // 2, i % 2].hist(data, bins=20, color='r', alpha=0.7) + + ax[i // 2, i % 2].set_xlim(-0.3, 0.5) # Adjust x limits for visibility + ax[i // 2, i % 2].tick_params(labelsize=20) + ax[i // 2, i % 2].grid() + else: + print(f"Skipping {key} due to insufficient data for KDE.") + +fig.tight_layout() +plt.show() + +# Second Plot +fig, ax = plt.subplots(figsize=(12, 8)) +plt.rcParams["axes.labelsize"] = 25 + +if any(len(data_dict[key]) > 1 for key in data_dict): + sns.kdeplot(data=pd.DataFrame(data_dict), linewidth=4, bw_adjust=0.5) + ax.tick_params(labelsize=20) + ax.set_xlim(-1, 1) # Adjust x limits for visibility + ax.set_ylim(0, 20) # Adjust x limits for visibility + ax.grid() + plt.setp(ax.get_legend().get_texts(), fontsize='20') # Reduce fontsize for better visibility + plt.setp(ax.get_legend().get_title(), fontsize='25') + plt.show() + +# %% +# plot hist of r: model fitting with gains +start_time = time.time() +mask = np.tril_indices(66,-1) +corr_sim ={} +r =[] +r1 =[] +for sub in Syn_ts_sim: + r.append(np.corrcoef(fc_true[mask], np.corrcoef(Syn_ts_sim[sub])[mask])[0,1]) + r1.append(np.corrcoef(fc_true[mask], np.corrcoef(Syn_ts_test[sub])[mask])[0,1]) +corr_sim['train'] = np.array(r) +corr_sim['test'] = np.array(r1) +fig, ax = plt.subplots(1,1, figsize=(12,8)) +plt.rcParams["axes.labelsize"] = 35 +sns.kdeplot(data = corr_sim, linewidth = 10) +#ax.set_title('hist of Pearlson correlation between the fitting and empirical BOLD', fontsize =20) +ax.tick_params(labelsize= 25) +ax.grid() + +plt.setp(ax.get_legend().get_texts(), fontsize='35') +plt.setp(ax.get_legend().get_title(), fontsize ='50') + +# %% +# plot hist of r: model fitting with gains +mask = np.tril_indices(66,-1) +corr_sim ={} +r =[] +r1 =[] +for sub in Syn_ts_sim: + r.append(np.corrcoef(fc_true[mask], np.corrcoef(Syn_ts_sim[sub][:,10:])[mask])[0,1]) + r1.append(np.corrcoef(fc_cross[mask], np.corrcoef(Syn_ts_test[sub])[mask])[0,1]) +corr_sim['train'] = np.array(r) +corr_sim['test'] = np.array(r1) +fig, ax = plt.subplots(1,1, figsize=(12,8)) +plt.rcParams["axes.labelsize"] = 35 +sns.kdeplot(data = corr_sim, linewidth = 10) +ax.tick_params(labelsize= 25) +ax.grid() + +plt.setp(ax.get_legend().get_texts(), fontsize='35') +plt.setp(ax.get_legend().get_title(), fontsize ='50') + + +# plot hist of r: model fitting with gains +# %% +mask = np.tril_indices(66,-1) +corr_sim ={} +r =[] +r1 =[] +for sub in Syn_ts_sim: + r.append(np.corrcoef(fc_true[mask], np.corrcoef(Syn_ts_sim[sub][:,10:])[mask])[0,1]) + r1.append(np.corrcoef(fc_true[mask], np.corrcoef(Syn_ts_test[sub])[mask])[0,1]) +corr_sim['train'] = np.array(r) +corr_sim['test'] = np.array(r1) +fig, ax = plt.subplots(1,1, figsize=(12,8)) +plt.rcParams["axes.labelsize"] = 35 +sns.kdeplot(data = corr_sim, linewidth = 10) + +ax.tick_params(labelsize= 25) +ax.grid() + +plt.setp(ax.get_legend().get_texts(), fontsize='35') +plt.setp(ax.get_legend().get_title(), fontsize ='50') +plt.show() + +start_time = time.time() +mask = np.tril_indices(66,-1) +subs = ['randSClmtd', + 'diagrandSClmtd', + 'randSC', + 'dist'] + +titles = ['RandNonzeros', + 'RandNonzeroes+diag', + 'Rand', + 'DistRecp'] + +sc_diff_file = base_dir + 'Syn_diff_sc.npy' +fc_diff_sim_file = base_dir + 'Syn_diff_fc_sim.npy' +sc_diff_mod_file = base_dir + 'Syn_diff_sc_mod.npy' + +Syn_diff_sc_data = np.load(sc_diff_file, allow_pickle=True) +Syn_diff_fc_sim_data = np.load(fc_diff_sim_file, allow_pickle=True) +Syn_diff_sc_mod_data = np.load(sc_diff_mod_file, allow_pickle=True) +Syn_diff_sc = Syn_diff_sc_data.item() +Syn_diff_fc_sim = Syn_diff_fc_sim_data.item() + +Syn_diff_sc_mod = Syn_diff_sc_mod_data.item() +vmax_ls = [0,0,0] +vmin_ls = [1,1,1] +delta = np.array([0,0,0,0,-0.001]) +i = 0 +for sub in Syn_diff_sc_mod: + if i < 4: + + if Syn_diff_sc[sub].max() > vmax_ls[0]: + vmax_ls[0] = Syn_diff_sc[sub].max() + if Syn_diff_sc_mod[sub].max() > vmax_ls[1]: + vmax_ls[1] = Syn_diff_sc_mod[sub].max() + if (Syn_diff_fc_sim[sub]- np.diag(np.diag(Syn_diff_fc_sim[sub]))).max() > vmax_ls[2]: + vmax_ls[2] = (Syn_diff_fc_sim[sub]- np.diag(np.diag(Syn_diff_fc_sim[sub]))).max() + + if Syn_diff_sc[sub].min() < vmin_ls[0]: + vmin_ls[0] = Syn_diff_sc[sub].min() + if Syn_diff_sc_mod[sub].min() < vmin_ls[1]: + vmin_ls[1] = Syn_diff_sc_mod[sub].min() + if (Syn_diff_fc_sim[sub]- np.diag(np.diag(Syn_diff_fc_sim[sub]))).min() < vmin_ls[2]: + vmin_ls[2] = (Syn_diff_fc_sim[sub]- np.diag(np.diag(Syn_diff_fc_sim[sub]))).min() + i += 1 +fig, ax = plt.subplots(4, 3, figsize=(24,24)) +i = 0 +for sub in Syn_diff_sc_mod: + if i < 4: + + im0 = ax[i,0].imshow(Syn_diff_sc[sub], cmap='bwr', vmax = vmax_ls[0], vmin=vmin_ls[0]) + im1 = ax[i,1].imshow(Syn_diff_sc_mod[sub], cmap='bwr', vmax = vmax_ls[1], vmin=vmin_ls[1]) + im2 = ax[i,2].imshow(Syn_diff_fc_sim[sub] - np.diag(np.diag(Syn_diff_fc_sim[sub])), cmap='bwr', vmax = vmax_ls[2], vmin=vmin_ls[2]) + + #ax[i,0].set_title(titles[i] +': initial SC', fontsize= 20) + ax[i,0].tick_params(labelsize= 15) + #ax[i,0].tick_params(axis='y', labelsize= 15) + #ax[i,1].set_title(titles[i] +': modified SC', fontsize= 20) + ax[i,1].tick_params(labelsize= 15) + #ax[i,1].tick_params(axis='y', labelsize= 15) + #ax[i,2].set_title(sub +': sim FC', fontsize= 20) + ax[i,2].tick_params(labelsize= 15)# + #ax[i,2].tick_params(axis='y', labelsize= 15) + ax[i,0].set_xticklabels([]) + ax[i,0].set_yticklabels([]) + ax[i,1].set_xticklabels([]) + ax[i,1].set_yticklabels([]) + ax[i,2].set_xticklabels([]) + ax[i,2].set_yticklabels([]) + i += 1 + else: + break +cbar0 = fig.colorbar(im0, ax = ax.T[0], location='bottom', aspect=10, fraction = 0.046, pad = 0.02) +cbar0.set_ticks(delta+np.round(np.linspace(vmin_ls[0], vmax_ls[0], 5), 3)) +cbar0.set_ticklabels(delta+ np.round(np.linspace(vmin_ls[0], vmax_ls[0], 5),3)) + +cbar1 = fig.colorbar(im1, ax = ax.T[1], location='bottom', aspect=10,fraction = 0.046, pad = 0.02) +cbar1.set_ticks(delta+np.round(np.linspace(vmin_ls[1], vmax_ls[1], 5), 3)) +cbar1.set_ticklabels(delta+np.round(np.linspace(vmin_ls[1], vmax_ls[1], 5), 3)) + +cbar2 = fig.colorbar(im2, ax = ax.T[2], location='bottom', aspect=10,fraction = 0.046, pad = 0.02) +cbar2.set_ticks(delta+np.round(np.linspace(vmin_ls[2], vmax_ls[2], 5), 3)) +cbar2.set_ticklabels(delta+np.round(np.linspace(vmin_ls[2], vmax_ls[2], 5),3)) +plt.show() + +# %% + +fig, ax = plt.subplots(2,2, figsize=(40,28)) + +mask = np.tril_indices(66,-1) +data_dict = {} +data_dict['g'] = [] +data_dict['gEE'] = [] +data_dict['gIE'] = [] +data_dict['gEI'] = [] + +para_name = ['g', + 'gEE', + 'gIE', + 'gEI'] + +for sub in Syn_paras: + + data_dict['g'].append(Syn_paras[sub][0].mean()) + data_dict['gEE'].append(Syn_paras[sub][1].mean()) + data_dict['gIE'].append(Syn_paras[sub][2].mean()) + data_dict['gEI'].append(Syn_paras[sub][3].mean()) + +data_dict['g'] = np.array(data_dict['g']) +data_dict['gEE'] = np.array(data_dict['gEE']) +data_dict['gIE'] = np.array(data_dict['gIE']) +data_dict['gEI'] = np.array(data_dict['gEI']) + +for i in range(4): + sns.histplot(ax=ax[i//2, i%2], data=data_dict[para_name[i]], kde=True, line_kws={'linewidth':10}) + ax[i//2, i%2].set_title('hist of ' + para_name[i] +' across subjects', fontsize='40') + ax[i//2, i%2].tick_params(labelsize=40) + ax[i//2, i%2].set_xlabel(para_name[i], fontsize=40) + +plt.show() + + +# %% +sc_m = Syn_sc_mod[sub] + +fig, ax = plt.subplots(1,3, figsize=(12,8), sharex=True, sharey=True) +ax[0].scatter(fc_true[mask], sc_m[mask], c='r') +#ax[0].set_title('The modified SC against FC in elements', fontsize=20) +ax[0].tick_params(labelsize = 20) +ax[0].set_xlabel('FC weights', fontsize = 20) +ax[0].set_ylabel('Modified SC weights', fontsize = 20) +ax[0].set_ylim(0,0.18) +ax[1].scatter(fc_true[mask], sc_true[mask]) +#ax[1].set_title('The SC against FC in elements', fontsize=20) +ax[1].set_ylim(0,0.18) +ax[1].tick_params(labelsize = 20) +ax[1].set_xlabel('FC weights', fontsize = 20) +ax[1].set_ylabel('SC weights', fontsize = 20) +ax[2].scatter(fc_true[mask], Syn_sc[sub][mask]) +#ax[2].set_title('The initial SC against FC in elements', fontsize=20) +ax[2].set_ylim(0,0.18) +ax[2].tick_params(labelsize = 20) +ax[2].set_xlabel('FC weights', fontsize = 20) +ax[2].set_ylabel('SC initial weights', fontsize = 20) +ax[0].grid() +ax[1].grid() +ax[2].grid() + +fig.tight_layout() +plt.show() + +# %% +sub = 'sub_25' + +sc = Syn_sc[sub] +sc_m = Syn_sc_mod[sub] +fig, ax = plt.subplots(1,1, figsize=(12,8)) +ax.scatter(fc_true[mask], np.corrcoef(Syn_ts_sim[sub])[mask], c='r') +#ax.set_title('The simulated FC against FC in elements',fontsize=20) +ax.tick_params(labelsize=30) +ax.set_xlabel('FC weights', fontsize=30) +ax.set_ylabel('simFC weights', fontsize=30) +ax.grid() + +# %% +## CS add in +import seaborn as sns +import matplotlib.pyplot as plt +import numpy as np + +sub = 'sub_25' + +# Extract data +empirical_fc = fc_true[mask] # (Synthetic) empirical FC values +simulated_fc = np.corrcoef(Syn_ts_sim[sub])[mask] # Fitted FC values + +# Create histogram plot +fig, ax = plt.subplots(figsize=(12, 8)) + +# Plot empirical FC in blue +sns.histplot(empirical_fc, bins=30, color='b', label="Empirical FC", alpha=0.6) + +# Plot simulated FC in red +sns.histplot(simulated_fc, bins=30, color='r', label="Fitted FC", alpha=0.6) + +# Customize the plot +ax.tick_params(labelsize=25) +ax.set_xlabel("FC weights", fontsize=30) +ax.set_ylabel("Count", fontsize=30) +ax.legend(fontsize=20) +ax.grid() + +# Show the plot +plt.show() + +# %% +# 4. Model fitting and training +# -------------------------------------------------- + +# 4.1: Analysis of the HCP data model-fitting +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# + +#Approximate run duration: 10-20 seconds + +#The model-fitting methods on the data from the Human Connectome Project (HCP) let us view how well the parameters fit the model. +#The HCP (a.k.a."original" or main HCP, HCP Young Adult, HCP-YA) maps the healthy human connectome by collecting and freely distributing neuroimaging and behavioral data on 1,200 normal young adults, aged 22-35. This notebook contains the post-analysis in which we determine how well the parameters fit the model using the model-fitting method. +#After fitting the model to the empirical BOLD signal, one may get a set of model parameters that fit the model well. We can use the forward model and the fitted model parameters to generate simulated BOLD data, fit the model again to the simulated BOLD signal with the fitted model parameters (true values), and the new model-fitted data are saved at the files with the `idt` suffix. +#In this analysis, we demonstrate how well they fit. +eg_dir = fetch_eggriffiths2022() + +# Set base_dir to the 'hcp/Outputs' subfolder +base_dir = os.path.join(eg_dir, 'hcp', 'Outputs') + +# Define file paths +fc_file = os.path.join(base_dir, 'HCP_ts.npy') # functional connectivity +sc_file = os.path.join(base_dir, 'HCP_sc.npy') # structural connectivity +fc_sim_file = os.path.join(base_dir, 'HCP_ts_sim.npy') # simulated FC +fc_sim_idt_file = os.path.join(base_dir, 'HCP_ts_sim_idt.npy') # simulated FC (identified params) +sc_mod_file = os.path.join(base_dir, 'HCP_sc_mod.npy') # modified SC +sc_mod_idt_file = os.path.join(base_dir, 'HCP_sc_mod_idt.npy') # modified SC (identified) +paras_file = os.path.join(base_dir, 'HCP_fitparas.npy') # parameters +paras_idt_file = os.path.join(base_dir, 'HCP_fitparas_idt.npy') # parameters (identified) +fc_test_file = os.path.join(base_dir, 'HCP_ts_test.npy') # FC test +fc_test_idt_file = os.path.join(base_dir, 'HCP_ts_test_idt.npy')# FC test (identified) + +HCP_fc_data = np.load(fc_file, allow_pickle=True) +HCP_sc_data = np.load(sc_file, allow_pickle=True) +HCP_fc_sim_data = np.load(fc_sim_file, allow_pickle=True) +HCP_fc_sim_idt_data = np.load(fc_sim_idt_file, allow_pickle=True) +HCP_sc_mod_data = np.load(sc_mod_file, allow_pickle=True) +HCP_sc_mod_idt_data = np.load(sc_mod_idt_file, allow_pickle=True) +HCP_para_data = np.load(paras_file, allow_pickle=True) +HCP_para_idt_data = np.load(paras_idt_file, allow_pickle=True) +HCP_fc_sim_data = np.load(fc_sim_file, allow_pickle=True) +HCP_fc_test_data = np.load(fc_test_file, allow_pickle=True) +HCP_fc_test_idt_data = np.load(fc_test_idt_file, allow_pickle=True) + +HCP_ts = HCP_fc_data.item() +HCP_sc = HCP_sc_data.item() +HCP_ts_sim = HCP_fc_sim_data.item() +HCP_ts_sim_idt = HCP_fc_sim_idt_data.item() +HCP_sc_mod = HCP_sc_mod_data.item() +HCP_sc_mod_idt = HCP_sc_mod_idt_data.item() +HCP_para = HCP_para_data.item() +HCP_para_idt = HCP_para_idt_data.item() +HCP_ts_test = HCP_fc_test_data.item() +HCP_ts_test_idt = HCP_fc_test_idt_data.item() + +i =0 +for sub in HCP_para: + if i < 43: + print(sub) + fig, ax = plt.subplots(1, 4, figsize=(12,8)) + im0 = ax[0].imshow(HCP_sc[sub], cmap='bwr') + sc_mod = HCP_sc_mod[sub] + sc = HCP_sc[sub] + w = (1 + np.tanh(sc_mod))*sc + w_n = 0.5*(w + w.T)/np.linalg.norm(0.5*(w + w.T)) + im1 = ax[1].imshow(w_n, cmap='bwr') + im2 = ax[2].imshow(np.corrcoef((HCP_ts[sub]).T)-np.diag(np.diag(np.corrcoef((HCP_ts[sub]).T))), cmap='bwr', vmax =1) #-HCP_ts[sub].mean(1))- np.diag(np.diag(np.corrcoef(HCP_ts[sub].T-HCP_ts[sub].mean(1)))) + im3 = ax[3].imshow(np.corrcoef(HCP_ts_test[sub])-np.diag(np.diag(np.corrcoef(HCP_ts_test[sub]))), cmap='bwr',vmax = 1) #-HCP_ts_test[sub].mean(0))- np.diag(np.diag(np.corrcoef(HCP_ts_test[sub]-HCP_ts_test[sub].mean(0)))) + plt.show() + i +=1 + +# %% +# 4.2 Plot the connectivity weight matrices +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# +fig, ax = plt.subplots(3, 4, figsize=(20,16)) + +vmax_ls = [0,-50,0, 0] +vmin_ls = [1,50,1, 1] +delta= np.array([0,0,0,0,-0.001]) +i = 0 +#for sub in ['562446', '257542', '154936']: +for sub in ['150423', '257542', '154936']: # CS choose + if i < 4: + if HCP_sc[sub].max() > vmax_ls[0]: + vmax_ls[0] = HCP_sc[sub].max() + if HCP_sc_mod[sub].max() > vmax_ls[1]: + + vmax_ls[1] = HCP_sc_mod[sub].max() + if (np.corrcoef(HCP_ts[sub].T)- np.diag(np.diag(np.corrcoef(HCP_ts[sub].T)))).max() > vmax_ls[2]: + vmax_ls[2] = (np.corrcoef(HCP_ts[sub].T)- np.diag(np.diag(np.corrcoef(HCP_ts[sub].T)))).max() + if (np.corrcoef(HCP_ts_test[sub])- np.diag(np.diag(np.corrcoef(HCP_ts_test[sub])))).max() > vmax_ls[3]: + vmax_ls[3] = (np.corrcoef(HCP_ts_test[sub])- np.diag(np.diag(np.corrcoef(HCP_ts_test[sub])))).max() + if HCP_sc[sub].min() < vmin_ls[0]: + vmin_ls[0] = HCP_sc[sub].min() + if HCP_sc_mod[sub].min() < vmin_ls[1]: + + vmin_ls[1] = HCP_sc_mod[sub].min() + if (np.corrcoef(HCP_ts[sub].T)- np.diag(np.diag(np.corrcoef(HCP_ts[sub].T)))).min() < vmin_ls[2]: + vmin_ls[2] = (np.corrcoef(HCP_ts[sub].T)- np.diag(np.diag(np.corrcoef(HCP_ts[sub].T)))).min() + if(np.corrcoef(HCP_ts_test[sub])- np.diag(np.diag(np.corrcoef(HCP_ts_test[sub])))).min() < vmin_ls[3]: + vmin_ls[3] = (np.corrcoef(HCP_ts_test[sub])- np.diag(np.diag(np.corrcoef(HCP_ts_test[sub])))).min() + i += 1 +i_st = 0 +i = 0 +delta = np.array([0.001,0,0,0,0.0]) +for sub in ['150423', '257542', '154936']: # CS choose +#for sub in ['562446', '257542', '154936'] : + if (i > i_st-1 ) and (i < i_st + 3): + im0 = ax[i-i_st,0].imshow(HCP_sc[sub], cmap='bwr') + sc_mod = HCP_sc_mod[sub] + sc = HCP_sc[sub] + + im1 = ax[i-i_st,1].imshow(sc_mod, cmap='bwr') + im2 = ax[i-i_st,2].imshow(np.corrcoef((HCP_ts[sub]).T)-np.diag(np.diag(np.corrcoef((HCP_ts[sub]).T))), cmap='bwr', vmax = 0.8, vmin=-0.2) #-HCP_ts_test[sub].mean(0))- np.diag(np.diag(np.corrcoef(HCP_ts_test[sub]-HCP_ts_test[sub].mean(0)))) + im3 = ax[i-i_st,3].imshow(np.corrcoef(HCP_ts_test[sub])-np.diag(np.diag(np.corrcoef(HCP_ts_test[sub]))), cmap='bwr', vmax= 0.8, vmin=-0.2)#-HCP_ts_test[sub].mean(0))- np.diag(np.diag(np.corrcoef(HCP_ts_test[sub]-HCP_ts_test[sub].mean(0))) + #ax[i-i_st,0].set_title('s'+sub +': empirical SC', fontsize=20) + ax[i-i_st,0].set_xticklabels([]) + ax[i-i_st,0].set_yticklabels([]) + #ax[i-i_st,1].set_title('s'+sub +': modified SC', fontsize=20) + ax[i-i_st,1].set_xticklabels([]) + ax[i-i_st,1].set_yticklabels([]) + #ax[i,0].tick_params(axis='y', labelsize=15) + #ax[i-i_st,2].set_title('s'+sub +': empirical FC', fontsize=20) + ax[i-i_st,2].set_xticklabels([]) + ax[i-i_st,2].set_yticklabels([]) + #ax[i,1].tick_params(axis='y', labelsize=15) + #ax[i-i_st,3].set_title('s'+sub +' simulated FC', fontsize=20) + ax[i-i_st,3].set_xticklabels([]) + ax[i-i_st,3].set_yticklabels([]) + #ax[i,2].tick_params(axis='y', labelsize= 15) + i += 1 + elif i < i_st: + i += 1 + else: + break + +cbar0 = fig.colorbar(im0, + ax=ax.T[0], + location='bottom', + aspect=10, + fraction=0.046, + pad=0.02) + +cbar1 = fig.colorbar(im1, + ax=ax.T[1], + location='bottom', + aspect=10, + fraction=0.046, + pad=0.02) + +cbar2 = fig.colorbar(im2, + ax=ax.T[2], + location='bottom', + aspect=10, + fraction=0.046, + pad=0.02) + +cbar3 = fig.colorbar(im3, + ax=ax.T[3], + location='bottom', + aspect=10, + fraction=0.046, + pad=0.02) + +cbar0.set_ticks(0.001*np.floor(1000*np.linspace(vmin_ls[0], vmax_ls[0], 5))) + +cbar1.set_ticks(0.001*np.floor(1000*np.linspace(vmin_ls[1], vmax_ls[1], 5))) +cbar2.set_ticks(0.001*np.floor(1000*np.linspace(-0.2, 0.8, 5))) +cbar3.set_ticks(0.001*np.floor(1000*np.linspace(-0.2, 0.8, 5))) +plt.show() +# fig.savefig(base_dir+'HCP_sc_fc_bestfit_3subjects.png') + +# %% +# 4.3 Model parameter fitting estimation +# ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +# +start_time = time.time() +mask = np.tril_indices(83,-1) + +data_dict = {} +data_max = {} +data_dict['g'] = [] +data_dict['gEE'] = [] +data_dict['gIE'] = [] +data_dict['gEI'] = [] + +data_max['g'] = [] +data_max['gEE'] = [] +data_max['gIE'] = [] +data_max['gEI'] = [] + +para_name = ['g', 'gEE', 'gIE', 'gEI'] + +""" +Across the HCP test subjects, take the maximum of the absolute value of the mean of the parameter +values identified for the model and the same for the parameter values. Append the two values and +add them to the `data_max` dictionary for the coupling strength parameters. +""" + +for sub in HCP_para: + + data_max['g'].append(max(np.abs(HCP_para_idt[sub][0].mean()), np.abs(HCP_para[sub][0].mean()))) + data_max['gEE'].append(max(np.abs(HCP_para_idt[sub][1].mean()), np.abs(HCP_para[sub][1].mean()))) + data_max['gIE'].append(max(np.abs(HCP_para_idt[sub][2].mean()), np.abs(HCP_para[sub][2].mean()))) + data_max['gEI'].append(max(np.abs(HCP_para_idt[sub][3].mean()), np.abs(HCP_para[sub][3].mean()))) + +""" +Take these parameter values and use these maximum values to normalize the range of parameter values. +""" + +for sub in HCP_para: + + data_dict['g'].append((HCP_para_idt[sub][0].mean() - HCP_para[sub][0].mean())/max(data_max['g'])) + data_dict['gEE'].append((HCP_para_idt[sub][1].mean() - HCP_para[sub][1].mean())/max(data_max['gEE'])) + data_dict['gIE'].append((HCP_para_idt[sub][2].mean() - HCP_para[sub][2].mean())/max(data_max['gIE'])) + data_dict['gEI'].append((HCP_para_idt[sub][3].mean() - HCP_para[sub][3].mean())/max(data_max['gEI'])) + +data_dict['g'] = np.array(data_dict['g']) +data_dict['gEE'] = np.array(data_dict['gEE']) +data_dict['gIE'] = np.array(data_dict['gIE']) +data_dict['gEI'] = np.array(data_dict['gEI']) +# %% +fig, ax = plt.subplots(2,2, figsize=(12,8), sharex=True) +plt.rcParams["axes.labelsize"] = 40 +for i in range(4): + sns.kdeplot(ax=ax[i//2, i%2], data=data_dict[para_name[i]], color='r', linewidth=10) + #ax[i//2, i%2].set_title('hist of error of ' + para_name[i] +' across subjects', fontsize = 20) + ax[i//2, i%2].set_xlim(-0.5,0.5) + ax[i//2, i%2].tick_params(labelsize = 30) + ax[i//2, i%2].grid() + if i%2 == 1: + ax[i//2, i%2].set_ylabel('') +fig.tight_layout() +plt.show() +# fig.savefig(base_dir+'HCP_hist_modelparaserror_normalized_pt.png') + +fig, ax = plt.subplots(1,1, figsize=(12,8)) +plt.rcParams["axes.labelsize"] = 25 +sns.kdeplot(data=pd.DataFrame(data_dict), linewidth=4) +#ax.set_title('hist of error of parameters across subjects', fontsize =25) +ax.tick_params(labelsize=20) +ax.set_xlim(-0.2,0.2) +ax.set_ylim(0,50) +ax.grid() +plt.setp(ax.get_legend().get_texts(), fontsize='35') +plt.setp(ax.get_legend().get_title(), fontsize='50') +plt.show() +# fig.savefig(base_dir+'HCP_hist_modelparaserror_normalized_in1fig_pt.png') +fig, ax = plt.subplots(2,2, figsize=(12,8)) +plt.rcParams["axes.labelsize"] = 40 +mask = np.tril_indices(83,-1) +data_dict = {} +data_dict['g'] = [] +data_dict['gEE'] = [] +data_dict['gIE'] = [] +data_dict['gEI'] = [] +data_dict['cA'] = [] +data_dict['cB'] = [] +data_dict['cC'] = [] +para_name = ['g', 'gEE', 'gIE', 'gEI'] #, 'cA', 'cB', 'cC'] + +""" +Get the mean of the parameters in the dictionary of data. +""" + +for sub in HCP_para: + data_dict['g'].append(HCP_para[sub][0].mean()) + data_dict['gEE'].append(HCP_para[sub][1].mean()) + data_dict['gIE'].append(HCP_para[sub][2].mean()) + data_dict['gEI'].append(HCP_para[sub][3].mean()) + """ + data_dict['cA'].append(HCP_para[sub][-10:,4].mean()) + data_dict['cB'].append(HCP_para[sub][-10:,5].mean()) + data_dict['cC'].append(HCP_para[sub][-10:,6].mean()) + """ +data_dict['g'] = np.array(data_dict['g']) +data_dict['gEE'] = np.array(data_dict['gEE']) +data_dict['gIE'] = np.array(data_dict['gIE']) +data_dict['gEI'] = np.array(data_dict['gEI']) + +""" +data_dict['cA'] = np.array(data_dict['cA']) +data_dict['cB'] = np.array(data_dict['cB']) +data_dict['cC'] = np.array(data_dict['cC']) +""" +for i in range(4): + sns.kdeplot(ax =ax[i//2, i%2], data= data_dict[para_name[i]], color='r', linewidth = 10) + #ax[i//2, i%2].set_title('hist of ' + para_name[i] +' across subjects', fontsize = 20) + ax[i//2, i%2].tick_params(labelsize = 20) + ax[i//2, i%2].set_xlabel(para_name[i], fontsize = 20) + ax[i//2, i%2].grid() + if i%2 == 1: + ax[i//2, i%2].set_ylabel('') +fig.tight_layout() +plt.show() +#fig.savefig(base_dir+'HCP_hist_modelparas.png') +fig, ax = plt.subplots(2,2, figsize=(12,8), sharex=True, sharey=True) +plt.rcParams["axes.labelsize"] = 20 +i = 0 +mask = np.tril_indices(83,-1) +for sub in ['198653', '257542', '154936', '280941']: + if i < 4: + data_dict = {} + data_dict['empFC'] = np.corrcoef(HCP_ts[sub].T)[mask] + data_dict['simFC'] = np.corrcoef(HCP_ts_sim[sub][:,10:])[mask] + data_dict['testFC'] = np.corrcoef(HCP_ts_test[sub])[mask] + sns.kdeplot(ax =ax[i//2, i%2], data=pd.DataFrame(data_dict), linewidth=10) + #ax[i//2, i%2].set_title('s'+sub, fontsize= 30) + ax[i//2, i%2].tick_params(labelsize=20) + ax[i//2, i%2].grid() + plt.setp(ax[i//2, i%2].get_legend().get_texts(), fontsize='25') + plt.setp(ax[i//2, i%2].get_legend().get_title(), fontsize ='40') + if i%2 == 1: + ax[i//2, i%2].set_ylabel('') + i += 1 +fig.tight_layout() +plt.show() +# fig.savefig(base_dir+'HCP_hist_FCs_pt.png') +fig, ax = plt.subplots(2,2,figsize=(12,8), sharex=True, sharey=True) +plt.rcParams["axes.labelsize"] = 20 +i = 0 +mask = np.tril_indices(83, -1) +for sub in ['198653', '257542', '154936', '280941']: + """ + The difference between original and modified structural connectivity weights. + """ + if i < 4: + data_dict = {} + data_dict['originalSC'] = HCP_sc[sub][mask] + sc_mod = HCP_sc_mod[sub] + sc = HCP_sc[sub] + w = (1 + np.tanh(sc_mod))*sc + w_n = 0.5*(w + w.T)/np.linalg.norm(0.5*(w + w.T)) + data_dict['modifiedSC'] = sc_mod[mask] + sns.kdeplot(ax=ax[i//2, i%2], data=pd.DataFrame(data_dict), linewidth=10) + #ax[i//2, i%2].set_title('s'+sub, fontsize = 30) + ax[i//2, i%2].tick_params(labelsize = 25) + ax[i//2, i%2].grid() + plt.setp(ax[i//2, i%2].get_legend().get_texts(), fontsize='25') + plt.setp(ax[i//2, i%2].get_legend().get_title(), fontsize='40') + + if i%2 == 1: + ax[i//2, i%2].set_ylabel('') + i += 1 +fig.tight_layout() +plt.show() +# fig.savefig(base_dir+'HCP_hist_SCs_pt.png') +corr_sim = {} +corr_sim['test'] = [] +corr_sim['train'] = [] +mask = np.tril_indices(83,-1) +for sub in HCP_ts: + corr_sim['train'].append(np.corrcoef(np.corrcoef(HCP_ts_sim[sub][:,10:])[mask], np.corrcoef(HCP_ts[sub].T)[mask])[0,1]) + corr_sim['test'].append(np.corrcoef(np.corrcoef(HCP_ts_test[sub])[mask], np.corrcoef(HCP_ts[sub].T)[mask])[0,1]) +fig, ax = plt.subplots(1,1, figsize=(12,8)) +plt.rcParams["axes.labelsize"] = 40 +sns.kdeplot(data=corr_sim, linewidth=10) +#ax.set_title('hist of Pearson correlation between the fitting and empirical BOLD', fontsize=20) +ax.tick_params(labelsize=25) +ax.grid() +plt.setp(ax.get_legend().get_texts(), fontsize='35') +plt.setp(ax.get_legend().get_title(), fontsize='50') + +# plt.savefig(base_dir+'HCP_hist_fittingcorrelation_pt.png') +# cross valid for identifiability... trained on the fist 2400 and tested on the last 2400 data points +corr_sim = {} +corr_sim['test'] = [] +corr_sim['train'] = [] +mask = np.tril_indices(83,-1) +for sub in HCP_ts_test: + corr_sim['train'].append(np.corrcoef(np.corrcoef(HCP_ts_sim_idt[sub][:,10:])[mask], np.corrcoef(HCP_ts_test[sub][:,:1200])[mask])[0,1]) + corr_sim['test'].append(np.corrcoef(np.corrcoef(HCP_ts_test_idt[sub])[mask], np.corrcoef(HCP_ts_test[sub][:,1200:])[mask])[0,1]) +fig, ax = plt.subplots(1,1, figsize=(12,8)) +plt.rcParams["axes.labelsize"] = 40 +sns.kdeplot(data=corr_sim, linewidth=10) +#ax.set_title('hist of Pearson correlation between the fitting and empirical BOLD', fontsize=20) +ax.tick_params(labelsize=25) +ax.grid() +plt.setp(ax.get_legend().get_texts(), fontsize='35') +plt.setp(ax.get_legend().get_title(), fontsize='50') + +fig, ax = plt.subplots(1,1, figsize=(12,8)) + +mask = np.tril_indices(83,-1) +data_dict = {} +data_dict['g'] = [] +data_dict['gEE'] = [] +data_dict['gIE'] = [] +data_dict['gEI'] = [] +""" +data_dict['cA'] = [] +data_dict['cB'] = [] +data_dict['cC'] = [] +""" +data_dict['$\mathregular{R_{train}}$'] = [] +data_dict['$\mathregular{R_{test}}$'] = [] +para_name = ['g', + 'gEE', + 'gIE', + 'gEI', + '$\mathregular{R_{train}}$', + '$\mathregular{R_{test}}$'] #'cA', 'cB', 'cC', + +for sub in HCP_ts: + data_dict['g'].append(HCP_para[sub][0].mean()) + data_dict['gEE'].append(HCP_para[sub][1].mean()) + data_dict['gIE'].append(HCP_para[sub][2].mean()) + data_dict['gEI'].append(HCP_para[sub][3].mean()) + """ + data_dict['cA'].append(HCP_para[sub][-10:,4].mean()) + data_dict['cB'].append(HCP_para[sub][-10:,5].mean()) + data_dict['cC'].append(HCP_para[sub][-10:,6].mean()) + """ + data_dict['$\mathregular{R_{train}}$'].append(np.corrcoef(np.corrcoef(HCP_ts[sub].T)[mask], np.corrcoef(HCP_ts_sim[sub][:,10:])[mask])[0,1]) + data_dict['$\mathregular{R_{test}}$'].append(np.corrcoef(np.corrcoef(HCP_ts[sub].T)[mask], np.corrcoef(HCP_ts_test[sub])[mask])[0,1]) +data_dict['g'] = np.array(data_dict['g']) +data_dict['gEE'] = np.array(data_dict['gEE']) +data_dict['gIE'] = np.array(data_dict['gIE']) +data_dict['gEI'] = np.array(data_dict['gEI']) +""" +data_dict['cA'] = np.array(data_dict['cA']) +data_dict['cB'] = np.array(data_dict['cB']) +data_dict['cC'] = np.array(data_dict['cC']) +""" +data_dict['$\mathregular{R_{train}}$'] = np.array(data_dict['$\mathregular{R_{train}}$']) +data_dict['$\mathregular{R_{test}}$'] = np.array(data_dict['$\mathregular{R_{test}}$']) + +corr_paras = np.zeros((len(data_dict), len(data_dict))) +for i in range(len(data_dict)): + for j in range(len(data_dict)): + corr_paras[i,j] = np.corrcoef(data_dict[para_name[i]], data_dict[para_name[j]])[0,1] + +im = ax.imshow(corr_paras, cmap='bwr') +ax.set_xticklabels(['0' ] + para_name, fontsize=20) +ax.set_yticklabels(['0'] + para_name, fontsize=20) +cbar = fig.colorbar(im, ax=ax) +cbar.ax.tick_params(labelsize=30) +cbar.set_ticks(np.linspace(-1,1,5)) +cbar.set_ticklabels(np.linspace(-1,1,5)) +plt.show() +# fig.savefig(base_dir + 'HCP_corr_modelparas+fit_pt.png') + +# %% +# 5. Generate the plots of the default mode network structure. +# --------------------------------------------------- +# + +l2k8_to_aparc_sort_idx = np.array([82, 30, 12, 8, 82, 20, 26, 24, 18, 28, 14, 22, 0, 23, 3, 29, 25, + 10, 5, 1, 4, 21, 15, 13, 9, 19, 11, 6, 7, 17, 31, 16, 2, 27, + 32, 33, 82, 71, 53, 49, 82, 61, 67, 65, 59, 69, 55, 63, 41, 64, 44, + 70, 66, 51, 46, 42, 45, 62, 56, 54, 50, 60, 52, 47, 48, 58, 72, 57, + 43, 68, 73, 74]) + +# Load the data from the .npy files. +l2k8_w = np.loadtxt('../data/HCP/hcp_grpavg_l2k8_sc33_weights.npy') # structural connectivity weights +l2k8_tl = np.loadtxt('../data/HCP/hcp_grpavg_l2k8_sc33_tractlengths.npy') # tract lengths between regions + +# Convert the files to pandas DataFrames. +df_l2k8_w = pd.DataFrame(l2k8_w) +df_l2k8_tl = pd.DataFrame(l2k8_tl) + +df_l2k8_w *= ((np.eye(83)*-1) + 1) # Multiply by a matrix of all 1s with 0s on the diagonal. +df_l2k8_wlog = np.log1p(df_l2k8_w) # Obtain log(1+x) for the weights x for the 83 x 83 matrix of regions. + +df_l2k8_wlog_divmax = df_l2k8_wlog.copy() # Copy the DataFrame. +df_l2k8_wlog_divmax /= df_l2k8_wlog.max() # Divide the DataFrame by its maximum amount to normalize it. + +_df_w = df_l2k8_wlog_divmax.copy() +_df_w[_df_w<0.2] = 0 +# Use these regions of the brain as focal nodes in representing the DMN. + +l2k8_lV1_idx = 21 # pericalcarine +l2k8_lM1_idx = 24 # precentral +l2k8_lA1_idx = 30 # superior temporal +l2k8_lS1_idx = 22 # postcentral +l2k8_lDLPFC_idx = 12 # lateral orbitofrontal +l2k8_lVLPFC_idx = 19 # pars orbitalis + +focal_nodes = [l2k8_lV1_idx, + l2k8_lA1_idx, + l2k8_lM1_idx, + l2k8_lS1_idx, + l2k8_lDLPFC_idx, + l2k8_lVLPFC_idx] + +focal_node_names = ['lV1', + 'lA1', + 'lM1', + 'lS1', + 'lDLPFC', + 'lVLPFC'] + +# Various surface representations for left and right hemisphere +# With the inflated surface representations, you can fully see the sulci. +lhi_file = '../data/HCP/fsav5_lh.inflated' +rhi_file = '../data/HCP/fsav5_rh.inflated' + +# The pial surface is used to calculate cortical gray matter volume and +# should accurately follow the boundary between the gray matter and the Cerebrospinal fluid (CSF). +lhp_file = '../data/HCP/fsav5_lh.pial' +rhp_file = '../data/HCP/fsav5_rh.pial' + +# Annotatation files for the parcellated regions +lh_annot_file = '../data/HCP/fsav5_lh.aparc.annot' +rh_annot_file = '../data/HCP/fsav5_rh.aparc.annot' + +# Return the vertex coordinates (rhi_vtx) and mesh triangles (rhi_tri) for the files +lhi_vtx, lhi_tri = nib.freesurfer.read_geometry(lhi_file) +rhi_vtx, rhi_tri = nib.freesurfer.read_geometry(rhi_file) + +lhp_vtx, lhp_tri = nib.freesurfer.read_geometry(lhp_file) +rhp_vtx, rhp_tri = nib.freesurfer.read_geometry(rhp_file) + +lh_annot = nib.freesurfer.read_annot(lh_annot_file) +rh_annot = nib.freesurfer.read_annot(rh_annot_file) + +vtx_l2k8_lh = lhp_vtx +tri_l2k8_lh = lhp_tri +#vtx_l2k8_lh = lhi_vtx +#tri_l2k8_lh = lhi_tri +rm_l2k8_lh = lh_annot[0] + +vtx_l2k8_rh = rhp_vtx +tri_l2k8_rh = rhp_tri +#vtx_l2k8_rh = rhi_vtx +#tri_l2k8_rh = rhi_tri +rm_l2k8_rh = rh_annot[0] + +# Concatenate them together for the regions. +vtx_l2k8_lhrh = np.concatenate([vtx_l2k8_lh,vtx_l2k8_rh], axis=0) +tri_l2k8_lhrh = np.concatenate([tri_l2k8_lh,tri_l2k8_rh+10242], axis=0) + +rm_l2k8_rh_p1 = rm_l2k8_rh + 36 +rm_l2k8_rh_p1[rm_l2k8_rh == 0] = 0 + +rm_l2k8_lhrh = np.concatenate([rm_l2k8_lh, rm_l2k8_rh_p1], axis=0) +rm_l2k8_lhrh + +hemi = np.ones_like(rm_l2k8_lhrh) +hemi[:10242] = 0 +hemi + +# %% +# functional connectivity values across regions for three different subjects +# The number coreresponds to the different columns of the FC (e.g., fc_54 has the 54th column). +fc_54 = np.array([0.07682509, -0.02128439, -0.03384062, 0.05320176, 0.05471503, + 0.13934857, 0.20621781, 0.17206639, 0.10512446, 0.10776393, + 0.29780534, 0.14799411, 0.30186789, 0.43196902, 0.33156828, + 0.25412462, 0.30922221, 0.12608491, 0.27343265, 0.41269201, + 0.25245327, 0.09809735, -0.11438562, 0.20858446, 0.07015369, + 0.08555236, 0.08659411, 0.04477133, 0.08317634, 0.11445271, + 0.12312393, 0.16196342, 0.06828829, 0.26020704, 0.12544554, + 0.05920674, 0.03225544, 0.01331936, -0.0215654 , 0.07715058, + -0.0389014 , 0.17158385, 0.0094781 , 0.02607858, 0.11923761, + 0.06269088, 0.10788386, 0.36012597, 0.22229945, 0.21825426, + 0.02970901, 0.24230718, 0.15896073, 0.40808317, 1. , + 0.33747972, 0.23850011, 0.30159616, 0.1130841 , 0.2969263 , + 0.47867023, 0.19959395, 0.06715806, -0.09413906, 0.20176679, + 0.10768835, 0.02944969, 0.04499624, -0.03408385, 0.17570769, + 0.16808477, 0.06259806, 0.26563371, 0.14068841, 0.27621275, + 0.1306482 , 0.09807702, 0.09850341, 0.05073915, 0.03550151, + 0.17014859, 0.06636479, 0.17870341]) + +fc_60 = np.array([0.15569094, -0.03885425, 0.04890683, 0.16171774, 0.04725465, + 0.07753457, 0.30875898, 0.35187976, 0.34662134, 0.09912189, + 0.29338462, 0.25739101, 0.17445902, 0.4670589 , 0.62694002, + 0.25998084, 0.25925487, 0.33025557, 0.62163083, 0.8369526 , + 0.52574254, 0.31510396, 0.06566058, 0.44301474, 0.27148304, + 0.30111694, 0.08374708, 0.11905168, 0.13817447, 0.35589242, + 0.20325078, 0.28584166, 0.13142756, 0.2040255 , 0.27938807, + 0.11800877, 0.03426762, 0.03599968, 0.01689927, 0.24260188, + -0.03815829, 0.19310196, 0.03499076, 0.10427147, 0.13123463, + 0.11359923, 0.2292006 , 0.51829568, 0.3945326 , 0.47134833, + 0.14234746, 0.17401509, 0.24885599, 0.29500217, 0.47867023, + 0.62890929, 0.34624069, 0.36416877, 0.3573325 , 0.69456097, + 1. , 0.45229606, 0.34654782, 0.12841021, 0.47623741, + 0.23457237, 0.06286334, 0.02768237, 0.00291156, 0.27546106, + 0.37390449, 0.29662638, 0.3605275 , 0.22705125, 0.26378275, + 0.22736693, 0.2490869 , 0.03850265, 0.0448428 , 0.07056435, + 0.13632041, 0.06124779, 0.13272672]) + +fc_48 = np.array([ 0.08348805, 0.25433363, 0.30718973, 0.21778039, -0.03452652, + -0.07191728, 0.04859293, 0.58070912, 0.3228884 , 0.03022416, + 0.096414 , 0.2067467 , 0.02934248, 0.11758467, 0.3395558 , + 0.2618392 , -0.04729804, 0.02403979, 0.33068409, 0.18542544, + 0.11157786, 0.24558548, 0.35921583, 0.14047045, 0.16072442, + 0.08725231, 0.04247384, 0.21984861, 0.15750137, 0.49704743, + 0.38107799, 0.18360744, 0.13072013, 0.03044336, 0.25038236, + 0.21535622, 0.04945266, 0.09180092, 0.07482232, 0.19016411, + -0.01817344, 0.30486544, 0.58977128, 0.33755094, 0.19382802, + 0.33391812, 0.38968696, 0.47185658, 1. , 0.74660754, + 0.25389967, 0.11339373, 0.30413138, 0.18101733, 0.22229945, + 0.48579973, 0.29879378, 0.14413303, 0.12710287, 0.61299693, + 0.3945326 , 0.03012108, 0.25349449, 0.36321952, 0.1842467 , + 0.09006558, 0.00976703, 0.01410019, 0.08918964, 0.24576006, + 0.72460302, 0.25165407, 0.31029377, 0.155622 , 0.07071528, + 0.26798855, 0.34361951, 0.15074631, 0.07711461, 0.13811093, + 0.20971265, 0.03274206, 0.08477825]) + +# Set the seed for the seed-based functional connectivity analysis. +seed = np.zeros_like(fc_48) +seed[60] = 1 +_cm = 'RdBu_r' + +# keywords +kws = {'edgecolors': 'k', + 'cmap': _cm, # 'Reds_r', #RdBu_r', # jet', + 'vmax': 0.5, + 'vmin': -.3, + 'alpha': None, + 'linewidth': 0.05} + +plot_surface_mpl_mv(vtx=vtx_l2k8_lhrh, + tri=tri_l2k8_lhrh, + hemi=hemi, + reorient='fs', + data=fc_60[l2k8_to_aparc_sort_idx][rm_l2k8_lhrh], + shade_kwargs=kws) +"""_cm = 'BuGn' +kws = {'edgecolors': 'k', 'cmap': _cm, #$ 'Reds_r', #RdBu_r', # jet', + 'vmax': 1.0,'vmin': 0.0, 'alpha': 0.15, 'linewidth': 0.05} +mask = np.array([False]*40960) +mask[:vtx_l2k8_lhrh.shape[0]]= seed[l2k8_to_aparc_sort_idx][rm_l2k8_lhrh] == 1 +mask1 = seed[l2k8_to_aparc_sort_idx][rm_l2k8_lhrh] == 1 +plot_surface_mpl_mv(vtx=vtx_l2k8_lhrh,tri=tri_l2k8_lhrh,hemi=hemi,reorient='fs',data=seed[l2k8_to_aparc_sort_idx][rm_l2k8_lhrh], + shade_kwargs=kws) """ + +fc_test_54 = np.array([4.04483037e-02, 7.32430791e-02, 7.76209493e-02, -4.46576574e-02, + 6.90320864e-02, 6.08216248e-02, 1.25866896e-01, 1.48259337e-01, + 1.46379082e-01, 1.60607085e-01, 1.99941445e-01, -4.73537989e-02, + 9.32955181e-02, 1.84146202e-01, 1.12089813e-01, 1.58828217e-01, + 1.37430776e-01, 3.38672182e-02, 8.48897343e-02, 1.80693377e-01, + 3.68792057e-02, -1.50098567e-02, 3.20046130e-02, 1.00762267e-02, + 1.56791288e-02, 2.76290084e-02, -3.96479495e-02, 2.80856785e-02, + 3.91929104e-02, 4.77790081e-02, 4.41821994e-02, 6.01223904e-02, + 5.10560110e-02, 7.76196128e-02, 1.12436866e-01, 2.73838335e-02, + -2.39388323e-02, 1.71982611e-02, 7.29186160e-02, 5.29034320e-02, + -1.99933514e-03, 7.06724821e-02, 5.93341413e-02, -4.13728002e-02, + 1.29869528e-04, 1.23697792e-01, 1.43169035e-01, 2.10767740e-01, + 3.11019690e-01, 1.38353247e-01, 1.86178990e-01, 2.70765146e-01, + -9.97468529e-03, 3.78946661e-01, 1.00000000e+00, 3.15421661e-01, + 1.52251462e-01, 1.33510163e-01, 1.12895918e-01, 1.14309818e-01, + 2.51496899e-01, 4.27383293e-02, 4.84373107e-02, 4.49691775e-02, + 6.54305540e-02, 6.93010133e-02, -3.49596325e-02, -1.16978796e-02, + -5.19326833e-02, 1.19297520e-01, 1.13479555e-01, 9.82750559e-02, + 1.37192609e-01, 1.10936182e-01, 1.36305410e-01, 1.35283865e-01, + 1.46130802e-01, 1.08109444e-01, 2.04031992e-02, -1.04834287e-02, + -5.35527088e-02, 3.55320361e-02, 5.05661223e-02]) + +fc_test_13 = np.array([ 0.02064446, -0.02125758, 0.06787312, -0.05500374, -0.03885027, + -0.01193319, 0.02590079, 0.07862877, 0.11140791, 0.12894951, + 0.29816439, -0.07257691, 0.39432721, 1. , 0.41285587, + 0.15905054, 0.02590495, 0.02414737, 0.16389443, 0.33678748, + 0.06684584, 0.02343256, 0.01415101, -0.00517585, 0.0300626 , + 0.14180168, 0.04087309, -0.02108286, 0.01249072, 0.02652301, + 0.0147603 , 0.05418977, 0.01977914, 0.03116899, 0.03394 , + 0.0515946 , -0.0622778 , -0.03677459, 0.02683279, 0.09717596, + 0.04846726, 0.0334019 , 0.06863676, -0.02671001, -0.04969303, + 0.04434678, 0.0240701 , 0.05507989, 0.09096978, 0.08247514, + 0.10845567, 0.22910861, -0.04563424, 0.13561862, 0.1841462 , + 0.16418864, 0.11336982, 0.11471379, 0.11565096, 0.05493549, + 0.27033029, 0.01242047, 0.03619441, 0.08456932, 0.05367578, + 0.10041786, -0.04090774, -0.01034624, -0.03077638, 0.08534603, + 0.06689338, 0.06317927, 0.07676879, 0.11862627, 0.05053353, + 0.10158778, 0.09397584, 0.05032458, 0.0745082 , 0.04716247, + 0.08466447, -0.05279958, 0.03736607]) + +fc_test_60 = np.array([3.51703550e-02, 5.18501633e-02, 3.00072240e-02, -3.22250993e-03, + 5.70771272e-03, 3.04967810e-02, 6.09221308e-02, 1.09301758e-01, + 1.10533317e-01, 1.28118958e-01, 2.10016752e-01, 1.09933921e-03, + 9.57229525e-02, 2.70330287e-01, 1.90774291e-01, 1.47761080e-01, + 8.97752877e-02, 1.11693599e-01, 1.95514360e-01, 2.72475220e-01, + 1.00020912e-01, 6.70112479e-02, 6.97584082e-02, 7.93095470e-02, + 8.93715930e-02, 7.42034261e-02, 3.39152899e-02, -4.49641276e-02, + -2.76792903e-02, 2.70008066e-02, -4.64044487e-03, 7.18338669e-02, + 2.04857838e-03, 8.14713997e-02, 1.17725733e-01, 7.42598533e-02, + 2.96991318e-03, -5.39451022e-02, 2.48211990e-02, 7.97761141e-02, + 1.85885648e-02, 1.56141833e-01, 3.34977809e-02, -3.00353457e-02, + -1.62362070e-02, 8.00848864e-02, 1.34135085e-01, 1.40980153e-01, + 1.77252780e-01, 1.22131104e-01, 2.10072865e-01, 2.10706170e-01, + -7.68060410e-04, 1.37570662e-01, 2.51496899e-01, 2.69339188e-01, + 2.11895067e-01, 1.73590178e-01, 2.42313753e-01, 2.07098915e-01, + 1.00000000e+00, 8.07081133e-02, 1.00150620e-01, 8.85262622e-02, + 9.18910872e-02, 1.83301361e-01, -1.60720267e-02, 3.92561061e-02, + 1.93154903e-03, 2.41607488e-01, 2.17235434e-01, 1.30866821e-01, + 2.47044743e-01, 2.56881848e-01, 1.50520001e-01, 2.11369632e-01, + 1.62973354e-01, 1.44820176e-01, 2.40272233e-02, 4.40490309e-02, + 3.92930927e-02, -4.07958819e-02, 4.43558957e-02]) + +fc_test_48 = np.array([ 0.04456317, 0.08430425, 0.10337379, -0.0751953 , 0.09872982, + 0.1291111 , 0.18314354, 0.20767963, 0.23303186, 0.19998463, + 0.14810478, -0.05109569, 0.05073159, 0.09096978, 0.08419494, + 0.17615717, 0.16614261, 0.0196335 , 0.09354732, 0.08995714, + 0.05359471, -0.00912246, 0.03514153, 0.02057001, 0.02612656, + 0.01342597, -0.0294798 , 0.00696833, 0.02931964, 0.03325134, + 0.07043438, 0.0502618 , 0.02825642, 0.18071165, 0.18340093, + 0.06189989, -0.00533773, 0.06097923, 0.10050126, 0.06412034, + 0.01178386, 0.17734297, 0.06325622, -0.02823774, -0.02826472, + 0.23305105, 0.33101713, 0.32282376, 1. , 0.31766683, + 0.28207603, 0.209288 , 0.00738726, 0.27140204, 0.31101969, + 0.22895553, 0.2234192 , 0.18360893, 0.16395947, 0.21453826, + 0.17725278, 0.0341468 , 0.03666224, 0.00832946, 0.05212808, + 0.07890769, 0.01180786, -0.01867186, -0.04819648, 0.11294647, + 0.16827467, 0.14320286, 0.16911658, 0.13737125, 0.19928486, + 0.21347781, 0.14709652, 0.16063771, 0.02247446, 0.04255669, + -0.00888587, 0.04339403, 0.07802365]) + +fc_test_41 = np.array([ 2.91574121e-02, 3.99662172e-02, 9.14379959e-02, -8.00334527e-03, + -2.50339264e-02, 3.63609413e-02, 1.06375335e-01, 1.06307706e-01, + 1.41449230e-01, 1.52185927e-01, 1.40082493e-01, 5.44236473e-03, + -2.84602123e-03, 3.34019030e-02, 7.24177217e-02, 1.29194284e-01, + 1.10487105e-01, 9.90362175e-02, 1.04988485e-01, 1.12635543e-01, + 1.02215036e-01, 3.63280972e-02, 9.13308055e-02, 6.25911165e-02, + 6.41443914e-02, 5.56257774e-02, 4.64617068e-02, -1.15959055e-02, + 9.60144121e-03, 5.03256381e-02, 3.93803611e-03, 5.98747396e-02, + -3.76220011e-02, 6.67289208e-02, 1.50632776e-01, 5.85915414e-02, + -5.79427852e-03, -2.67941385e-02, 2.60525511e-02, 8.43785417e-02, + 6.25090275e-02, 1.00000000e+00, 8.24955180e-02, -3.07733096e-02, + 3.37654868e-02, 1.69057126e-01, 1.98151129e-01, 2.99425215e-01, + 1.77342974e-01, 1.43791454e-01, 1.87007384e-01, 1.55162276e-01, + 2.97501823e-03, 1.24536132e-01, 7.06724821e-02, 1.09327160e-01, + 2.14722698e-01, 1.89543923e-01, 1.55433207e-01, 2.17070757e-01, + 1.56141833e-01, 7.86344152e-02, 8.39345856e-02, 9.96740939e-02, + 7.45054762e-02, 1.31522853e-01, 4.30021043e-02, -1.37114672e-02, + -3.84809861e-04, 1.87042770e-01, 2.42837501e-01, 1.73869281e-01, + 2.42828040e-01, 1.65126795e-01, 3.06438319e-01, 2.48172135e-01, + 2.74101333e-01, 2.87125048e-01, 3.19166679e-03, 1.65886750e-02, + 3.04314795e-02, 3.21617761e-02, 4.86736314e-03]) + +seed = np.zeros_like(fc_48) +seed[48] = 1 +_cm = 'RdBu_r' +kws = {'edgecolors': 'k', + 'cmap': _cm, # 'Reds_r', #RdBu_r', # jet', + 'vmax': 0.35, + 'vmin': -.2, + 'alpha': None, + 'linewidth': 0.05} + +plot_surface_mpl_mv(vtx=vtx_l2k8_lhrh, + tri=tri_l2k8_lhrh, + hemi=hemi, + reorient='fs', + data=fc_test_48[l2k8_to_aparc_sort_idx][rm_l2k8_lhrh], + shade_kwargs=kws) +"""_cm = 'BuGn' +kws = {'edgecolors': 'k', 'cmap': _cm, #$ 'Reds_r', #RdBu_r', # jet', + 'vmax': 1.0,'vmin': 0.0, 'alpha': 0.15, 'linewidth': 0.05} +mask = np.array([False]*40960) +mask[:vtx_l2k8_lhrh.shape[0]]= seed[l2k8_to_aparc_sort_idx][rm_l2k8_lhrh] == 1 +mask1 = seed[l2k8_to_aparc_sort_idx][rm_l2k8_lhrh] == 1 +plot_surface_mpl_mv(vtx=vtx_l2k8_lhrh,tri=tri_l2k8_lhrh,hemi=hemi,reorient='fs',data=seed[l2k8_to_aparc_sort_idx][rm_l2k8_lhrh], + shade_kwargs=kws)""" +# %% +# References +# --------------------------------------------------- +#Deep Learning-Based Parameter Estimation for Neurophysiological Models of Neuroimaging Data. (2022). In Health & Medicine Week (pp. 1830-). NewsRX LLC.