Add Griffiths2022 update depr files

whobpyt/depr/griffiths2022/bifurcation_analysis.py

@@ -0,0 +1,151 @@
+# -*- coding: utf-8 -*-
+"""bifurcation_analysis.ipynb
+
+Automatically generated by Colab.
+
+Original file is located at
+    https://colab.research.google.com/drive/1-N-GjmkkJdycabTr4S95Z5unyT-EyvtD
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.optimize import fsolve
+
+def h_tf(a, b, d, x):
+    return (a*x-b)/(1.0000 -np.exp(-d*(a*x-b)))
+
+def dh_tf(a, b, d, x):
+    tmp_e = np.exp(-d*(a*x-b))
+    tmp_d = 1. - np.exp(-d*(a*x-b))
+    slope_E = (a*tmp_d - (a*x-b)*d*a*tmp_e) / tmp_d**2
+    return slope_E
+
+def smooth_normalize(x):
+    return max(x, 0.000001)
+
+def derivative_orig(x, param):
+    E = x.reshape((2, 50))[0]
+    I = x.reshape((2, 50))[1]
+
+    IE = param["W_E"]*param["I_0"] + param["g_EE"]*E - param["g_IE"]*I
+    II = param["W_I"]*param["I_0"] + param["g_EI"]*E - I
+
+    rE = h_tf(param["aE"], param["bE"], param["dE"], IE)
+    rI = h_tf(param["aI"], param["bI"], param["dI"], II)
+
+    ddE = -E / param["tau_E"] + param["gamma_E"] * (1. - E) * rE
+    ddI = -I / param["tau_I"] + param["gamma_I"] * rI
+
+    return np.array([ddE, ddI]).ravel()
+
+def derivative(x, param):
+    E = x[0]
+    I = x[1]
+
+    IE = param["W_E"]*param["I_0"] + param["g_EE"]*E - param["g_IE"]*I
+    II = param["W_I"]*param["I_0"] + param["g_EI"]*E - I
+
+    rE = h_tf(param["aE"], param["bE"], param["dE"], IE)
+    rI = h_tf(param["aI"], param["bI"], param["dI"], II)
+
+    ddE = -E / param["tau_E"] + param["gamma_E"] * (1. - E) * rE
+    ddI = -I / param["tau_I"] + param["gamma_I"] * rI
+
+    return 10000.0 * np.array([ddE, ddI])
+
+def get_eig_sys(E, I, param):
+    IE = param["W_E"]*param["I_0"] + param["g_EE"]*E - param["g_IE"]*I
+    II = param["W_I"]*param["I_0"] + param["g_EI"]*E - I
+
+    rE = h_tf(param["aE"], param["bE"], param["dE"], IE)
+    rI = h_tf(param["aI"], param["bI"], param["dI"], II)
+    drEdIE = dh_tf(param["aE"], param["bE"], param["dE"], IE)
+    drIIdII = dh_tf(param["aI"], param["bI"], param["dI"], II)
+
+    A = np.zeros((2, 2))
+    A[0, 0] = -1 / param["tau_E"] - param["gamma_E"] * rE + (1 - E) * param["gamma_E"] * drEdIE * param["g_EE"]
+    A[0, 1] = -(1 - E) * param["gamma_E"] * drEdIE * param["g_IE"]
+    A[1, 0] = param["gamma_I"] * drIIdII * param["g_EI"]
+    A[1, 1] = -param["gamma_I"] * drIIdII
+
+    A = np.nan_to_num(A)
+    d, _ = np.linalg.eig(A)
+    return d
+
+def regime_search_I0(I0_rng, gEE_rng, gIE_rng, gEI_rng, param):
+    num_param = 3
+    num_trials = len(gEE_rng)
+
+    c_I0 = []
+    for I0 in I0_rng:
+        param["I_0"] = I0
+        c = []
+        for i in range(num_trials ** num_param):
+            ind_0 = i // (num_trials ** (num_param - 1))
+            ind_1 = (i % (num_trials ** (num_param - 1))) // num_trials
+            ind_2 = i % num_trials
+
+            param["g_EE"] = gEE_rng[ind_0]
+            param["g_IE"] = gIE_rng[ind_1]
+            param["g_EI"] = gEI_rng[ind_2]
+
+            initial = np.random.uniform(0., 2, [2, 50])
+            solns = []
+            for j in range(initial.shape[1]):
+                x0 = initial[:, j]
+                x0 = np.round(fsolve(lambda x: derivative(x, param), x0), decimals=4)
+                if (np.abs(derivative(x0, param)) > 1.0).sum() == 0:
+                    solns.append(tuple(x0))
+
+            good_sols = []
+            for sol in set(solns):
+                sol_good = True
+                for g_sol in good_sols:
+                    if np.sqrt(((np.array(g_sol)-np.array(sol))**2).mean()) < 1e-3:
+                        sol_good = False
+                        break
+                if sol_good:
+                    good_sols.append(sol)
+            c.append(len(good_sols))
+        c_I0.append(max(c))
+    return c_I0
+
+def regime_search_gEE(I0_rng, gEE_rng, gIE_rng, gEI_rng, param):
+    n_I0 = len(I0_rng)
+    n_gEE = len(gEE_rng)
+    n_gIE = len(gIE_rng)
+    n_gEI = len(gEI_rng)
+
+    c_gEE = []
+    for gEE in gEE_rng:
+        param["g_EE"] = gEE
+        c = []
+        for i in range(n_I0 * n_gIE * n_gEI):
+            ind_0 = i // (n_gIE * n_gEI)
+            ind_1 = (i % (n_gIE * n_gEI)) // n_gEI
+            ind_2 = i % n_gEI
+
+            param["I_0"] = I0_rng[ind_0]
+            param["g_IE"] = gIE_rng[ind_1]
+            param["g_EI"] = gEI_rng[ind_2]
+
+            initial = np.random.uniform(0., 2, [2, 50])
+            solns = []
+            for j in range(initial.shape[1]):
+                x0 = initial[:, j]
+                x0 = np.round(fsolve(lambda x: derivative(x, param), x0), decimals=4)
+                if (np.abs(derivative(x0, param)) > 1.0).sum() == 0:
+                    solns.append(tuple(x0))
+
+            good_sols = []
+            for sol in set(solns):
+                sol_good = True
+                for g_sol in good_sols:
+                    if np.sqrt(((np.array(g_sol)-np.array(sol))**2).mean()) < 1e-3:
+                        sol_good = False
+                        break
+                if sol_good:
+                    good_sols.append(sol)
+            c.append(len(good_sols))
+        c_gEE.append(max(c))
+    return c_gEE