diff --git a/multiarea_model/simulation.py b/multiarea_model/simulation.py
index a3f1e68a20272beaf7a05d2019972865a8e3c255..4d474babfc12033dd265f47a6cd41b8f72733a5c 100644
--- a/multiarea_model/simulation.py
+++ b/multiarea_model/simulation.py
@@ -622,6 +622,7 @@ def connect(simulation,
     source_area : Area instance
         Source area of the projection
     """
+
     network = simulation.network
     synapses = extract_area_dict(network.synapses,
                                  network.structure,
@@ -635,40 +636,77 @@ def connect(simulation,
                              network.structure,
                              target_area.name,
                              source_area.name)
-    for target in target_area.populations:
-        for source in source_area.populations:
-            conn_spec = {'rule': 'fixed_total_number',
-                         'N': int(synapses[target][source])}
-
-            syn_weight = {'distribution': 'normal_clipped',
-                          'mu': W[target][source],
-                          'sigma': W_sd[target][source]}
-            if target_area == source_area:
-                if 'E' in source:
-                    syn_weight.update({'low': 0.})
-                    mean_delay = network.params['delay_params']['delay_e']
-                elif 'I' in source:
-                    syn_weight.update({'high': 0.})
-                    mean_delay = network.params['delay_params']['delay_i']
-            else:
-                v = network.params['delay_params']['interarea_speed']
-                s = network.distances[target_area.name][source_area.name]
-                mean_delay = s / v
-
-            syn_delay = {'distribution': 'normal_clipped',
-                         'low': simulation.params['dt'],
-                         'mu': mean_delay,
-                         'sigma': mean_delay * network.params['delay_params']['delay_rel']}
-            syn_spec = {'weight': syn_weight,
-                        'delay': syn_delay,
-                        'model': 'static_synapse'}
-
-            if network.params['USING_NEST_3']:
+
+    if network.params['USING_NEST_3']:
+        for target in target_area.populations:
+            for source in source_area.populations:
+                conn_spec = {'rule': 'fixed_total_number',
+                             'N': int(synapses[target][source])}
+
+                if target_area == source_area:
+                    if 'E' in source:
+                        w_min = 0.
+                        w_max = np.Inf
+                        mean_delay = network.params['delay_params']['delay_e']
+                    elif 'I' in source:
+                        w_min = np.NINF
+                        w_max = 0.
+                        mean_delay = network.params['delay_params']['delay_i']
+                else:
+                    w_min = 0.
+                    w_max = np.Inf
+                    v = network.params['delay_params']['interarea_speed']
+                    s = network.distances[target_area.name][source_area.name]
+                    mean_delay = s / v
+
+                syn_spec = {
+                    'synapse_model': 'static_synapse',
+                    'weight': nest.math.redraw(
+                        nest.random.normal(mean=W[target][source],
+                                           std=W_sd[target][source]),
+                                           min=w_min,
+                                           max=w_max),
+                    'delay': nest.math.redraw(
+                        nest.random.normal(
+                            mean=mean_delay,
+                            std=(mean_delay *
+                                 network.params['delay_params']['delay_rel'])),
+                            min=simulation.params['dt'],
+                            max=np.Inf)
+                    }
+
                 nest.Connect(source_area.gids[source],
                              target_area.gids[target],
                              conn_spec,
                              syn_spec)
-            else:
+    else:
+        for target in target_area.populations:
+            for source in source_area.populations:
+                conn_spec = {'rule': 'fixed_total_number',
+                             'N': int(synapses[target][source])}
+
+                syn_weight = {'distribution': 'normal_clipped',
+                              'mu': W[target][source],
+                              'sigma': W_sd[target][source]}
+                if target_area == source_area:
+                    if 'E' in source:
+                        syn_weight.update({'low': 0.})
+                        mean_delay = network.params['delay_params']['delay_e']
+                    elif 'I' in source:
+                        syn_weight.update({'high': 0.})
+                        mean_delay = network.params['delay_params']['delay_i']
+                else:
+                    v = network.params['delay_params']['interarea_speed']
+                    s = network.distances[target_area.name][source_area.name]
+                    mean_delay = s / v
+                syn_delay = {'distribution': 'normal_clipped',
+                             'low': simulation.params['dt'],
+                             'mu': mean_delay,
+                             'sigma': mean_delay * network.params['delay_params']['delay_rel']}
+                syn_spec = {'weight': syn_weight,
+                            'delay': syn_delay,
+                            'model': 'static_synapse'}
+
                 nest.Connect(tuple(range(source_area.gids[source][0],
                                          source_area.gids[source][1] + 1)),
                              tuple(range(target_area.gids[target][0],