Repository: Hjorthmedh/Snudda (language: python)

Issue:
Pr connectivity
Added preliminary docstring for save_network_activity.py


Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index ca01071a3..6c0c325e2 100644
--- a/README.md
+++ b/README.md
@@ -7,7 +7,7 @@ Human Brain Project
 hjorth@kth.se
 
 ## Funding
-Horizon 2020 Framework Programme (785907, HBP SGA2); Horizon 2020 Framework Programme (945539, HBP SGA3); Vetenskapsrådet (VR-M-2017-02806, VR-M-2020-01652); Swedish e-science Research Center (SeRC); KTH Digital Futures.
+Horizon 2020 Framework Programme (785907, HBP SGA2); Horizon 2020 Framework Programme (945539, HBP SGA3); Vetenskapsrådet (VR-M-2017-02806, VR-M-2020-01652); Swedish e-science Research Center (SeRC); KTH Digital Futures. The computations are enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at PDC KTH partially funded by the Swedish Research Council through grant agreement no. 2018-05973. We acknowledge the use of Fenix Infrastructure resources, which are partially funded from the European Union's Horizon 2020 research and innovation programme through the ICEI project under the grant agreement No. 800858.
 
 ## Citation
 Please cite the first paper for the general Snudda network creation and simulation methods, and the second paper for the Striatal microcircutiry model.
diff --git a/requirements.txt b/requirements.txt
index 3d5ea387b..5d9f0d606 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,20 +1,20 @@
 bluepyopt >=1.11.7
-h5py >=3.1.0
+h5py >=3.2.1
 ipyparallel>=6.3.0
 matplotlib>=3.3.4
 mpi4py>=3.0.3
-numpy>=1.20.1
-scipy>=1.6.1
+numpy>=1.20.2
+scipy>=1.6.3
 sonata>=0.0.2
 pyzmq>=22.0.3
-numexpr
+numexpr>=2.7.3
 # allensdk  # Optional, only needed for downloading meshes from Allen Institute
 argparse
 snudda
-NEURON
-pyswarms
+NEURON>=7.8.2
+pyswarms>=1.3.0
 setuptools
-psutil
+psutil >= 5.8.0
 cython  # Needed for compiling NEURON
-numba>=0.53.1  # Optimisation
+numba>=0.55.1  # Optimisation
 
