diff --git a/multiarea_model/data_multiarea/VisualCortex_Data.py b/multiarea_model/data_multiarea/VisualCortex_Data.py
index eb7b17b7abccaf5165ad119e73a8ca66ea5e6dfd..bc84ee5c63c2740dda3e77f87f08f5baa4a4dcab 100644
--- a/multiarea_model/data_multiarea/VisualCortex_Data.py
+++ b/multiarea_model/data_multiarea/VisualCortex_Data.py
@@ -38,8 +38,6 @@ will be stored in the corresponding dictionaries:
 6. CoCoMac data about the existence and patterns of connections
    between areas
    ----> cocomac
-6. CoCoMac data about the existence of connections between
-   areas and their laminar patterns
 7. FLN data about extrinsic connections to three areas (V1,V2,V4)
    from Markov et al. (2011)
    ---> FLN_Data
@@ -310,7 +308,6 @@ def process_raw_data():
 
     for source in dat:
         for target in dat[source]:
-            # import pdb
             source_pattern = dat[source][target][0]
             target_pattern = dat[source][target][1]
 
@@ -1054,7 +1051,6 @@ def process_raw_data():
                 FV91_source = re.sub(
                     "FVE.", "", re.sub("FVE_all.", "", FV91_key))
                 if FV91_source in area_list:
-                    # if norm > 0.0 and FV91_source in area_list :
                     if FV91_source in FLN_Data_FV91_mapped[target]:
                         FLN_Data_FV91_mapped[target][FV91_source] += overlap['all'][
                             source_key][FV91_key] / 100. * FLN_Data_FV91[target][source]
@@ -1138,7 +1134,7 @@ def process_raw_data():
                     distances_FV91 = np.append(distances_FV91, median_distance_data[
                                                target_area][source_area])
 
-    # Linear Fit to log values"
+    # Linear Fit to log values
     gradient, intercept, r_value, p_value, std_err = stats.linregress(
         distances_FV91, np.log(FLN_values_FV91))
     EDR_params = [intercept, gradient]
@@ -1282,7 +1278,6 @@ def process_raw_data():
                 FV91_source = re.sub(
                     "FVE.", "", re.sub("FVE_all.", "", FV91_key))
                 if FV91_source in area_list:
-                    # if norm > 0.0 and FV91_source in area_list :
                     if FV91_source in SLN_Data_FV91_mapped[target]:
                         SLN_Data_FV91_mapped[target][FV91_source]['S'] += (overlap['all'][
                             source_key][FV91_key] / 100. * SLN_Data_FV91[
diff --git a/multiarea_model/stabilize.py b/multiarea_model/stabilize.py
index 79fef2e7b7d55d013948e5a1f6fa54f1890e4e9e..cf891553838289bcdc2aa2c8c1cc340ec9ec379b 100644
--- a/multiarea_model/stabilize.py
+++ b/multiarea_model/stabilize.py
@@ -54,7 +54,7 @@ def stabilize(theo, theo_prime, fixed_point, a='fac_nu_ext_5E_6E', b='indegree')
     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
@@ -70,8 +70,6 @@ def stabilize(theo, theo_prime, fixed_point, a='fac_nu_ext_5E_6E', b='indegree')
     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
@@ -106,9 +104,6 @@ def S_T(theo, fixed_point):
                   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
@@ -126,27 +121,6 @@ def fixed_point_shift(a, theo, theo_prime, fixed_point):
         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
diff --git a/multiarea_model/theory.py b/multiarea_model/theory.py
index 200019923dc61d64ed6319465d967fcfa81d3c1c..ba9ec58012b6e41fe75e84fe5feab469f2dfb321 100644
--- a/multiarea_model/theory.py
+++ b/multiarea_model/theory.py
@@ -294,9 +294,6 @@ class Theory:
             rates = np.hstack((rates, self.network.params['input_params']['rate_ext']))
         else:
             rates = np.hstack((rates, np.zeros(1)))
-        # if dist:
-        #     # due to distributed weights with std = 0.1
-        #     J2[:, :7] += 0.01 * J[:, :7] * J[:, :7]
         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
diff --git a/multiarea_model/theory_helpers.py b/multiarea_model/theory_helpers.py
index 77cecdaf913354bd9350663f0c04de673467e63c..ac3832f4d508afb12fab6132f4b2e01d12db73c4 100644
--- a/multiarea_model/theory_helpers.py
+++ b/multiarea_model/theory_helpers.py
@@ -86,7 +86,6 @@ def nu_0(tau_m, tau_r, V_th, V_r, mu, sigma):
         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)
     else:
         return siegert2(tau_m, tau_r, V_th, V_r, mu, sigma)