-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathoptimize_SC_JR.py
121 lines (91 loc) · 4.86 KB
/
optimize_SC_JR.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
# -*- coding: utf-8 -*-
"""
Created on Wen Jan 06 12:23:19 2021
@author: Carlos Coronel
Code for SC optimization based in the method proposed by Deco et al. (2019) [1]. First, the model was fitted to the
empirical Functionl Connectivity matrix using the original Structural Connectivity matrix. Then, Structural
Connectivity was updated iteratively employing the point-to-point difference between empirical and simulated
Functional Connectivity matrices. Negative values within the optimized Structural Connectivity matrix were discarded.
The code started with an alpha = 1.37 employing the human Structural Connectivity matrix used in [2]
[1] Deco, G., Cruzat, J., Cabral, J., Tagliazucchi, E., Laufs, H., Logothetis, N. K.,
& Kringelbach, M. L. (2019). Awakening: Predicting external stimulation to
force transitions between different brain states. Proceedings of the National
Academy of Sciences, 116(36), 18088-18097.
[2] Luppi, A. I., & Stamatakis, E. A. (2021). Combining network topology and
information theory to construct representative brain networks. Network Neuroscience,
5(1), 96-124.
"""
import numpy as np
from scipy import signal
import BOLDModel as BD
import JansenRitModel as JR
from scipy import stats
import graph_utils
import networkx as nx
#Uppter triangular of the empirical Functional Connectivity matrix
dist_placebo_rest = graph_utils.get_uptri(np.load('FC_placebo_rest.npy'))
#Simulation parameters
JR.dt = 1E-3 #Integration step
JR.teq = 60 #Simulation time for stabilizing the system
JR.tmax = 1200 + JR.teq * 2 #Simulation time
JR.downsamp = 10 #Downsampling to reduce the number of points
Neq = int(JR.teq / JR.dt / JR.downsamp)
Nmax = int(JR.tmax / JR.dt / JR.downsamp)
Ntotal = Neq + Nmax
ttotal = JR.teq + JR.tmax #Total simulation time
#Structural connectivity
nnodes = 100 #number of nodes
JR.nnodes = nnodes
JR.M = nx.to_numpy_array(nx.watts_strogatz_graph(nnodes, 8, 0.05, 1)) #Toy matrix. Replace with the one used in [2]
JR.norm = np.mean(np.sum(JR.M,0)) #Normalization factor
#Node parameters
JR.alpha = 1.37 * np.ones(nnodes) #Long-range pyramidal-pyramidal coupling
JR.C4 = (0.3 + JR.alpha * 0.3 / 0.5) * 135 #Connectivity between inhibitory pop. and pyramidal pop.
#Current copy of the SC matrix
C = np.copy(JR.M)
JR.norm = np.mean(np.sum(JR.M,0))
seeds = 50 #Number of FCs computed in each iterations
iters = 200 #Number of iterations
fitting = np.zeros((4,iters)) #Matrix to save the results: fitting-related metrics across iterations
all_SCs = np.zeros((nnodes,nnodes,iters)) #Array to save all SCs across iterations
epsilon = 0.01 #Convergence rate
#JR update
JR.update()
for i in range(0,iters):
FC_BOLD = np.zeros((nnodes,nnodes))
for ss in range(0,seeds):
all_SCs[:,:,i] = np.copy(C)
#np.random.seed(seed)
JR.seed = ss
y, time_vector = JR.Sim(verbose = False)
pyrm = JR.C2 * y[:,1] - JR.C4 * y[:,2] + JR.C * JR.alpha * y[:,3] #EEG-like output of the model
#Firing rates
rE = JR.s(pyrm,JR.r0)
#Simulating BOLD
BOLD_signals = BD.Sim(rE, nnodes, JR.dt * JR.downsamp)
BOLD_signals = BOLD_signals[Neq:,:]
#Same dt as empirical signals (RT = 1.5 seconds)
BOLD_dt = 1.5
BOLD_signals = signal.decimate(BOLD_signals, n = 3,
q = int(BOLD_dt / JR.dt / JR.downsamp), axis = 0)
#Filter the BOLD-like signal between 0.01 and 0.1 Hz
Fmin, Fmax = 0.01, 0.08
a0, b0 = signal.bessel(3, [2 * BOLD_dt * Fmin, 2 * BOLD_dt * Fmax], btype = 'bandpass')
BOLD_filt = signal.filtfilt(a0, b0, BOLD_signals[:,:], axis = 0)
cut0, cut1 = int(Neq / (BOLD_dt / JR.dt / JR.downsamp)), int((Nmax - Neq) / (BOLD_dt / JR.dt / JR.downsamp))
BOLD_filt = BOLD_filt[cut0:cut1,:]
FC_BOLD += np.corrcoef((BOLD_filt.T))
print([ss,i])
FC_BOLD /= seeds #Averaged FC matrix across random seeds
dist_sim = graph_utils.get_uptri((FC_BOLD)) #Upper triangular of the simulated FC matrix
#Fitting metrics
fitting[0,i] = np.mean(dist_placebo_rest) - np.mean(dist_sim) #Difference in mean Overall FC (global correlations)
fitting[1,i] = stats.ks_2samp(dist_placebo_rest, dist_sim)[0] #Difference in weights' distributions
fitting[2,i] = np.linalg.norm(dist_placebo_rest - dist_sim) #Euclidean distance
fitting[3,i] = stats.pearsonr(dist_placebo_rest,dist_sim)[0] #Pearson Correlation
#For following up the optimization
print(ss,fitting[:,i],np.sum(C), np.sum(C>0))
#Updating the C matrix
C += graph_utils.matrix_recon(epsilon*(dist_placebo_rest - dist_sim))
C[C < 0] = 0 #Discarding negative values
JR.M = C #Update the SC within the model