diff --git a/setup.py b/setup.py
index da1162381..79283a320 100644
--- a/setup.py
+++ b/setup.py
@@ -20,19 +20,20 @@
 else:
     print(f"READTHEDOCS = {os.environ.get('READTHEDOCS')}") 
     install_requires = [
-        "bluepyopt>=1.9.126",
-        "h5py>=3.1.0",
+        "bluepyopt>=1.11.7",
+        "h5py>=3.2.1",
         "ipyparallel>=6.3.0",
         "matplotlib>=3.3.4",
         "mpi4py>=3.0.3",
-        "numpy>=1.20.1",
-        "scipy>=1.6.1",
+        "numpy>=1.20.2",
+        "scipy>=1.6.3",
         "sonata>=0.0.2",
         "pyzmq>=22.0.3",
         "setuptools",
         "psutil",
         "argparse",
-        "numexpr"
+        "numexpr>=2.7.3",
+        "numba>=0.53.1"
     ]
     
 setuptools.setup(
diff --git a/snudda/__init__.py b/snudda/__init__.py
index 01e654492..e13852f16 100644
--- a/snudda/__init__.py
+++ b/snudda/__init__.py
@@ -1,6 +1,6 @@
 from .core import Snudda
 
-__version__ = "1.2.16"
+__version__ = "1.2.17"
 
 from .init import SnuddaInit
 from .place import SnuddaPlace
diff --git a/snudda/cli.py b/snudda/cli.py
index 0b06992fd..508377091 100644
--- a/snudda/cli.py
+++ b/snudda/cli.py
@@ -125,6 +125,7 @@ def snudda_cli():
     simulate_parser.add_argument("--profile", help="Run python cProfile", action="store_true")
     simulate_parser.add_argument("--verbose", action="store_true")
     simulate_parser.add_argument("--exportCoreNeuron", action="store_true")
+    simulate_parser.add_argument("--recordALL", dest="record_all", type=str, default=None)
 
     export_parser = sub_parsers.add_parser("export")
     export_parser.add_argument("path", help="Location of network")
diff --git a/snudda/core.py b/snudda/core.py
index d9c2023d5..5b5f309ad 100755
--- a/snudda/core.py
+++ b/snudda/core.py
@@ -47,6 +47,7 @@
 
 import pkg_resources
 import json
+import numpy as np
 
 from snudda.utils import snudda_path
 from snudda.utils.snudda_path import snudda_isfile
@@ -562,6 +563,7 @@ def simulate(self, args):
 
                 sim = SnuddaSimulateNeuromodulation(network_file=network_file,
                                                     input_file=input_file,
+                                                    output_file=output_file,
                                                     disable_gap_junctions=disable_gj,
                                                     log_file=log_file,
                                                     verbose=args.verbose)
@@ -576,6 +578,7 @@ def simulate(self, args):
 
                 sim = SnuddaSimulateNeuromodulationSynapse(network_file=network_file,
                                                            input_file=input_file,
+                                                           output_file=output_file,
                                                            disable_gap_junctions=disable_gj,
                                                            log_file=log_file,
                                                            neuromodulator_description=neuromod_dict)
@@ -590,6 +593,7 @@ def simulate(self, args):
             # Simulate is deterministic, no random seed.
             sim = SnuddaSimulate(network_file=network_file,
                                  input_file=input_file,
+                                 output_file=output_file,
                                  disable_gap_junctions=disable_gj,
                                  log_file=log_file,
                                  verbose=args.verbose)
@@ -599,9 +603,15 @@ def simulate(self, args):
         sim.check_memory_status()
 
         if args.record_volt:
-            sim.add_recording(side_len=None)  # Side len let you record from a subset
+            # sim.add_volt_recording_all()
+            sim.add_volt_recording_soma()
             # sim.addRecordingOfType("dSPN",5) # Side len let you record from a subset
 
+        if args.record_all:
+            record_cell_id = np.array([int(x) for x in args.record_all.split(",")])
+            sim.add_volt_recording_all(cell_id=record_cell_id)
+            sim.add_synapse_current_recording_all(record_cell_id)
+
         t_sim = args.time * 1000  # Convert from s to ms for Neuron simulator
 
         if args.exportCoreNeuron:
@@ -613,7 +623,7 @@ def simulate(self, args):
         sim.run(t_sim)  # In milliseconds
 
         print("Simulation done, saving output")
-        sim.write_output(output_file)
+        sim.write_output()
 
         stop = timeit.default_timer()
         if sim.pc.id() == 0:
diff --git a/snudda/detect/detect.py b/snudda/detect/detect.py
index 2e42bcc6f..461b7026e 100644
--- a/snudda/detect/detect.py
+++ b/snudda/detect/detect.py
@@ -28,6 +28,7 @@
 
 from numba import jit
 
+from snudda.detect.projection_detection import ProjectionDetection
 from snudda.neurons.neuron_morphology import NeuronMorphology
 from snudda.neurons.neuron_prototype import NeuronPrototype
 
@@ -209,7 +210,7 @@ def __init__(self,
         # 5: hyperVoxelID, 6: channelModelID,
         # 7: sourceAxonSomaDist (not SI scaled 1e6, micrometers),
         # 8: destDendSomaDist (not SI scaled 1e6, micrometers)
-        # 9: destSegID, 10: destSegX (int 0 - 1000, SONATA wants float 0.0-1.0)
+        # 9: destSecID, 10: destSecX (int 0 - 1000, SONATA wants float 0.0-1.0)
         # 11: conductance (int, not SI scaled 1e12, in pS)
         # 12: parameterID
         #
@@ -221,7 +222,7 @@ def __init__(self,
         self.hyper_voxel_synapses = None
 
         # Columns in hyperVoxelGapJunctions
-        # 0: sourceCellID, 1: destCellID, 2: sourceSegID, 3: destSegID,
+        # 0: sourceCellID, 1: destCellID, 2: sourceSecID, 3: destSecID,
         # 4: sourceSegX, 5: destSegX, 6: voxelX, 7: voxelY, 8: voxelZ,
         # 9: hyperVoxelID, 10: conductance (integer, in pS)
         self.hyper_voxel_gap_junctions = None
@@ -250,16 +251,6 @@ def __init__(self,
         self.max_synapses = 2000000
         self.max_gap_junctions = 100000
 
-        # We have to dynamically create this lookup
-        # self.synapseTypeLookup = { 1 : "GABA",
-        #                           2 : "AMPA_NMDA",
-        #                           3 : "GapJunction",
-        #                           4 : "ACh",
-        #                           5 : "NO"}
-        #
-        # self.synapseTypeReverseLookup = \
-        #    {v: k for k, v in self.synapseTypeLookup.items()}
-
         self.connectivity_distributions = dict([])
         # self.connectivityDistributionsGJ = dict([])
         self.next_channel_model_id = 10
@@ -277,6 +268,8 @@ def __init__(self,
         # Read positions
         self.read_neuron_positions(position_file)
 
+        self.projection_detection = None  # Helper class for handling projections between structures
+
     def detect(self, restart_detection_flag=True, rc=None):
 
         """
@@ -328,6 +321,10 @@ def detect(self, restart_detection_flag=True, rc=None):
             # (can be more than one)
             self.distribute_neurons_parallel(d_view=d_view)
 
+            # We also need to start the projection code
+            self.projection_detection = ProjectionDetection(snudda_detect=self, role=self.role, rc=self.rc)
+            self.projection_detection.find_neurons_projections_in_hyper_voxels()
+
             if d_view is not None:
                 self.parallel_process_hyper_voxels(rc=self.rc, d_view=d_view)
 
@@ -974,7 +971,7 @@ def setup_hyper_voxel(self, hyper_voxel_origo, hyper_voxel_id):
         random_seed = self.hyper_voxels[hyper_voxel_id]["randomSeed"]
         self.hyper_voxel_rng = np.random.default_rng(random_seed)
 
-        self.hyper_voxel_coords[hyper_voxel_id] = hyper_voxel_origo  # Used???
+        self.hyper_voxel_coords[hyper_voxel_id] = hyper_voxel_origo
 
         self.hyper_voxel_origo = hyper_voxel_origo
         self.hyper_voxel_id = hyper_voxel_id
@@ -1097,7 +1094,7 @@ def detect_synapses(self):
             # Maybe make dendrite loop outer, since it has more variables?
             # speedup??
             for (ax_id, ax_dist) in zip(axon_id_list, axon_dist):
-                for (d_id, d_seg_id, d_seg_x, d_dist) \
+                for (d_id, d_sec_id, d_sec_x, d_dist) \
                         in zip(dend_id_list, dend_sec_id, dend_sec_x, dend_dist):
 
                     if ax_id == d_id:
@@ -1122,37 +1119,83 @@ def detect_synapses(self):
                             mean_synapse_cond, std_synapse_cond = con_dict[con_type]["conductance"]
                             channel_model_id = con_dict[con_type]["channelModelID"]
 
-                            # We can not do pruning at this stage, since we only see
-                            # synapses within hyper voxel, and pruning depends on
-                            # all synapses between two connected cells.
+                            # Should we add just one synapse, or a cluster of synapses
+                            cluster_size = con_dict[con_type]["clusterSize"]
+                            if cluster_size > 1:
+                                cluster_spread = con_dict[con_type]["clusterSpread"]
+
+                                # This uses clone in neuron_prototype which should be cached
+                                neuron = self.load_neuron(self.neurons[d_id])
+
+                                try:
+                                    cluster_sec_x, syn_coords, soma_dist \
+                                        = neuron.cluster_synapses(sec_id=d_sec_id, sec_x=d_sec_x,
+                                                                  count=cluster_size, distance=cluster_spread,
+                                                                  rng=self.hyper_voxel_rng)
+                                except:
+                                    import traceback
+                                    tstr = traceback.format_exc()
+                                    print(tstr)
+                                    import pdb
+                                    pdb.set_trace()
+
+                                if self.hyper_voxel_synapse_ctr + cluster_size >= self.max_synapses:
+                                    self.resize_hyper_voxel_synapses_matrix()
+
+                                cluster_cond = self.hyper_voxel_rng.lognormal(synapse_mu, synapse_sigma, cluster_size)
+                                cluster_cond = np.maximum(cluster_cond, mean_synapse_cond * 0.1)
+                                cluster_param_id = self.hyper_voxel_rng.integers(1000000, size=cluster_size)
+
+                                # We need to convert coords to hyper voxel coords, to fit with other coords
+                                coords_all = np.round((syn_coords - self.hyper_voxel_origo) / self.voxel_size).astype(int)
+
+                                for d_sec_x, x, y, z, d_dist, cond, param_id \
+                                    in zip(cluster_sec_x, coords_all[:, 0], coords_all[:, 1], coords_all[:, 2],
+                                           soma_dist, cluster_cond, cluster_param_id):
+
+                                    assert cond > 0, f"Conductance should be larger than 0. cond = {cond}"
+
+                                    self.hyper_voxel_synapses[self.hyper_voxel_synapse_ctr, :] = \
+                                        [ax_id, d_id, x, y, z, self.hyper_voxel_id, channel_model_id,
+                                         ax_dist, d_dist, d_sec_id, d_sec_x * 1000, cond * 1e12, param_id]
 
-                            # Do we have enough space allocated?
-                            if self.hyper_voxel_synapse_ctr >= self.max_synapses:
-                                self.resize_hyper_voxel_synapses_matrix()
+                                    # !!! OBS, dSegX is a value between 0 and 1, multiplied by 1000
+                                    # need to divide by 1000 later
 
-                            # Synapse conductance varies between synapses
-                            # cond = self.hyper_voxel_rng.normal(mean_synapse_cond, std_synapse_cond)
+                                    self.hyper_voxel_synapse_ctr += 1
 
-                            # lognormal distribution -- https://www.nature.com/articles/nrn3687
-                            # https://en.wikipedia.org/wiki/Log-normal_distribution
-                            cond = self.hyper_voxel_rng.lognormal(synapse_mu, synapse_sigma)
+                            else:
+                                # We can not do pruning at this stage, since we only see
+                                # synapses within hyper voxel, and pruning depends on
+                                # all synapses between two connected cells.
 
-                            # Need to make sure the conductance is not negative,
-                            # set lower cap at 10% of mean value
-                            cond = np.maximum(cond, mean_synapse_cond * 0.1)
-                            assert cond > 0, f"Conductance should be larger than 0. cond = {cond}"
+                                # Do we have enough space allocated?
+                                if self.hyper_voxel_synapse_ctr >= self.max_synapses:
+                                    self.resize_hyper_voxel_synapses_matrix()
 
-                            param_id = self.hyper_voxel_rng.integers(1000000)
+                                # Synapse conductance varies between synapses
+                                # cond = self.hyper_voxel_rng.normal(mean_synapse_cond, std_synapse_cond)
 
-                            # Add synapse
-                            self.hyper_voxel_synapses[self.hyper_voxel_synapse_ctr, :] = \
-                                [ax_id, d_id, x, y, z, self.hyper_voxel_id, channel_model_id,
-                                 ax_dist, d_dist, d_seg_id, d_seg_x * 1000, cond * 1e12, param_id]
+                                # lognormal distribution -- https://www.nature.com/articles/nrn3687
+                                # https://en.wikipedia.org/wiki/Log-normal_distribution
+                                cond = self.hyper_voxel_rng.lognormal(synapse_mu, synapse_sigma)
 
-                            # !!! OBS, dSegX is a value between 0 and 1, multiplied by 1000
-                            # need to divide by 1000 later
+                                # Need to make sure the conductance is not negative,
+                                # set lower cap at 10% of mean value
+                                cond = np.maximum(cond, mean_synapse_cond * 0.1)
+                                assert cond > 0, f"Conductance should be larger than 0. cond = {cond}"
 
-                            self.hyper_voxel_synapse_ctr += 1
+                                param_id = self.hyper_voxel_rng.integers(1000000)
+
+                                # Add synapse
+                                self.hyper_voxel_synapses[self.hyper_voxel_synapse_ctr, :] = \
+                                    [ax_id, d_id, x, y, z, self.hyper_voxel_id, channel_model_id,
+                                     ax_dist, d_dist, d_sec_id, d_sec_x * 1000, cond * 1e12, param_id]
+
+                                # !!! OBS, dSegX is a value between 0 and 1, multiplied by 1000
+                                # need to divide by 1000 later
+
+                                self.hyper_voxel_synapse_ctr += 1
 
         # Sort the synapses (note sortIdx will not contain the empty rows
         # at the end.
@@ -1163,8 +1206,7 @@ def detect_synapses(self):
         # basically how many voxel steps do we need to take to go from
         # simulationOrigo to hyperVoxelOrigo (those were not included, so add them)
         hyper_voxel_offset = np.round((self.hyper_voxel_origo - self.simulation_origo)
-                                      / self.hyper_voxel_width).astype(int) \
-                             * self.hyper_voxel_size
+                                      / self.hyper_voxel_width).astype(int) * self.hyper_voxel_size
 
         # Just a double check...
         assert self.hyper_voxel_id_lookup[int(np.round(hyper_voxel_offset[0] / self.hyper_voxel_size))][
@@ -1236,6 +1278,10 @@ def place_synapses_no_axon(self, hyper_id, voxel_space, voxel_space_ctr,
             # No neurons without axons
             return
 
+        # TODO: We need to update the code here to handle projection touch detection also
+        #       for neurons with "probability cloud axons" that should not be rotated
+        #       with the neuron rotation.
+
         for na_neuron in no_axon_neurons:
 
             # There are two types of axon density specified
@@ -1326,7 +1372,7 @@ def place_synapses_no_axon(self, hyper_id, voxel_space, voxel_space_ctr,
 
         end_time = timeit.default_timer()
 
-        self.write_log(f"place_synapses_no_axon_sphere: {end_time - start_time:.1f} s, hyper_id: {hyper_id}")
+        self.write_log(f"place_synapses_no_axon: {end_time - start_time:.1f} s, hyper_id: {hyper_id}")
 
     ############################################################################
 
@@ -1956,6 +2002,12 @@ def read_prototypes(self, config_file=None, axon_stump_id_flag=False):
                 sigma = np.sqrt(np.log(1 + std_cond ** 2 / mean_cond ** 2))
                 con_def[key]["lognormal_mu_sigma"] = [mu, sigma]
 
+                if "clusterSize" not in con_def[key]:
+                    con_def[key]["clusterSize"] = 1
+
+                if "clusterSpread" not in con_def[key]:
+                    con_def[key]["clusterSpread"] = 20e-3
+
             self.connectivity_distributions[pre_type, post_type] = con_def
 
     ############################################################################
@@ -2399,6 +2451,8 @@ def distribute_neurons(self, neuron_idx=None, distribution_seeds=None, min_coord
 
                 elif neuron.axon_density_type == "xyz":
 
+                    # TODO: Maybe replace random points by a grid for this test step?
+
                     rng = np.random.default_rng(d_seed)
 
                     # Estimate how many points we need to randomly place
@@ -2441,6 +2495,9 @@ def distribute_neurons(self, neuron_idx=None, distribution_seeds=None, min_coord
                                    f"{neuron.axon_density_type}", is_error=True)
                     assert False, f"No axon for {neuron.name}"
 
+                # TODO: We need to add the neurons that have touch detection projection also here
+                #       to the list of hyper voxels the neuron belongs to
+
                 # Find unique hyper voxel coordinates
                 h_loc = np.unique(np.concatenate([axon_loc, dend_loc]), axis=0).astype(int)
 
@@ -2497,10 +2554,10 @@ def setup_parallel(self, d_view=None):
         """ Prepares workers for parallel execution if d_view is not None. """
 
         assert self.role == "master", \
-            "setupParallel: Should only be called by master node"
+            "setup_parallel: Should only be called by master node"
 
         if d_view is None:
-            self.write_log("setupParallel called without dView, aborting.")
+            self.write_log("setup_parallel called without dView, aborting.")
             return
 
         if self.workers_initialised:
@@ -3207,6 +3264,18 @@ def process_hyper_voxel(self, hyper_id):
             # This detects the synapses where we use a density distribution for axons
             # self.detectSynapsesNoAxonSLOW (hyperID) # --replaced by placeSynapseNoAxon
 
+            # Finally this adds axon voxels for projections comming from other structures using projection maps
+            try:
+                self.projection_detection.voxelise_projections()
+            except:
+                # !!! TODO remove this bit of logging code
+                import traceback
+                t_str = traceback.format_exc()
+                self.write_log(t_str, is_error=True)
+                print(t_str)
+                import pdb
+                pdb.set_trace()
+
             # The normal voxel synapse detection
             self.detect_synapses()
 
diff --git a/snudda/detect/project.py b/snudda/detect/project.py
index 08f72d413..2d517ae2f 100644
--- a/snudda/detect/project.py
+++ b/snudda/detect/project.py
@@ -158,7 +158,7 @@ def connect_projection_helper(self, pre_neuron_type, post_neuron_type, connectio
 
         for connection_type, con_info in connection_info.items():
 
-            if "projectionFile" not in con_info:
+            if "projectionFile" not in con_info or "projectionName" not in con_info:
                 # Not a projection, skipping
                 continue
 
diff --git a/snudda/detect/projection_detection.py b/snudda/detect/projection_detection.py
new file mode 100644
index 000000000..b5220de4d
--- /dev/null
+++ b/snudda/detect/projection_detection.py
@@ -0,0 +1,372 @@
+import numpy as np
+import json
+from collections import OrderedDict
+from scipy.interpolate import griddata
+
+
+class ProjectionDetection:
+
+    def __init__(self, snudda_detect, rc=None, role=None):
+
+        """
+
+            Args:
+                snudda_detect : SnuddaDetect object
+                rc (default None): Ipyparallel remote client
+
+        """
+
+        self.snudda_detect = snudda_detect
+        self.rc = rc
+        self.d_view = None
+        self.workers_initialised = False
+
+        self.hyper_voxel_projections = dict()
+        self.projections = dict()
+
+        if not role:
+            self.role = "master"
+        else:
+            self.role = role
+
+        self.parse_config()
+
+        if self.rc is not None:
+            self.d_view = self.rc.direct_view(targets='all')
+            self.setup_parallel()
+
+    def setup_parallel(self):
+
+        if self.d_view is None:
+            return
+
+        """ Prepares workers for parallel execution if d_view is not None. """
+
+        assert self.role == "master", \
+            "setup_parallel: Should only be called by master node"
+
+        if self.workers_initialised:
+            self.write_log("Workers already initialised.")
+            return
+
+        with self.d_view.sync_imports():
+            from snudda.detect.projection_detection import ProjectionDetection
+
+        cmd_str = "spd = ProjectionDetection(snudda_detect=sd, role='worker')"
+        self.d_view.execute(cmd_str, block=True)
+        cmd_str2 = "sd.projection_detection = spd"
+        self.d_view.execute(cmd_str2, block=True)
+
+        self.write_log(f"Workers set up.")
+
+    def sync_projection_info(self):
+
+        if self.role == "master" and self.d_view is not None:
+
+            self.write_log("Synchronising projection info to workers.")
+            for var_name in ["projections", "hyper_voxel_projections"]:
+                self.d_view.push({"tmp_val": eval(f"self.{var_name}")}, block=True)
+                cmd_str = f"spd.{var_name} = tmp_val"
+                self.d_view.execute(cmd_str, block=True)
+
+    def write_log(self, *args, **kvargs):
+        # Reuse the log files for snudda_detect
+        self.snudda_detect.write_log(*args, **kvargs)
+
+    def find_neurons_projections_in_hyper_voxels(self):
+
+        """ For each hyper voxel, list which neurons project to that hypervoxel"""
+
+        if self.role != "master":
+            # Master node sets this up
+            return
+
+        try:
+            # In case the user has been naughty and added extra projections in the code.
+            self.sync_projection_info()
+
+            ss = np.random.SeedSequence(self.snudda_detect.random_seed + 202222)
+            all_seeds = ss.generate_state(len(self.snudda_detect.neurons))
+            all_neuron_id = sorted([n["neuronID"] for n in self.snudda_detect.neurons])
+
+            seed_lookup = dict()
+            for s, n in zip(all_seeds, all_neuron_id):
+                seed_lookup[n] = s
+
+            proj_list = []
+
+            self.write_log("!! about to loop projections")
+
+            for neuron_type in self.projections:
+
+                neuron_id_list = self.get_neurons_of_type(neuron_type)
+                for neuron_id in neuron_id_list:
+                    proj_list.append((neuron_id, neuron_type, seed_lookup[neuron_id]))
+
+            self.hyper_voxel_projections = dict()
+
+            if self.d_view is not None:
+                self.d_view.scatter("proj_list", proj_list, block=True)
+
+                cmd_str = "neuron_hv_list = spd.find_hyper_voxel_helper_parallel(proj_info=proj_list)"
+                self.d_view.execute(cmd_str, block=True)
+
+                neuron_hv_list = self.d_view.gather("neuron_hv_list", block=True)  # List of list of tuples
+
+                for n_hv in neuron_hv_list:
+                    neuron_id = n_hv[0]
+                    neuron_hv = n_hv[1]
+
+                    for hid in neuron_hv:
+                        if hid not in self.hyper_voxel_projections:
+                            self.hyper_voxel_projections[hid] = set([neuron_id])
+                        else:
+                            self.hyper_voxel_projections[hid].update([neuron_id])
+
+            else:
+                neuron_hv_list = self.find_hyper_voxel_helper_parallel(proj_list)
+
+                for n_hv in neuron_hv_list:
+                    neuron_id = n_hv[0]
+                    neuron_hv = n_hv[1]
+
+                    for hid in neuron_hv:
+                        if hid not in self.hyper_voxel_projections:
+                            self.hyper_voxel_projections[hid] = set([neuron_id])
+                        else:
+                            self.hyper_voxel_projections[hid].update([neuron_id])
+
+            # Make sure the workers get the updated info
+            self.sync_projection_info()
+
+        except:
+            # TODO: Remove this logging code
+            import traceback
+            t_str = traceback.format_exc()
+            self.write_log(t_str, is_error=True)
+            print(t_str)
+            import pdb
+            pdb.set_trace()
+
+    def find_hyper_voxel_helper_parallel(self, proj_info):
+
+        neuron_hv = []
+
+        for neuron_id, neuron_type, random_seed in proj_info:
+            n_id, hyper_voxels = self.find_hyper_voxel_helper(neuron_id=neuron_id, neuron_type=neuron_type,
+                                                              random_seed=random_seed)
+
+            assert n_id == neuron_id
+            neuron_hv.append((n_id, hyper_voxels))
+
+        return neuron_hv
+
+    def find_hyper_voxel_helper(self, neuron_id, neuron_type, random_seed):
+
+        """ Returns hyper voxels that neuron_id (of neuron_type) projects to. """
+
+        try:
+            rng = np.random.default_rng(random_seed)
+
+            hyper_voxels = set()
+
+            if neuron_type in self.projections:
+                for proj in self.projections[neuron_type]:
+                    target_pos, target_rotation, axon_dist = proj["target_info"][neuron_id]
+                    num_points = proj["num_points"]
+
+                    if isinstance(proj["radius"], (np.ndarray, list)):
+                        rx, ry, rz = proj["radius"]
+                    else:
+                        rx = ry = rz = proj["radius"]
+
+                    # Find better way to calculate intersection between ellipsoid and hyper voxel cubes?
+                    pos = self.ellipsoid_coordinates(target_pos, rx, ry, rz, target_rotation, num_points, rng)
+
+                    hv_idx = ((pos - self.snudda_detect.simulation_origo)
+                              / (self.snudda_detect.hyper_voxel_size * self.snudda_detect.voxel_size)).astype(int)
+
+                    hv_id = map(lambda x, y, z: self.snudda_detect.hyper_voxel_id_lookup[x, y, z],
+                                hv_idx[:, 0], hv_idx[:, 1], hv_idx[:, 2])
+
+                    hyper_voxels.update(hv_id)
+
+        except:
+            import traceback
+            t_str = traceback.format_exc()
+            self.write_log(t_str, is_error=True)
+
+        return neuron_id, hyper_voxels
+
+    @staticmethod
+    def ellipsoid_coordinates(target_pos, rx, ry, rz, rotation, num_points, rng):
+
+        """
+            Return num points centred around target_pos, within an ellipsoid with axis rx, ry, rz rotated by rotation
+
+            Args:
+                target_pos : 3 x 1 numpy array
+                rx, ry, rz : float values, axis of ellipsoid (before rotation)
+                rotation : 3x3 rotation matrix, or None
+                num_points : number of points to place
+                rng : numpy random stream
+
+            Returns:
+                coordinates in meters
+        """
+
+        rnd_values = rng.random((num_points, 3))  # Use one call for all random values
+        theta = 2 * np.pi * rnd_values[:, 0]
+        phi = np.arccos(2 * rnd_values[:, 1] - 1)
+        r_scale = rnd_values[:, 2] ** (1 / 3)
+
+        x_coord = np.multiply(rx * r_scale, np.multiply(np.sin(phi), np.cos(theta)))
+        y_coord = np.multiply(ry * r_scale, np.multiply(np.sin(phi), np.sin(theta)))
+        z_coord = np.multiply(rz * r_scale, np.cos(phi))
+
+        pos_offset = np.vstack((x_coord, y_coord, z_coord))
+
+        if rotation is not None:
+            pos_offset = np.matmul(rotation, pos_offset)
+
+        pos = target_pos + pos_offset.T
+
+        return pos
+
+    def voxelise_projections(self):
+
+        """ Add the projection of each neuron in pre_neuron_list to the axon space in currently active hyper voxel
+         """
+
+        if self.snudda_detect.hyper_voxel_id not in self.hyper_voxel_projections:
+            # No projections in this hyper voxel, skip it
+            return
+
+        rng = self.snudda_detect.hyper_voxel_rng
+
+        pre_neuron_list = self.hyper_voxel_projections[self.snudda_detect.hyper_voxel_id]
+
+        for neuron_id in pre_neuron_list:
+            neuron_type = self.snudda_detect.neurons[neuron_id]["type"]
+
+            for proj in self.projections[neuron_type]:
+                target_pos, target_rotation, axon_dist = proj["target_info"][neuron_id]
+                num_points = proj["num_points"]
+                if isinstance(proj["radius"], (np.ndarray, list)):
+                    rx, ry, rz = proj["radius"]
+                else:
+                    rx = ry = rz = proj["radius"]
+
+                pos = self.ellipsoid_coordinates(target_pos, rx, ry, rz, target_rotation, num_points, rng)
+
+                # Convert to coordinates in the hyper voxel
+                voxel_coords = np.round((pos - self.snudda_detect.hyper_voxel_origo) / self.snudda_detect.voxel_size).astype(int)
+
+                inside_idx = np.where(np.sum(np.logical_and(0 <= voxel_coords,
+                                                            voxel_coords < self.snudda_detect.hyper_voxel_size),
+                                             axis=1) == 3)[0]
+
+                for x, y, z in voxel_coords[inside_idx, :]:
+                    if self.snudda_detect.axon_voxel_ctr[x, y, z] < self.snudda_detect.max_axon:
+                        self.snudda_detect.axon_voxels[x, y, z, self.snudda_detect.axon_voxel_ctr[x, y, z]] = neuron_id
+                        self.snudda_detect.axon_voxel_ctr[x, y, z] += 1
+                        self.snudda_detect.axon_soma_dist[x, y, z] = axon_dist  # This is an underestimation
+                    else:
+                        self.snudda_detect.voxel_overflow_ctr += 1
+
+    def get_neurons_of_type(self, neuron_type):
+        neuron_id = np.array([x["neuronID"] for x in self.snudda_detect.neurons if x["type"] == neuron_type])
+        return neuron_id
+
+    def get_neuron_type(self, neuron_id):
+        return self.snudda_detect.neurons[neuron_id]["type"]
+
+    def get_neuron_positions(self, neuron_id):
+        return self.snudda_detect.neuron_positions[neuron_id, :]
+
+    def delete_projection(self, projection_name, pre_neuron_type):
+
+        projection_list = self.projections[pre_neuron_type]
+        new_projection_list = []
+        for proj in projection_list:
+            if proj["name"] != projection_name:
+                new_projection_list.append(proj)
+
+        self.projections[pre_neuron_type] = new_projection_list
+
+    def parse_config(self):
+
+        for con_name, con_info in self.snudda_detect.config["Connectivity"].items():
+            for con_type, con_config in con_info.items():
+                if "projectionFile" in con_config and "projectionName" not in con_config:
+                    pre_neuron_type = con_name.split(",")[0]
+                    self.add_projection(projection_name=f"{con_name},{con_type}", pre_neuron_type=pre_neuron_type,
+                                        projection_file=con_config["projectionFile"])
+
+    def add_projection(self, projection_name, pre_neuron_type, projection_file):
+
+        """
+
+            Args:
+                projection_name (str) : Name of projection
+                pre_neuron_type (str) : Neuron type of presynaptic neuron
+                projection_file (str) : Path to projection file, used by griddata to find target region coordinates
+
+            Projection mapping and rotations are specified in the projection_file.
+
+        """
+
+        with open(projection_file, "r") as f:
+            projection_data = json.load(f, object_pairs_hook=OrderedDict)
+
+        if projection_name and projection_name in projection_data:
+            proj_file_info = projection_data[projection_name]
+        else:
+            proj_file_info = projection_data
+
+        assert "source" in proj_file_info, f"'source' must exist in {projection_file}"
+        assert "destination" in proj_file_info, f"'destination' must exist in {projection_file}"
+        assert "radius" in proj_file_info, f"'radius' must exist in {projection_file}"
+        assert "numPoints" in proj_file_info, f"'numPoints' must exist in {projection_file}"
+
+        proj_info = dict()
+        proj_info["name"] = projection_name
+        proj_info["source"] = np.array(proj_file_info["source"]) * 1e-6
+        proj_info["dest"] = np.array(proj_file_info["destination"]) * 1e-6
+
+        proj_info["radius"] = proj_file_info["radius"]  # Either single radius, or (rx, ry, rz)
+        proj_info["num_points"] = proj_file_info["numPoints"]  # Number of axon points placed
+
+        if "rotation" in projection_data[projection_name]:
+            proj_info["rotation"] = np.array(proj_file_info["rotation"])  # n x 9 matrix
+
+        # Pre-compute target centres for each neuron of pre_neuron_type
+        neuron_id = self.get_neurons_of_type(pre_neuron_type)
+        pre_positions = self.get_neuron_positions(neuron_id)
+        target_centres = griddata(points=proj_info["source"], values=proj_info["dest"],
+                                  xi=pre_positions, method="linear")
+
+        axon_dist = np.linalg.norm(target_centres - pre_positions, axis=1)
+
+        if "rotation" in proj_info:
+            target_rotation = griddata(points=proj_info["dest"], values=proj_info["rotation"],
+                                       xi=target_centres, method="linear")
+        else:
+            target_rotation = [None for x in neuron_id]
+
+        proj_info["target_info"] = dict()
+        for nid, pos, rot, ad in zip(neuron_id, target_centres, target_rotation, axon_dist):
+            if rot is not None:
+                rot_mat = rot.reshape((3, 3))
+                assert np.abs(np.linalg.det(rot_mat) - 1) < 1e-6, \
+                    f"Invalid rotation matrix for {projection_name}: {rot_mat}"
+                proj_info["target_info"][nid] = pos, rot_mat, ad
+            else:
+                proj_info["target_info"][nid] = pos, None, ad
+
+        if pre_neuron_type not in self.projections:
+            self.projections[pre_neuron_type] = [proj_info]
+        else:
+            self.projections[pre_neuron_type].append(proj_info)
+
+
diff --git a/snudda/init/init.py b/snudda/init/init.py
index 0517f9cbd..e471b3943 100644
--- a/snudda/init/init.py
+++ b/snudda/init/init.py
@@ -231,7 +231,11 @@ def add_neuron_target(self, neuron_name, target_name, connection_type,
                           cluster_synapses=False,
                           mod_file=None,
                           parameter_file=None,
-                          channel_param_dictionary=None):
+                          channel_param_dictionary=None,
+                          projection_file=None):
+
+        # OBS, projection file is only needed if you want to create a projection between structures.
+        # For normal touch detection within a volume it is not needed.
 
         if conductance is None:
             conductance = [1.0e-9, 0]
@@ -312,6 +316,9 @@ def add_neuron_target(self, neuron_name, target_name, connection_type,
             # Different pruning rules for within and between neuron units
             con_info["pruningOther"] = pruning_info_other
 
+        if projection_file:
+            con_info["projectionFile"] = projection_file
+
         # Json did not like tuples in keys, so we separate by comma
         nt_key = f"{neuron_name},{target_name}"
         if nt_key not in self.network_data["Connectivity"]:
diff --git a/snudda/input/input.py b/snudda/input/input.py
index 5e4720f7f..4db5b23c2 100644
--- a/snudda/input/input.py
+++ b/snudda/input/input.py
@@ -12,6 +12,9 @@
 
 # Smith, Galvan, ..., Bolam 2014 -- Bra info om thalamic inputs, CM/PF
 #
+
+# TODO: Randomise conductance for the inputs and store it, use it later when adding external synapses in simulate.py
+
 import sys
 from collections import OrderedDict
 
@@ -430,6 +433,7 @@ def make_neuron_input_parallel(self):
         parameter_file_list = []
         parameter_list_list = []
         cluster_size_list = []
+        cluster_spread_list = []
 
         dendrite_location_override_list = []
 
@@ -608,7 +612,13 @@ def make_neuron_input_parallel(self):
                     else:
                         cluster_size = None
 
+                    if "clusterSpread" in input_inf:
+                        cluster_spread = input_inf["clusterSpread"]
+                    else:
+                        cluster_spread = 20e-6
+
                     cluster_size_list.append(cluster_size)
+                    cluster_spread_list.append(cluster_spread)
 
                     if "dendriteLocation" in input_info:
                         assert "morphologyKey" in input_info, \
@@ -673,6 +683,7 @@ def make_neuron_input_parallel(self):
                                   parameter_list_list,
                                   seed_list,
                                   cluster_size_list,
+                                  cluster_spread_list,
                                   dendrite_location_override_list))
 
             self.d_view.scatter("input_list", input_list, block=True)
@@ -713,6 +724,7 @@ def make_neuron_input_parallel(self):
                       parameter_list_list,
                       seed_list,
                       cluster_size_list,
+                      cluster_spread_list,
                       dendrite_location_override_list)
 
         # Gather the spikes that were generated in parallel
@@ -1144,6 +1156,7 @@ def dendrite_input_locations(self,
                                  synapse_density=None,
                                  num_spike_trains=None,
                                  cluster_size=None,
+                                 cluster_spread=30e-6,
                                  input_type=None):
 
         """
@@ -1155,6 +1168,7 @@ def dendrite_input_locations(self,
             synapse_density (str): Distance function f(d)
             num_spike_trains (int): Number of spike trains
             cluster_size (int): Size of each synaptic cluster (None = No clustering)
+            cluster_spread (float): Spread of cluster along dendrite (in meters)
             input_type (str): Type of input, eg. "Cortical" or "Thalamic"
         """
 
@@ -1220,7 +1234,9 @@ def dendrite_input_locations(self,
 
         return morphology.dendrite_input_locations(synapse_density=synapse_density,
                                                    num_locations=num_spike_trains,
-                                                   rng=rng, cluster_size=cluster_size)
+                                                   rng=rng,
+                                                   cluster_size=cluster_size,
+                                                   cluster_spread=cluster_spread)
 
     ############################################################################
 
@@ -1342,7 +1358,8 @@ def make_input_helper_parallel(self, args):
 
             neuron_id, input_type, freq, start, end, synapse_density, num_spike_trains, p_keep, \
                 population_unit_spikes, jitter_dt, population_unit_id, conductance, correlation, mod_file, \
-                parameter_file, parameter_list, random_seed, cluster_size, dendrite_location_override = args
+                parameter_file, parameter_list, random_seed, cluster_size, cluster_spread, \
+                dendrite_location_override = args
 
             return self.make_input_helper_serial(neuron_id=neuron_id,
                                                  input_type=input_type,
@@ -1362,6 +1379,7 @@ def make_input_helper_parallel(self, args):
                                                  parameter_list=parameter_list,
                                                  random_seed=random_seed,
                                                  cluster_size=cluster_size,
+                                                 cluster_spread=cluster_spread,
                                                  dendrite_location=dendrite_location_override)
 
         except:
@@ -1398,6 +1416,7 @@ def make_input_helper_serial(self,
                                  parameter_list,
                                  random_seed,
                                  cluster_size=None,
+                                 cluster_spread=None,
                                  dendrite_location=None):
 
         """
@@ -1421,6 +1440,7 @@ def make_input_helper_serial(self,
             parameter_list: Parameter list (to inline parameters, instead of reading from file)
             random_seed: Random seed.
             cluster_size: Input synapse cluster size
+            cluster_spread: Spread of cluster along dendrite (in meters)
             dendrite_location: Override location of dendrites, list of (sec_id, sec_x) tuples.
             """
 
@@ -1467,6 +1487,7 @@ def make_input_helper_serial(self,
                                                           num_spike_trains=num_spike_trains,
                                                           rng=rng,
                                                           cluster_size=cluster_size,
+                                                          cluster_spread=cluster_spread,
                                                           input_type=input_type)
 
             num_inputs = input_loc[0].shape[0]
diff --git a/snudda/input/input_tuning.py b/snudda/input/input_tuning.py
index 2ae832a93..119e5365c 100644
--- a/snudda/input/input_tuning.py
+++ b/snudda/input/input_tuning.py
@@ -905,7 +905,7 @@ def simulate(self, mech_dir=None):
         sim.add_external_input()
         sim.check_memory_status()
 
-        sim.add_recording()
+        sim.add_volt_recording()
 
         t_sim = self.max_time * 1000  # Convert from s to ms for Neuron simulator
 
diff --git a/snudda/neurons/neuron_morphology.py b/snudda/neurons/neuron_morphology.py
index f73dc4834..54a545f78 100644
--- a/snudda/neurons/neuron_morphology.py
+++ b/snudda/neurons/neuron_morphology.py
@@ -11,8 +11,10 @@
 import numexpr
 from snudda.utils.snudda_path import snudda_parse_path
 
+# Use np.random.choice for selecting
 
-class NeuronMorphology(object):
+
+class NeuronMorphology:
 
     """ Neuron morphology object. Also see NeuronPrototype class which handles multiple morphology variations. """
     # axonStumpIDFlag should be True if running Network_simulate.py
@@ -44,12 +46,12 @@ def __init__(self,
                  virtual_neuron=False,
                  axon_stump_id_flag=False):
 
-        self.cache_version = 0.9
+        self.cache_version = 0.99
 
-        self.position = np.copy(np.array(position))
+        self.position = np.array(position)
 
         if rotation is not None:
-            self.rotation = np.copy(np.array(rotation))
+            self.rotation = np.array(rotation)
         else:
             self.rotation = None
 
@@ -91,7 +93,7 @@ def __init__(self,
         self.verbose = verbose
         self.use_cache = use_cache
         self.pickle_version = pickle_version
-        self.logFile = logfile
+        self.log_file = logfile
         self.virtual_neuron = virtual_neuron
 
         self.rotated_flag = False
@@ -119,6 +121,13 @@ def __init__(self,
         self.dend_sec_id = np.array((1,))
         self.dend_sec_x = np.array((2, 0))
 
+        # Lookup for section length for each section, based on sec_id
+        self.sec_id_to_len = None
+
+        # Lookup for index into self.dend_links based on sec_id
+        self.sec_id_links = dict()
+        self.sec_id_links_x = dict()
+
         if colour is None:
             self.colour = np.random.random((3,))
         else:
@@ -247,6 +256,9 @@ def clone(self,
             new_neuron.dend_links = self.dend_links
             new_neuron.dend_sec_x = self.dend_sec_x
             new_neuron.dend_sec_id = self.dend_sec_id
+            new_neuron.sec_id_links = self.sec_id_links
+            new_neuron.sec_id_links_x = self.sec_id_links_x
+            new_neuron.sec_id_to_len = self.sec_id_to_len
 
             new_neuron.axon_stump_id_flag = self.axon_stump_id_flag
 
@@ -274,8 +286,8 @@ def clone(self,
     def write_log(self, text, is_error=False):
 
         """ Write text to log file. Prints on screen if self.verbose or is_error """
-        if self.logFile is not None:
-            self.logFile.write(f"{text}\n")
+        if self.log_file is not None:
+            self.log_file.write(f"{text}\n")
 
         if self.verbose or is_error:
             print(text)
@@ -421,6 +433,11 @@ def save_cache(self, cache_file=None):
         morph["maxDendRadius"] = self.max_dend_radius
         morph["dendDensity"] = self.dend_density
         morph["axonDensity"] = self.axon_density
+
+        morph["secIdToSecLen"] = self.sec_id_to_len
+        morph["secIdLinks"] = self.sec_id_links
+        morph["secIdLinksX"] = self.sec_id_links_x
+
         morph["version"] = self.cache_version
 
         assert (cache_file != self.swc_filename)
@@ -511,6 +528,11 @@ def load_cache(self, cache_file=None):
         else:
             self.axon_density = None
 
+        self.sec_id_to_len = morph["secIdToSecLen"]
+
+        self.sec_id_links = morph["secIdLinks"]
+        self.sec_id_links_x = morph["secIdLinksX"]
+
         # Place neuron -- Do not place neuron, loadNeuronMorphology does that
         # self.place()
 
@@ -586,7 +608,7 @@ def load_swc(self, swc_file=None):
                 continue
 
             links[link_idx, 0:2] = [id0, id1]
-            links[link_idx, 5] = points[idx, 5]     # !!! TODO: Should take from points[idx, 6]
+            links[link_idx, 5] = points[idx, 6]
 
             link_idx += 1
 
@@ -607,7 +629,7 @@ def load_swc(self, swc_file=None):
                 import pdb
                 pdb.set_trace()
 
-        # Make sure soma has a child count > 1 --- no we dont want some as node
+        # Make sure soma has a child count > 1 --- no we dont want soma as node
         # if(points[0,9] == 0):
         #   points[0,9] = 100
 
@@ -787,6 +809,27 @@ def load_swc(self, swc_file=None):
             self.dend_sec_id = np.zeros((0,))
             self.dend_sec_x = np.zeros((0, 2))
 
+        # Lookup for section length based on section id
+        self.sec_id_to_len = np.ones((int(1+np.max(points[:, 10])),))
+        for sec_id, sec_len in zip(points[:, 10].astype(int), points[:, 11]):
+            if sec_id >= 0:
+                self.sec_id_to_len[sec_id] = sec_len
+
+        # Lookup table to find self.dend_links with a given section id
+        tmp_id_lookup = dict()
+        tmp_x_lookup = dict()
+        for idx, (sec_id, sec_x) in enumerate(zip(self.dend_sec_id, self.dend_sec_x)):
+            if sec_id not in tmp_id_lookup:
+                tmp_id_lookup[sec_id] = [idx]
+                tmp_x_lookup[sec_id] = [sec_x[1]]
+            else:
+                tmp_id_lookup[sec_id].append(idx)
+                tmp_x_lookup[sec_id].append(sec_x[1])
+
+        for k in tmp_id_lookup.keys():
+            self.sec_id_links[k] = np.array(tmp_id_lookup[k])
+            self.sec_id_links_x[k] = np.array(tmp_x_lookup[k])
+
         # self.dendriteDensity() # -- depricated
         self.find_radius()
         self.place()
@@ -1048,7 +1091,7 @@ def compartment_length(self, comp_type="dend"):
     # TODO: Update the code so that it gives exactly num_locations positions (currently it varies)
 
     def dendrite_input_locations(self, synapse_density, rng, num_locations=None,
-                                 cluster_size=None):
+                                 cluster_size=None, cluster_spread=20e-6):
 
         """
         Randomises input locations on dendrites.
@@ -1065,7 +1108,8 @@ def dendrite_input_locations(self, synapse_density, rng, num_locations=None,
             return self.dendrite_input_locations_helper(synapse_density=synapse_density,
                                                         rng=rng,
                                                         num_locations=num_locations,
-                                                        cluster_size=cluster_size)
+                                                        cluster_size=cluster_size,
+                                                        cluster_spread=cluster_spread)
 
         expected_synapses = self.get_expected_synapses_per_compartment(synapse_density=synapse_density)
 
@@ -1155,11 +1199,60 @@ def get_expected_synapses_per_compartment(self, synapse_density):
 
     ############################################################################
 
+    def cluster_synapses(self, sec_id, sec_x, count, distance, rng):
+
+        """
+            Randomize sec_x for cluster of 'count' synapses with centre placed at 'sec_id', 'sec_x'
+            spread over 'distance' (but constrained to current section extent).
+
+            Args:
+                sec_id : Section id of cluster centre
+                sec_x : Section x of cluster centre
+                count : Number of synapses in cluster
+                distance : Maximal spread of cluster along dendrite
+                rng : Numpy random stream
+
+            Returns:
+                cluster_sec_x : Section x for cluster synapse
+                coords : Coordinates for synapse (in meters)
+                soma_dist : Distance to soma (in meters)
+        """
+
+        sec_len = self.sec_id_to_len[sec_id]
+        min_sec_x = np.maximum(1e-3, sec_x - 0.5 * distance / sec_len)
+        max_sec_x = np.minimum(1 - 1e-3, sec_x + 0.5 * distance / sec_len)
+
+        cluster_sec_x = rng.uniform(low=min_sec_x, high=max_sec_x, size=count)
+
+        syn_coords = np.zeros((count, 3))
+        soma_dist = np.zeros((count,))
+
+        for i, sx in enumerate(cluster_sec_x):
+            # We also need to calculate the x,y,z coordinates and distance to soma
+            idx = len(self.sec_id_links_x[sec_id]) - np.sum(sx < self.sec_id_links_x[sec_id])
+            link_idx = self.sec_id_links[sec_id][idx]
+            assert self.dend_sec_id[link_idx] == sec_id
+
+            x_start, x_end = self.dend_sec_x[link_idx]
+            comp_x = (x_start - sx) / (x_start - x_end)
+            start_info = self.dend[self.dend_links[link_idx, 0], :]  # x,y,z,radie,soma dist
+            end_info = self.dend[self.dend_links[link_idx, 1], :]
+            syn_info = start_info * (1 - comp_x) + end_info * comp_x
+
+            syn_coords[i, :] = syn_info[:3]
+            soma_dist[i] = syn_info[4]
+
+        # OBS, syn_coords in meters, and soma dist in meters also
+        return cluster_sec_x, syn_coords, soma_dist
+
+    ############################################################################
+
     def dendrite_input_locations_helper(self,
                                         synapse_density,
                                         rng,
                                         num_locations,
-                                        cluster_size=1):
+                                        cluster_size=1,
+                                        cluster_spread=20e-6):
 
         """
         Helper function. Places num_location clusters, each containing cluster_size synapses. Density is relative, scaled.
@@ -1180,21 +1273,28 @@ def dendrite_input_locations_helper(self,
             cluster_size = 1
 
         num_synapses = num_locations*cluster_size
-        input_loc = np.zeros((num_synapses, 5))
-        dist_syn_soma = np.zeros((num_synapses,))
+        input_loc = np.zeros((num_synapses, 6))
 
         d = self.dend[:, 4]
         syn_ctr = 0
 
-        for r in rand_vals:
-            comp_idx = len(p_cum) - np.sum(r < p_cum)
+        assert np.all(expected_synapses >= 0), f"Check your synapse density {synapse_density}, found negative values."
+        exp_sum = np.sum(expected_synapses)
+        assert exp_sum > 0, f"Check your synapse density {synapse_density}, all elements zero."
 
-            for j in range(cluster_size):
+        syn_idx = rng.choice(len(expected_synapses), size=num_locations, replace=True,
+                             p=expected_synapses/exp_sum)
 
-                comp_x = rng.random()
+        for comp_idx in syn_idx:
+
+            # Cluster centre sec_id:
+            sec_id = self.dend_sec_id[comp_idx]
+            # We place cluster centre in the middle of the 5 micrometer morphology piece
+            sec_x = (self.dend_sec_x[comp_idx, 0] + self.dend_sec_x[comp_idx, 1]) / 2
 
-                dist_syn_soma[syn_ctr] = (d[self.dend_links[comp_idx, 0]] * (1 - comp_x)
-                                          + d[self.dend_links[comp_idx, 1]] * comp_x)
+            if cluster_size is None or cluster_size == 1:
+
+                comp_x = rng.random()
 
                 coords = (self.dend[self.dend_links[comp_idx, 0], :3] * (1 - comp_x)
                           + self.dend[self.dend_links[comp_idx, 1], :3] * comp_x)
@@ -1203,13 +1303,36 @@ def dendrite_input_locations_helper(self,
                 input_loc[syn_ctr, 3] = self.dend_sec_id[comp_idx]
                 input_loc[syn_ctr, 4] = (self.dend_sec_x[comp_idx, 0] * (1 - comp_x)
                                          + comp_x * self.dend_sec_x[comp_idx, 1])
+                input_loc[syn_ctr, 5] = (d[self.dend_links[comp_idx, 0]] * (1 - comp_x)
+                                         + d[self.dend_links[comp_idx, 1]] * comp_x)
 
                 syn_ctr += 1
 
+            else:
+                syn_info = self.cluster_synapses(sec_id=sec_id, sec_x=sec_x, count=cluster_size,
+                                                 distance=cluster_spread, rng=rng)
+                input_loc[syn_ctr:syn_ctr+cluster_size, :3] = syn_info[1]  # Synapse x,y,z -- in global coordinate system(!)
+                input_loc[syn_ctr:syn_ctr + cluster_size, 3] = sec_id      # Sec id
+                input_loc[syn_ctr:syn_ctr + cluster_size, 4] = syn_info[0] # Sec x
+                input_loc[syn_ctr:syn_ctr + cluster_size, 5] = syn_info[2] # Soma distance
+
+                syn_ctr += cluster_size
+
         assert syn_ctr == input_loc.shape[0], f"Not all input_loc was set. Rows {input_loc.shape[0]}, syn_ctr={syn_ctr}"
 
         # Return xyz,secID,secX, dist_to_soma (update: now also added distance synapse to soma)
-        return input_loc[:, :3], input_loc[:, 3], input_loc[:, 4], dist_syn_soma
+        return input_loc[:, :3], input_loc[:, 3], input_loc[:, 4], input_loc[:, 5]
+
+    ############################################################################
+
+    def place_cluster_helper(self, cluster_xyz, cluster_radius, cluster_size):
+
+        raise NotImplementedError("This is not yet implemented")
+
+        # 1. Identify points on neuron within the cluster radius (length)
+        # 2. Pick n synapses from the positions available, restrict choice to segment id of
+        #    parent/child/grandparent/grandchild segments
+        # help: lookup matrix --> 1 x 5
 
     ############################################################################
 
diff --git a/snudda/place/projection_map_finder.py b/snudda/place/projection_map_finder.py
new file mode 100755
index 000000000..f399eaff2
--- /dev/null
+++ b/snudda/place/projection_map_finder.py
@@ -0,0 +1,150 @@
+#!/usr/bin/env python3
+
+import numpy as np
+from scipy.spatial.distance import pdist, squareform
+
+
+class PointCluster:
+
+    def __init__(self):
+        self.coords = None
+        self.dist = None
+        self.num = 0
+
+    def generate_cube(self, centre, side, num):
+
+        self.coords = side * np.random.uniform(size=(num, 3))
+        self.coords[:, 0] += centre[0]
+        self.coords[:, 1] += centre[1]
+        self.coords[:, 2] += centre[2]
+
+        self.num = num
+        self.calc_distance()
+
+    def generate_disk(self, centre, side, num):
+
+        self.coords = side * np.random.uniform(size=(num, 3))
+        self.coords[:, 0] += centre[0]
+        self.coords[:, 1] += centre[1]
+        self.coords[:, 2] = centre[2]
+
+        self.num = num
+        self.calc_distance()
+
+    def calc_distance(self):
+        self.dist = squareform(pdist(self.coords))
+
+    def get_colour(self, idx):
+        min_coord = np.min(self.coords)
+        max_coord = np.max(self.coords)
+
+        return np.divide(self.coords[idx, :] - min_coord, max_coord-min_coord)
+
+
+class PointMap:
+
+    def __init__(self, source, dest):
+
+        self.source = source
+        self.dest = dest
+
+        self.map = np.arange(source.num)
+
+        assert source.num == dest.num, f"source and dest should have equal number of points"
+
+        self.start_energy = None
+
+    def delta_energy_WRONG(self, a, b):
+
+        # For each pair of neurons in the source structure, find their distance, and the distance of their projections.
+        # The sum of the product of these distances is the total energy.
+        # So the question is how does this sum change if we would swap the neurons.
+
+        a_dest = self.map[a]
+        b_dest = self.map[b]
+
+        test_map = self.map.copy()
+        test_map[a] = self.map[b]
+        test_map[b] = self.map[a]
+
+        energy_before = np.sum(np.multiply(self.source.dist[a, :], self.dest.dist[a_dest, :][test_map])) \
+            + np.sum(np.multiply(self.source.dist[b, :], self.dest.dist[b_dest, :][test_map]))
+
+        energy_after = np.sum(np.multiply(self.source.dist[a, :], self.dest.dist[b_dest, :][test_map])) \
+            + np.sum(np.multiply(self.source.dist[b, :], self.dest.dist[a_dest, :][test_map]))
+
+        return energy_after - energy_before
+
+    def delta_energy(self, a, b):
+
+        energy_before = self.energy()
+        alt_map = self.map.copy()
+        alt_map[a] = self.map[b]
+        alt_map[b] = self.map[a]
+        energy_after = self.energy(alt_map)
+
+        return energy_after - energy_before
+
+    def energy(self, alt_map=None):
+
+        if alt_map is not None:
+            return np.sum(np.power(self.source.dist - self.dest.dist[:, alt_map][alt_map, :], 2))
+
+        return np.sum(np.power(self.source.dist - self.dest.dist[:, self.map][self.map, :], 2))
+
+    def iterate(self, n_iter):
+
+        if self.start_energy is None:
+            self.start_energy = self.energy()
+
+        # Pick a pair of points in source, check if their dest points should be swapped
+        ab = np.random.randint(low=0, high=self.source.coords.shape[0], size=(n_iter, 2))
+        p_val = np.random.uniform(size=n_iter)
+
+        for idx, ((a, b), px) in enumerate(zip(ab, p_val)):
+            energy_change = self.delta_energy(a, b)
+            p_swap = 1 / (1 + np.exp(4 * energy_change / self.start_energy * 100))
+            print(f"{idx}/{n_iter}  energy = {self.energy()}, p_swap = {p_swap}")
+            if px < p_swap:
+
+                # Swap the two if energy
+                old_a = self.map[a]
+                self.map[a] = self.map[b]
+                self.map[b] = old_a
+
+            if idx % 100000 == 0:
+                self.plot_map()
+
+        print(f"Energy = {self.energy()}")
+        self.plot_map()
+
+    def plot_map(self):
+
+        import matplotlib.pyplot as plt
+        fig = plt.figure()
+        ax = fig.add_subplot(111, projection='3d')
+
+        for idx in range(self.source.num):
+            col = self.source.get_colour(idx)
+            sp = self.source.coords[idx, :]
+            dp = self.dest.coords[self.map[idx], :]
+            ax.plot(xs=sp[0], ys=sp[1], zs=sp[2], color=col, marker='o')
+            ax.plot(xs=dp[0], ys=dp[1], zs=dp[2], color=col, marker='o')
+
+        fig.show()
+        plt.pause(0.1)
+
+
+if __name__ == "__main__":
+    pc_a = PointCluster()
+    pc_b = PointCluster()
+
+    pc_a.generate_disk(centre=np.array([0, 0, 0]), side=1000e-6, num=100)
+    pc_b.generate_disk(centre=np.array([0, 0, 5e-3]), side=1000e-6, num=100)
+
+    pm = PointMap(source=pc_a, dest=pc_b)
+
+    pm.iterate(10000)
+
+    import pdb
+    pdb.set_trace()
\ No newline at end of file
diff --git a/snudda/place/region_mesh.py b/snudda/place/region_mesh.py
index 85e962c49..ed3f2d3f2 100644
--- a/snudda/place/region_mesh.py
+++ b/snudda/place/region_mesh.py
@@ -209,7 +209,7 @@ def mark_borders(self):
                     ic = np.floor((x_point - self.min_coord) / self.bin_width)
                     self.voxel_mask_border[int(ic[0]), int(ic[1]), int(ic[2])] = 1
 
-        self.write_log(f"maxN = {max_num}")
+        self.write_log(f"max_num = {max_num}")
 
     def __del__(self):
         if self.close_log_file:
diff --git a/snudda/place/volumetric_mapping.py b/snudda/place/volumetric_mapping.py
new file mode 100644
index 000000000..1925a6b00
--- /dev/null
+++ b/snudda/place/volumetric_mapping.py
@@ -0,0 +1,181 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Jan 28 11:56:04 2022
+
+@author: wst
+"""
+
+import os
+import numpy as np
+import json
+from scipy.spatial.distance import pdist, squareform
+from snudda.place.region_mesh import RegionMesh
+
+class VolumetricMap:
+    
+    def __init__(self, meshfile_a, meshfile_b, n = 100, inversion = False, initial = True, predefined_points = None):
+        self.meshfile_a = meshfile_a
+        self.meshfile_b = meshfile_b
+        self.n = n
+        self.inversion = inversion
+        self.initial = initial
+        
+        if predefined_points:                                    #user defined nested list of points to map between                    
+            self.a_points = np.array(predefined_points[0])
+            self.b_points = np.array(predefined_points[1])
+        else:                                                   #highjack RegionMesh to sample the volume
+            self.a = RegionMesh(meshfile_a)
+            self.b = RegionMesh(meshfile_b)
+            self.a_points = self.a.place_neurons(self.n)
+            self.b_points = self.b.place_neurons(self.n)
+            
+    def distance(self, coords):
+        
+        return squareform(pdist(coords))
+        
+    def initial_map(self):
+        
+        bcopy = self.b_points.copy()
+        if self.inversion:                                              #invert ventrodorsal 
+            bcopy[:,2]*= -1
+            
+        mapping = []
+        for i in range(len(self.a_points)):
+            j= np.linalg.norm(self.a_points[i] - bcopy, axis = 1).argmin()
+            k = np.argwhere(self.b_points == bcopy[j])[0,0]
+            mapping.append(k)
+            bcopy = np.delete(bcopy, j, axis = 0)
+            
+        return mapping
+    
+    def energy(self, mapping):
+        
+        a_dist = self.distance(self.a_points)
+        if self.inversion:  #invert ventrodorsal topography
+            inverted_points = self.b_points.copy() 
+            inverted_points[:,2]*= -1
+            # inverted_points[:,1]*= 1
+            b_dist = self.distance(inverted_points)
+        else:
+            b_dist = self.distance(self.b_points)
+
+        return np.sum(np.power(a_dist - b_dist[:, mapping][mapping, :], 2))
+
+    def optimize_map(self, n_iter):
+        
+        if self.initial:
+            self.mapping = self.initial_map()
+        else:
+            self.mapping = np.arange(0,len(self.a_points) ,1)
+            
+        history = []                                                        #track optimization
+        e1 = self.energy(self.mapping)
+        e = e1
+        history.append(e)
+        
+        m = np.random.randint(low=0, high=self.n, size=(n_iter, 2)) 
+        p_val = np.random.uniform(size=n_iter)
+
+        for i, ((m1,m2), px) in enumerate(zip(m, p_val)):
+            
+            new = self.mapping.copy()
+            new[m1] = self.mapping[m2]
+            new[m2] = self.mapping[m1]
+            e_new = self.energy(new)
+            e_change = e_new - e
+            
+            p_swap = 1 / (1 + np.exp(50 * e_change / e1 * 100))
+            
+            if px < p_swap:
+                self.mapping = new
+                e = e_new
+            history.append(e)
+            
+        self.opt_history = np.array(history)
+        self.a_points = self.a_points
+        self.b_points = self.b_points[self.mapping]
+        
+        return 
+    
+    def plot_energy(self):
+        
+        import matplotlib.pyplot as plt
+        if self.opt_history is None:
+            self.optimize_map(1000)
+            
+        fig, ax = plt.subplots()
+        ax.plot(self.opt_history)
+        ax.set_xlabel('Iterations')
+        ax.set_ylabel('Energy')
+        title = '{} points, {} iterations'.format(self.n, len(self.opt_history) -1)
+        ax.set_title(title)
+        
+        return
+    
+    def get_colour(self, coords, idx):
+        min_coord = np.min(coords)
+        max_coord = np.max(coords)
+
+        return np.divide(coords[idx, :] - min_coord, max_coord-min_coord)
+
+    def visualize_map(self, lines = True, points = True):
+            
+        import trimesh
+        from vedo import embedWindow
+        embedWindow(backend = False) 
+        from vedo import Mesh
+        from vedo import Line
+        import brainrender
+        brainrender.settings.SHADER_STYLE = 'plastic'
+        brainrender.settings.SHOW_AXES = False
+        from brainrender.scene import Scene
+        
+        a_tri = trimesh.load_mesh(self.meshfile_a)
+        b_tri = trimesh.load_mesh(self.meshfile_b)
+        
+        scene = Scene(root = True)
+        scene.jupyter = True
+        
+        scene.add(Mesh(a_tri, c = 'red', alpha = 0.1))
+        scene.add(Mesh(b_tri, c = 'blue', alpha = 0.1))
+
+        if lines:
+            for i in range(self.n):
+                col = self.get_colour(self.a_points[:], i)
+                scene.add(Line(self.a_points[i,:]*1e6,self.b_points[i,:]*1e6, lw= 2.5, c = col))
+        if points:
+            for i in range(self.n):
+                col = self.get_colour(self.a_points[:], i)
+                scene.add(Point(self.a_points[i,:]*1e6, radius= 25, color = col, alpha = 0.4))
+                scene.add(Point(self.b_points[i,:]*1e6, radius= 25, color = col, alpha = 0.4))
+
+        scene.slice(scene.atlas.get_plane(norm = [0, 0, 1], sx=15000, sy=15000), close_actors=False)
+        scene.render(camera = 'top')
+        
+        return
+    
+    def write_json(self, projection_name, outfile):
+        
+        mapping ={projection_name: {'source': self.a_points.tolist(), 'destination': self.b_points.tolist()}}
+        
+        with open(outfile, 'w+') as f:
+            json.dump(mapping, f, indent = 4)
+            
+        return
+
+
+if __name__ == "__main__":
+
+    source_file = os.path.join("$SNUDDA_DATA", "mesh", "striatum-d.obj")
+    target_file = os.path.join("$SNUDDA_DATA", "mesh", "SNr.obj")
+    outfile = os.path.join("$SNUDDA_DATA", "projections", "direct.json")
+    projection_name = 'DirectPathway'
+    
+    vm = VolumetricMap(source_file, target_file, n = 50, inversion = True)
+    vm.optimize_map(10000)
+    vm.plot_energy()
+    vm.visualize_map()
+    vm.write_json(projection_name = projection_name, outfile = outfile )
+
+
diff --git a/snudda/simulate/model_current_injections.py b/snudda/simulate/model_current_injections.py
index 41e161124..b8e1f4f0b 100644
--- a/snudda/simulate/model_current_injections.py
+++ b/snudda/simulate/model_current_injections.py
@@ -305,7 +305,7 @@ def _simulate_network_helper(self, sim_name, sim_type, stim_id, measure_id):
         # Also add voltage recording for debugging reasons
         save_voltage = True  # False #True
         if save_voltage:
-            self.snudda_sim.add_recording(cell_id=stim_id)
+            self.snudda_sim.add_volt_recording(cell_id=stim_id)
 
         self.set_gaba_rev(self.gaba_rev)
 
diff --git a/snudda/simulate/network_pair_pulse_simulation.py b/snudda/simulate/network_pair_pulse_simulation.py
index 6a52b22b3..c2457307f 100644
--- a/snudda/simulate/network_pair_pulse_simulation.py
+++ b/snudda/simulate/network_pair_pulse_simulation.py
@@ -280,7 +280,7 @@ def run_sim(self, gaba_rev, pre_id=None):
         # record voltage from all neurons
 
         if self.post_type == "ALL":
-            self.snudda_sim.add_recording()
+            self.snudda_sim.add_volt_recording()
         else:
             # Record from all the potential post synaptic neurons
             self.snudda_sim.add_recording_of_type(self.post_type)
diff --git a/snudda/simulate/nrn_simulator_parallel.py b/snudda/simulate/nrn_simulator_parallel.py
index da4a20b00..05af6295b 100644
--- a/snudda/simulate/nrn_simulator_parallel.py
+++ b/snudda/simulate/nrn_simulator_parallel.py
@@ -68,7 +68,7 @@ def run(
         else:
             self.neuron.h.dt = dt
             self.neuron.h.steps_per_ms = 1.0 / dt
-            print(f'Running Neuron simulator {tstop:.6g} ms, with dt={dt}%r')
+            print(f'Running Neuron simulator {tstop:.6g} ms, with dt={dt}')
 
         if random123_globalindex is None:
             random123_globalindex = self.random123_globalindex
diff --git a/snudda/simulate/pair_recording.py b/snudda/simulate/pair_recording.py
index 485a32767..584ec05a0 100644
--- a/snudda/simulate/pair_recording.py
+++ b/snudda/simulate/pair_recording.py
@@ -22,7 +22,7 @@
 import numpy as np
 import json
 from snudda.simulate import SnuddaSimulate
-from snudda.utils.save_network_activity import SnuddaSaveNetworkActivity
+from snudda.utils.save_network_recording import SnuddaSaveNetworkRecordings
 
 import neuron
 neuron.h.load_file("stdrun.hoc")
@@ -124,7 +124,7 @@ def parse_experiment_config(self):
                                         end_times=stim_end_time, amplitudes=stim_amplitude)
 
         # Add voltage recordings to neurons
-        self.add_recording()
+        self.add_volt_recording()
 
     @staticmethod
     def to_list(val, new_list_len=1):
@@ -279,14 +279,14 @@ def run(self):
 
         # Write results to disk
         try:
-            save = SnuddaSaveNetworkActivity(output_file=self.output_file_name)
+            save = SnuddaSaveNetworkRecordings(output_file=self.output_file_name)
             save.write(t_save=self.t_save, v_save=self.v_save, v_key=self.v_key,
                        t_spikes=self.t_spikes, id_spikes=self.id_spikes)
 
             pre_id = np.array([x[0] for x in self.synapse_currents])
             post_id = np.array([x[1] for x in self.synapse_currents])
             cur = [np.array(x[2]) for x in self.synapse_currents]
-            save.write_currents(t_save=self.t_save, i_save=cur, pre_id=pre_id, post_id=post_id)
+            save.write_currents_OLD(t_save=self.t_save, i_save=cur, pre_id=pre_id, post_id=post_id)
 
         except:
             import traceback
diff --git a/snudda/simulate/simulate.py b/snudda/simulate/simulate.py
index dfaaf0403..41f63e058 100644
--- a/snudda/simulate/simulate.py
+++ b/snudda/simulate/simulate.py
@@ -26,7 +26,7 @@
 import json
 import timeit
 
-from snudda.utils.save_network_activity import SnuddaSaveNetworkActivity
+from snudda.utils.save_network_recording import SnuddaSaveNetworkRecordings
 from snudda.neurons.neuron_model_extended import NeuronModel
 # from Network_place_neurons import NetworkPlaceNeurons
 import numpy as np
@@ -56,6 +56,7 @@ def __init__(self,
                  network_path=None,
                  network_file=None,
                  input_file=None,
+                 output_file=None,
                  verbose=False,
                  log_file=None,
                  disable_gap_junctions=False,
@@ -68,6 +69,7 @@ def __init__(self,
             network_path (str): Path to network directory
             network_file (str, optional): Path to network file, default network-synapses.hdf5 in network_path
             input_file (str, optional): Path to input file, default input-spikes.hdf5 in network_path
+            output_file (str, optional): Path to output file, default simulation/output.hdf5
             verbose (bool): Extra printouts
             log_file (str): Path to logfile
             disable_gap_junctions (bool): Disable gap junctions, default False
@@ -101,6 +103,11 @@ def __init__(self,
         else:
             self.input_file = input_file
 
+        if output_file:
+            self.output_file = output_file
+        else:
+            self.output_file = os.path.join(self.network_path, "simulation", "output.hdf5")
+
         # Init
         self.snudda_loader = None
         self.network_info = None
@@ -162,15 +169,16 @@ def __init__(self,
 
         self.virtual_neurons = {}
 
-        self.net_con_list = []  # Avoid premature garbage collection
-        self.synapse_list = []
+        self.net_con_list = []  # Avoid premature garbage collection -- todo: THIS WILL BE REMOVED
+        self.synapse_list = []  # todo: THIS WILL BE REMOVED, replaced by synapse_dict
+        self.synapse_dict = dict()
         self.i_stim = []
         self.v_clamp_list = []
         self.gap_junction_list = []
         self.external_stim = dict([])
         self.t_save = []
-        self.v_save = []
-        self.v_key = []
+        # self.v_save = []
+        # self.v_key = []
         self.i_save = []
         self.i_key = []
 
@@ -178,8 +186,7 @@ def __init__(self,
 
         self.gap_junction_next_gid = 0  # Are these gids separate from cell gids?
 
-        self.t_spikes = h.Vector()
-        self.id_spikes = h.Vector()
+        self.check_id_recordings = [] # Prevent segmentation fault due to garbage collection of spike id
 
         # Make sure the output dir exists, so we don't fail at end because we
         # cant write file
@@ -194,6 +201,13 @@ def __init__(self,
         # We need to initialise random streams, see Lytton el at 2016 (p2072)
 
         self.load_network_info(self.network_file)
+        
+        self.record = SnuddaSaveNetworkRecordings(output_file=self.output_file, network_data=self.network_info)
+        self.record.add_unit(data_type="voltage", target_unit="V", conversion_factor=1e-3)
+        self.record.add_unit(data_type="synaptic_current", target_unit="A", conversion_factor=1e-9)
+        self.record.add_unit(data_type="spikes", target_unit="s", conversion_factor=1e-3)
+        self.record.add_unit(data_type="time", target_unit="s", conversion_factor=1e-3)
+        # TODO: Add more units as needed https://www.neuron.yale.edu/neuron/static/docs/units/unitchart.html
 
     def setup(self):
 
@@ -508,7 +522,16 @@ def setup_neurons(self):
                 # self.netConList.append(nc) -- Not needed according to Lytton et al 2016
 
                 # Record all spikes
-                self.pc.spike_record(ID, self.t_spikes, self.id_spikes)
+                t_spikes = h.Vector()
+                id_spikes = h.Vector()
+
+                self.pc.spike_record(ID, t_spikes, id_spikes)
+                self.check_id_recordings.append((ID, id_spikes))
+                self.record.register_compartment_data("spikes", neuron_id=ID, data=t_spikes,
+                                                      sec_id=-1, sec_x=0.5)
+
+
+                # self.pc.spike_record(ID, self.t_spikes, self.id_spikes)
 
     ############################################################################
 
@@ -692,13 +715,15 @@ def connect_neuron_synapses(self, start_row, end_row):
         source_id_list, dest_id, dend_sections, sec_id, sec_x, synapse_type_id, \
             axon_distance, conductance, parameter_id = self.get_synapse_info(start_row=start_row, end_row=end_row)
 
-        for (src_id, section, section_x, s_type_id, axon_dist, cond, p_id) \
-                in zip(source_id_list, dend_sections, sec_x, synapse_type_id,
+        for (src_id, section, section_id, section_x, s_type_id, axon_dist, cond, p_id) \
+                in zip(source_id_list, dend_sections, sec_id, sec_x, synapse_type_id,
                        axon_distance, conductance, parameter_id):
 
             try:
                 self.add_synapse(cell_id_source=src_id,
+                                 cell_id_dest=dest_id,
                                  dend_compartment=section,
+                                 section_id=section_id,
                                  section_dist=section_x,
                                  synapse_type_id=s_type_id,
                                  axon_dist=axon_dist,
@@ -798,13 +823,7 @@ def connect_network_gap_junctions_local(self):
 
         (neuron_id, compartment, seg_x, gj_gid_src, gj_gid_dest, cond) = self.find_local_gap_junctions()
 
-        # import pdb
-        # pdb.set_trace()
-
         try:
-            # WHY??!
-            # ValueError: too many values to unpack (expected 6)
-
             for nid, comp, s_x, gid_src, gid_dest, g \
                     in zip(neuron_id, compartment, seg_x, gj_gid_src, gj_gid_dest, cond):
                 self.add_gap_junction(section=comp,
@@ -827,7 +846,7 @@ def get_synapse(channel_module, dend_compartment, section_dist):
         """ Helper method, returns channel_module(dend_compartment(section_dist)) """
         return channel_module(dend_compartment(section_dist))
 
-    def add_synapse(self, cell_id_source, dend_compartment, section_dist, conductance,
+    def add_synapse(self, cell_id_source, cell_id_dest, dend_compartment, section_id, section_dist, conductance,
                     parameter_id, synapse_type_id, axon_dist=None):
 
         """
@@ -892,52 +911,31 @@ def add_synapse(self, cell_id_source, dend_compartment, section_dist, conductanc
                     import pdb
                     pdb.set_trace()
 
-        # Just create a default expsyn for test, will need to create proper GABA
-        # synapses later
-        # if(synapseType == 'ExpSyn'):
-        #  syn = self.sim.neuron.h.ExpSyn(dendCompartment(sectionDist))
-        # elif(synapseType == 'GABA'):
-        #  syn = self.sim.neuron.h.tmGabaA(dendCompartment(sectionDist))
-        # elif(synapseType == "AMPA_NMDA"):
-        #  syn = self.sim.neuron.h.tmGlut(dendCompartment(sectionDist))
-        # else:
-        #  self.writeLog("Synapse type not implemented: ", synapseType)
-        #  import pdb
-        #  pdb.set_trace()
-
         if axon_dist is not None:
             # axon dist is in micrometer, want delay in ms
             synapse_delay = (1e3 * 1e-6 * axon_dist) / self.axon_speed + self.synapse_delay
         else:
             synapse_delay = self.synapse_delay
 
-        #    self.write_log(f"Synapse delay: {synapse_delay} ms")
-
-        # What do we do if the GID does not exist?
-        # print("GID exists:" + str(self.pc.gid_exists(cellIDsource)))
-
         if self.is_virtual_neuron[cell_id_source]:
             # Source is a virtual neuron, need to read and connect input
             src_name = self.network_info["neurons"][cell_id_source]["name"]
 
-            nc = self.pc.gid_connect(cell_id_source, syn)
-            nc.weight[0] = conductance
-            nc.delay = synapse_delay
-            nc.threshold = self.spike_threshold
-
-            # Prevent garbage collection in python
-            self.net_con_list.append(nc)
-            self.synapse_list.append(syn)
+        # Prevent garbage collection in python
+        if (cell_id_source, cell_id_dest) not in self.synapse_dict:
+            self.synapse_dict[cell_id_source, cell_id_dest] = []
 
-        else:
+        nc = self.pc.gid_connect(cell_id_source, syn)
+        nc.weight[0] = conductance
+        nc.delay = synapse_delay
+        nc.threshold = self.spike_threshold
 
-            nc = self.pc.gid_connect(cell_id_source, syn)
-            nc.weight[0] = conductance
-            nc.delay = synapse_delay
-            nc.threshold = self.spike_threshold
+        # This prevents garbage collection of syn and nc
+        self.synapse_dict[cell_id_source, cell_id_dest].append((syn, nc, synapse_type_id, section_id))
 
-            self.net_con_list.append(nc)
-            self.synapse_list.append(syn)
+        # TODO: Johanna promised to remove these when she is ready for it.
+        self.synapse_list.append(syn)
+        self.net_con_list.append(nc)
 
         return syn
 
@@ -1018,8 +1016,8 @@ def add_external_input(self, input_file=None):
                 self.write_log(f"Warning - No input specified for {name}", is_error=True)
                 continue
 
-            for inputType in self.input_data["input"][str(neuron_id)]:
-                neuron_input = self.input_data["input"][str(neuron_id)][inputType]
+            for input_type in self.input_data["input"][str(neuron_id)]:
+                neuron_input = self.input_data["input"][str(neuron_id)][input_type]
 
                 loc_type = 1 * np.ones((neuron_input["sectionID"].shape[0],))  # Axon-Dend
 
@@ -1033,17 +1031,16 @@ def add_external_input(self, input_file=None):
                 eval_str = f"self.sim.neuron.h.{mod_file}"
                 channel_module = eval(eval_str)
 
-                for inputID, (section, section_x, paramID, nSpikes) \
+                for input_id, (section, section_x, param_id, n_spikes) \
                         in enumerate(zip(sections,
                                          neuron_input["sectionX"],
                                          neuron_input["parameterID"],
                                          neuron_input["nSpikes"])):
                     # We need to find cellID (int) from neuronID (string, eg. MSD1_3)
 
-                    idx = inputID
-                    spikes = neuron_input["spikes"][inputID, :nSpikes] * 1e3  # Neuron uses ms
-                    assert (spikes >= 0).all(), \
-                        "Negative spike times for neuron " + str(neuron_id) + " " + inputType
+                    idx = input_id
+                    spikes = neuron_input["spikes"][input_id, :n_spikes] * 1e3  # Neuron uses ms
+                    assert (spikes >= 0).all(), f"Negative spike times for neuron {neuron_id} {input_type}"
 
                     # Creating NEURON VecStim and vector
                     # https://www.neuron.yale.edu/phpBB/viewtopic.php?t=3125
@@ -1068,21 +1065,17 @@ def add_external_input(self, input_file=None):
                     elif section_x == 1.0:
                         section_x = 0.99
 
-                    # !!! Parameters for the tmGlut should be possible to set in the
-                    # input specification !!!
-                    # syn = self.sim.neuron.h.tmGlut(section(sectionX))
                     syn = self.get_external_input_synapse(channel_module, section, section_x)
                     nc = h.NetCon(vs, syn)
 
                     nc.delay = 0.0
-                    # Should weight be between 0 and 1, or in microsiemens?
-                    nc.weight[0] = neuron_input["conductance"][()] * 1e6  # !! what is unit? microsiemens?
+                    nc.weight[0] = neuron_input["conductance"][()] * 1e6  # Neurons needs microsiemens
                     nc.threshold = 0.1
 
                     # Get the modifications of synapse parameters, specific to
                     # this synapse
                     if param_list is not None and len(param_list) > 0:
-                        syn_params = param_list[paramID % len(param_list)]["synapse"]
+                        syn_params = param_list[param_id % len(param_list)]["synapse"]
 
                         for par in syn_params:
                             if par == "expdata":
@@ -1103,21 +1096,10 @@ def add_external_input(self, input_file=None):
                                      f"but expected >= 0.01 and < 10000")
 
                             setattr(syn, par, par_value)
-                            # eval_str = "syn." + par + "=" + str(syn_params[par])
-                            # self.writeLog("Updating synapse: " + evalStr)
-                            # !!! Can we avoid an eval here, it is soooo SLOW
-                            # exec(eval_str)
-
-                    # !!! Set parameters in synParams
 
                     # Need to save references, otherwise they will be freed
-                    # So sorry, but that is how neuron is
                     self.external_stim[neuron_id].append((v, vs, nc, syn, spikes))
 
-                    # ps = h.PatternStim()
-
-                    # HOW DO WE USE PATTERNSTIM?
-
     ############################################################################
 
     def set_resting_voltage(self, neuron_id, rest_volt=None):
@@ -1219,7 +1201,7 @@ def add_recording_of_type(self, neuron_type, num_neurons=None):
         cell_id = self.snudda_loader.get_neuron_id_of_type(neuron_type=neuron_type,
                                                            num_neurons=num_neurons)
 
-        self.add_recording(cell_id)
+        self.add_volt_recording(cell_id)
 
     ############################################################################
 
@@ -1284,7 +1266,82 @@ def add_voltage_clamp(self, cell_id, voltage, duration, res=1e-9, save_i_flag=Fa
 
     ############################################################################
 
-    def add_recording(self, cell_id=None, side_len=None):
+    def add_volt_recording_all(self, cell_id=None):
+
+        if cell_id is None:
+            cell_id = self.neuron_id
+
+        if isinstance(cell_id, (int, np.integer)):
+            cell_id = [cell_id]
+
+        for cid in cell_id:
+            sec_id = [0]
+            sec_x = [0.5]
+
+            for sid, sec in enumerate(self.neurons[cid].icell.dend):
+                for seg in sec.allseg():
+                    sec_id.append(sid+1)  # NEURON indexes from 0, Snudda has soma as 0, and first dendrite is 1
+                    sec_x.append(seg.x)
+
+            self.add_volt_recording(cell_id=cid, sec_id=sec_id, sec_x=sec_x)
+
+    def add_volt_recording_soma(self, cell_id=None):
+
+        if cell_id is None:
+            cell_id = self.neuron_id
+
+        for cid in cell_id:
+            self.add_volt_recording(cid, [0], [0.5])
+
+    def add_volt_recording(self, cell_id : int, sec_id, sec_x):
+
+        # cell_id cant be a list, but sec_id and sec_x must be iterable
+
+        sections = self.neurons[cell_id].map_id_to_compartment(sec_id)
+
+        v = self.sim.neuron.h.Vector()
+
+        self.pc.threshold(cell_id, self.spike_threshold)  # TODO: Set individual spike thresholds based on parameters
+        for s, sx, sid in zip(sections, sec_x, sec_id):
+            v.record(getattr(s(sx), '_ref_v'))
+
+            # From the Snudda synapse matrix. sec_id 0 is soma, sec_id >= 1 is dendrite, sec_id <= -1 is axon
+            self.record.register_compartment_data(neuron_id=cell_id, data_type="voltage", data=v,
+                                                  sec_id=sid, sec_x=sx)
+
+        if self.record.time is None:
+            t_save = self.sim.neuron.h.Vector()
+            t_save.record(self.sim.neuron.h._ref_t)
+            self.record.register_time(time=t_save)
+
+    def add_synapse_current_recording_all(self, dest_id):
+
+        if isinstance(dest_id, (int, np.integer)):
+            dest_id = [dest_id]
+
+        for d_id in dest_id:
+            for sid, did in self.synapse_dict.keys():
+                if did == d_id:
+                    self.add_synapse_current_recording(sid, did)
+
+    def add_synapse_current_recording(self, source_id, dest_id):
+
+        synapse_info_list = self.synapse_dict[source_id, dest_id]
+
+        for syn, nc, synapse_type_id, sec_id in synapse_info_list:
+
+            data = self.sim.neuron.h.Vector()
+            data.record(syn._ref_i)
+            seg = syn.get_segment()
+
+            self.record.register_synapse_data(neuron_id=dest_id, data_type="synaptic_current", data=data,
+                                              synapse_type=synapse_type_id,
+                                              presynaptic_id=source_id,
+                                              sec_id=sec_id,
+                                              sec_x=seg.x,
+                                              cond=nc.weight[0])
+
+    def add_recording_OLD2(self, cell_id=None, side_len=None):
 
         """
         Adds somatic voltage recording to neuron cell_id.
@@ -1319,6 +1376,7 @@ def add_recording(self, cell_id=None, side_len=None):
             sys.exit(-1)
 
         if len(self.t_save) == 0 or self.t_save is None:
+            # TODO: !!! This needs to be registered in save_network_activity
             self.t_save = self.sim.neuron.h.Vector()
             self.t_save.record(self.sim.neuron.h._ref_t)
 
@@ -1330,6 +1388,7 @@ def add_recording(self, cell_id=None, side_len=None):
 
                 self.pc.threshold(cell_key, self.spike_threshold)    # TODO: Set individual spike thresholds based on parameters
                 v.record(getattr(cell.icell.soma[0](0.5), '_ref_v'))
+                # TODO: !!! This needs to be registered in save_network_activity, and ability to save voltage in compartments and currents
 
                 self.v_save.append(v)
                 self.v_key.append(cell_key)
@@ -1376,6 +1435,8 @@ def run(self, t=1000.0, hold_v=None):
 
     def plot(self):
 
+        assert False, "plot is deprecated now"
+
         """ Save a plot of voltage trace. """
 
         import matplotlib.pyplot as pyplot
@@ -1558,18 +1619,9 @@ def verify_synapse_placement(self, sec_list, sec_x_list, dest_id, voxel_coords):
 
     ############################################################################
 
-    def write_output(self, output_file=None):
-
-        """ Save neuron voltage to HDF5 file """
+    def write_output(self):
 
-        if not output_file:
-            output_file = os.path.join(self.network_path, "simulation", "network-output.hdf5")
-        elif os.path.sep not in output_file:
-            output_file = os.path.join(self.network_path, "simulation", output_file)
-            
-        sv = SnuddaSaveNetworkActivity(output_file=output_file, network_data=self.network_info)
-        sv.write(t_save=self.t_save, v_save=self.v_save, v_key=self.v_key,
-                 t_spikes=self.t_spikes, id_spikes=self.id_spikes)
+        self.record.write()
 
     ############################################################################
 
@@ -1578,42 +1630,6 @@ def write_output(self, output_file=None):
     # cellID,v0,v1,v2,v3, ... (voltage for cell #ID)
     # repeat
 
-    def write_voltage_OLD(self,
-                          output_file=None,
-                          down_sampling=10):
-
-        if not output_file:
-            output_file = os.path.join(self.network_path, "simulation", "network-voltage.txt")
-
-        if not os.path.exists(os.path.dirname(output_file)):
-            os.mkdir(os.path.dirname(output_file))
-
-        """ Writes voltage to output_file, with the option to down sample data to save space. """
-
-        for i in range(int(self.pc.nhost())):
-            self.pc.barrier()
-
-            if i == int(self.pc.id()):
-                if i == 0:
-                    mode = 'w'
-                else:
-                    mode = 'a'
-
-                with open(output_file, mode) as voltage_file:
-                    if mode == 'w':
-                        voltage_file.write('-1')  # Indicate that first column is time
-
-                        for t_idx in range(0, len(self.t_save), down_sampling):
-                            voltage_file.write(',%.4f' % self.t_save[t_idx])
-
-                    for v_id, voltage in zip(self.v_key, self.v_save):
-                        voltage_file.write('\n%d' % v_id)
-
-                        for v_idx in range(0, len(voltage), down_sampling):
-                            voltage_file.write(',%.4f' % voltage[v_idx])
-
-            self.pc.barrier()
-
     ############################################################################
 
     # File format for csv current file:
@@ -1625,6 +1641,8 @@ def write_current(self,
                       output_file="save/traces/network-current",
                       down_sampling=20):
 
+        print("This function will be deprecated in the future")
+
         """ Writes current to output_file, with option to down sample data to save space. """
 
         for i in range(int(self.pc.nhost())):
@@ -1786,9 +1804,6 @@ def convert_to_natural_units(self, param_name, param_value):
 
         return val
 
-
-
-
     ############################################################################
 
     # Use event handling
@@ -1922,7 +1937,7 @@ def check_memory_status(self, threshold=0.1):
     sim.check_memory_status()
 
     if args.record_volt:
-        sim.add_recording(side_len=None)  # Side len let you record from a subset
+        sim.add_volt_recording(side_len=None)  # Side len let you record from a subset
         # sim.addRecordingOfType("dSPN",5) # Side len let you record from a subset
         # sim.addRecordingOfType("dSPN",2)
         # sim.addRecordingOfType("iSPN",2)
diff --git a/snudda/utils/export_connection_matrix.py b/snudda/utils/export_connection_matrix.py
index 359eee774..88ee4210a 100755
--- a/snudda/utils/export_connection_matrix.py
+++ b/snudda/utils/export_connection_matrix.py
@@ -29,6 +29,7 @@ def __init__(self, in_file, out_file, save_sparse=True, save_on_init=True):
         self.neuron_type = [x["type"] for x in data["neurons"]]
         self.neuron_name = [x["name"] for x in data["neurons"]]
         self.neuron_morph = [x["morphology"] for x in data["neurons"]]
+        self.population_unit = data["populationUnit"]
 
         self.pos = data["neuronPositions"]
 
@@ -56,8 +57,9 @@ def save(self):
 
         print("Writing " + self.out_file_meta)
         with open(self.out_file_meta, "w") as f_out_meta:
-            for i, (nt, nn, p, mf) in enumerate(zip(self.neuron_type, self.neuron_name, self.pos, self.neuron_morph)):
-                s = "%d,%s,%s,%f,%f,%f,%s\n" % (i, nt, nn, p[0], p[1], p[2], mf)
+            for i, (nt, nn, p, mf, pu) in enumerate(zip(self.neuron_type, self.neuron_name, self.pos, self.neuron_morph,
+                                                        self.population_unit)):
+                s = "%d,%s,%s,%f,%f,%f,%s, %d\n" % (i, nt, nn, p[0], p[1], p[2], mf, pu)
                 f_out_meta.write(s)
             f_out_meta.close()
 
diff --git a/snudda/utils/pip_upgrade.py b/snudda/utils/pip_upgrade.py
new file mode 100644
index 000000000..6fc0b2680
--- /dev/null
+++ b/snudda/utils/pip_upgrade.py
@@ -0,0 +1,5 @@
+import pkg_resources
+from subprocess import call
+
+packages = [dist.project_name for dist in pkg_resources.working_set]
+call("pip3 install --upgrade " + ' '.join(packages), shell=True)
diff --git a/snudda/utils/save_network_activity.py b/snudda/utils/save_network_activity.py
deleted file mode 100644
index 060df1d35..000000000
--- a/snudda/utils/save_network_activity.py
+++ /dev/null
@@ -1,228 +0,0 @@
-import os.path
-
-import h5py
-import numpy as np
-from mpi4py import MPI  # This must be imported before neuron, to run parallel
-from neuron import h  # , gui
-
-
-class SnuddaSaveNetworkActivity:
-
-    def __init__(self, output_file, network_data=None):
-
-        self.output_file = output_file
-        self.network_data = network_data
-
-        self.pc = h.ParallelContext()
-
-    def spike_sort(self, t_spikes, id_spikes):
-
-        spike_count = dict()
-
-        for t, idx in zip(t_spikes, id_spikes):
-            if idx in spike_count:
-                spike_count[idx] += 1
-            else:
-                spike_count[idx] = 1
-
-        spikes = dict()
-        spike_ctr = dict()
-
-        for idx, num in spike_count.items():
-            spikes[idx] = np.full((num,), np.nan)  # Initialise to NaN to catch any unassigned values
-            spike_ctr[idx] = 0
-
-        for t, idx in zip(t_spikes, id_spikes):
-            spikes[idx][spike_ctr[idx]] = t
-            spike_ctr[idx] += 1
-
-        # Internal consistency
-        for idx in spikes:
-            assert len(spikes[idx]) == 0 or not np.isnan(spikes[idx][-1])
-            assert (idx in spike_ctr and spikes[idx].shape[0] == spike_ctr[idx]) \
-                or (idx not in spike_ctr and len(spikes[idx]) == 0)
-
-        return spikes
-
-    def write(self, t_save, v_save, v_key, t_spikes, id_spikes, output_file=None):
-
-        """ Write spike data and voltage data to output_file
-
-        Args:
-            t_save : array with time
-            v_save : dictionary with arrays of voltage
-            v_key : neuron_id of voltage data
-            t_spikes : spike times
-            id_spikes : neuron_id of spike times
-            """
-
-        if not output_file:
-            output_file = self.output_file
-
-        self.pc.barrier()
-
-        if int(self.pc.id()) == 0:
-
-            if not os.path.isdir(os.path.dirname(output_file)):
-                os.mkdir(os.path.dirname(output_file))
-
-            print(f"Writing network output to {output_file}")
-            out_file = h5py.File(output_file, "w")
-
-            meta_data = out_file.create_group("metaData")
-            out_file.create_group("voltData")
-            out_file.create_group("spikeData")
-
-            if self.network_data:
-                neuron_id = np.array([x["neuronID"] for x in self.network_data["neurons"]])
-                meta_data.create_dataset("ID", data=neuron_id)
-
-                neuron_names = [x["name"] for x in self.network_data["neurons"]]
-                str_type = 'S' + str(max(1, max([len(x) for x in neuron_names])))
-                meta_data.create_dataset("name", (len(neuron_names),), str_type, neuron_names, compression="gzip")
-
-                neuron_types = [x["type"] for x in self.network_data["neurons"]]
-                str_type = 'S' + str(max(1, max([len(x) for x in neuron_types])))
-                meta_data.create_dataset("type", (len(neuron_names),), str_type, neuron_names, compression="gzip")
-
-                swc_list = [n["morphology"] for n in self.network_data["neurons"]]
-                max_swc_len = max([len(x) for x in swc_list])
-                meta_data.create_dataset("morphology", (len(swc_list),), f"S{max_swc_len}",
-                                         swc_list, compression="gzip")
-
-                parameter_keys = [n["parameterKey"] for n in self.network_data["neurons"]]
-                parameter_key_length = max([len(x) for x in parameter_keys])
-                meta_data.create_dataset("parameterKey", (len(parameter_keys),), f"S{parameter_key_length}",
-                                         parameter_keys, compression="gzip")
-
-                morphology_keys = [n["morphologyKey"] for n in self.network_data["neurons"]]
-                morphology_key_length = max(1, max([len(x) for x in morphology_keys]))
-                meta_data.create_dataset("morphologyKey", (len(morphology_keys),), f"S{morphology_key_length}",
-                                         morphology_keys, compression="gzip")
-
-                modulation_keys = [n["modulationKey"] for n in self.network_data["neurons"]]
-                modulation_key_length = max(1, max([len(x) for x in modulation_keys]))
-                meta_data.create_dataset("modulationKey", (len(modulation_keys),), f"S{modulation_key_length}",
-                                         modulation_keys, compression="gzip")
-
-                meta_data.create_dataset("populationUnit", data=self.network_data["populationUnit"], compression="gzip")
-                meta_data.create_dataset("position", data=self.network_data["neuronPositions"], compression="gzip")
-
-            out_file.close()
-
-        if not t_save or not v_save or not v_key:
-            print("No voltage data saved.")
-        else:
-            print("Saving voltage data...")
-            for i in range(int(self.pc.nhost())):
-
-                if i == int(self.pc.id()):
-                    out_file = h5py.File(output_file, "a")
-
-                    if i == 0:
-                        out_file["voltData"].create_dataset("time", data=t_save * 1e-3, compression="gzip")
-
-                    for neuron_id, voltage in zip(v_key, v_save):
-                        out_file["voltData"].create_dataset(str(neuron_id), data=voltage*1e-3, compression="gzip")
-
-                    out_file.close()
-
-                self.pc.barrier()
-
-        # Write spike data
-        print("Sorting spikes")
-        spikes = self.spike_sort(t_spikes=t_spikes, id_spikes=id_spikes)
-
-        print("Saving spike data...")
-
-        for i in range(int(self.pc.nhost())):
-            if i == int(self.pc.id()):
-                out_file = h5py.File(output_file, "a")
-
-                for idx, spike_times in spikes.items():
-                    out_file["spikeData"].create_dataset(f"{idx:.0f}", data=spike_times*1e-3, compression="gzip")
-
-                out_file.close()
-
-            self.pc.barrier()
-
-    def write_currents(self, t_save, i_save, pre_id, post_id, section_id=None, section_x=None, output_file=None):
-
-        """ This adds currents to an already existing hdf5 output file.
-
-        Args:
-            t_save : array with times
-            i_save : list of arrays with current
-            pre_id : array with neuron id of presynaptic neuron
-            post_id : array with neuron id of postsynaptic neuron
-            section_id : section id of synapse current is taken from (optional)
-            section_x : section_x of synapse current is taken from (optional)
-            output_file : output file (optional)
-
-        """
-
-        if not output_file:
-            output_file = self.output_file
-
-        assert os.path.isfile(output_file), f"write_current only appends data to exist file, please use write() first."
-
-        self.pc.barrier()
-
-        # First figure out how much space we need to allocate in the hdf5 file
-        n_cur = len(i_save)
-        n_cur_all = self.pc.py_gather(n_cur, 0)
-
-        if int(self.pc.id()) == 0:
-            n_cur_total = np.sum(n_cur_all)
-
-            out_file = h5py.File(output_file, "a")
-            cur_data = out_file.create_group("currentData")
-            cur_data.create_dataset("current", shape=(len(t_save), n_cur_total), dtype=np.float32, compression="gzip")
-            cur_data.create_dataset("preID", shape=(n_cur_total,), dtype=np.int32)
-            cur_data.create_dataset("postID", shape=(n_cur_total,), dtype=np.int32)
-
-            if section_id is not None:
-                cur_data.create_dataset("sectionID", shape=(n_cur_total,), dtype=np.int32)
-
-            if section_x is not None:
-                cur_data.create_dataset("sectionX", shape=(n_cur_total,), dtype=np.float32)
-
-            out_file.close()
-
-        # Next problem, figure out which column each worker has allocated.
-        pre_id_all = np.concatenate(self.pc.py_allgather(pre_id))
-        post_id_all = np.concatenate(self.pc.py_allgather(post_id))
-        node_all = np.concatenate(self.pc.py_allgather(np.full(pre_id.shape, self.pc.id(), dtype=int)))
-        idx_all = np.concatenate(self.pc.py_allgather(np.arange(0, len(pre_id))))
-
-        pre_post_all = np.vstack([pre_id_all, post_id_all]).T
-        assert pre_post_all.shape[1] == 2, f"Incorrect shape. {pre_post_all.shape}, expected n x 2"
-        sort_idx = np.lexsort(pre_post_all.T)
-
-        my_cols = np.where(node_all[sort_idx] == self.pc.id())[0]
-        my_idx = idx_all[sort_idx][my_cols]
-
-        for i in range(int(self.pc.nhost())):
-            if i == int(self.pc.id()):
-                out_file = h5py.File(output_file, "a")
-
-                cur_data = out_file["currentData"]
-                my_worker_id = self.pc.id()
-
-                for idx, col_id in zip(my_idx, my_cols):
-                    assert node_all[idx] == my_worker_id, \
-                        f"Problem with sorting. Worker {my_worker_id} found data for worker {node_all[idx]}."
-
-                    cur_data["current"][:, col_id] = i_save[idx]
-                    cur_data["preID"][col_id] = pre_id[idx]
-                    cur_data["postID"][col_id] = post_id[idx]
-
-                    if section_id is not None:
-                        cur_data["sectionID"][col_id] = section_id[idx]
-
-                    if section_x is not None:
-                        cur_data["sectionX"][col_id] = section_x[idx]
-
-                out_file.close()
-
-            self.pc.barrier()
diff --git a/snudda/utils/save_network_recording.py b/snudda/utils/save_network_recording.py
new file mode 100644
index 000000000..ba68d601d
--- /dev/null
+++ b/snudda/utils/save_network_recording.py
@@ -0,0 +1,413 @@
+import os.path
+
+import h5py
+import numpy as np
+from mpi4py import MPI  # This must be imported before neuron, to run parallel
+from neuron import h  # , gui
+
+# TODO:
+
+# Check if we can get synapse row number stored somehow...
+
+# 4. Hur ska spikes sparas? (cell id, tid)
+# 5. Gör en generell write funktion som tittar på vilken data som finns och anropar rätt sub-funktioner
+# -- Skriv tutorial hur man lägger till en egen recording
+#
+# Förslag på namn: TimeSeriesData
+#                  EventData
+
+# NeuronActivity håller information om alla olika mätningar för en given neuron,
+# Simulate har en lista med NeuronActivity
+#
+# I simulate måste vi:
+# För varje neuron som ska mäta något på, skapa en NeuronActivity och lägg i en lista
+# För varje measurement vi vill spara (ström, voltage, etc), använd add_measurements
+# och sedan
+
+
+class NeuronRecordings:
+    """
+        Container class for all recordings associated with a neuron.
+    """
+    def __init__(self, neuron_id):
+        """
+        Constructor
+
+        Args:
+            neuron_id  (int): Neuron ID
+        """
+
+        self.neuron_id = neuron_id
+
+        self.data = dict()
+
+    def register_compartment_data(self, data, data_type, sec_id, sec_x):
+        """
+        Adds a new measurement.
+            
+        Args:
+            data (neuron.h.Vector): NEURON vector holding recording.
+            data_type (str): Name of the tracked data.
+            sec_id (int): Section ID
+            sec_x (float): Section X (segment location)
+
+        """
+
+        if data_type not in self.data:
+            self.data[data_type] = CompartmentData(neuron_id=self.neuron_id, data_type=data_type)
+
+        self.data[data_type].append(data=data, sec_id=sec_id, sec_x=sec_x)
+
+    def register_synapse_data(self, data, data_type, synapse_type, presynaptic_id, sec_id, sec_x, cond):
+
+        if data_type not in self.data:
+            self.data[data_type] = SynapseData(neuron_id=self.neuron_id, data_type=data_type)
+
+        self.data[data_type].append(data=data, synapse_type=synapse_type, presynaptic_id=presynaptic_id,
+                                    sec_id=sec_id, sec_x=sec_x, cond=cond)
+
+
+class SynapseData:
+
+    def __init__(self, neuron_id, data_type):
+
+        self.neuron_id = neuron_id
+        self.data = []
+        self.data_type = data_type
+        self.sec_id = []
+        self.sec_x = []
+        self.synapse_type = []
+        self.presynaptic_id = []
+        self.cond = []
+
+    def append(self, data, synapse_type, presynaptic_id, sec_id, sec_x, cond):
+        self.data.append(data)
+        self.synapse_type.append(synapse_type)
+        self.sec_id.append(sec_id)
+        self.sec_x.append(sec_x)
+        self.presynaptic_id.append(presynaptic_id)
+        self.cond.append(cond)
+
+    def to_numpy(self):
+        """
+            Returns:
+                (np.ndarray): Data represented as np.ndarrays
+        """
+        return np.vstack([np.array(d) for d in self.data])
+
+
+class CompartmentData:
+    """
+        Container class for recordings from a compartment. 
+    """
+
+    def __init__(self, neuron_id, data_type):
+        """
+        Constructor
+
+        Args:
+            neuron_id (int): Neuron ID
+            data_type (str): Name of the tracked data.
+        """
+        self.neuron_id = neuron_id
+        self.data_type = data_type
+        self.data = []
+        self.sec_id = []
+        self.sec_x = []
+
+    def append(self, data, sec_id, sec_x):
+
+        # !!! Issue (?): The same compartment can hold several recordings now.
+
+        """
+        Appends a recording for a given compartment.
+
+        Args:
+            data (neuron.h.Vector): NEURON vector for the recording.
+            sec_id (int): Section ID.
+            sec_x (float): Section X (segment location). 
+
+        """
+        self.data.append(data)
+        
+        self.sec_id.append(sec_id)
+        self.sec_x.append(sec_x)
+
+    def to_numpy(self):  # Misnomer? (original data is preserved). (Instead "as_ndarray" (like in NEURON)?)
+        """
+            Returns:
+                (np.ndarray): Data represented as np.ndarrays 
+        """
+        return np.vstack([np.array(d) for d in self.data])
+
+
+class SnuddaSaveNetworkRecordings:
+
+    def __init__(self, output_file, network_data=None):
+
+        self.output_file = output_file
+        self.network_data = network_data
+        self.header_exists = False
+        self.neuron_activities = dict()
+        self.time = None
+
+        self.units = dict()
+
+        self.pc = h.ParallelContext()
+
+    def add_unit(self, data_type, target_unit, conversion_factor):
+        # Units reference chart: https://www.neuron.yale.edu/neuron/static/docs/units/unitchart.html
+        self.units[data_type] = {"unit": target_unit, "conversion_factor": conversion_factor}
+
+    def get_conversion(self, data_type):
+        if data_type in self.units:
+            return self.units[data_type]["conversion_factor"]
+        else:
+            return 1.0
+
+    def register_compartment_data(self, data_type, neuron_id, data, sec_id, sec_x):
+        if neuron_id not in self.neuron_activities:
+            self.neuron_activities[neuron_id] = NeuronRecordings(neuron_id)
+
+        self.neuron_activities[neuron_id].register_compartment_data(data=data, data_type=data_type,
+                                                                    sec_id=sec_id, sec_x=sec_x)
+
+    def register_synapse_data(self, data_type, neuron_id, data, synapse_type, presynaptic_id, sec_id, sec_x, cond):
+        if neuron_id not in self.neuron_activities:
+            self.neuron_activities[neuron_id] = NeuronRecordings(neuron_id)
+
+        self.neuron_activities[neuron_id].register_synapse_data(data=data, data_type=data_type,
+                                                                synapse_type=synapse_type,
+                                                                presynaptic_id=presynaptic_id,
+                                                                sec_id=sec_id, sec_x=sec_x, cond=cond)
+
+    def register_time(self, time):
+        self.time = time
+
+    @staticmethod
+    def spike_sort(t_spikes, id_spikes):
+
+        spike_count = dict()
+
+        for t, idx in zip(t_spikes, id_spikes):
+            if idx in spike_count:
+                spike_count[idx] += 1
+            else:
+                spike_count[idx] = 1
+
+        spikes = dict()
+        spike_ctr = dict()
+
+        for idx, num in spike_count.items():
+            spikes[idx] = np.full((num,), np.nan)  # Initialise to NaN to catch any unassigned values
+            spike_ctr[idx] = 0
+
+        for t, idx in zip(t_spikes, id_spikes):
+            spikes[idx][spike_ctr[idx]] = t
+            spike_ctr[idx] += 1
+
+        # Internal consistency
+        for idx in spikes:
+            assert len(spikes[idx]) == 0 or not np.isnan(spikes[idx][-1])
+            assert (idx in spike_ctr and spikes[idx].shape[0] == spike_ctr[idx]) \
+                or (idx not in spike_ctr and len(spikes[idx]) == 0)
+
+        return spikes
+
+    # TODO: save_network_activity
+    # 1. write_header -- node 0 skriver metadata till filen
+    # 2. write_spikes -- alla noder skriver sin spikdata
+    # 3. write_soma_voltage -- alla noder skriver soma voltage data
+    # 4. write_compartment_voltage -- alla noder skriver segment voltage
+
+    def write_string_meta_data(self, group, name):
+        string_data = [x[name] for x in self.network_data["neurons"]]
+        max_len = max(1, max([len(x) for x in string_data]))
+        str_type = f"S{max_len}"
+        group.create_dataset(name, (len(string_data),), str_type, string_data, compression="gzip")
+
+    def write_header(self):
+
+        if self.header_exists:
+            return
+
+        self.pc.barrier()
+
+        if int(self.pc.id()) == 0:
+
+            if not os.path.isdir(os.path.dirname(self.output_file)):
+                os.mkdir(os.path.dirname(self.output_file))
+
+            print(f"Writing network output to {self.output_file}")
+            out_file = h5py.File(self.output_file, "w")
+
+            meta_data = out_file.create_group("metaData")
+            out_file.create_group("neurons")
+
+            if self.network_data:
+                neuron_id = np.array([x["neuronID"] for x in self.network_data["neurons"]])
+                meta_data.create_dataset("ID", data=neuron_id)
+
+                for name in ["name", "type", "morphology", "parameterKey", "morphologyKey", "modulationKey"]:
+                    self.write_string_meta_data(group=meta_data, name=name)
+
+                meta_data.create_dataset("populationUnit", data=self.network_data["populationUnit"], compression="gzip")
+                meta_data.create_dataset("position", data=self.network_data["neuronPositions"], compression="gzip")
+
+            out_file.close()
+
+        self.header_exists = True
+
+    def write(self):
+
+        self.write_header()
+
+        if int(self.pc.id()) == 0:
+
+            out_file = h5py.File(self.output_file, "a")
+            if "time" not in out_file:
+                out_file.create_dataset("time", data=self.time * self.get_conversion("time"))
+                out_file.close()
+
+        for i in range(int(self.pc.nhost())):
+            out_file = h5py.File(self.output_file, "a")
+
+            for na in self.neuron_activities.values():
+                neuron_id_str = str(na.neuron_id)
+                if neuron_id_str not in out_file["neurons"]:
+                    out_file["neurons"].create_group(neuron_id_str)
+
+                for m in na.data.values():
+
+                    conversion_factor = self.get_conversion(m.data_type)
+
+                    out_file["neurons"][neuron_id_str].create_group(m.data_type)
+                    data_group = out_file["neurons"][neuron_id_str][m.data_type]
+                    data_group.create_dataset("data", data=m.to_numpy() * conversion_factor, compression="gzip")
+                    data_group.create_dataset("sec_id", data=np.array(m.sec_id), compression="gzip")
+                    data_group.create_dataset("sec_x", data=np.array(m.sec_x), compression="gzip")
+
+                    if isinstance(m, SynapseData):
+                        data_group.create_dataset("synapse_type", data=np.array(m.synapse_type), compression="gzip")
+                        data_group.create_dataset("presynaptic_id", data=np.array(m.presynaptic_id), compression="gzip")
+                        data_group.create_dataset("cond", data=np.array(m.cond), compression="gzip")
+
+            out_file.close()
+
+        self.pc.barrier()
+
+    def write_OLD(self, t_save, v_save, v_key, t_spikes, id_spikes, output_file=None):
+
+        """ Write spike data and voltage data to output_file
+
+        Args:
+            t_save : array with time
+            v_save : dictionary with arrays of voltage
+            v_key : neuron_id of voltage data
+            t_spikes : spike times
+            id_spikes : neuron_id of spike times
+        """
+
+        self.write_header()
+
+        if not t_save or not v_save or not v_key:
+            print("No voltage data saved.")
+        else:
+            print("Saving voltage data...")
+            for i in range(int(self.pc.nhost())):
+
+                if i == int(self.pc.id()):
+                    out_file = h5py.File(output_file, "a")
+
+                    if i == 0:
+                        out_file["voltData"].create_dataset("time", data=t_save * 1e-3, compression="gzip")
+
+                    for neuron_id, voltage in zip(v_key, v_save):
+                        out_file["voltData"].create_dataset(str(neuron_id), data=voltage*1e-3, compression="gzip")
+
+                    out_file.close()
+
+                self.pc.barrier()
+
+        self.write_spikes_OLD(t_spikes, id_spikes)
+
+    def write_currents_OLD(self, t_save, i_save, pre_id, post_id, section_id=None, section_x=None, output_file=None):
+
+        """ This adds currents to an already existing hdf5 output file.
+
+        Args:
+            t_save : array with times
+            i_save : list of arrays with current
+            pre_id : array with neuron id of presynaptic neuron
+            post_id : array with neuron id of postsynaptic neuron
+            section_id : section id of synapse current is taken from (optional)
+            section_x : section_x of synapse current is taken from (optional)
+            output_file : output file (optional)
+
+        """
+
+        if not output_file:
+            output_file = self.output_file
+
+        assert os.path.isfile(output_file), f"write_current only appends data to exist file, please use write() first."
+
+        self.pc.barrier()
+
+        # First figure out how much space we need to allocate in the hdf5 file
+        n_cur = len(i_save)
+        n_cur_all = self.pc.py_gather(n_cur, 0)
+
+        if int(self.pc.id()) == 0:
+            n_cur_total = np.sum(n_cur_all)
+
+            out_file = h5py.File(output_file, "a")
+            cur_data = out_file.create_group("currentData")
+            cur_data.create_dataset("current", shape=(len(t_save), n_cur_total), dtype=np.float32, compression="gzip")
+            cur_data.create_dataset("preID", shape=(n_cur_total,), dtype=np.int32)
+            cur_data.create_dataset("postID", shape=(n_cur_total,), dtype=np.int32)
+
+            if section_id is not None:
+                cur_data.create_dataset("sectionID", shape=(n_cur_total,), dtype=np.int32)
+
+            if section_x is not None:
+                cur_data.create_dataset("sectionX", shape=(n_cur_total,), dtype=np.float32)
+
+            out_file.close()
+
+        # Next problem, figure out which column each worker has allocated.
+        pre_id_all = np.concatenate(self.pc.py_allgather(pre_id))
+        post_id_all = np.concatenate(self.pc.py_allgather(post_id))
+        node_all = np.concatenate(self.pc.py_allgather(np.full(pre_id.shape, self.pc.id(), dtype=int)))
+        idx_all = np.concatenate(self.pc.py_allgather(np.arange(0, len(pre_id))))
+
+        pre_post_all = np.vstack([pre_id_all, post_id_all]).T
+        assert pre_post_all.shape[1] == 2, f"Incorrect shape. {pre_post_all.shape}, expected n x 2"
+        sort_idx = np.lexsort(pre_post_all.T)
+
+        my_cols = np.where(node_all[sort_idx] == self.pc.id())[0]
+        my_idx = idx_all[sort_idx][my_cols]
+
+        for i in range(int(self.pc.nhost())):
+            if i == int(self.pc.id()):
+                out_file = h5py.File(output_file, "a")
+
+                cur_data = out_file["currentData"]
+                my_worker_id = self.pc.id()
+
+                for idx, col_id in zip(my_idx, my_cols):
+                    assert node_all[idx] == my_worker_id, \
+                        f"Problem with sorting. Worker {my_worker_id} found data for worker {node_all[idx]}."
+
+                    cur_data["current"][:, col_id] = i_save[idx]
+                    cur_data["preID"][col_id] = pre_id[idx]
+                    cur_data["postID"][col_id] = post_id[idx]
+
+                    if section_id is not None:
+                        cur_data["sectionID"][col_id] = section_id[idx]
+
+                    if section_x is not None:
+                        cur_data["sectionX"][col_id] = section_x[idx]
+
+                out_file.close()
+
+            self.pc.barrier()


Test patch:
diff --git a/tests/data/ExampleProjectionDetection.json b/tests/data/ExampleProjectionDetection.json
new file mode 100644
index 000000000..1164928a9
--- /dev/null
+++ b/tests/data/ExampleProjectionDetection.json
@@ -0,0 +1,32 @@
+{
+    "dSPN,iSPN,GABA" : {
+	"source" : [[4850, -150, -150],
+                    [4850, -150, 150],
+                    [4850, 150, -150],
+                    [4850, 150, 150],
+                    [5150, -150, -150],
+                    [5150, -150, 150],
+                    [5150, 150, -150],
+                    [5150, 150, 150]],
+	
+	"destination" : [[-5150, -150, -150],
+                         [-5150, -150, 150],
+                         [-5150, 150, -150],
+                         [-5150, 150, 150],
+                         [-4850, -150, -150],
+                         [-4850, -150, 150],
+                         [-4850, 150, -150],
+                         [-4850, 150, 150]],
+	"radius" : 30e-6,
+	"numPoints" : 10000,
+	"rotation" : [[1, 0, 0, 0, 1, 0, 0, 0, 1],
+		      [1, 0, 0, 0, 1, 0, 0, 0, 1],
+		      [1, 0, 0, 0, 1, 0, 0, 0, 1],
+		      [1, 0, 0, 0, 1, 0, 0, 0, 1],
+		      [1, 0, 0, 0, 1, 0, 0, 0, 1],
+		      [1, 0, 0, 0, 1, 0, 0, 0, 1],
+		      [1, 0, 0, 0, 1, 0, 0, 0, 1],
+		      [1, 0, 0, 0, 1, 0, 0, 0, 1]]
+    }
+    
+}
diff --git a/tests/networks/network_testing_detect_cluster/network-config.json b/tests/networks/network_testing_detect_cluster/network-config.json
new file mode 100644
index 000000000..4a4555316
--- /dev/null
+++ b/tests/networks/network_testing_detect_cluster/network-config.json
@@ -0,0 +1,49 @@
+{
+    "RandomSeed": {
+        "masterseed": 1234,
+        "init": 2906597030,
+        "place": 1602421836,
+        "detect": 216676975,
+        "prune": 2698621808,
+        "input": 507409703,
+        "simulate": 2825158027
+    },
+    "Volume": {
+	"TestStriatum": {
+	    "type": "mesh",
+	    "dMin": 1e-5,
+	    "meshFile": "mesh/simple_mesh.obj",
+	    "meshBinWidth": 5e-5
+	}
+    },
+    "Connectivity": {
+	"ballandstick,ballandstick": {
+	    "GABA": {
+		"conductance": [1.1e-9, 1.5e-9],
+		"clusterSize": 5,
+		"clusterSpread": 25e-6,
+                "pruning": {
+                    "f1": 1,
+                    "softMax": 5,
+                    "mu2": null,
+                    "a3": 1,
+                    "distPruning": null
+                }
+	    }
+	}
+    },
+    "Neurons": {
+	"ballandstick_0": {
+            "morphology": "validation/ballanddoublestick/double.swc",
+            "parameters": "",
+            "mechanisms": "",
+            "modulation": "",
+            "num": 8,
+            "hoc": null,
+            "neuronType": "neuron",
+            "rotationMode": "",
+            "volumeID": "TestStriatum"
+
+	}
+    }
+}
diff --git a/tests/networks/network_testing_init_cube/network-config.json b/tests/networks/network_testing_init_cube/network-config.json
index 3499a31c1..7f49fea24 100644
--- a/tests/networks/network_testing_init_cube/network-config.json
+++ b/tests/networks/network_testing_init_cube/network-config.json
@@ -1,12 +1,12 @@
 {
     "RandomSeed": {
-        "init": 247808882,
-        "place": 799680898,
-        "detect": 1560220428,
-        "project": 1228648387,
-        "prune": 299996022,
-        "input": 1052838105,
-        "simulate": 2972898657
+        "init": 3714952339,
+        "place": 1974606074,
+        "detect": 4039039081,
+        "project": 458631388,
+        "prune": 675094707,
+        "input": 1495983570,
+        "simulate": 3200261559
     },
     "Volume": {
         "Striatum": {
@@ -154,7 +154,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 0.38,
+                    "f1": 0.28500000000000003,
                     "softMax": 3,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -162,7 +162,7 @@
                     "cluster": false
                 },
                 "pruningOther": {
-                    "f1": 0.38,
+                    "f1": 0.28500000000000003,
                     "softMax": 3,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -191,7 +191,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 0.2,
+                    "f1": 0.164,
                     "softMax": 3,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -199,7 +199,7 @@
                     "cluster": false
                 },
                 "pruningOther": {
-                    "f1": 0.2,
+                    "f1": 0.164,
                     "softMax": 3,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -248,7 +248,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 0.3,
+                    "f1": 0.279,
                     "softMax": 4,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -256,7 +256,7 @@
                     "cluster": false
                 },
                 "pruningOther": {
-                    "f1": 0.3,
+                    "f1": 0.279,
                     "softMax": 4,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -340,7 +340,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 1.0,
+                    "f1": 0.3,
                     "softMax": 15,
                     "mu2": 3,
                     "a3": 0.3,
@@ -367,7 +367,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 1.0,
+                    "f1": 0.3,
                     "softMax": 15,
                     "mu2": 3,
                     "a3": 0.3,
@@ -398,59 +398,54 @@
     },
     "Neurons": {
         "FS_0": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20210210/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20210210/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20210210/mechanisms.json",
-            "num": 260,
+            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20190301/MTC251001A-IDB-cor-rep.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20190301/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20190301/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20190301/modulation.json",
+            "num": 325,
             "hoc": null,
             "neuronType": "neuron",
             "rotationMode": "random",
             "volumeID": "Striatum"
         },
         "FS_1": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20210210/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20210210/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20210210/mechanisms.json",
-            "num": 260,
+            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20190312/MTC180800A-IDB-cor-rep.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20190312/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20190312/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20190312/modulation.json",
+            "num": 325,
             "hoc": null,
             "neuronType": "neuron",
             "rotationMode": "random",
             "volumeID": "Striatum"
         },
         "FS_2": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20210211/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20210211/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20210211/mechanisms.json",
-            "num": 260,
+            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20190225/DR-rat-Mar-13-08-1-536-R-cor-rep.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20190225/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20190225/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20190225/modulation.json",
+            "num": 325,
             "hoc": null,
             "neuronType": "neuron",
             "rotationMode": "random",
             "volumeID": "Striatum"
         },
         "FS_3": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20210210/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20210210/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20210210/mechanisms.json",
-            "num": 260,
-            "hoc": null,
-            "neuronType": "neuron",
-            "rotationMode": "random",
-            "volumeID": "Striatum"
-        },
-        "FS_4": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mBE104E-v20210209/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mBE104E-v20210209/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mBE104E-v20210209/mechanisms.json",
-            "num": 260,
+            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20190226/MTC180800A-IDB-cor-rep.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20190226/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20190226/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20190226/modulation.json",
+            "num": 325,
             "hoc": null,
             "neuronType": "neuron",
             "rotationMode": "random",
             "volumeID": "Striatum"
         },
         "dSPN_0": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_D1-mWT-1215MSN03-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_D1-mWT-1215MSN03-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_D1-mWT-1215MSN03-v20211026/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20190508/WT-0728MSN01-cor-rep-ax.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20190508/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20190508/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20190508/modulation.json",
             "num": 11875,
             "hoc": null,
             "neuronType": "neuron",
@@ -458,9 +453,10 @@
             "volumeID": "Striatum"
         },
         "dSPN_1": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20211028/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20211028/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20211028/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20190521/WT-P270-20-15ak-cor.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20190521/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20190521/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20190521/modulation.json",
             "num": 11875,
             "hoc": null,
             "neuronType": "neuron",
@@ -468,9 +464,10 @@
             "volumeID": "Striatum"
         },
         "dSPN_2": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20211026/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_d1-mWT-1215MSN03-v20190521/WT-1215MSN03-cor-rep-ax2.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_d1-mWT-1215MSN03-v20190521/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_d1-mWT-1215MSN03-v20190521/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_d1-mWT-1215MSN03-v20190521/modulation.json",
             "num": 11875,
             "hoc": null,
             "neuronType": "neuron",
@@ -478,9 +475,10 @@
             "volumeID": "Striatum"
         },
         "dSPN_3": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20211026/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20190503/21-6-DE-cor-rep-ax.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20190503/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20190503/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20190503/modulation.json",
             "num": 11875,
             "hoc": null,
             "neuronType": "neuron",
@@ -488,9 +486,10 @@
             "volumeID": "Striatum"
         },
         "iSPN_0": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20211026/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20190603/WT-MSN1-cor-rep-ax.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20190603/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20190603/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20190603/modulation.json",
             "num": 11875,
             "hoc": null,
             "neuronType": "neuron",
@@ -498,9 +497,10 @@
             "volumeID": "Striatum"
         },
         "iSPN_1": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20211026/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20190529/46-3-DE-cor-rep-ax.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20190529/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20190529/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20190529/modulation.json",
             "num": 11875,
             "hoc": null,
             "neuronType": "neuron",
@@ -508,9 +508,10 @@
             "volumeID": "Striatum"
         },
         "iSPN_2": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20211026/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20190527/WT-P270-09-15ak-cor.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20190527/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20190527/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20190527/modulation.json",
             "num": 11875,
             "hoc": null,
             "neuronType": "neuron",
@@ -518,9 +519,10 @@
             "volumeID": "Striatum"
         },
         "iSPN_3": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20211026/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20190611/51-5-DE-cor-rep-ax.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20190611/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20190611/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20190611/modulation.json",
             "num": 11875,
             "hoc": null,
             "neuronType": "neuron",
diff --git a/tests/networks/network_testing_init_full/network-config.json b/tests/networks/network_testing_init_full/network-config.json
index 16a2f2c5a..43c7c3ee5 100644
--- a/tests/networks/network_testing_init_full/network-config.json
+++ b/tests/networks/network_testing_init_full/network-config.json
@@ -172,7 +172,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 0.38,
+                    "f1": 0.28500000000000003,
                     "softMax": 3,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -180,7 +180,7 @@
                     "cluster": false
                 },
                 "pruningOther": {
-                    "f1": 0.38,
+                    "f1": 0.28500000000000003,
                     "softMax": 3,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -209,7 +209,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 0.2,
+                    "f1": 0.164,
                     "softMax": 3,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -217,7 +217,7 @@
                     "cluster": false
                 },
                 "pruningOther": {
-                    "f1": 0.2,
+                    "f1": 0.164,
                     "softMax": 3,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -266,7 +266,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 0.3,
+                    "f1": 0.279,
                     "softMax": 4,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -274,7 +274,7 @@
                     "cluster": false
                 },
                 "pruningOther": {
-                    "f1": 0.3,
+                    "f1": 0.279,
                     "softMax": 4,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -358,7 +358,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 1.0,
+                    "f1": 0.3,
                     "softMax": 15,
                     "mu2": 3,
                     "a3": 0.3,
@@ -385,7 +385,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 1.0,
+                    "f1": 0.3,
                     "softMax": 15,
                     "mu2": 3,
                     "a3": 0.3,
@@ -416,69 +416,65 @@
     },
     "Neurons": {
         "FS_0": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20210210/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20210210/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20210210/mechanisms.json",
-            "num": 4426,
+            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20190301/MTC251001A-IDB-cor-rep.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20190301/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20190301/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20190301/modulation.json",
+            "num": 5532,
             "hoc": null,
             "neuronType": "neuron",
             "rotationMode": "random",
             "volumeID": "Striatum"
         },
         "FS_1": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20210210/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20210210/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20210210/mechanisms.json",
-            "num": 4426,
+            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20190312/MTC180800A-IDB-cor-rep.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20190312/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20190312/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20190312/modulation.json",
+            "num": 5533,
             "hoc": null,
             "neuronType": "neuron",
             "rotationMode": "random",
             "volumeID": "Striatum"
         },
         "FS_2": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20210211/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20210211/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20210211/mechanisms.json",
-            "num": 4426,
+            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20190225/DR-rat-Mar-13-08-1-536-R-cor-rep.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20190225/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20190225/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20190225/modulation.json",
+            "num": 5533,
             "hoc": null,
             "neuronType": "neuron",
             "rotationMode": "random",
             "volumeID": "Striatum"
         },
         "FS_3": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20210210/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20210210/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20210210/mechanisms.json",
-            "num": 4426,
-            "hoc": null,
-            "neuronType": "neuron",
-            "rotationMode": "random",
-            "volumeID": "Striatum"
-        },
-        "FS_4": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mBE104E-v20210209/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mBE104E-v20210209/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mBE104E-v20210209/mechanisms.json",
-            "num": 4426,
+            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20190226/MTC180800A-IDB-cor-rep.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20190226/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20190226/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20190226/modulation.json",
+            "num": 5532,
             "hoc": null,
             "neuronType": "neuron",
             "rotationMode": "random",
             "volumeID": "Striatum"
         },
         "dSPN_0": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_D1-mWT-1215MSN03-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_D1-mWT-1215MSN03-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_D1-mWT-1215MSN03-v20211026/mechanisms.json",
-            "num": 202153,
+            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20190508/WT-0728MSN01-cor-rep-ax.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20190508/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20190508/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20190508/modulation.json",
+            "num": 202154,
             "hoc": null,
             "neuronType": "neuron",
             "rotationMode": "random",
             "volumeID": "Striatum"
         },
         "dSPN_1": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20211028/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20211028/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20211028/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20190521/WT-P270-20-15ak-cor.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20190521/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20190521/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20190521/modulation.json",
             "num": 202153,
             "hoc": null,
             "neuronType": "neuron",
@@ -486,19 +482,21 @@
             "volumeID": "Striatum"
         },
         "dSPN_2": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20211026/mechanisms.json",
