forked from dagush/WholeBrain
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathDecoEtAl2018_Fig3E-Right-Simulated.py
129 lines (110 loc) · 5.63 KB
/
DecoEtAl2018_Fig3E-Right-Simulated.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
122
123
124
125
126
127
128
129
# ==========================================================================
# ==========================================================================
# Computes the Functional Connectivity Dynamics (FCD)
#
# From the original code:
# --------------------------------------------------------------------------
#
# Computes simulations with the Dynamic Mean Field Model (DMF) using
# Feedback Inhibitory Control (FIC) and Regional Drug Receptor Modulation (RDRM):
#
# - the optimal coupling (we=2.1) for fitting the placebo condition
# - the optimal neuromodulator gain for fitting the LSD condition (wge=0.2)
#
# Taken from the code (Code_Figure3.m) from:
#
# [DecoEtAl_2018] Deco,G., Cruzat,J., Cabral, J., Knudsen,G.M., Carhart-Harris,R.L., Whybrow,P.C.,
# Whole-brain multimodal neuroimaging model using serotonin receptor maps explain non-linear functional effects of LSD
# Logothetis,N.K. & Kringelbach,M.L. (2018) Current Biology
# https://www.cell.com/current-biology/pdfExtended/S0960-9822(18)31045-5
#
# Code written by Gustavo Deco [email protected] 2017
# Reviewed by Josephine Cruzat and Joana Cabral
#
# Translated to Python by Gustavo Patow
# ==========================================================================
# ==========================================================================
from pathlib import Path
import time
# --------------------------------------------------------------------------
# Begin setup...
# --------------------------------------------------------------------------
from DecoEtAl2018_Setup import *
import functions.simulateFCD as simulateFCD
simulateFCD.integrator = integrator
simulateFCD.simModel = simulateBOLD
# --------------------------------------------------------------------------
# End setup...
# --------------------------------------------------------------------------
def my_hist(x, bin_centers):
bin_edges = np.r_[-np.Inf, 0.5 * (bin_centers[:-1] + bin_centers[1:]), np.Inf]
counts, edges = np.histogram(x, bin_edges)
return [counts, bin_centers]
# ============================================================================
# ============= Compute the J values for Balance conditions ==================
# ============================================================================
# Define optimal parameters
# ==== J is calculated this only once, then saved
baseName = "Data_Produced/SC90/J_Balance_we2.1.mat"
balancedG = BalanceFIC.Balance_J9(2.1, C, False, baseName)
balancedG['J'] = balancedG['J'].flatten()
balancedG['we'] = balancedG['we'].flatten()
serotonin2A.setParms(balancedG)
# serotonin2A.wgaini = 0. # Placebo conditions, to calibrate the J's...
# serotonin2A.wgaine = 0.
serotonin2A.setParms({'S_E':0., 'S_I':0.})
recompileSignatures()
initRandom()
# ============================================================================
# ============= Simulate Placebo =============================================
# ============================================================================
pla_path = 'Data_Produced/SC90/FCD_values_placebo.mat'
if not Path(pla_path).is_file():
# SIMULATION OF OPTIMAL PLACEBO
print("\n\nSIMULATION OF OPTIMAL PLACEBO")
# serotonin2A.wgaini = 0.
# serotonin2A.wgaine = 0. # 0 for placebo, 0.2 for LSD
serotonin2A.setParms({'S_E':0., 'S_I':0.})
recompileSignatures()
start_time = time.clock()
cotsampling_pla_s = simulateFCD.simulate(NumSubjects, C)
print("\n\n--- TOTAL TIME: {} seconds ---\n\n".format(time.clock() - start_time))
sio.savemat(pla_path, {'cotsampling_pla_s': cotsampling_pla_s}) # save FCD_values_placebo cotsampling_pla_s
else:
print("LOADING OPTIMAL PLACEBO")
cotsampling_pla_s = sio.loadmat(pla_path)['cotsampling_pla_s']
# ============================================================================
# ============= Simulate LSD =================================================
# ============================================================================
lsd_path = 'Data_Produced/SC90/FCD_values_lsd.mat'
if True: #not Path(lsd_path).is_file():
# SIMULATION OF OPTIMAL LSD fit
print("\n\nSIMULATION OF OPTIMAL LSD fit ")
# serotonin2A.wgaini = 0.
# serotonin2A.wgaine = 0.2 # 0 for placebo, 0.2 for LSD
serotonin2A.setParms({'S_E':0.2, 'S_I':0.})
recompileSignatures()
start_time = time.clock()
cotsampling_lsd_s = simulateFCD.simulate(NumSubjects, C)
print("\n\n--- TOTAL TIME: {} seconds ---\n\n".format(time.clock() - start_time))
sio.savemat(lsd_path, {'cotsampling_lsd_s': cotsampling_lsd_s}) # save FCD_values_lsd cotsampling_lsd_s
else:
print("LOADING OPTIMAL LSD fit")
cotsampling_lsd_s = sio.loadmat(lsd_path)['cotsampling_lsd_s']
# ============================================================================
# Plot
# ============================================================================
[h_pla, x1] = my_hist(cotsampling_pla_s[:].T.flatten(), np.arange(-.1, 1.025, .025))
[h_lsd, x] = my_hist(cotsampling_lsd_s[:].T.flatten(), np.arange(-.1, 1.025, .025))
import matplotlib.pyplot as plt
width=0.01
plaBar = plt.bar(x, h_pla, width=width, color="red", label="Placebo")
lsdBar = plt.bar(x+width, h_lsd, width=width, color="blue", label="LSD")
plt.xlabel('FCD values')
plt.ylabel('Count')
plt.legend(handles=[plaBar, lsdBar], loc='upper right')
plt.title('Simulated data')
plt.show()
# ================================================================================================================
# ================================================================================================================
# ================================================================================================================EOF