work on 2D network of Figure 2

figures/SchueckerSchmidt2017/Fig2_EE_example.py

@@ -342,9 +342,12 @@ ax.set_xticks([])
 ax.set_yticks([])
 
 """
-Panel D: Phase space
+Panel D and G: Phase space with flux and nullclines
 """
 ax = axes['D']
+axG = axes['G']
+axG.set_ylabel(r'Excitatory rate 1 $(10^{-2}/\mathrm{s})$')
+axG.set_xlabel(r'Excitatory rate 2 $(10^{-2}/\mathrm{s})$')
 
 rates_init = np.arange(0., 200., 10.)
 
@@ -362,14 +365,15 @@ colors = ['k', myblue, myred]
 for i, netp in enumerate([network_params_base,
                           network_params_inc,
                           network_params_stab]):
-    # For the stabilized network, compute fixed point of network with
-    # increased ext. input and adapt indegree to stabilize the network
-    if i == 2:
-        fp = net.fsolve(([18., 18.]))['rates'][0][0]
-    deltaK = -1. * network_params['K'] * (161. - 160.) / fp
-    network_params_stab.update({'K': network_params['K'] + deltaK})
     net = network2D(netp)
 
+    # For the stabilized network, compute fixed point of base network
+    # and adapt indegree to stabilize the network
+    if i == 0:
+        fp = net.fsolve(([18., 18.]))['rates'][0][0]
+        deltaK = -1. * network_params['K'] * (161. - 160.) / fp
+        network_params_stab.update({'K_stable': network_params['K'] + deltaK})
+
     fp_list = []
     for init in rate_init:
         res = net.fsolve([init, init])