-            "num": 202154,
+            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_d1-mWT-1215MSN03-v20190521/WT-1215MSN03-cor-rep-ax2.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_d1-mWT-1215MSN03-v20190521/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_d1-mWT-1215MSN03-v20190521/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_d1-mWT-1215MSN03-v20190521/modulation.json",
+            "num": 202153,
             "hoc": null,
             "neuronType": "neuron",
             "rotationMode": "random",
             "volumeID": "Striatum"
         },
         "dSPN_3": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20211026/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20190503/21-6-DE-cor-rep-ax.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20190503/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20190503/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20190503/modulation.json",
             "num": 202154,
             "hoc": null,
             "neuronType": "neuron",
@@ -506,19 +504,21 @@
             "volumeID": "Striatum"
         },
         "iSPN_0": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20211026/mechanisms.json",
-            "num": 202153,
+            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20190603/WT-MSN1-cor-rep-ax.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20190603/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20190603/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20190603/modulation.json",
+            "num": 202154,
             "hoc": null,
             "neuronType": "neuron",
             "rotationMode": "random",
             "volumeID": "Striatum"
         },
         "iSPN_1": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20211026/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20190529/46-3-DE-cor-rep-ax.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20190529/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20190529/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20190529/modulation.json",
             "num": 202153,
             "hoc": null,
             "neuronType": "neuron",
