From 3fd2a76d0cbb07cad01512198b57f6dbe241768a Mon Sep 17 00:00:00 2001
From: Maximilian Schmidt <max.schmidt@fz-juelich.de>
Date: Fri, 11 May 2018 16:22:26 +0900
Subject: [PATCH] Implement stabilization method, adapt model classes and
parameters accordingly
---
figures/SchueckerSchmidt2017/stabilization.py | 113 ++++++++++++
multiarea_model/default_params.py | 2 -
multiarea_model/multiarea_model.py | 2 +-
multiarea_model/simulation.py | 4 +-
multiarea_model/stabilize.py | 171 ++++++++++++++++++
multiarea_model/theory.py | 65 ++++---
multiarea_model/theory_helpers.py | 4 +-
7 files changed, 320 insertions(+), 41 deletions(-)
create mode 100644 figures/SchueckerSchmidt2017/stabilization.py
create mode 100644 multiarea_model/stabilize.py
diff --git a/figures/SchueckerSchmidt2017/stabilization.py b/figures/SchueckerSchmidt2017/stabilization.py
new file mode 100644
index 0000000..d9cb70e
--- /dev/null
+++ b/figures/SchueckerSchmidt2017/stabilization.py
@@ -0,0 +1,113 @@
+import copy
+import pylab as pl
+import numpy as np
+from multiarea_model import MultiAreaModel
+from multiarea_model.multiarea_helpers import create_vector_mask
+from multiarea_model.stabilize import stabilize
+import utils
+import seaborn as sns
+cp = sns.color_palette()
+
+"""
+Initialization
+"""
+conn_params = {'g': -16.,
+ 'av_indegree_V1': 3950.,
+ 'fac_nu_ext_TH': 1.2}
+input_params = {'rate_ext': 8.}
+
+network_params = {'connection_params': conn_params,
+ 'input_params': input_params}
+theory_params = {'dt': 0.01,
+ 'T': 30.}
+time = np.arange(0., theory_params['T'], theory_params['dt'])
+
+M_base = MultiAreaModel(network_params, theory=True, theory_spec=theory_params)
+
+
+c_target = copy.deepcopy(conn_params)
+c_target.update({'fac_nu_ext_5E': 1.2,
+ 'fac_nu_ext_6E': 10/3.*1.2-7/3.})
+network_params_target = {'connection_params': c_target,
+ 'input_params': input_params}
+M_target = MultiAreaModel(network_params_target, theory=True,
+ theory_spec=theory_params)
+
+THREADS = 4
+load_list = [1, 2, 3, 4]
+
+# This list defines which of the detected minima of the velocity
+# vector is identified as the unstable fixed point. It has to be
+# created manually.
+ind_list = [1, 1, 0, 1]
+
+"""
+Main loop
+"""
+data = {}
+for iteration in [1, 2, 3, 4, 5]:
+ print("Iteration {}".format(iteration))
+ if iteration == 1:
+ K_stable = None
+ else:
+ K_stable = 'iteration_{}/K_prime.npy'.format(iteration - 1)
+ conn_params = {'g': -16.,
+ 'av_indegree_V1': 3950.,
+ 'fac_nu_ext_TH': 1.2,
+ 'K_stable': K_stable}
+ network_params = {'connection_params': conn_params,
+ 'input_params': input_params}
+
+ if iteration in load_list:
+ data[iteration] = utils.load_iteration(iteration)
+ else:
+ fac_nu_ext_5E_list = np.append(np.arange(1., 1.2, 0.01), np.array([1.125]))
+
+ # Prepare base instance of the network
+ M_base = MultiAreaModel(network_params, theory=True, theory_spec=theory_params)
+ # Scan parameter space to find a good approximation of the
+ # critical parameter value where the model crosses the
+ # separatrix for the initial condition of zero rates
+ if iteration < 5: # For iteration 5, we just analyze the behavior without performing the stabilization
+ data[iteration] = utils.compute_iteration(7, fac_nu_ext_5E_list,
+ theory_params, M_base, threads=THREADS)
+ else:
+ data[iteration] = utils.compute_iteration(1, fac_nu_ext_5E_list,
+ theory_params, M_base, threads=THREADS)
+
+ if iteration != 5:
+ # Determine the transition parameter and the minima of the
+ # velocity of the trajectories (see Fig. 4 of Schuecker, Schmidt et al.)
+ (par_transition, r_low, r_high,
+ minima_low, minima_high) = utils.determine_velocity_minima(time, data[iteration])
+ # Retrieve unstable fixed point for both trajectories
+ unstable_low = r_low[:, minima_low[ind_list[iteration - 1]]]
+ unstable_high = r_high[:, minima_high[ind_list[iteration - 1]]]
+ # Simple check to guarantee that the unstable fixed point is
+ # approximately equal for both trajectories
+ assert(np.allclose(unstable_low, unstable_high, rtol=2e-1))
+
+ c = copy.deepcopy(conn_params)
+ c.update({'fac_nu_ext_5E': par_transition,
+ 'fac_nu_ext_6E': 10/3.*par_transition-7/3.})
+ network_params = {'connection_params': c,
+ 'input_params': input_params}
+ M = MultiAreaModel(network_params, theory=True,
+ theory_spec=theory_params)
+
+ K_prime = stabilize(M.theory,
+ M_target.theory,
+ unstable_low,
+ a='fac_nu_ext_5E_6E', b='indegree')
+ data[iteration]['K_prime'] = K_prime[:, :-1]
+ utils.save_iteration(iteration, data[iteration])
+
+
+fig = pl.figure()
+
+mask = create_vector_mask(M.structure, pops=['5E', '6E'])
+for iteration in [1, 2, 3, 4]:
+ pl.plot(data[iteration]['parameters'],
+ np.mean(data[iteration]['results'][:, mask, -1], axis=1), '.-')
+pl.yscale('Log')
+pl.show()
diff --git a/multiarea_model/default_params.py b/multiarea_model/default_params.py
index cf96621..dbd61a3 100644
--- a/multiarea_model/default_params.py
+++ b/multiarea_model/default_params.py
@@ -228,7 +228,6 @@ input_params = {
'dc_stimulus': False,
}
-sim_params.update({'input_params': input_params})
network_params.update({'input_params': input_params})
"""
@@ -266,7 +265,6 @@ Theory params
"""
theory_params = {'neuron_params': neuron_params,
- 'input_params': input_params,
# Initial rates can be None (start integration at
# zero rates), a numpy.ndarray defining the initial
# rates or 'random_uniform' which leads to randomly
diff --git a/multiarea_model/multiarea_model.py b/multiarea_model/multiarea_model.py
index 6eb06a4..d947c48 100644
--- a/multiarea_model/multiarea_model.py
+++ b/multiarea_model/multiarea_model.py
@@ -139,7 +139,7 @@ class MultiAreaModel:
raise NotImplementedError('Stabilization procedure has '
'to be integrated.')
elif isinstance(self.params['connection_params']['K_stable'], np.ndarray):
- raise ValueError("Not supported. Please store the "
+ raise NotImplementedError("Not supported. Please store the "
"matrix in a file and define the path to the file as "
"the parameter value.")
else: # Assume that the parameter defines a filename containing the matrix
diff --git a/multiarea_model/simulation.py b/multiarea_model/simulation.py
index a2b6276..7ff5a4d 100644
--- a/multiarea_model/simulation.py
+++ b/multiarea_model/simulation.py
@@ -220,7 +220,7 @@ class Simulation:
elif replace_non_simulated_areas == 'hom_poisson_stat':
non_simulated_cc_input = {source_area_name:
{source_pop:
- self.params['input_params']['rate_ext']
+ self.network.params['input_params']['rate_ext']
for source_pop in
self.network.structure[source_area_name]}
for source_area_name in self.network.area_list}
@@ -238,7 +238,7 @@ class Simulation:
elif replace_cc == 'hom_poisson_stat':
cc_input = {source_area_name:
{source_pop:
- self.params['input_params']['rate_ext']
+ self.network.params['input_params']['rate_ext']
for source_pop in
self.network.structure[source_area_name]}
for source_area_name in self.network.area_list}
diff --git a/multiarea_model/stabilize.py b/multiarea_model/stabilize.py
new file mode 100644
index 0000000..79fef2e
--- /dev/null
+++ b/multiarea_model/stabilize.py
@@ -0,0 +1,171 @@
+import numpy as np
+from .multiarea_helpers import create_mask, create_vector_mask
+from copy import deepcopy
+from .theory_helpers import d_nu_d_mu_fb_numeric, d_nu_d_sigma_fb_numeric
+import copy
+
+"""
+Implementation of the stabilization method of [1].
+The notation follows Eqs. (6-13) of [1].
+
+1. Schuecker J, Schmidt M, van Albada SJ, Diesmann M & Helias M (2017)
+ Fundamental Activity Constraints Lead to Specific Interpretations of the Connectome.
+ PLOS Computational Biology, 13(2).
+ [https://doi.org/10.1371/journal.pcbi.1005179](https://doi.org/10.1371/journal.pcbi.1005179)
+"""
+
+
+def stabilize(theo, theo_prime, fixed_point, a='fac_nu_ext_5E_6E', b='indegree'):
+ """
+ Implementation of the stabilization algorithm.
+
+ Parameters
+ ----------
+ theo : Instance of Theory class
+ Unperturbed network.
+ theo_prime : Instance of Theory class
+ Network perturbed by a change in the a parameter
+ fixed_point : numpy.ndarray
+ Unstable fixed point that we want to preserve.
+ a : str
+ The first parameter to be changed. Defaults to
+ 'fac_nu_ext_5E_6E' which is the relative change of the
+ external indegree onto populations 5E and 6E.
+ b : str
+ The second parameter to be changed in order to preserve the
+ location of the separatrix. Defaults to the indegrees.
+ """
+ if b != 'indegree':
+ raise NotImplementedError("Stabilizing using b = {} is not implemented.".format(b))
+
+ """
+ First calculate the change of the fixed point that, to first
+ order, is described by Eq. 6 of [1], using Eq. 8.
+ """
+ S_vector, S, T_vector, T, M = S_T(theo, fixed_point)
+ delta_bar_nu_star = fixed_point_shift(a, theo, theo_prime, fixed_point)
+ delta_nu_star = np.dot(np.linalg.inv(np.identity(M.shape[0]) - M), delta_bar_nu_star)
+
+ """
+ Next, determine the change of the parameter b that is
+ necessary to revert the change (Eq. 9).
+
+ Calculate eigen decomposition of the effective connectivity
+ matrix M
+ """
+ lambda_ev, u, v = eigen_decomp_M(M)
+
+ a_hat = np.dot(v, delta_bar_nu_star)
+ v_hat = np.dot(v, fixed_point)
+ epsilon = - 1. * a_hat / v_hat
+
+ # Calculate the single terms of the sum in Eq. (13)
+ eta_tilde = []
+ d = np.dot(delta_nu_star, delta_nu_star)
+ for l in range(epsilon.size):
+ eta_tilde.append(-1. * a_hat[l] / (1 - lambda_ev[l]) * np.dot(u[:, l], delta_nu_star) / d)
+
+ """
+ Calculate perturbation of beta (Eq. 11)
+ Only take the most critical eigendirection into account.
+ """
+ eigen_proj = np.outer(u[:, 0], v[0])
+ # fac = (theo.NP['tau_syn'] /
+ # theo.network.params['neuron_params']['single_neuron_dict']['C_m'])
+ fac = 1.
+ denom = (S * theo.network.J_matrix[:, :-1] +
+ T * theo.network.J_matrix[:, :-1]**2) * fac * theo.NP['tau_m'] * 1.e-3
+ delta_K = epsilon[0] * eigen_proj / denom
+
+ """
+ Apply three constraints:
+ 1. No inhibitory cortico-cortical connections
+ 2. No cortico-cortical connections from population 4E
+ 3. Indegree have to be > 0 -> Negative entries are set to zero.
+ """
+ index = np.zeros_like(delta_K, dtype=np.bool)
+ for area in theo.network.area_list:
+ for area2 in theo.network.area_list:
+ if area2 != area:
+ mask = create_mask(theo.network.structure,
+ target_areas=[area],
+ source_areas=[area2],
+ source_pops=['23I', '4E', '4I', '5I', '6I'])
+ index = np.logical_or(index, mask[:, :-1])
+ delta_K[index] = 0.
+ K_prime = copy.copy(theo.network.K_matrix)
+ K_prime[:, :-1] += np.real(delta_K)
+ K_prime[np.where(K_prime < 0.0)] = 0.0
+ return K_prime
+
+
+def S_T(theo, fixed_point):
+ mu, sigma = theo.mu_sigma(fixed_point)
+ S_vector, T_vector = theo.d_nu(mu, sigma)
+ S = np.array([S_vector[i] * np.ones(theo.network.K_matrix.shape[0])
+ for i in range(theo.network.K_matrix.shape[0])])
+ T = np.array([T_vector[i] * np.ones(theo.network.K_matrix.shape[0])
+ for i in range(theo.network.K_matrix.shape[0])])
+ # fac = (theo.NP['tau_syn'] /
+ # theo.network.params['neuron_params']['single_neuron_dict']['C_m']) * 1.e3
+ # import pdb; pdb.set_trace()
+ fac = 1.
+ W = theo.network.K_matrix[:, :-1] * theo.network.J_matrix[:, :-1]
+ W2 = theo.network.K_matrix[:, :-1] * theo.network.J_matrix[:, :-1]**2
+ M = (S * W * fac * theo.NP['tau_m'] * 1.e-3 +
+ T * W2 * fac ** 2 * theo.NP['tau_m'] * 1.e-3)
+ return S_vector, S, T_vector, T, M
+
+
+def fixed_point_shift(a, theo, theo_prime, fixed_point):
+ S_vector, S, T_vector, T, SJ_TJ2 = S_T(theo, fixed_point)
+ if a in ['fac_nu_ext_5E_6E']:
+ W_ext = deepcopy(theo.network.J_matrix[:, -1])
+
+ K_ext = deepcopy(theo.network.K_matrix[:, -1])
+ K_ext_prime = theo_prime.network.K_matrix[:, -1]
+ delta_Kext = K_ext_prime - K_ext
+
+ # if a == 'fac_nu_ext_5E_6E':
+ # mask = create_vector_mask(theo.network.structure, pops=['5E'])
+ # K_ext[mask] /= theo.network.params['connection_params']['fac_nu_ext_5E']
+ # delta_param = np.zeros_like(K_ext)
+ # delta_a = (np.array(theo_prime.network.params['connection_params'][
+ # 'fac_nu_ext_5E']) -
+ # np.array(theo.network.params['connection_params']['fac_nu_ext_5E']))
+ # delta_param[mask] = delta_a * theo.network.params['input_params']['rate_ext']
+
+ # mask = create_vector_mask(theo.network.structure, pops=['6E'])
+ # # in fact we realize a change in nu_ext via a change in K_ext. Here
+ # # we again shift this change to a change in the external rate.
+ # # Therefore we need to divide the indegree by the factor here.
+ # delta_a = (np.array(theo_prime.network.params['connection_params'][
+ # 'fac_nu_ext_6E']) -
+ # np.array(theo.network.params['connection_params']['fac_nu_ext_6E']))
+ # K_ext[mask] /= theo.network.params['connection_params']['fac_nu_ext_6E']
+ # delta_param[mask] = delta_a * theo.network.params['input_params']['rate_ext']
+
+ # fac = (theo.NP['tau_syn'] /
+ # theo.network.params['neuron_params']['single_neuron_dict']['C_m']) * 1.e3
+ fac = 1.
+ rate_ext = theo.network.params['input_params']['rate_ext']
+ v_mu = fac * theo.NP['tau_m'] * 1.e-3 * S_vector * delta_Kext * W_ext * rate_ext
+ v_sigma = fac ** 2 * theo.NP['tau_m'] * 1.e-3 * T_vector * delta_Kext * W_ext**2 * rate_ext
+ v = v_mu + v_sigma
+ else:
+ raise NotImplementedError('a = {} not implemented.'.format(a))
+ return v
+
+
+def eigen_decomp_M(M):
+ eig = np.linalg.eig(M)
+ evec_left = np.linalg.inv(eig[1])
+ evec_right = eig[1]
+ evals = eig[0]
+
+ index = np.argsort(np.real(evals))[::-1]
+ evals_sorted = evals[index]
+ evec_right_sorted = evec_right[:, index]
+ evec_left_sorted = evec_left[index]
+
+ return evals_sorted, evec_right_sorted, evec_left_sorted
diff --git a/multiarea_model/theory.py b/multiarea_model/theory.py
index c0a0808..a9aab52 100644
--- a/multiarea_model/theory.py
+++ b/multiarea_model/theory.py
@@ -18,7 +18,6 @@ provedure.
import json
import pprint
-import multiprocessing
import nest
import numpy as np
@@ -26,13 +25,11 @@ from copy import copy
from .default_params import nested_update, theory_params
from .default_params import check_custom_params
from dicthash import dicthash
-from functools import partial
from .multiarea_helpers import create_mask, create_vector_mask, dict_to_vector
from .theory_helpers import d_nu_d_mu_fb_numeric, d_nu_d_sigma_fb_numeric
-from .theory_helpers import nu0_fb
-class Theory():
+class Theory:
def __init__(self, network, theory_spec):
self.params = copy(theory_params)
check_custom_params(theory_spec, self.params)
@@ -48,7 +45,6 @@ class Theory():
'tau_syn': self.params['neuron_params']['single_neuron_dict']['tau_syn_ex'],
't_ref': self.params['neuron_params']['single_neuron_dict']['t_ref'],
'tau': 1.}
-
self.label = dicthash.generate_hash_from_dict({'params': self.params,
'network_label': self.network.label})
@@ -71,13 +67,13 @@ class Theory():
"""
dt = self.params['dt']
T = self.params['T']
- rate_ext = self.params['input_params']['rate_ext']
+ rate_ext = self.network.params['input_params']['rate_ext']
K = copy(self.network.K_matrix)
J = copy(self.network.J_matrix)
tau = self.NP['tau_m'] * 1e-3
dim = np.shape(K)[0]
nest.ResetKernel()
- nest.set_verbosity('M_ERROR')
+ nest.set_verbosity('M_FATAL')
nest.SetKernelStatus({'resolution': dt,
'use_wfr': False,
'print_time': False,
@@ -85,6 +81,7 @@ class Theory():
# create neurons for external drive
drive = nest.Create(
'siegert_neuron', 1, params={'rate': rate_ext, 'mean': rate_ext})
+
# create neurons representing populations
neurons = nest.Create(
'siegert_neuron', dim, params=self.NP)
@@ -152,14 +149,13 @@ class Theory():
# Loop over all iterations of different initial conditions
for nsim in range(num_iter):
- print("Iteration {}".format(nsim))
initial_rates = next(gen)
for ii in range(dim):
nest.SetStatus([neurons[ii]], {'rate': initial_rates[ii]})
# create recording device
multimeter = nest.Create('multimeter', params={'record_from':
- ['rate'], 'interval': 1.,
+ ['rate'], 'interval': dt,
'to_screen': False,
'to_file': False,
'to_memory': True})
@@ -202,7 +198,7 @@ class Theory():
self.network.structure)
elif self.network.params['connection_params']['replace_cc'] == 'hom_poisson_stat':
self.cc_input_rates = (np.ones(self.network.K_matrix.shape[0]) *
- self.network.params['input_params']['rate_ext'])
+ self.network.network.params['input_params']['rate_ext'])
for area in self.network.area_list:
area_dim = len(self.network.structure[area])
mask = create_mask(self.network.structure,
@@ -289,14 +285,14 @@ class Theory():
sigma2_CC = np.zeros_like(rates)
KJ = K * J
J2 = J * J
+ KJ2 = K * J2
if external:
- rates = np.hstack((rates, self.params['input_params']['rate_ext']))
+ rates = np.hstack((rates, self.network.params['input_params']['rate_ext']))
else:
rates = np.hstack((rates, np.zeros(self.dim_ext)))
# if dist:
# # due to distributed weights with std = 0.1
# J2[:, :7] += 0.01 * J[:, :7] * J[:, :7]
- KJ2 = K * J2
C_m = self.network.params['neuron_params']['single_neuron_dict']['C_m']
mu = self.NP['tau_m'] * 1e-3 * np.dot(KJ, rates) + mu_CC + self.NP[
'tau_m'] / C_m * self.network.add_DC_drive
@@ -305,6 +301,24 @@ class Theory():
return mu, sigma
+ def d_nu(self, mu, sigma):
+ d_nu_d_mu = np.array([d_nu_d_mu_fb_numeric(1.e-3*self.NP['tau_m'],
+ 1.e-3*self.NP['tau_syn'],
+ 1.e-3*self.NP['t_ref'],
+ self.NP['theta'],
+ self.NP['V_reset'],
+ mu[ii], sigma[ii]) for ii in range(len(mu))])
+ # Unit: 1/(mV)**2
+ d_nu_d_sigma = np.array([d_nu_d_sigma_fb_numeric(1.e-3*self.NP['tau_m'],
+ 1.e-3*self.NP['tau_syn'],
+ 1.e-3*self.NP['t_ref'],
+ self.NP['theta'],
+ self.NP['V_reset'],
+ mu[ii], sigma[ii])*1/(
+ 2. * sigma[ii])
+ for ii in range(len(mu))])
+ return d_nu_d_mu, d_nu_d_sigma
+
def stability_matrix(self, rates, matrix_filter=None,
vector_filter=None, full_output=False):
"""
@@ -350,35 +364,20 @@ class Theory():
if np.any(vector_filter is not None):
mu = mu[vector_filter]
sigma = sigma[vector_filter]
- slope = np.array([d_nu_d_mu_fb_numeric(1.e-3*self.NP['tau_m'],
- 1.e-3*self.NP['tau_syn'],
- 1.e-3*self.NP['t_ref'],
- self.NP['V_th'],
- self.NP['V_reset'],
- mu[ii], sigma[ii]) for ii in range(N.size)])
+ d_nu_d_mu, d_nu_d_sigma = self.d_nu(mu, sigma)
- # Unit: 1/(mV)**2
- slope_sigma = np.array([d_nu_d_sigma_fb_numeric(1.e-3*self.NP['tau_m'],
- 1.e-3*self.NP['tau_syn'],
- 1.e-3*self.NP['t_ref'],
- self.NP['V_th'],
- self.NP['V_reset'],
- mu[ii], sigma[ii])*1/(
- 2. * sigma[ii])
- for ii in range(N.size)])
-
- slope_matrix = np.zeros_like(J)
+ slope_mu_matrix = np.zeros_like(J)
slope_sigma_matrix = np.zeros_like(J)
for ii in range(N.size):
- slope_matrix[:, ii] = slope
- slope_sigma_matrix[:, ii] = slope_sigma
+ slope_mu_matrix[:, ii] = d_nu_d_mu
+ slope_sigma_matrix[:, ii] = d_nu_d_sigma
V = C*(1-C)
- G = (self.NP['tau_m'] * 1e-3)**2 * (slope_matrix*J +
+ G = (self.NP['tau_m'] * 1e-3)**2 * (slope_mu_matrix*J +
slope_sigma_matrix*J**2)**2
G_N = N_pre * G
M = G_N * V
if full_output:
- return M, slope, slope_sigma, M, C, V, G_N, J, N_pre
+ return M, d_nu_d_mu, d_nu_d_sigma, M, C, V, G_N, J, N_pre
else:
return M
diff --git a/multiarea_model/theory_helpers.py b/multiarea_model/theory_helpers.py
index 5f53d99..77cecda 100644
--- a/multiarea_model/theory_helpers.py
+++ b/multiarea_model/theory_helpers.py
@@ -30,7 +30,7 @@ def nu0_fb(mu, sigma, tau_m, tau_s, tau_r, V_th, V_r):
"""
Compute the stationary firing rate of a neuron with synaptic
filter of time constant tau_s driven by Gaussian white noise, from
- Fourcaud & Brunel 2002. #
+ Fourcaud & Brunel 2002.
Parameters
----------
@@ -56,7 +56,6 @@ def nu0_fb(mu, sigma, tau_m, tau_s, tau_r, V_th, V_r):
# effective reset
V_r1 = V_r + sigma * alpha / 2. * np.sqrt(tau_s / tau_m)
-
# use standard Siegert with modified threshold and reset
return nu_0(tau_m, tau_r, V_th1, V_r1, mu, sigma)
@@ -86,7 +85,6 @@ def nu_0(tau_m, tau_r, V_th, V_r, mu, sigma):
sigma : float
Variance of the input current to the neurons in mV
"""
-
if mu <= V_th + (0.95 * abs(V_th) - abs(V_th)):
# if mu <= 0.95*V_th:
return siegert1(tau_m, tau_r, V_th, V_r, mu, sigma)
--
GitLab