Last attempts of het

DynamicsAndNeuralSystems · Feb 25, 2021 · 7b1d823 · 7b1d823
1 parent 32c4cfe
commit 7b1d823
Show file tree

Hide file tree

Showing 2 changed files with 198 additions and 5 deletions.
diff --git a/Loop_sim.py b/Loop_sim.py
@@ -89,13 +89,180 @@
 ################################################################################################################################
 
 # i is PBS_ARRAY_INDEX - Allows for creation of multiple jobs 
-i = int(sys.argv[1])
-ParamsDict["G"] = np.array([i*0.05]) 
-ParamsDict["BOLD"] = True
+I = int(sys.argv[1])
 
 ParamsDict["Simul_length"] = 1.2e5
 #ParamsDict["Snip"] = 100 # Gets multiplied by 100 so this is 1 e4 (or 1%)
 
+# Obtain Het Data
+df = pd.read_csv("CortexDensitiesAlter.csv",delimiter=",")
+E_pop = df.excitatory.values
+I_pop = df.inhibitory.values
+E_mean = np.mean(E_pop)
+I_mean = np.mean(I_pop)
+
+# E_normalised is (when excluding region 7) -0.28 to 0.54
+E_normalised = (E_pop-E_mean)/E_mean
+# I_normalised is (when excluding region 7) -0.45 to 1.44
+I_normalised = (I_pop-I_mean)/I_mean
+
+
+# Gradient - Version 2 
+E_prop = E_pop / (E_pop + I_pop)
+I_prop = I_pop / (E_pop + I_pop)
+
+Mean_e_prop = np.mean(E_prop)
+Mean_i_prop = np.mean(I_prop)
+
+# -0.1 to 0.08 roughly. 
+E_prop_norm = (E_prop-Mean_e_prop)/Mean_e_prop
+# -0.4 to 0.7
+I_prop_norm = (I_prop-Mean_i_prop)/Mean_i_prop
+
+E_normalised = E_prop_norm
+I_normalised = I_prop_norm 
+
+# Simulate LCycle Het 1.2 e5 no G sweep and finer sigma. 
+# LCycle
+ParamsDict["MODEL_c_ee"] = np.array([11.0])
+ParamsDict["MODEL_c_ei"] = np.array([10.0])
+ParamsDict["MODEL_c_ie"] = np.array([10.0])
+ParamsDict["MODEL_c_ii"] = np.array([1.0])
+ParamsDict["B_e"] = np.array([1.5]) 
+ParamsDict["G"] = np.array([0.75]) 
+# Homogeneous Coupling constants
+h_ee = ParamsDict["MODEL_c_ee"] 
+h_ei = ParamsDict["MODEL_c_ei"] 
+h_ie = ParamsDict["MODEL_c_ie"] 
+h_ii = ParamsDict["MODEL_c_ii"] 
+
+# Simualte V2 het (so using proportions.)
+
+ParamsDict["Sigma_e"] =I*0.1*5
+sigma_e = ParamsDict["Sigma_e"] 
+
+ParamsDict["Sigma_i"] = I*0.1
+sigma_i = ParamsDict["Sigma_i"] 
+
+        # Heterogeneous Coupling Constants (array)
+ParamsDict["MODEL_c_ie"] = h_ie * (1 + sigma_e * E_normalised) 
+ParamsDict["MODEL_c_ee"] = h_ee  * (1 + sigma_e * E_normalised) 
+ParamsDict["MODEL_c_ii"] = h_ii  * (1 + sigma_i * I_normalised) 
+ParamsDict["MODEL_c_ei"] = h_ei  * (1 + sigma_i * I_normalised) 
+
+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=ParamsDict["B_e"],b_i=numpy.array([2.8]),tau_e=numpy.array([10.0]),
+                                    tau_i=numpy.array([65.0])) 
+ParamsDict["tag"] = "LCycle_HetV2_G" + str(ParamsDict["G"]) + "_b_e" + str(ParamsDict["B_e"]) + "_sig_e" + str(sigma_e) + "_sig_i" + str(sigma_i)
+Simul_Pipeline(ParamsDict=ParamsDict)
+
+# Simualte V3 het (so using proportions and adjusting a_e, a_i)
+ParamsDict["MODEL_c_ee"] = np.array([11.0])
+ParamsDict["MODEL_c_ei"] = np.array([10.0])
+ParamsDict["MODEL_c_ie"] = np.array([10.0])
+ParamsDict["MODEL_c_ii"] = np.array([1.0])
+
+ParamsDict["Sigma_e"] =I*0.1*5
+sigma_e = ParamsDict["Sigma_e"] 
+
+ParamsDict["Sigma_i"] = I*0.1
+sigma_i = ParamsDict["Sigma_i"] 
+
+        # Heterogeneous Coupling Constants (array)
+ParamsDict["MODEL_a_e"] = numpy.array([1.0]) * (1 + sigma_e * E_normalised) 
+ParamsDict["MODEL_a_i"] = numpy.array([1.0])  * (1 + sigma_i * I_normalised) 
+
+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=ParamsDict["MODEL_a_e"],a_i=ParamsDict["MODEL_a_i"],b_e=ParamsDict["B_e"],b_i=numpy.array([2.8]),tau_e=numpy.array([10.0]),
+                                    tau_i=numpy.array([65.0])) 
+ParamsDict["tag"] = "LCycle_HetV3_G" + str(ParamsDict["G"]) + "_b_e" + str(ParamsDict["B_e"]) + "_sig_e" + str(sigma_e) + "_sig_i" + str(sigma_i)
+Simul_Pipeline(ParamsDict=ParamsDict)
+
+# Finally try Het when I remove Nodes 26,19,3 (of the original order)
+
+ParamsDict["REMOVE"] = [26,19,7,3]
+
+df = pd.read_csv("CortexDensitiesAlter.csv",delimiter=",")
+E_pop = df.excitatory.values
+I_pop = df.inhibitory.values
+
+E_pop = np.delete(E_pop,[26,19,3])
+I_pop = np.delete(I_pop,[26,19,3])
+E_mean = np.mean(E_pop)
+I_mean = np.mean(I_pop)
+
+# E_normalised is (when excluding region 7,26,19,3) -0.27 -> 0.41
+E_normalised = (E_pop-E_mean)/E_mean
+# I_normalised is (when excluding region 7,26,19,3) -0.4 -> 0.63
+I_normalised = (I_pop-I_mean)/I_mean
+
+# With Alt LCycle B_e = 1.5
+# Simulate LCycle Het 1.2 e5 no G sweep and finer sigma. 
+# LCycle
+ParamsDict["MODEL_c_ee"] = np.array([11.0])
+ParamsDict["MODEL_c_ei"] = np.array([10.0])
+ParamsDict["MODEL_c_ie"] = np.array([10.0])
+ParamsDict["MODEL_c_ii"] = np.array([1.0])
+ParamsDict["B_e"] = np.array([1.5]) 
+ParamsDict["G"] = np.array([0.75]) 
+# Homogeneous Coupling constants
+h_ee = ParamsDict["MODEL_c_ee"] 
+h_ei = ParamsDict["MODEL_c_ei"] 
+h_ie = ParamsDict["MODEL_c_ie"] 
+h_ii = ParamsDict["MODEL_c_ii"] 
+
+ParamsDict["Sigma_e"] =I*0.1
+sigma_e = ParamsDict["Sigma_e"] 
+
+ParamsDict["Sigma_i"] = I*0.1
+sigma_i = ParamsDict["Sigma_i"] 
+
+        # Heterogeneous Coupling Constants (array)
+ParamsDict["MODEL_c_ie"] = h_ie * (1 + sigma_e * E_normalised) 
+ParamsDict["MODEL_c_ee"] = h_ee  * (1 + sigma_e * E_normalised) 
+ParamsDict["MODEL_c_ii"] = h_ii  * (1 + sigma_i * I_normalised) 
+ParamsDict["MODEL_c_ei"] = h_ei  * (1 + sigma_i * I_normalised) 
+
+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=ParamsDict["B_e"],b_i=numpy.array([2.8]),tau_e=numpy.array([10.0]),
+                                    tau_i=numpy.array([65.0])) 
+ParamsDict["tag"] = "LCycleAlt_Het_Removed_G" + str(ParamsDict["G"]) + "_b_e" + str(ParamsDict["B_e"]) + "_sig_e" + str(sigma_e) + "_sig_i" + str(sigma_i)
+Simul_Pipeline(ParamsDict=ParamsDict)
+
+# With Best LCycle B_e = 0.8
+# Simulate LCycle Het 1.2 e5 no G sweep and finer sigma. 
+# LCycle
+ParamsDict["MODEL_c_ee"] = np.array([11.0])
+ParamsDict["MODEL_c_ei"] = np.array([10.0])
+ParamsDict["MODEL_c_ie"] = np.array([10.0])
+ParamsDict["MODEL_c_ii"] = np.array([1.0])
+ParamsDict["B_e"] = np.array([0.8]) 
+ParamsDict["G"] = np.array([0.5]) 
+# Homogeneous Coupling constants
+h_ee = ParamsDict["MODEL_c_ee"] 
+h_ei = ParamsDict["MODEL_c_ei"] 
+h_ie = ParamsDict["MODEL_c_ie"] 
+h_ii = ParamsDict["MODEL_c_ii"] 
+
+ParamsDict["Sigma_e"] =I*0.1
+sigma_e = ParamsDict["Sigma_e"] 
+
+ParamsDict["Sigma_i"] = I*0.1
+sigma_i = ParamsDict["Sigma_i"] 
+
+        # Heterogeneous Coupling Constants (array)
+ParamsDict["MODEL_c_ie"] = h_ie * (1 + sigma_e * E_normalised) 
+ParamsDict["MODEL_c_ee"] = h_ee  * (1 + sigma_e * E_normalised) 
+ParamsDict["MODEL_c_ii"] = h_ii  * (1 + sigma_i * I_normalised) 
+ParamsDict["MODEL_c_ei"] = h_ei  * (1 + sigma_i * I_normalised) 
+
+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=ParamsDict["B_e"],b_i=numpy.array([2.8]),tau_e=numpy.array([10.0]),
+                                    tau_i=numpy.array([65.0])) 
+ParamsDict["tag"] = "LCycleB_Het_Removed_G" + str(ParamsDict["G"]) + "_b_e" + str(ParamsDict["B_e"]) + "_sig_e" + str(sigma_e) + "_sig_i" + str(sigma_i)
+Simul_Pipeline(ParamsDict=ParamsDict)
+
+'''
 # B_e: 0 -> 3
 for b_e in np.arange(start=0,stop=3,step=0.1): 
     # LCycle
@@ -139,7 +306,7 @@
                                             tau_i=numpy.array([10.0])) 
     ParamsDict["tag"] = "Hysteresis_BOLD_G"+ str(ParamsDict["G"]) + "_b_e" + str(ParamsDict["B_e"]) 
     Simul_Pipeline(ParamsDict=ParamsDict)
-
+'''
 
 '''
 # LCycle-RHopf = LCycleCut
@@ -155,7 +322,6 @@
 ParamsDict["tag"] = "LCycleCut_i" + str(i) + "Length_e5"
 Simul_Pipeline(ParamsDict=ParamsDict)
 '''
-
 '''
 # Gradient - Version 2 
 E_prop = E_pop / (E_pop + I_pop)

diff --git a/loop-yoss.csh b/loop-yoss.csh
@@ -0,0 +1,27 @@
+# This is the job scripts file to run jobs on the super cluster.  For Array job on Yossarian. 
+
+#!/bin/csh
+#PBS -N loop-yoss
+#PBS -q yossarian
+#PBS -l select=1:ncpus=1:mem=4GB
+# Minimum acceptable walltime: day-hours:mins:secs
+#PBS -l walltime=48:00:00
+# Send e-mail if job ends or aborts
+#PBS -m bea
+#PBS -M psiu5120@uni.sydney.edu.au
+#PBS -V
+
+# -J Specifies Job Array with Indices 
+#PBS -J 0-10
+
+# Show the host on which the job ran
+hostname
+module load Anaconda3-5.1.0
+
+# Load TVB virtual env.
+source ~/TVB/bin/activate.csh
+
+# Working directory is where I ran qsub
+cd "$PBS_O_WORKDIR"
+python Loop_sim.py $PBS_ARRAY_INDEX
+exit