@@ -526,9 +526,10 @@
             "volumeID": "Striatum"
         },
         "iSPN_2": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20211026/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20190527/WT-P270-09-15ak-cor.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20190527/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20190527/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20190527/modulation.json",
             "num": 202154,
             "hoc": null,
             "neuronType": "neuron",
@@ -536,19 +537,21 @@
             "volumeID": "Striatum"
         },
         "iSPN_3": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20211026/mechanisms.json",
-            "num": 202154,
+            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20190611/51-5-DE-cor-rep-ax.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20190611/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20190611/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20190611/modulation.json",
+            "num": 202153,
             "hoc": null,
             "neuronType": "neuron",
             "rotationMode": "random",
             "volumeID": "Striatum"
         },
         "ChIN_0": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/chin/str-chin-e170614_cell6-m17JUL301751_170614_no6_MD_cell_1_x63-v20190710/morphology",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/chin/str-chin-e170614_cell6-m17JUL301751_170614_no6_MD_cell_1_x63-v20190710/optim_chin_morph_renamed2019-11-08.swc",
             "parameters": "$SNUDDA_DATA/neurons/striatum/chin/str-chin-e170614_cell6-m17JUL301751_170614_no6_MD_cell_1_x63-v20190710/parameters.json",
             "mechanisms": "$SNUDDA_DATA/neurons/striatum/chin/str-chin-e170614_cell6-m17JUL301751_170614_no6_MD_cell_1_x63-v20190710/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/chin/str-chin-e170614_cell6-m17JUL301751_170614_no6_MD_cell_1_x63-v20190710/modulation.json",
             "num": 18726,
             "hoc": null,
             "neuronType": "neuron",