@@ -378,23 +382,69 @@ for i, netp in enumerate([network_params_base,
     print(fp_list)
     for fp in fp_list:
         if fp[0] < 3.:
-            ax.plot(fp[0], fp[1], 'D', color=colors[i], markersize=2)
+            if i == 1:
+                ax.plot(fp[0], fp[1], 's', color=colors[i], markersize=2 + 3)
+                axG.plot(fp[0] * 1e2, fp[1] * 1e2, 's',
+                         color=colors[i], markersize=2 + 3)
+            else:
+                ax.plot(fp[0], fp[1], 'D', color=colors[i], markersize=2)
+                axG.plot(fp[0] * 1e2, fp[1] * 1e2, 'D',
+                         color=colors[i], markersize=2)
         elif fp[0] > 3. and fp[0] < 28.8:
             ax.plot(fp[0], fp[1], '.', color=colors[i], markersize=5)
+            axG.plot(fp[0] * 1e2, fp[1] * 1e2, '.',
+                     color=colors[i], markersize=5)
+            separatrix = ([0, 2 * fp[0]], [2 * fp[1], 0])
+            if i == 0:
+                ax.plot(separatrix[0], separatrix[1], lw=5, color=colors[i])
+            else:
+                ax.plot(separatrix[0], separatrix[1], color=colors[i])
+
         elif fp[0] > 28.8:
-            ax.plot(fp[0], fp[1], '+', color=colors[i], markersize=3)
+            if i == 0:
+                ax.plot(fp[0], fp[1], '+', color=colors[i], markersize=3)
+                axG.plot(fp[0] * 1e2, fp[1] * 1e2, '+',
+                         color=colors[i], markersize=3)
 
 
 ax.set_ylabel(r'Excitatory rate 1 $(1/\mathrm{s})$')
 ax.set_xlabel(r'Excitatory rate 2 $(1/\mathrm{s})$')
-
-
-"""
-Panel G: Flux space
-"""
-ax = axes['G']
-ax.set_ylabel(r'Excitatory rate 1 $(10^{-2}/\mathrm{s})$')
-ax.set_xlabel(r'Excitatory rate 2 $(10^{-2}/\mathrm{s})$')
+y0 = 0.
+x1 = 0.015
+y1 = 0.015
+
+# vector fiels
+netp = network_params_base
+net = network2D(netp)
+range_vector = np.arange(0., 51., 10.)
+x, y, vx, vy = net.vector_field(range_vector, range_vector)
+ax.quiver(x, y, vx, vy, angles='xy', scale_units='xy', scale=8)
+
+range_vector = np.linspace(0.001, 0.015, 4)
+x, y, vx, vy = net.vector_field(range_vector, range_vector)
+axG.quiver(x * 1e2, y * 1e2, vx * 1e2, vy * 1e2,
+           angles='xy', scale_units='xy', scale=10)
+
+netp = network_params_stab
+net = network2D(netp)
+range_vector = np.linspace(0.001, 0.015, 4)
+x, y, vx, vy = net.vector_field(range_vector, range_vector)
+axG.quiver(x * 1e2, y * 1e2, vx * 1e2, vy * 1e2,
+           angles='xy', scale_units='xy', scale=10, color=myred)
+
+# nullclines
+netp = network_params_base
+net = network2D(netp)
+x0_vec = np.arange(0, 50, 0.1)
+nullcline_x0 = net.nullclines_x0(x0_vec)
+ax.plot(nullcline_x0, x0_vec, '--', color='black', label='x0')
+ax.plot(x0_vec, nullcline_x0, '--', color='black', label='x0')
+
+# set limes
+axG.set_xlim((x0 * 1e2, x1 * 1e2))
+axG.set_ylim((y0 * 1e2, y1 * 1e2))
+ax.set_xlim([-5, 50])
+ax.set_ylim([-5, 50])
 
 
 """

figures/SchueckerSchmidt2017/Fig2_EE_network.py

@@ -70,7 +70,7 @@ class network2D:
     def __init__(self, params):
         self.label = '2D'
         self.params = {'input_params': params['input_params'],
-                       'neuron_params': {'single_neuron_dict': single_neuron_dict},
+                       'neuron_params': {'single_neuron_dict': copy(single_neuron_dict)},
                        'connection_params': {'replace_cc': None,
                                              'replace_cc_input_source': None}
                        }
@@ -79,8 +79,13 @@ class network2D:
         self.structure = {'A': {'E1', 'E2'}}
         self.structure_vec = ['A-E1', 'A-E2']
         self.area_list = ['A']
-        self.K_matrix = np.array(
-            [[params['K'] / 2., params['K'] / 2., params['K']]])
+        if 'K_stable' in params.keys():
+            self.K_matrix = np.array(
+                [[params['K_stable'] / 2., params['K_stable'] / 2., params['K']]])
+        else:
+            self.K_matrix = np.array(
+                [[params['K'] / 2., params['K'] / 2., params['K']]])
+
         self.W_matrix = np.array([[params['W'], params['W'], params['W']]])
         self.J_matrix = convert_syn_weight(self.W_matrix,
                                            self.params['neuron_params']['single_neuron_dict'])
@@ -105,3 +110,47 @@ class network2D:
         result_dic = {'rates': np.array([result[0]]), 'mus': np.array(
             [mu]), 'sigmas': np.array([sigma]), 'eps': result[-1], 'time': np.array([0])}
         return result_dic
+
+    def vector_field(self, x_vec, y_vec):
+        NP = self.params['neuron_params']['single_neuron_dict']
+        vector_matrix_x = np.zeros((len(y_vec), len(x_vec)))
+        vector_matrix_y = np.zeros((len(y_vec), len(x_vec)))
+        for i, x in enumerate(y_vec):
+            for j, y in enumerate(x_vec):
+                mu, sigma = self.theory.mu_sigma([x, y])
+                new_rates = np.array(
+                    list(map(lambda mu, sigma: nu0_fb(mu, sigma,
+                                                      1.e-3 * NP['tau_m'],
+                                                      1.e-3 * NP['tau_syn_ex'],
+                                                      1.e-3 * NP['t_ref'],
+                                                      NP['V_th'] - NP['E_L'],
+                                                      NP['V_reset'] - NP['E_L']),
+                             mu, sigma)))
+                vector_matrix_x[i, j] = (new_rates[1] - y)
+                vector_matrix_y[i, j] = (new_rates[0] - x)
+        x, y = np.meshgrid(x_vec, y_vec)
+        return x, y, vector_matrix_x, vector_matrix_y
+
+    def nullclines_x0(self, x0_vec):
+        NP = self.params['neuron_params']['single_neuron_dict']
+
+        def nullcline(x0, x1):
+            rates = np.zeros(2)
+            rates[0] = x0
+            rates[1] = x1
+            mu, sigma = self.theory.mu_sigma(rates)
+            new_rates = np.array(
+                list(map(lambda mu, sigma: nu0_fb(mu, sigma,
+                                                  1.e-3 * NP['tau_m'],
+                                                  1.e-3 * NP['tau_syn_ex'],
+                                                  1.e-3 * NP['t_ref'],
+                                                  NP['V_th'] - NP['E_L'],
+                                                  NP['V_reset'] - NP['E_L']), mu, sigma)))[0]
+            return new_rates - x0
+
+        nullcline_x0 = []
+        for x0 in x0_vec:
+            result = optimize.fsolve(
+                lambda x: nullcline(x0, x), 0, full_output=1)
+            nullcline_x0.append(result[0][0])
+        return nullcline_x0

multiarea_model/theory.py

@@ -298,7 +298,6 @@ class Theory:
             'tau_m'] / C_m * self.network.add_DC_drive
         sigma2 = self.NP['tau_m'] * 1e-3 * np.dot(KJ2, rates) + sigma2_CC
         sigma = np.sqrt(sigma2)
-
         return mu, sigma
 
     def d_nu(self, mu, sigma):