diff --git a/figures/SchueckerSchmidt2017/stabilization.py b/figures/SchueckerSchmidt2017/stabilization.py
new file mode 100644
index 0000000000000000000000000000000000000000..d9cb70e4f03b4da58a50dd25faffa20923c8335e
--- /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 cf966213c4cb343287093ae683172d574041a153..dbd61a39f9a262a6935b4654bd3dccbc86ff6b1b 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 6eb06a4a009aa0234d673028da4541b37c3afe81..d947c48ba6f220eb5bc864233601903b2e559fad 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 a2b62762fd866564a847f887232eec0fff262cb3..7ff5a4ddd01e8b770a77486f9f39f60d1a9cc784 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 0000000000000000000000000000000000000000..79fef2e7b7d55d013948e5a1f6fa54f1890e4e9e
--- /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 c0a08088165583f6d6315ea211c3cba6191dc20d..a9aab52b78f240a77bc1da27620331cfc9772f08 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 5f53d99572a63eb3b787035a32be533222b998ff..77cecdaf913354bd9350663f0c04de673467e63c 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)