@@ -561,9 +564,10 @@
             ]
         },
         "LTS_0": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/lts/LTS_180118_morp_9862_updated_20210301/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/lts/LTS_180118_morp_9862_updated_20210301/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/lts/LTS_180118_morp_9862_updated_20210301/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/lts/LTS_Experiment-9862_20181211/lts_morp_2019-11-07_centered_no_axon.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/lts/LTS_Experiment-9862_20181211/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/lts/LTS_Experiment-9862_20181211/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/lts/LTS_Experiment-9862_20181211/modulation.json",
             "num": 11916,
             "hoc": null,
             "neuronType": "neuron",
diff --git a/tests/networks/network_testing_init_slice/network-config.json b/tests/networks/network_testing_init_slice/network-config.json
index a706adcfb..5e4efc325 100644
--- a/tests/networks/network_testing_init_slice/network-config.json
+++ b/tests/networks/network_testing_init_slice/network-config.json
@@ -1,12 +1,12 @@
 {
     "RandomSeed": {
-        "init": 2350074332,
-        "place": 2771748326,
-        "detect": 1845204186,
-        "project": 1679767388,
-        "prune": 3303766660,
-        "input": 4237944138,
-        "simulate": 4081984847
+        "init": 2832289324,
+        "place": 556803843,
+        "detect": 643836513,
+        "project": 3602177942,
+        "prune": 4265081994,
+        "input": 3814079397,
+        "simulate": 2826350743
     },
     "Volume": {
         "Striatum": {
@@ -154,7 +154,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 0.38,
+                    "f1": 0.28500000000000003,
                     "softMax": 3,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -162,7 +162,7 @@
                     "cluster": false
                 },
                 "pruningOther": {
-                    "f1": 0.38,
+                    "f1": 0.28500000000000003,
                     "softMax": 3,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -191,7 +191,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 0.2,
+                    "f1": 0.164,
                     "softMax": 3,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -199,7 +199,7 @@
                     "cluster": false
                 },
                 "pruningOther": {
-                    "f1": 0.2,
+                    "f1": 0.164,
                     "softMax": 3,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -248,7 +248,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 0.3,
+                    "f1": 0.279,
                     "softMax": 4,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -256,7 +256,7 @@
                     "cluster": false
                 },
                 "pruningOther": {
-                    "f1": 0.3,
+                    "f1": 0.279,
                     "softMax": 4,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -340,7 +340,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 1.0,
+                    "f1": 0.3,
                     "softMax": 15,
                     "mu2": 3,
                     "a3": 0.3,
@@ -367,7 +367,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 1.0,
+                    "f1": 0.3,
                     "softMax": 15,
                     "mu2": 3,
                     "a3": 0.3,
@@ -398,59 +398,54 @@
     },
     "Neurons": {
         "FS_0": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20210210/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20210210/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20210210/mechanisms.json",
-            "num": 260,
+            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20190301/MTC251001A-IDB-cor-rep.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20190301/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20190301/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20190301/modulation.json",
+            "num": 325,
             "hoc": null,
             "neuronType": "neuron",
             "rotationMode": "random",
             "volumeID": "Striatum"
         },
         "FS_1": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20210210/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20210210/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20210210/mechanisms.json",
-            "num": 260,
+            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20190312/MTC180800A-IDB-cor-rep.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20190312/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20190312/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20190312/modulation.json",
+            "num": 325,
             "hoc": null,
             "neuronType": "neuron",
             "rotationMode": "random",
             "volumeID": "Striatum"
         },
         "FS_2": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20210211/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20210211/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e180418_FS5-mMTC251001A-IDB-v20210211/mechanisms.json",
-            "num": 260,
+            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20190225/DR-rat-Mar-13-08-1-536-R-cor-rep.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20190225/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20190225/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161024_FS16-mDR-rat-Mar-13-08-1-536-R-v20190225/modulation.json",
+            "num": 325,
             "hoc": null,
             "neuronType": "neuron",
             "rotationMode": "random",
             "volumeID": "Striatum"
         },
         "FS_3": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20210210/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20210210/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mMTC180800A-IDB-v20210210/mechanisms.json",
