forked from Spokhim/MouseBrainModelling
-
Notifications
You must be signed in to change notification settings - Fork 1
/
SimulationPipeline.py
202 lines (165 loc) · 9.55 KB
/
SimulationPipeline.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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
# Simulation Pipeline as Python Script
# Intended for use in SuperComputerCluster Automation
# Outputs Time Series Matrix, Parameters, FCM and SCM vs FCM Spearson Correlation
# First Import the packages.
import numpy as np
import pandas as pd
import pylab
import matplotlib.pyplot as plt
from scipy import stats
import inspect
import os
import csv
import time
from tvb.simulator.lab import *
from tvb.simulator.plot.tools import *
def Simul_Pipeline(ParamsDict):
'''
Simulation Pipeline, Does the simulations and returns Scorr of FC vs SC. Saves some csvs in do-no-track folder. Takes in inputs as ParamsDict
'''
# Define the model if it is not yet defined.
if "MODEL" not in ParamsDict:
ParamsDict["MODEL"] = models.WilsonCowan(c_ee=ParamsDict["MODEL_c_ee"],c_ei=ParamsDict["MODEL_c_ei"],c_ie=ParamsDict["MODEL_c_ie"] ,c_ii=ParamsDict["MODEL_c_ii"],
a_e=numpy.array([1.0]),a_i=numpy.array([1.0]),b_e=numpy.array([2.8]),b_i=numpy.array([2.8]),tau_e=numpy.array([10.0]),
tau_i=numpy.array([65.0]),)
# Load the connectivity data from a zip file.
con = connectivity.Connectivity.from_file(os.getcwd() +"/Connectomes/" + ParamsDict["name"] + ".zip")
# Now need to prepare the connectivity data accordingly. Unfortuantely doesn't load eveyrthing in properly. May need to adjust con.undirected in future.
# Remove the ith row and column in centres, tract_lengths and weights. i.e. the specified region(s)
con.centres = np.delete(con.centres,ParamsDict["REMOVE"])
con.weights = np.delete(con.weights,obj=ParamsDict["REMOVE"],axis=0)
con.weights = np.delete(con.weights,obj=ParamsDict["REMOVE"],axis=1)
con.tract_lengths = np.delete(con.tract_lengths,obj=ParamsDict["REMOVE"],axis=0)
con.tract_lengths = np.delete(con.tract_lengths,obj=ParamsDict["REMOVE"],axis=1)
# Number of regions
con.number_of_regions = con.weights.shape[0]
# Change to Connectome to Binary if desired:
if ParamsDict["BINARY"]==True:
con.weights = con.weights!=0
# Set the parameter of the resting state simulation
# Bold Simulation, Initial Conditions Provided
if (ParamsDict["BOLD"] == True and "Init_Cons" in ParamsDict ):
sim = simulator.Simulator(model=ParamsDict["MODEL"],
connectivity=con,
coupling=coupling.Linear(a=ParamsDict["G"]),
integrator=integrators.EulerStochastic(dt=ParamsDict["dt"],noise=noise.Additive(nsig=ParamsDict["noise"],
random_stream=np.random.RandomState(ParamsDict["RandState"]))),
monitors=(monitors.Bold(period=1e3),
monitors.TemporalAverage(period=1e3)),
simulation_length=ParamsDict["Simul_length"],
initial_conditions=ParamsDict["Init_Cons"],
).configure()
# Run the resting state simulation
(bold_time, bold_data), _ = sim.run()
# Bold Simulation, Initial Conditions Not Provided
elif (ParamsDict["BOLD"] == True and "Init_Cons" not in ParamsDict ):
sim = simulator.Simulator(model=ParamsDict["MODEL"],
connectivity=con,
coupling=coupling.Linear(a=ParamsDict["G"]),
integrator=integrators.EulerStochastic(dt=ParamsDict["dt"],noise=noise.Additive(nsig=ParamsDict["noise"],
random_stream=np.random.RandomState(ParamsDict["RandState"]))),
monitors=(monitors.Bold(period=1e3),
monitors.TemporalAverage(period=1e3)),
simulation_length=ParamsDict["Simul_length"],
#initial_conditions=0.5 + numpy.zeros((con.number_of_regions*con.number_of_regions,2,con.number_of_regions,1)),
).configure()
# Run the resting state simulation
(bold_time, bold_data), _ = sim.run()
# No Monitors, Initial Conditions Provided
elif (ParamsDict["BOLD"] == False and "Init_Cons" in ParamsDict ):
sim = simulator.Simulator(model=ParamsDict["MODEL"],
connectivity=con,
coupling=coupling.Linear(a=ParamsDict["G"]),
integrator=integrators.EulerStochastic(dt=ParamsDict["dt"],noise=noise.Additive(nsig=ParamsDict["noise"],
random_stream=np.random.RandomState(ParamsDict["RandState"]))),
simulation_length=ParamsDict["Simul_length"],
initial_conditions=ParamsDict["Init_Cons"],
).configure()
# Run the resting state simulation
awer = sim.run()
bold_time = awer[0][0]
bold_data = awer[0][1]
# No Monitors, Initial Conditions Not Provided
else:
sim = simulator.Simulator(model=ParamsDict["MODEL"],
connectivity=con,
coupling=coupling.Linear(a=ParamsDict["G"]),
integrator=integrators.EulerStochastic(dt=ParamsDict["dt"],noise=noise.Additive(nsig=ParamsDict["noise"],
random_stream=np.random.RandomState(ParamsDict["RandState"]))),
simulation_length=ParamsDict["Simul_length"],
#initial_conditions=0.5 + numpy.zeros((con.number_of_regions*con.number_of_regions,2,con.number_of_regions,1)),
).configure()
# Run the resting state simulation
awer = sim.run()
bold_time = awer[0][0]
bold_data = awer[0][1]
# Got lazy
Snip = ParamsDict["Snip"]
# Note, if BOLD = False, Snip gets multiplied by 100.
if ParamsDict["BOLD"] == False:
Snip = 100 * Snip
# Functional Connectivity Matrix.
# We note that this is a static analysis. More advanced version would be a Dynamic version with windowing.
# Convert Simulation output into a form usable by Numpy.
TSeriesMatrix = np.empty((bold_data.shape[2], bold_data.shape[0]-Snip))
for i in range(bold_data.shape[2]):
TSeriesMatrix[i] = bold_data[Snip:,0,i].flatten()
# Functional Conenctivity MAtrix = Pearson Correlation.
FCM = np.corrcoef(TSeriesMatrix)
# Set diagonals to NaN
FCM1 = FCM
np.fill_diagonal(FCM1,np.nan)
# Comparing SC vs FC with Spearman Corr (Known as SCFC)
# Check if SCM is symmetric:
SCM = con.weights
Sym_check = numpy.allclose(SCM, SCM.T,equal_nan=True)
if Sym_check == True:
#It is a symmetric SCM, so only use upper triangles
# Grab Upper triangles
FCM_Upper = FCM[np.triu_indices(FCM.shape[0], k = 1)]
SCM_Upper = con.weights[np.triu_indices(con.weights.shape[0], k = 1)]
elif Sym_check == False:
# If SCM is not symmetric, need to calcualte spearman corr for entire matrix.
# Set Diagonal to Nans
np.fill_diagonal(SCM,np.nan)
# Remove all Nans for SCM and FCM
SCM_Upper = SCM[~numpy.isnan(SCM)]
FCM_Upper = FCM1[~numpy.isnan(FCM1)]
# Spearman Correlation
Scorr = stats.spearmanr(a=FCM_Upper,b=SCM_Upper)
print(Scorr)
# Calculate FC-FC score for mouse if requested.
if ParamsDict["FCFC"] == True:
# FCM_exp
FCM_exp = np.genfromtxt('FCM_MouseExperimental.csv',delimiter = "\t")
# Set diagonals to NaN
np.fill_diagonal(FCM_exp,np.nan)
# Remove the ith row and column in FCM (i.e. the specified region)
FCM_exp = np.delete(FCM_exp,obj=ParamsDict["REMOVE"],axis=0)
FCM_exp = np.delete(FCM_exp,obj=ParamsDict["REMOVE"],axis=1)
# Comparing FC_experimental Vs FC_Simulation with Spearman Correlation
FCM_Exp_U = FCM_exp[np.triu_indices(FCM_exp.shape[0], k = 1)]
FCM_Upper = FCM[np.triu_indices(FCM.shape[0], k = 1)]
# FC-FC Spearman Correlation
FCFC = stats.spearmanr(a=FCM_Exp_U,b=FCM_Upper)
print(FCFC)
# Concancatanate to end of Scorr output file.
Scorr = Scorr + FCFC
# Export the simulation information
time_now = time.strftime("%Y%m%d-%H%M%S")
date = time.strftime("%Y_%m_%d/")
# Create new directory which is the date.
os.makedirs("do-not-track/" + date,exist_ok=True)
# Params Dictionary - Note how we sort the dictionary.
with open("do-not-track/" + date + ParamsDict["tag"] + "_" + ParamsDict["name"] + "_Params_" + time_now + "_.csv", "w") as outfile:
writer = csv.writer(outfile)
for key, val in sorted(ParamsDict.items()):
writer.writerow([key, val])
np.savetxt("do-not-track/" + date + ParamsDict["tag"] + "_" + ParamsDict["name"] + "_FCM_" + time_now + "_.csv", FCM, delimiter = "\t")
np.savetxt("do-not-track/" + date + ParamsDict["tag"] + "_" + ParamsDict["name"] + "_Scorr_" + time_now + "_.csv", Scorr, delimiter = "\t")
# Export the time series simulation (which is a large file)
if ParamsDict["ExportSim"] == True:
# Create Time Series and save.
TSeries = np.concatenate((bold_time[Snip:].reshape(1,len(bold_time[Snip:])),TSeriesMatrix))
np.savetxt("do-not-track/" + date + ParamsDict["tag"] + "_" + ParamsDict["name"] + "_Tseries_" + time_now + "_.csv", TSeries, delimiter="\t")
return Scorr