-            "num": 260,
-            "hoc": null,
-            "neuronType": "neuron",
-            "rotationMode": "random",
-            "volumeID": "Striatum"
-        },
-        "FS_4": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mBE104E-v20210209/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mBE104E-v20210209/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e161205_FS1-mBE104E-v20210209/mechanisms.json",
-            "num": 260,
+            "morphology": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20190226/MTC180800A-IDB-cor-rep.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20190226/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20190226/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/fs/str-fs-e160628_FS2-mMTC180800A-IDB-v20190226/modulation.json",
+            "num": 325,
             "hoc": null,
             "neuronType": "neuron",
             "rotationMode": "random",
             "volumeID": "Striatum"
         },
         "dSPN_0": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_D1-mWT-1215MSN03-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_D1-mWT-1215MSN03-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_D1-mWT-1215MSN03-v20211026/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20190508/WT-0728MSN01-cor-rep-ax.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20190508/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20190508/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20190508/modulation.json",
             "num": 11875,
             "hoc": null,
             "neuronType": "neuron",
@@ -458,9 +453,10 @@
             "volumeID": "Striatum"
         },
         "dSPN_1": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20211028/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20211028/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20211028/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20190521/WT-P270-20-15ak-cor.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20190521/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20190521/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20190521/modulation.json",
             "num": 11875,
             "hoc": null,
             "neuronType": "neuron",
@@ -468,9 +464,10 @@
             "volumeID": "Striatum"
         },
         "dSPN_2": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150602_c1_D1-mWT-0728MSN01-v20211026/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_d1-mWT-1215MSN03-v20190521/WT-1215MSN03-cor-rep-ax2.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_d1-mWT-1215MSN03-v20190521/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_d1-mWT-1215MSN03-v20190521/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c9_d1-mWT-1215MSN03-v20190521/modulation.json",
             "num": 11875,
             "hoc": null,
             "neuronType": "neuron",
@@ -478,9 +475,10 @@
             "volumeID": "Striatum"
         },
         "dSPN_3": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c10_D1-mWT-P270-20-v20211026/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20190503/21-6-DE-cor-rep-ax.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20190503/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20190503/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/dspn/str-dspn-e150917_c6_D1-m21-6-DE-v20190503/modulation.json",
             "num": 11875,
             "hoc": null,
             "neuronType": "neuron",
@@ -488,9 +486,10 @@
             "volumeID": "Striatum"
         },
         "iSPN_0": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20211026/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20190603/WT-MSN1-cor-rep-ax.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20190603/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20190603/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20190603/modulation.json",
             "num": 11875,
             "hoc": null,
             "neuronType": "neuron",
@@ -498,9 +497,10 @@
             "volumeID": "Striatum"
         },
         "iSPN_1": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20211026/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20190529/46-3-DE-cor-rep-ax.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20190529/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20190529/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e160118_c10_D2-m46-3-DE-v20190529/modulation.json",
             "num": 11875,
             "hoc": null,
             "neuronType": "neuron",
@@ -508,9 +508,10 @@
             "volumeID": "Striatum"
         },
         "iSPN_2": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150917_c11_D2-mWT-MSN1-v20211026/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20190527/WT-P270-09-15ak-cor.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20190527/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20190527/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e151123_c1_D2-mWT-P270-09-v20190527/modulation.json",
             "num": 11875,
             "hoc": null,
             "neuronType": "neuron",
@@ -518,9 +519,10 @@
             "volumeID": "Striatum"
         },
         "iSPN_3": {
-            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20211026/morphology",
-            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20211026/parameters.json",
-            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20211026/mechanisms.json",
+            "morphology": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20190611/51-5-DE-cor-rep-ax.swc",
+            "parameters": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20190611/parameters.json",
+            "mechanisms": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20190611/mechanisms.json",
+            "modulation": "$SNUDDA_DATA/neurons/striatum/ispn/str-ispn-e150908_c4_D2-m51-5-DE-v20190611/modulation.json",
             "num": 11875,
             "hoc": null,
             "neuronType": "neuron",
diff --git a/tests/networks/network_testing_input/input-test-1.json b/tests/networks/network_testing_input/input-test-1.json
index b945a082c..70d809281 100644
--- a/tests/networks/network_testing_input/input-test-1.json
+++ b/tests/networks/network_testing_input/input-test-1.json
@@ -19,6 +19,7 @@
 	    "synapseDensity" : "0.45*0.05/(1+exp(-(d-30e-6)/5e-6))",
 	    "nInputs": 20,
 	    "clusterSize" : 3,
+	    "clusterSpread" : 30e-6,
 	    "frequency" : [1, 2, 3],
 	    "populationUnitCorrelation" : 0.0,
 	    "jitter" : 0.002,
diff --git a/tests/networks/network_testing_place/network-config.json b/tests/networks/network_testing_place/network-config.json
index 0bf5ff352..972b5ef38 100644
--- a/tests/networks/network_testing_place/network-config.json
+++ b/tests/networks/network_testing_place/network-config.json
@@ -155,7 +155,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 0.38,
+                    "f1": 0.28500000000000003,
                     "softMax": 3,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -163,7 +163,7 @@
                     "cluster": false
                 },
                 "pruningOther": {
-                    "f1": 0.38,
+                    "f1": 0.28500000000000003,
                     "softMax": 3,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -192,7 +192,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 0.2,
+                    "f1": 0.164,
                     "softMax": 3,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -200,7 +200,7 @@
                     "cluster": false
                 },
                 "pruningOther": {
-                    "f1": 0.2,
+                    "f1": 0.164,
                     "softMax": 3,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -249,7 +249,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 0.3,
+                    "f1": 0.279,
                     "softMax": 4,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -257,7 +257,7 @@
                     "cluster": false
                 },
                 "pruningOther": {
-                    "f1": 0.3,
+                    "f1": 0.279,
                     "softMax": 4,
                     "mu2": 2.4,
                     "a3": 1.0,
@@ -341,7 +341,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 1.0,
+                    "f1": 0.3,
                     "softMax": 15,
                     "mu2": 3,
                     "a3": 0.3,
@@ -368,7 +368,7 @@
                     "modFile": "tmGabaA"
                 },
                 "pruning": {
-                    "f1": 1.0,
+                    "f1": 0.3,
                     "softMax": 15,
                     "mu2": 3,
                     "a3": 0.3,
diff --git a/tests/test_detect.py b/tests/test_detect.py
index 5df2bf9d3..5823d9649 100644
--- a/tests/test_detect.py
+++ b/tests/test_detect.py
@@ -29,9 +29,7 @@ def setUp(self):
 
         #  TODO: If d_view is None code run sin serial, add test parallel
         sp = SnuddaPlace(config_file=config_file, d_view=None, verbose=True)
-
-        sp.parse_config()
-        sp.write_data(position_file)
+        sp.place()
 
         # We want to load in the ball and stick neuron that has 20 micrometer soma diameter, and axon (along y-axis),
         # and dendrite along (x-axis) out to 100 micrometer distance from centre of soma.
@@ -117,16 +115,16 @@ def test_detect(self):
 
             self.assertEqual(self.sd.hyper_voxel_synapse_ctr, 101)
 
-            for pre_id in range(0, 10):
-                for post_id in range(10, 20):
+            for post_id in range(0, 10):
+                for pre_id in range(10, 20):
                     self.assertEqual(self.check_neuron_pair_has_synapse(pre_id, post_id), 1)
 
-            for pre_id in range(0, 10):
-                for post_id in range(0, 10):
+            for post_id in range(0, 10):
+                for pre_id in range(0, 10):
                     self.assertFalse(self.check_neuron_pair_has_synapse(pre_id, post_id))
 
-            for pre_id in range(10, 20):
-                for post_id in range(10, 20):
+            for post_id in range(10, 20):
+                for pre_id in range(10, 20):
                     self.assertFalse(self.check_neuron_pair_has_synapse(pre_id, post_id))
 
         with self.subTest(stage="synapse_sorting_check"):
@@ -184,9 +182,9 @@ def check_neuron_pair_has_synapse(self, pre_neuron, post_neuron):
 
         connections = dict()
 
-        for synapse_row in self.sd.hyper_voxel_synapses[0:self.sd.hyper_voxel_synapse_ctr,:]:
+        for synapse_row in self.sd.hyper_voxel_synapses[0:self.sd.hyper_voxel_synapse_ctr, :]:
 
-            loc = (synapse_row[1], synapse_row[0])
+            loc = (synapse_row[0], synapse_row[1])
             if loc in connections:
                 connections[loc] += 1
             else:
@@ -283,8 +281,8 @@ def test_detect_lines(self):
             self.assertTrue(self.sd.hyper_voxel_synapse_ctr >= 10)
 
             # Verify that all pairs are connected
-            for pre_id in range(0, 10):
-                post_id = pre_id + 10
+            for post_id in range(0, 10):
+                pre_id = post_id + 10
                 self.assertTrue(1 <= self.check_neuron_pair_has_synapse(pre_id, post_id) <= 2)
 
 
diff --git a/tests/test_detect_cluster.py b/tests/test_detect_cluster.py
new file mode 100644
index 000000000..1f1695e7d
--- /dev/null
+++ b/tests/test_detect_cluster.py
@@ -0,0 +1,149 @@
+import unittest
+import os
+import sys
+import numpy as np
+
+from snudda.place.create_cube_mesh import create_cube_mesh
+from snudda.detect.detect import SnuddaDetect
+from snudda.place.place import SnuddaPlace
+
+from scipy.spatial.distance import pdist
+
+sys.path.insert(0, os.path.join(os.path.dirname(__file__), ".."))
+
+
+class TestDetectSynapseCluster(unittest.TestCase):
+
+    def setUp(self):
+
+        if os.path.dirname(__file__):
+            os.chdir(os.path.dirname(__file__))
+
+        network_path = os.path.join(os.path.dirname(__file__), "networks", "network_testing_detect_cluster")
+
+        create_cube_mesh(file_name=os.path.join(network_path, "mesh", "simple_mesh.obj"),
+                         centre_point=(0, 0, 0),
+                         side_len=500e-6)
+
+        config_file = os.path.join(network_path, "network-config.json")
+        position_file = os.path.join(network_path, "network-neuron-positions.hdf5")
+        save_file = os.path.join(network_path, "voxels", "network-putative-synapses.hdf5")
+
+        #  TODO: If d_view is None code run sin serial, add test parallel
+        sp = SnuddaPlace(config_file=config_file, d_view=None, verbose=True)
+        sp.place()
+
+        # We want to load in the ball and stick neuron that has 20 micrometer soma diameter, and axon (along y-axis),
+        # and dendrite along (x-axis) out to 100 micrometer distance from centre of soma.
+
+        self.sd = SnuddaDetect(config_file=config_file, position_file=position_file,
+                               save_file=save_file, rc=None,
+                               hyper_voxel_size=150, verbose=True)
+
+
+    def test_clustering(self):
+
+        neuron_positions = np.array([[0, 20, 0],    # Postsynaptiska
+                                     [0, 80, 0],
+                                     [0, 140, 0],
+                                     [0, 200, 0],
+                                     [20, 0, 0],    # Presynaptiska
+                                     [80, 0, 0],
+                                     [140, 0, 0],
+                                     [200, 0, 0]
+                                     ])*1e-6
+
+        for idx, pos in enumerate(neuron_positions):
+            self.sd.neurons[idx]["position"] = pos
+
+        ang = -np.pi/2
+        R_x = np.array([[1, 0, 0],
+                        [0, np.cos(ang), -np.sin(ang)],
+                        [0, np.sin(ang), np.cos(ang)]])
+
+        ang = np.pi/2
+        R_y = np.array([[np.cos(ang), 0, np.sin(ang)],
+                        [0, 1, 0],
+                        [-np.sin(ang), 0, np.cos(ang)]])
+
+        for idx in range(0, 4):    # Post-synaptic neurons
+            self.sd.neurons[idx]["rotation"] = R_x
+
+        for idx in range(4, 8):   # Pre-synaptic neurons
+            self.sd.neurons[idx]["rotation"] = R_y
+
+        self.sd.detect(restart_detection_flag=True)
+
+        synapse_voxel_loc = self.sd.hyper_voxel_synapses[:self.sd.hyper_voxel_synapse_ctr, 2:5]
+        synapse_coords = synapse_voxel_loc * self.sd.voxel_size + self.sd.hyper_voxel_origo
+
+        fig_path = os.path.join("networks", "network_testing_detect_cluster", "figures")
+
+        if not os.path.exists(fig_path):
+            os.mkdir(fig_path)
+
+        if False:   # Set to True to include plot
+            self.sd.plot_hyper_voxel(plot_neurons=True, fig_file_name="touch-detection-clusters-validation")
+
+        # Check that each neuron pair
+        with self.subTest(stage="cluster-size"):
+
+            connection_list = [(4, 0, 5), (4, 1, 5), (4, 2, 5), (4, 3, 5),
+                               (5, 0, 10), (5, 1, 10), (5, 2, 10), (5, 3, 10),
+                               (6, 0, 10), (6, 1, 10), (6, 2, 10), (6, 3, 10),
+                               (7, 0, 10), (7, 1, 10), (7, 2, 10), (7, 3, 10)]
+
+            for pre_id, post_id, n_synapses in connection_list:
+                self.assertTrue(self.check_neuron_pair_has_synapse(pre_id, post_id) == n_synapses)
+
+        with self.subTest(stage="cluster-spread"):
+
+            self.assertTrue(0 <= self.get_synapse_spread(4, 0) <= 25e-6)
+            self.assertTrue(0 <= self.get_synapse_spread(4, 1) <= 25e-6)
+            self.assertTrue(0 <= self.get_synapse_spread(4, 2) <= 25e-6)
+            self.assertTrue(0 <= self.get_synapse_spread(4, 3) <= 25e-6)
+
+        # Lägg till tester som kollar positionerna på synapserna
+        # import pdb
+        # pdb.set_trace()
+
+    def check_neuron_pair_has_synapse(self, pre_neuron, post_neuron):
+
+        connections = dict()
+
+        for synapse_row in self.sd.hyper_voxel_synapses[0:self.sd.hyper_voxel_synapse_ctr, :]:
+
+            loc = (synapse_row[0], synapse_row[1])
+            if loc in connections:
+                connections[loc] += 1
+            else:
+                connections[loc] = 1
+
+        if (pre_neuron, post_neuron) in connections:
+            # print(f"pre: {pre_neuron}, post: {post_neuron}, connections: {connections[(pre_neuron, post_neuron)]}")
+            return connections[(pre_neuron, post_neuron)]
+        else:
+            # print(f"No connection between pre {pre_neuron} and post {post_neuron}")
+            return False
+
+    def get_synapse_spread(self, pre_neuron, post_neuron):
+
+        coords = self.get_synapse_coordinates(pre_neuron, post_neuron)
+        return np.max(pdist(coords))
+
+    def get_synapse_coordinates(self, pre_neuron, post_neuron):
+
+        coord_list = []
+
+        for synapse_row in self.sd.hyper_voxel_synapses[0:self.sd.hyper_voxel_synapse_ctr, :]:
+            if synapse_row[0] == pre_neuron and synapse_row[1] == post_neuron:
+                coord_list.append(self.convert_to_coordinates(synapse_row[2:5]))
+
+        return np.vstack(coord_list)
+
+    def convert_to_coordinates(self, voxel_xyz):
+        return np.array(voxel_xyz) * self.sd.voxel_size + self.sd.hyper_voxel_origo
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/tests/test_input.py b/tests/test_input.py
index b49054386..0519a311f 100644
--- a/tests/test_input.py
+++ b/tests/test_input.py
@@ -94,6 +94,11 @@ def test_input_1(self):
                     print(f"Checking number of inputs is {config_data[neuron_type][input_type]['nInputs']} * {cluster_size}")
                     self.assertEqual(config_data[neuron_type][input_type]["nInputs"]*cluster_size, n_traces)
 
+                    if cluster_size > 1:
+                        # Verify that all the clusters have the right size
+                        for ctr in range(0, cluster_size-1):
+                            self.assertTrue(np.all(np.diff(input_info["sectionID"])[ctr::cluster_size] == 0))
+
                 max_len = 1
                 if type(start_time) is np.ndarray:
                     max_len = np.maximum(max_len, len(start_time))
diff --git a/tests/test_neuron_morphology.py b/tests/test_neuron_morphology.py
index d16a33de7..43687b323 100644
--- a/tests/test_neuron_morphology.py
+++ b/tests/test_neuron_morphology.py
@@ -1,6 +1,7 @@
 import os
 import unittest
 import numpy as np
+import scipy
 
 from snudda.neurons.neuron_morphology import NeuronMorphology
 
@@ -13,7 +14,7 @@ def setUp(self):
 
         swc_file = os.path.join(os.path.dirname(__file__), "validation", "striatum", "fs",
                                 "str-fs-e161205_FS1-mMTC180800A-IDB-v20190312", "MTC180800A-IDB-cor-rep.swc")
-        self.nm = NeuronMorphology(swc_filename=swc_file, load_morphology=True)
+        self.nm = NeuronMorphology(swc_filename=swc_file, load_morphology=True, use_cache=False)
 
     def test_input_location(self, stage="dist_to_soma"):
 
@@ -71,7 +72,41 @@ def test_clone(self, stage="clone"):
         self.assertTrue((np.abs(self.nm.axon[:, 2] + new_nm.axon[:, 2]) < 1e-6).all())
 
         new_nm.place(position=np.array([1, 2, 3]))
-        self.assertTrue((np.abs(new_nm.soma[0,:3] - np.array([1, 2, 3])) < 1e-6).all())
+        self.assertTrue((np.abs(new_nm.soma[0, :3] - np.array([1, 2, 3])) < 1e-6).all())
+
+    def test_cluster_synapses(self, stage="cluster_synapses"):
+
+        cluster_spread = 30e-6
+        n_synapses = 10
+
+        cluster_sec_x, syn_coords, soma_dist = self.nm.cluster_synapses(sec_id=30, sec_x=0.5,
+                                                                        count=n_synapses, distance=cluster_spread,
+                                                                        rng=np.random.default_rng(20220125))
+
+        self.assertEqual(len(cluster_sec_x), n_synapses)
+        self.assertTrue(np.max(soma_dist) - np.min(soma_dist) <= cluster_spread)
+
+        d = scipy.spatial.distance.pdist(syn_coords)
+        self.assertTrue(np.all(d <= cluster_spread))
+        self.assertTrue(np.any(d >= 10e-6))
+
+        # TODO: Verify self.nm.sec_id_to_len[0] value ... ie segment length of soma?
+        # xx = self.nm.compartment_length()
+        # np.sum(self.nm.sec_id_to_len[1:])
+        # ska vara samma som np.sum(xx)
+        # och self.nm.sec_id_to_len[0] deals with the soma.
+
+        if False:
+            import matplotlib.pyplot as plt
+            fig = plt.figure()
+            ax = fig.add_subplot(111, projection='3d')
+            ax = self.nm.plot_neuron(plot_axon=False, plot_dendrite=True, axis=ax)
+            ax.scatter(syn_coords[:, 0], syn_coords[:, 1], syn_coords[:, 2], s=5, c="red")
+            fig.show()
+            plt.pause(3)
+
+            import pdb
+            pdb.set_trace()
 
 
 if __name__ == '__main__':
diff --git a/tests/test_projection_detection.py b/tests/test_projection_detection.py
new file mode 100644
index 000000000..a9df1c079
--- /dev/null
+++ b/tests/test_projection_detection.py
@@ -0,0 +1,205 @@
+import os
+import time
+import unittest
+
+
+class TestProjectionDetection(unittest.TestCase):
+
+    def setUp(self):
+        from snudda.place.create_cube_mesh import create_cube_mesh
+
+        # Create cube meshes
+        self.network_path = os.path.join("networks", "network_testing_projection_detection")
+        mesh_file_a = os.path.join(self.network_path, "mesh", "volume_A.obj")
+        mesh_file_b = os.path.join(self.network_path, "mesh", "volume_B.obj")
+
+        create_cube_mesh(mesh_file_a, [5e-3, 0, 0], 300e-6, "Volume A - connect structures example")
+        create_cube_mesh(mesh_file_b, [-5e-3, 0, 0], 300e-6, "Volume B - connect structures example")
+
+        # Define network
+        from snudda.init.init import SnuddaInit
+
+        cnc = SnuddaInit(network_path=self.network_path, random_seed=123)
+        cnc.define_structure(struct_name="VolumeA", struct_mesh=mesh_file_a, d_min=15e-6, mesh_bin_width=50e-6)
+        cnc.define_structure(struct_name="VolumeB", struct_mesh=mesh_file_b, d_min=15e-6, mesh_bin_width=50e-6)
+
+        cnc.add_neurons(name="dSPN", num_neurons=20, volume_id="VolumeA",
+                        neuron_dir=os.path.join("$DATA", "neurons", "striatum", "dspn"))
+        cnc.add_neurons(name="iSPN", num_neurons=20, volume_id="VolumeB",
+                        neuron_dir=os.path.join("$DATA", "neurons", "striatum", "ispn"))
+
+        # Add the projection we want to test dSPN->iSPN
+        proj_file = os.path.join("data", "ExampleProjectionDetection.json")
+
+        cnc.add_neuron_target(neuron_name="dSPN",
+                              target_name="iSPN",
+                              connection_type="GABA",
+                              dist_pruning=None,
+                              f1=None,
+                              soft_max=None,
+                              mu2=None,
+                              a3=None,
+                              conductance=5e-9,
+                              cluster_synapses=False,
+                              mod_file="tmGlutA",
+                              parameter_file=None,
+                              channel_param_dictionary=None,
+                              projection_file=proj_file)
+
+        # Also add dSPN-dSPN and iSPN-iSPN synapses
+        # Note we do NOT add dSPN-iSPN again this way, as that would overwrite the above connections
+        # (The above neuron_projection will also do normal touch detection)
+
+        SPN2SPNdistDepPruning = "1-exp(-(0.4*d/60e-6)**2)"
+
+        MSD1gGABA = [0.24e-9, 0.1e-9]
+        MSD2gGABA = [0.24e-9, 0.1e-9]
+
+        MSD1GABAfailRate = 0.7  # Taverna 2008, figure 2
+        MSD2GABAfailRate = 0.4  # Taverna 2008, 2mM
+
+        pfdSPNdSPN = os.path.join("$DATA", "synapses", "striatum", "PlanertFitting-DD-tmgaba-fit.json")
+        pfdSPNiSPN = os.path.join("$DATA", "synapses", "striatum", "PlanertFitting-DI-tmgaba-fit.json")
+        pfiSPNdSPN = os.path.join("$DATA", "synapses", "striatum", "PlanertFitting-ID-tmgaba-fit.json")
+        pfiSPNiSPN = os.path.join("$DATA", "synapses", "striatum", "PlanertFitting-II-tmgaba-fit.json")
+
+        cnc.add_neuron_target(neuron_name="dSPN",
+                              target_name="dSPN",
+                              connection_type="GABA",
+                              dist_pruning=SPN2SPNdistDepPruning,
+                              f1=0.38, soft_max=3, mu2=2.4,
+                              a3=1.0,
+                              conductance=MSD1gGABA,
+                              parameter_file=pfdSPNdSPN,
+                              mod_file="tmGabaA",
+                              channel_param_dictionary={"tau1": (1.3e-3, 1e3),
+                                                        "tau2": (12.4e-3, 1e3),
+                                                        "failRate": MSD1GABAfailRate})
+
+        cnc.add_neuron_target(neuron_name="iSPN",
+                              target_name="iSPN",
+                              connection_type="GABA",
+                              dist_pruning=SPN2SPNdistDepPruning,
+                              f1=0.55, soft_max=4, mu2=2.4,
+                              a3=1.0,
+                              conductance=MSD2gGABA,
+                              parameter_file=pfiSPNiSPN,
+                              mod_file="tmGabaA",
+                              channel_param_dictionary={"tau1": (1.3e-3, 1e3),
+                                                        "tau2": (12.4e-3, 1e3),
+                                                        "failRate": MSD2GABAfailRate})
+
+        cnc.write_json()
+
+        # Place neurons, then detect, project and prune
+
+        from snudda.place.place import SnuddaPlace
+        sp = SnuddaPlace(network_path=self.network_path, verbose=True)
+        sp.parse_config()
+        sp.write_data()
+
+        from snudda.detect.detect import SnuddaDetect
+        sd = SnuddaDetect(network_path=self.network_path, hyper_voxel_size=100, verbose=True)
+        sd.detect()
+
+        # Old version of projection
+        from snudda.detect.project import SnuddaProject
+        sp = SnuddaProject(network_path=self.network_path)
+        sp.project()
+
+        from snudda.detect.prune import SnuddaPrune
+        sp = SnuddaPrune(network_path=self.network_path, verbose=True)
+        sp.prune()
+
+    def test_project(self):
+
+        # Are there connections dSPN->iSPN
+        from snudda.utils.load import SnuddaLoad
+        network_file = os.path.join(self.network_path, "network-synapses.hdf5")
+        sl = SnuddaLoad(network_file)
+
+        dspn_id_list = sl.get_neuron_id_of_type("dSPN")
+        ispn_id_list = sl.get_neuron_id_of_type("iSPN")
+
+        tot_proj_ctr = 0
+
+        for dspn_id in dspn_id_list:
+            for ispn_id in ispn_id_list:
+
+                synapses, synapse_coords = sl.find_synapses(pre_id=dspn_id, post_id=ispn_id)
+                if synapses is not None:
+                    tot_proj_ctr += synapses.shape[0]
+
+        with self.subTest(stage="projection_exists"):
+            # There should be projection synapses between dSPN and iSPN in this toy example
+            self.assertTrue(tot_proj_ctr > 0)
+
+        tot_dd_syn_ctr = 0
+        for dspn_id in dspn_id_list:
+            for dspn_id2 in dspn_id_list:
+
+                synapses, synapse_coords = sl.find_synapses(pre_id=dspn_id, post_id=dspn_id2)
+                if synapses is not None:
+                    tot_dd_syn_ctr += synapses.shape[0]
+
+        tot_ii_syn_ctr = 0
+        for ispn_id in ispn_id_list:
+            for ispn_id2 in ispn_id_list:
+
+                synapses, synapse_coords = sl.find_synapses(pre_id=ispn_id, post_id=ispn_id2)
+                if synapses is not None:
+                    tot_ii_syn_ctr += synapses.shape[0]
+
+        with self.subTest(stage="normal_synapses_exist"):
+            # In this toy example neurons are quite sparsely placed, but we should have at least some
+            # synapses
+            self.assertTrue(tot_dd_syn_ctr > 0)
+            self.assertTrue(tot_ii_syn_ctr > 0)
+
+        # We need to run in parallel also to verify we get same result (same random seed)
+
+        serial_synapses = sl.data["synapses"].copy()
+        del sl  # Close old file so we can overwrite it
+
+        os.environ["IPYTHONDIR"] = os.path.join(os.path.abspath(os.getcwd()), ".ipython")
+        os.environ["IPYTHON_PROFILE"] = "default"
+        os.system("ipcluster start -n 4 --profile=$IPYTHON_PROFILE --ip=127.0.0.1&")
+        time.sleep(10)
+
+        # Run place, detect and prune in parallel by passing rc
+        from ipyparallel import Client
+        u_file = os.path.join(".ipython", "profile_default", "security", "ipcontroller-client.json")
+        rc = Client(url_file=u_file, timeout=120, debug=False)
+        d_view = rc.direct_view(targets='all')  # rc[:] # Direct view into clients
+
+        from snudda.detect.detect import SnuddaDetect
+        sd = SnuddaDetect(network_path=self.network_path, hyper_voxel_size=100, rc=rc, verbose=True)
+        sd.detect()
+
+        # Old version of project, the new version with touch detection is in detect.
+        from snudda.detect.project import SnuddaProject
+        # TODO: Currently SnuddaProject only runs in serial
+        sp = SnuddaProject(network_path=self.network_path)
+        sp.project()
+
+        from snudda.detect.prune import SnuddaPrune
+        # Prune has different methods for serial and parallel execution, important to test it!
+        sp = SnuddaPrune(network_path=self.network_path, rc=rc, verbose=True)
+        sp.prune()
+
+        with self.subTest(stage="check-parallel-identical"):
+            sl2 = SnuddaLoad(network_file)
+            parallel_synapses = sl2.data["synapses"].copy()
+
+            # ParameterID, sec_X etc are randomised in hyper voxel, so you need to use same
+            # hypervoxel size for reproducability between serial and parallel execution
+
+            # All synapses should be identical regardless of serial or parallel execution path
+            self.assertTrue(serial_synapses.shape == parallel_synapses.shape)
+            self.assertTrue((serial_synapses == parallel_synapses).all())
+
+        os.system("ipcluster stop")
+
+
+if __name__ == '__main__':
+    unittest.main()