From dc371e5d15fa6acc4b778c12df1f5c6e4adaadc2 Mon Sep 17 00:00:00 2001
From: Nora Abi Akar <nora.abiakar@gmail.com>
Date: Fri, 22 Oct 2021 11:23:05 +0200
Subject: [PATCH] Gap Junction mechanisms (#1682)

Adds support for user-defined gap junction mechanisms. Fixes #1600.

- API Change: Wrap `mechanism_desc` in one of `density`, `synapse`, `junction` to be able to differentiate the kind of placed mechanism between point and junction mechanisms.
- API Change: The `ggap` parameter of a `gap_junction_connection` renamed to `weight` and is not a unit-less parameter.
- Mechanism ABI change: New `peer_index` vector which stores the indices of the peer sites of a gap-junction-connection. New mechanism kind.
- Internal API change: `shared_state` no longer controls the gap junction current contribution, that is now handled in the mechanism callbacks.
- Modcc change: Add `JUNCTION_PROCESS` and `v_peer` keywords to nmodl, the first to identify a gap-junction mechanism, the second to access the peer voltage of a gap-junction connection. Modify modcc to be able to compile nmodl and generate correct code.
- Behavior change: gap-junctions now contribute to both the current and conductance of a CV as opposed to just the current.
- Add new "gj" mechanism, to replace the built-in constant-conductance based linear gap-junction implementation.
- Fix unit tests.
- Fix python wrapper.
- Fix docs.
---
 arbor/backends/gpu/gpu_store_types.hpp        |   1 -
 arbor/backends/gpu/shared_state.cpp           |  17 +-
 arbor/backends/gpu/shared_state.cu            |  24 -
 arbor/backends/gpu/shared_state.hpp           |   6 -
 arbor/backends/multicore/multicore_common.hpp |   1 -
 arbor/backends/multicore/shared_state.cpp     |  24 +-
 arbor/backends/multicore/shared_state.hpp     |   6 -
 arbor/cable_cell.cpp                          |  18 +-
 arbor/fvm_layout.cpp                          | 178 +++-
 arbor/fvm_layout.hpp                          |  27 +-
 arbor/fvm_lowered_cell_impl.hpp               |  76 +-
 arbor/include/arbor/cable_cell.hpp            |  20 +-
 arbor/include/arbor/cable_cell_param.hpp      |  42 +-
 arbor/include/arbor/fvm_types.hpp             |  13 +-
 arbor/include/arbor/mechanism.hpp             |   3 +
 arbor/include/arbor/mechanism_abi.h           |   2 +
 arbor/include/arbor/mechinfo.hpp              |   1 -
 arbor/include/arbor/recipe.hpp                |   5 +-
 arbor/mechinfo.cpp                            |   2 +-
 arborio/cableio.cpp                           |  35 +-
 doc/concepts/cable_cell.rst                   |   4 +-
 doc/concepts/decor.rst                        |  75 +-
 doc/concepts/interconnectivity.rst            |   9 +-
 doc/concepts/mechanisms.rst                   |  41 +-
 doc/concepts/recipe.rst                       |   4 +-
 doc/cpp/cable_cell.rst                        |  57 +-
 doc/cpp/interconnectivity.rst                 |  13 +-
 doc/cpp/morphology.rst                        |   2 +-
 doc/fileformat/cable_cell.rst                 |  36 +-
 doc/fileformat/nmodl.rst                      |   5 +
 doc/internals/mechanism_abi.rst               |   5 +
 doc/python/cable_cell.rst                     |   6 +-
 doc/python/decor.rst                          |  56 +-
 doc/python/interconnectivity.rst              |  11 +-
 doc/python/mechanisms.rst                     | 205 +++-
 doc/tutorial/single_cell_detailed.rst         |   6 +-
 doc/tutorial/single_cell_model.rst            |   4 +-
 example/dryrun/branch_cell.hpp                |   8 +-
 example/gap_junctions/gap_junctions.cpp       |  14 +-
 example/generators/generators.cpp             |   4 +-
 example/lfp/lfp.cpp                           |   6 +-
 example/probe-demo/probe-demo.cpp             |   2 +-
 example/ring/branch_cell.hpp                  |   8 +-
 example/single/single.cpp                     |   6 +-
 mechanisms/CMakeLists.txt                     |   2 +-
 mechanisms/default/gj.mod                     |  14 +
 modcc/identifier.hpp                          |   6 +-
 modcc/module.cpp                              |  17 +-
 modcc/parser.cpp                              |   5 +-
 modcc/printer/cprinter.cpp                    |   1 +
 modcc/printer/gpuprinter.cpp                  |   3 +-
 modcc/printer/printerutil.cpp                 |   5 +
 modcc/printer/printerutil.hpp                 |   7 +-
 modcc/token.cpp                               |   6 +-
 modcc/token.hpp                               |   2 +-
 python/cells.cpp                              |  81 +-
 python/example/network_ring.py                |   6 +-
 python/example/network_ring_mpi.py            |   6 +-
 python/example/single_cell_cable.py           |   2 +-
 python/example/single_cell_detailed.py        |   8 +-
 python/example/single_cell_detailed_recipe.py |   8 +-
 .../single_cell_extracellular_potentials.py   |   2 +-
 python/example/single_cell_model.py           |   2 +-
 python/example/single_cell_nml.py             |   8 +-
 python/example/single_cell_recipe.py          |   2 +-
 python/example/single_cell_stdp.py            |  11 +-
 python/example/single_cell_swc.py             |   6 +-
 python/recipe.cpp                             |  12 +-
 python/test/fixtures.py                       |   2 +-
 python/test/unit/test_cable_probes.py         |   6 +-
 python/test/unit/test_catalogues.py           |   2 +-
 test/common_cells.cpp                         |  10 +-
 test/unit-distributed/test_communicator.cpp   |   8 +-
 test/unit/CMakeLists.txt                      |   2 +
 test/unit/mod/gj0.mod                         |  14 +
 test/unit/mod/gj1.mod                         |  15 +
 test/unit/test_abi.cpp                        |  12 +-
 test/unit/test_cable_cell.cpp                 |  10 +-
 test/unit/test_domain_decomposition.cpp       |   2 +-
 test/unit/test_event_delivery.cpp             |   6 +-
 test/unit/test_fvm_layout.cpp                 | 945 ++++++++++++++++--
 test/unit/test_fvm_lowered.cpp                | 377 +------
 test/unit/test_kinetic_linear.cpp             |   3 +-
 test/unit/test_mc_cell_group.cpp              |   4 +-
 test/unit/test_mc_cell_group_gpu.cpp          |   4 +-
 test/unit/test_mech_temp_diam.cpp             |   6 +-
 test/unit/test_mechcat.cpp                    |   2 +-
 test/unit/test_probe.cpp                      |  44 +-
 test/unit/test_recipe.cpp                     |   4 +-
 test/unit/test_s_expr.cpp                     |  84 +-
 test/unit/test_spikes.cpp                     |   2 +-
 test/unit/test_synapses.cpp                   |  20 +-
 test/unit/unit_test_catalogue.cpp             |   4 +
 93 files changed, 1876 insertions(+), 1022 deletions(-)
 create mode 100644 mechanisms/default/gj.mod
 create mode 100644 test/unit/mod/gj0.mod
 create mode 100644 test/unit/mod/gj1.mod

diff --git a/arbor/backends/gpu/gpu_store_types.hpp b/arbor/backends/gpu/gpu_store_types.hpp
index 708233a9..ac215ed4 100644
--- a/arbor/backends/gpu/gpu_store_types.hpp
+++ b/arbor/backends/gpu/gpu_store_types.hpp
@@ -17,7 +17,6 @@ namespace gpu {
 
 using array  = memory::device_vector<fvm_value_type>;
 using iarray = memory::device_vector<fvm_index_type>;
-using gjarray = memory::device_vector<fvm_gap_junction>;
 
 using deliverable_event_stream = arb::gpu::multi_event_stream<deliverable_event>;
 using sample_event_stream = arb::gpu::multi_event_stream<sample_event>;
diff --git a/arbor/backends/gpu/shared_state.cpp b/arbor/backends/gpu/shared_state.cpp
index 566c2822..44acd5e4 100644
--- a/arbor/backends/gpu/shared_state.cpp
+++ b/arbor/backends/gpu/shared_state.cpp
@@ -38,9 +38,6 @@ void set_dt_impl(
     fvm_size_type nintdom, fvm_size_type ncomp, fvm_value_type* dt_intdom, fvm_value_type* dt_comp,
     const fvm_value_type* time_to, const fvm_value_type* time, const fvm_index_type* cv_to_intdom);
 
-void add_gj_current_impl(
-    fvm_size_type n_gj, const fvm_gap_junction* gj, const fvm_value_type* v, fvm_value_type* i);
-
 void take_samples_impl(
     const multi_event_stream_state<raw_probe_info>& s,
     const fvm_value_type* time, fvm_value_type* sample_time, fvm_value_type* sample_value);
@@ -179,7 +176,6 @@ shared_state::shared_state(
     fvm_size_type n_detector,
     const std::vector<fvm_index_type>& cv_to_intdom_vec,
     const std::vector<fvm_index_type>& cv_to_cell_vec,
-    const std::vector<fvm_gap_junction>& gj_vec,
     const std::vector<fvm_value_type>& init_membrane_potential,
     const std::vector<fvm_value_type>& temperature_K,
     const std::vector<fvm_value_type>& diam,
@@ -189,10 +185,8 @@ shared_state::shared_state(
     n_intdom(n_intdom),
     n_detector(n_detector),
     n_cv(cv_to_intdom_vec.size()),
-    n_gj(gj_vec.size()),
     cv_to_intdom(make_const_view(cv_to_intdom_vec)),
     cv_to_cell(make_const_view(cv_to_cell_vec)),
-    gap_junctions(make_const_view(gj_vec)),
     time(n_intdom),
     time_to(n_intdom),
     dt_intdom(n_intdom),
@@ -280,6 +274,7 @@ void shared_state::instantiate(mechanism& m, unsigned id, const mechanism_overri
     using util::value_by_key;
 
     bool mult_in_place = !pos_data.multiplicity.empty();
+    bool peer_indices = !pos_data.peer_cv.empty();
 
     // Set internal variables
     m.time_ptr_ptr   = &time_ptr;
@@ -358,7 +353,7 @@ void shared_state::instantiate(mechanism& m, unsigned id, const mechanism_overri
     // Allocate and initialize index vectors, viz. node_index_ and any ion indices.
     {
         // Allocate bulk storage
-        std::size_t count = mult_in_place + m.mech_.n_ions + 1;
+        std::size_t count = mult_in_place + peer_indices + m.mech_.n_ions + 1;
         store.indices_ = iarray(count*width_padded);
         chunk_writer writer(store.indices_.data(), width);
 
@@ -379,6 +374,10 @@ void shared_state::instantiate(mechanism& m, unsigned id, const mechanism_overri
         }
 
         m.ppack_.multiplicity = mult_in_place? writer.append(pos_data.multiplicity): nullptr;
+        // `peer_index` holds the peer CV of each CV in node_index.
+        // Peer CVs are only filled for gap junction mechanisms. They are used
+        // to index the voltage at the other side of a gap-junction connection.
+        m.ppack_.peer_index = peer_indices? writer.append(pos_data.peer_cv): nullptr;
     }
 
     // Shift data to GPU, set up pointers
@@ -443,10 +442,6 @@ void shared_state::set_dt() {
     set_dt_impl(n_intdom, n_cv, dt_intdom.data(), dt_cv.data(), time_to.data(), time.data(), cv_to_intdom.data());
 }
 
-void shared_state::add_gj_current() {
-    add_gj_current_impl(n_gj, gap_junctions.data(), voltage.data(), current_density.data());
-}
-
 void shared_state::add_stimulus_current() {
     stim_data.add_current(time, cv_to_intdom, current_density);
 }
diff --git a/arbor/backends/gpu/shared_state.cu b/arbor/backends/gpu/shared_state.cu
index 8fdd161b..5f6aba70 100644
--- a/arbor/backends/gpu/shared_state.cu
+++ b/arbor/backends/gpu/shared_state.cu
@@ -26,20 +26,6 @@ __global__ void update_time_to_impl(unsigned n,
     }
 }
 
-template <typename T, typename I>
-__global__ void add_gj_current_impl(unsigned n,
-                                    const T* __restrict__ const gj_info,
-                                    const I* __restrict__ const voltage,
-                                    I* __restrict__ const current_density) {
-    unsigned i = threadIdx.x+blockIdx.x*blockDim.x;
-    if (i<n) {
-        auto gj = gj_info[i];
-        auto curr = gj.weight * (voltage[gj.loc.second] - voltage[gj.loc.first]); // nA
-
-        gpu_atomic_sub(current_density + gj.loc.first, curr);
-    }
-}
-
 // Vector/scalar addition: x[i] += v ∀i
 template <typename T>
 __global__ void add_scalar(unsigned n,
@@ -118,16 +104,6 @@ void set_dt_impl(
     kernel::set_dt_impl<<<nblock, block_dim>>>(dt_intdom, time_to, time, ncomp, dt_comp, cv_to_intdom);
 }
 
-void add_gj_current_impl(
-    fvm_size_type n_gj, const fvm_gap_junction* gj_info, const fvm_value_type* voltage, fvm_value_type* current_density)
-{
-    if (!n_gj) return;
-
-    constexpr int block_dim = 128;
-    int nblock = block_count(n_gj, block_dim);
-    kernel::add_gj_current_impl<<<nblock, block_dim>>>(n_gj, gj_info, voltage, current_density);
-}
-
 void take_samples_impl(
     const multi_event_stream_state<raw_probe_info>& s,
     const fvm_value_type* time, fvm_value_type* sample_time, fvm_value_type* sample_value)
diff --git a/arbor/backends/gpu/shared_state.hpp b/arbor/backends/gpu/shared_state.hpp
index d74029f0..eb3301ab 100644
--- a/arbor/backends/gpu/shared_state.hpp
+++ b/arbor/backends/gpu/shared_state.hpp
@@ -117,11 +117,9 @@ struct shared_state {
     fvm_size_type n_intdom = 0;   // Number of distinct integration domains.
     fvm_size_type n_detector = 0; // Max number of detectors on all cells.
     fvm_size_type n_cv = 0;       // Total number of CVs.
-    fvm_size_type n_gj = 0;       // Total number of GJs.
 
     iarray cv_to_intdom;     // Maps CV index to intdom index.
     iarray cv_to_cell;       // Maps CV index to cell index.
-    gjarray gap_junctions;   // Stores gap_junction info.
     array time;              // Maps intdom index to integration start time [ms].
     array time_to;           // Maps intdom index to integration stop time [ms].
     array dt_intdom;         // Maps intdom index to (stop time) - (start time) [ms].
@@ -152,7 +150,6 @@ struct shared_state {
         fvm_size_type n_detector,
         const std::vector<fvm_index_type>& cv_to_intdom_vec,
         const std::vector<fvm_index_type>& cv_to_cell_vec,
-        const std::vector<fvm_gap_junction>& gj_vec,
         const std::vector<fvm_value_type>& init_membrane_potential,
         const std::vector<fvm_value_type>& temperature_K,
         const std::vector<fvm_value_type>& diam,
@@ -185,9 +182,6 @@ struct shared_state {
     // Set the per-intdom and per-compartment dt from time_to - time.
     void set_dt();
 
-    // Update gap_junction state
-    void add_gj_current();
-
     // Update stimulus state and add current contributions.
     void add_stimulus_current();
 
diff --git a/arbor/backends/multicore/multicore_common.hpp b/arbor/backends/multicore/multicore_common.hpp
index 4a3ee9bc..b0b3ed4e 100644
--- a/arbor/backends/multicore/multicore_common.hpp
+++ b/arbor/backends/multicore/multicore_common.hpp
@@ -23,7 +23,6 @@ using padded_vector = std::vector<V, util::padded_allocator<V>>;
 
 using array  = padded_vector<fvm_value_type>;
 using iarray = padded_vector<fvm_index_type>;
-using gjarray = padded_vector<fvm_gap_junction>;
 
 using deliverable_event_stream = arb::multicore::multi_event_stream<deliverable_event>;
 using sample_event_stream = arb::multicore::multi_event_stream<sample_event>;
diff --git a/arbor/backends/multicore/shared_state.cpp b/arbor/backends/multicore/shared_state.cpp
index 34a9853d..b282186f 100644
--- a/arbor/backends/multicore/shared_state.cpp
+++ b/arbor/backends/multicore/shared_state.cpp
@@ -191,7 +191,6 @@ shared_state::shared_state(
     fvm_size_type n_detector,
     const std::vector<fvm_index_type>& cv_to_intdom_vec,
     const std::vector<fvm_index_type>& cv_to_cell_vec,
-    const std::vector<fvm_gap_junction>& gj_vec,
     const std::vector<fvm_value_type>& init_membrane_potential,
     const std::vector<fvm_value_type>& temperature_K,
     const std::vector<fvm_value_type>& diam,
@@ -203,10 +202,8 @@ shared_state::shared_state(
     n_intdom(n_intdom),
     n_detector(n_detector),
     n_cv(cv_to_intdom_vec.size()),
-    n_gj(gj_vec.size()),
     cv_to_intdom(math::round_up(n_cv, alignment), pad(alignment)),
     cv_to_cell(math::round_up(cv_to_cell_vec.size(), alignment), pad(alignment)),
-    gap_junctions(math::round_up(n_gj, alignment), pad(alignment)),
     time(n_intdom, pad(alignment)),
     time_to(n_intdom, pad(alignment)),
     dt_intdom(n_intdom, pad(alignment)),
@@ -232,10 +229,6 @@ shared_state::shared_state(
         std::copy(cv_to_cell_vec.begin(), cv_to_cell_vec.end(), cv_to_cell.begin());
         std::fill(cv_to_cell.begin() + n_cv, cv_to_cell.end(), cv_to_cell_vec.back());
     }
-    if (n_gj>0) {
-        std::copy(gj_vec.begin(), gj_vec.end(), gap_junctions.begin());
-        std::fill(gap_junctions.begin()+n_gj, gap_junctions.end(), gj_vec.back());
-    }
 
     util::fill(time_since_spike, -1.0);
     for (unsigned i = 0; i<n_cv; ++i) {
@@ -323,16 +316,6 @@ void shared_state::set_dt() {
     }
 }
 
-void shared_state::add_gj_current() {
-    for (unsigned i = 0; i < n_gj; i++) {
-        auto gj = gap_junctions[i];
-        auto curr = gj.weight *
-                    (voltage[gj.loc.second] - voltage[gj.loc.first]); // nA
-
-        current_density[gj.loc.first] -= curr;
-    }
-}
-
 void shared_state::add_stimulus_current() {
      stim_data.add_current(time, cv_to_intdom, current_density);
 }
@@ -504,6 +487,7 @@ void shared_state::instantiate(arb::mechanism& m, unsigned id, const mechanism_o
     m.ppack_.vec_t            = nullptr;
 
     bool mult_in_place = !pos_data.multiplicity.empty();
+    bool peer_indices = !pos_data.peer_cv.empty();
 
     if (storage.find(id) != storage.end()) throw arb::arbor_internal_error("Duplicate mech id in shared state");
     auto& store = storage[id];
@@ -563,7 +547,7 @@ void shared_state::instantiate(arb::mechanism& m, unsigned id, const mechanism_o
     {
         // Allocate bulk storage
         std::size_t index_width_padded = extend_width<arb_index_type>(m, pos_data.cv.size());
-        std::size_t count = mult_in_place + m.mech_.n_ions + 1;
+        std::size_t count = mult_in_place + peer_indices + m.mech_.n_ions + 1;
         store.indices_ = iarray(count*index_width_padded, 0, pad);
         chunk_writer writer(store.indices_.data(), index_width_padded);
         // Setup node indices
@@ -597,6 +581,10 @@ void shared_state::instantiate(arb::mechanism& m, unsigned id, const mechanism_o
             arb_assert(compatible_index_constraints(node_index, util::range_n(m.ppack_.ion_states[idx].index, index_width_padded), m.iface_.partition_width));
         }
         if (mult_in_place) m.ppack_.multiplicity = writer.append(pos_data.multiplicity, 0);
+        // `peer_index` holds the peer CV of each CV in node_index.
+        // Peer CVs are only filled for gap junction mechanisms. They are used
+        // to index the voltage at the other side of a gap-junction connection.
+        if (peer_indices)  m.ppack_.peer_index   = writer.append(pos_data.peer_cv, pos_data.peer_cv.back());
     }
 }
 
diff --git a/arbor/backends/multicore/shared_state.hpp b/arbor/backends/multicore/shared_state.hpp
index d913ca9d..4d12e0f4 100644
--- a/arbor/backends/multicore/shared_state.hpp
+++ b/arbor/backends/multicore/shared_state.hpp
@@ -122,11 +122,9 @@ struct shared_state {
     fvm_size_type n_intdom = 0; // Number of integration domains.
     fvm_size_type n_detector = 0; // Max number of detectors on all cells.
     fvm_size_type n_cv = 0;   // Total number of CVs.
-    fvm_size_type n_gj = 0;   // Total number of GJs.
 
     iarray cv_to_intdom;      // Maps CV index to integration domain index.
     iarray cv_to_cell;        // Maps CV index to the first spike
-    gjarray gap_junctions;   // Stores gap_junction info.
     array time;               // Maps intdom index to integration start time [ms].
     array time_to;            // Maps intdom index to integration stop time [ms].
     array dt_intdom;          // Maps  index to (stop time) - (start time) [ms].
@@ -157,7 +155,6 @@ struct shared_state {
         fvm_size_type n_detector,
         const std::vector<fvm_index_type>& cv_to_intdom_vec,
         const std::vector<fvm_index_type>& cv_to_cell_vec,
-        const std::vector<fvm_gap_junction>& gj_vec,
         const std::vector<fvm_value_type>& init_membrane_potential,
         const std::vector<fvm_value_type>& temperature_K,
         const std::vector<fvm_value_type>& diam,
@@ -190,9 +187,6 @@ struct shared_state {
     // Set the per-integration domain and per-compartment dt from time_to - time.
     void set_dt();
 
-    // Update gap_junction state
-    void add_gj_current();
-
     // Update stimulus state and add current contributions.
     void add_stimulus_current();
 
diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp
index 7cd703b4..1432b6a1 100644
--- a/arbor/cable_cell.cpp
+++ b/arbor/cable_cell.cpp
@@ -74,8 +74,12 @@ struct cable_cell_impl {
         return location_map.get<T>();
     }
 
-    mlocation_map<mechanism_desc>& get_location_map(const mechanism_desc& desc) {
-        return location_map.get<mechanism_desc>()[desc.name()];
+    mlocation_map<synapse>& get_location_map(const synapse& desc) {
+        return location_map.get<synapse>()[desc.mech.name()];
+    }
+
+    mlocation_map<junction>& get_location_map(const junction& desc) {
+        return location_map.get<junction>()[desc.mech.name()];
     }
 
     template <typename Item>
@@ -98,8 +102,8 @@ struct cable_cell_impl {
         return region_map.get<T>();
     }
 
-    mcable_map<mechanism_desc>& get_region_map(const mechanism_desc& desc) {
-        return region_map.get<mechanism_desc>()[desc.name()];
+    mcable_map<density>& get_region_map(const density& desc) {
+        return region_map.get<density>()[desc.mech.name()];
     }
 
     mcable_map<init_int_concentration>& get_region_map(const init_int_concentration& init) {
@@ -210,11 +214,11 @@ const std::unordered_multimap<cell_tag_type, lid_range>& cable_cell::detector_ra
 }
 
 const std::unordered_multimap<cell_tag_type, lid_range>& cable_cell::synapse_ranges() const {
-    return impl_->labeled_lid_ranges.get<mechanism_desc>();
+    return impl_->labeled_lid_ranges.get<synapse>();
 }
 
-const std::unordered_multimap<cell_tag_type, lid_range>& cable_cell::gap_junction_ranges() const {
-    return impl_->labeled_lid_ranges.get<gap_junction_site>();
+const std::unordered_multimap<cell_tag_type, lid_range>& cable_cell::junction_ranges() const {
+    return impl_->labeled_lid_ranges.get<junction>();
 }
 
 } // namespace arb
diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index 12ed2a34..438f3313 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -738,8 +738,10 @@ fvm_mechanism_data& append(fvm_mechanism_data& left, const fvm_mechanism_data& r
 
             L.kind = R.kind;
             append(L.cv, R.cv);
+            append(L.peer_cv, R.peer_cv);
             append(L.multiplicity, R.multiplicity);
             append(L.norm_area, R.norm_area);
+            append(L.local_weight, R.local_weight);
             append_offset(L.target, target_offset, R.target);
 
             arb_assert(util::equal(L.param_values, R.param_values,
@@ -769,15 +771,70 @@ fvm_mechanism_data& append(fvm_mechanism_data& left, const fvm_mechanism_data& r
     return left;
 }
 
-fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties& gprop,
-    const cable_cell& cell, const fvm_cv_discretization& D, fvm_size_type cell_idx);
+std::unordered_map<cell_member_type, fvm_size_type> fvm_build_gap_junction_cv_map(
+    const std::vector<cable_cell>& cells,
+    const std::vector<cell_gid_type>& gids,
+    const fvm_cv_discretization& D)
+{
+    arb_assert(cells.size() == gids.size());
+    std::unordered_map<cell_member_type, fvm_size_type> gj_cvs;
+    for (auto cell_idx: util::make_span(0, cells.size())) {
+        for (const auto& mech : cells[cell_idx].junctions()) {
+            for (const auto& gj: mech.second) {
+                gj_cvs.insert({cell_member_type{gids[cell_idx], gj.lid}, D.geometry.location_cv(cell_idx, gj.loc, cv_prefer::cv_nonempty)});
+            }
+        }
+    }
+    return gj_cvs;
+}
 
-fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties& gprop,
-    const std::vector<cable_cell>& cells, const fvm_cv_discretization& D, const execution_context& ctx)
+std::unordered_map<cell_gid_type, std::vector<fvm_gap_junction>> fvm_resolve_gj_connections(
+    const std::vector<cell_gid_type>& gids,
+    const cell_label_range& gj_data,
+    const std::unordered_map<cell_member_type, fvm_size_type>& gj_cvs,
+    const recipe& rec)
+{
+    // Construct and resolve all gj_connections.
+    std::unordered_map<cell_gid_type, std::vector<fvm_gap_junction>> gj_conns;
+    label_resolution_map resolution_map({gj_data, gids});
+    auto gj_resolver = resolver(&resolution_map);
+    for (const auto& gid: gids) {
+        std::vector<fvm_gap_junction> local_conns;
+        for (const auto& conn: rec.gap_junctions_on(gid)) {
+            auto local_idx = gj_resolver.resolve({gid, conn.local});
+            auto peer_idx  = gj_resolver.resolve(conn.peer);
+
+            auto local_cv = gj_cvs.at({gid, local_idx});
+            auto peer_cv  = gj_cvs.at({conn.peer.gid, peer_idx});
+
+            local_conns.push_back({local_idx, local_cv, peer_cv, conn.weight});
+        }
+        // Sort local_conns by local_cv.
+        util::sort(local_conns);
+        gj_conns[gid] = std::move(local_conns);
+    }
+    return gj_conns;
+}
+
+fvm_mechanism_data fvm_build_mechanism_data(
+    const cable_cell_global_properties& gprop,
+    const cable_cell& cell,
+    const std::vector<fvm_gap_junction>& gj_conns,
+    const fvm_cv_discretization& D,
+    fvm_size_type cell_idx);
+
+fvm_mechanism_data fvm_build_mechanism_data(
+    const cable_cell_global_properties& gprop,
+    const std::vector<cable_cell>& cells,
+    const std::vector<cell_gid_type>& gids,
+    const std::unordered_map<cell_gid_type, std::vector<fvm_gap_junction>>& gj_conns,
+    const fvm_cv_discretization& D,
+    const execution_context& ctx)
 {
     std::vector<fvm_mechanism_data> cell_mech(cells.size());
-    threading::parallel_for::apply(0, cells.size(), ctx.thread_pool.get(),
-          [&] (int i) { cell_mech[i]=fvm_build_mechanism_data(gprop, cells[i], D, i);});
+    threading::parallel_for::apply(0, cells.size(), ctx.thread_pool.get(), [&] (int i) {
+        cell_mech[i] = fvm_build_mechanism_data(gprop, cells[i], gj_conns.at(gids[i]), D, i);
+    });
 
     fvm_mechanism_data combined;
     for (auto cell_idx: count_along(cells)) {
@@ -788,8 +845,12 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
 
 // Construct FVM mechanism data for a single cell.
 
-fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties& gprop,
-    const cable_cell& cell, const fvm_cv_discretization& D, fvm_size_type cell_idx)
+fvm_mechanism_data fvm_build_mechanism_data(
+    const cable_cell_global_properties& gprop,
+    const cable_cell& cell,
+    const std::vector<fvm_gap_junction>& gj_conns,
+    const fvm_cv_discretization& D,
+    fvm_size_type cell_idx)
 {
     using size_type = fvm_size_type;
     using index_type = fvm_index_type;
@@ -854,7 +915,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
 
     // Density mechanisms:
 
-    for (const auto& entry: cell.region_assignments().get<mechanism_desc>()) {
+    for (const auto& entry: cell.region_assignments().get<density>()) {
         const std::string& name = entry.first;
         mechanism_info info = catalogue[name];
 
@@ -862,7 +923,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
         if (info.kind != arb_mechanism_kind_density) {
             throw cable_cell_error("expected density mechanism, got " +name +" which has " +arb_mechsnism_kind_str(info.kind));
         }
-        config.kind = info.kind;
+        config.kind = arb_mechanism_kind_density;
 
         std::vector<std::string> param_names;
         assign(param_names, util::keys(info.parameters));
@@ -885,9 +946,10 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
         param_maps.resize(n_param);
 
         for (auto& on_cable: entry.second) {
-            verify_mechanism(info, on_cable.second);
+            const auto& mech = on_cable.second.mech;
+            verify_mechanism(info, mech);
             mcable cable = on_cable.first;
-            const auto& set_params = on_cable.second.values();
+            const auto& set_params = mech.values();
 
             support.insert(cable, 1.);
             for (std::size_t i = 0; i<n_param; ++i) {
@@ -995,10 +1057,11 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
         arb_assert(ix==n_param);
 
         std::size_t offset = 0;
-        for (const placed<mechanism_desc>& pm: entry.second) {
-            verify_mechanism(info, pm.item);
+        for (const placed<synapse>& pm: entry.second) {
+            const auto& mech = pm.item.mech;
+            verify_mechanism(info, mech);
 
-            synapse_instance in;
+            synapse_instance in{};
 
             in.param_values_offset = offset;
             offset += n_param;
@@ -1007,7 +1070,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
             double* in_param = all_param_values.data()+in.param_values_offset;
             std::copy(default_param_value.begin(), default_param_value.end(), in_param);
 
-            for (const auto& kv: pm.item.values()) {
+            for (const auto& kv: mech.values()) {
                 in_param[param_index.at(kv.first)] = kv.second;
             }
 
@@ -1051,7 +1114,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
         bool coalesce = catalogue[name].linear && gprop.coalesce_synapses;
 
         fvm_mechanism_config config;
-        config.kind = info.kind;
+        config.kind = arb_mechanism_kind_point;
         for (auto& kv: info.parameters) {
             config.param_values.emplace_back(kv.first, std::vector<value_type>{});
             if (!coalesce) {
@@ -1089,6 +1152,85 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
     }
     M.post_events = post_events;
 
+    // Gap junctions:
+
+    struct junction_desc {
+        std::string name;                     // mechanism name.
+        std::vector<value_type> param_values; // overridden parameter values.
+    };
+
+    // Gap-junction mechanisms are handled differently from point mechanisms.
+    // There is a separate mechanism instance at the local site of every gap-junction connection,
+    // meaning there can be multiple gap-junction mechanism instances of the same type (name) per
+    // lid.
+    // As a result, building fvm_mechanism_config per junction mechanism is split into 2 phases.
+    // (1) For every type (name) of gap-junction mechanism used on the cell, an fvm_mechanism_config
+    //     object is constructed and only the kind and parameter names are set. The object is
+    //     stored in the `junction_configs` map. Another map `lid_junction_desc` containing the
+    //     name and parameter values of the mechanism per lid is stored, needed to complete the
+    //     description of the fvm_mechanism_config object in the next step.
+    // (2) For every gap-junction connection, the cv, peer_cv, local_weight and parameter values
+    //     of the mechanism present on the local lid of the connection are added to the
+    //     fvm_mechanism_config of that mechanism. This completes the fvm_mechanism_config
+    //     description for each gap-junction mechanism.
+
+    std::unordered_map<std::string, fvm_mechanism_config> junction_configs;
+    std::unordered_map<cell_lid_type, junction_desc> lid_junction_desc;
+    for (const auto& [name, placements]: cell.junctions()) {
+        mechanism_info info = catalogue[name];
+        if (info.kind != arb_mechanism_kind_gap_junction) {
+            throw cable_cell_error("expected gap_junction mechanism, got " +name +" which has " +arb_mechsnism_kind_str(info.kind));
+        }
+
+        fvm_mechanism_config config;
+        config.kind = arb_mechanism_kind_gap_junction;
+
+        std::vector<std::string> param_names;
+        std::vector<double> param_dflt;
+
+        assign(param_names, util::keys(info.parameters));
+        std::size_t n_param = param_names.size();
+
+        param_dflt.reserve(n_param);
+        for (const auto& p: param_names) {
+            config.param_values.emplace_back(p, std::vector<value_type>{});
+            param_dflt.push_back(info.parameters.at(p).default_value);
+        }
+
+        for (const placed<junction>& pm: placements) {
+            const auto& mech = pm.item.mech;
+            verify_mechanism(info, mech);
+            const auto& set_params = mech.values();
+
+            junction_desc per_lid;
+            per_lid.name = name;
+            for (std::size_t i = 0; i<n_param; ++i) {
+                per_lid.param_values.push_back(value_by_key(set_params, param_names[i]).value_or(param_dflt[i]));
+            }
+            lid_junction_desc.insert({pm.lid, std::move(per_lid)});
+        }
+        junction_configs[name] = std::move(config);
+    }
+
+    // Iterate over the gj_conns local to the cell, and complete the fvm_mechanism_config.
+    // The gj_conns are expected to be sorted by local CV index.
+    for (const auto& conn: gj_conns) {
+        auto local_junction_desc = lid_junction_desc[conn.local_idx];
+        auto& config = junction_configs[local_junction_desc.name];
+
+        config.cv.push_back(conn.local_cv);
+        config.peer_cv.push_back(conn.peer_cv);
+        config.local_weight.push_back(conn.weight);
+        for (unsigned i = 0; i < local_junction_desc.param_values.size(); ++i) {
+            config.param_values[i].second.push_back(local_junction_desc.param_values[i]);
+        }
+    }
+
+    // Add non-empty fvm_mechanism_config to the fvm_mechanism_data
+    for (auto [name, config]: junction_configs) {
+        if (!config.cv.empty()) M.mechanisms[name] = std::move(config);
+    }
+
     // Stimuli:
 
     if (!cell.stimuli().empty()) {
@@ -1259,7 +1401,7 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
                 }
                 else {
                     fvm_mechanism_config config;
-                    config.kind = info.kind;
+                    config.kind = arb_mechanism_kind_reversal_potential;
                     config.cv = M.ions[ion].cv;
                     config.norm_area.assign(config.cv.size(), 1.);
 
diff --git a/arbor/fvm_layout.hpp b/arbor/fvm_layout.hpp
index 6d39d93e..dc9ffa82 100644
--- a/arbor/fvm_layout.hpp
+++ b/arbor/fvm_layout.hpp
@@ -242,6 +242,12 @@ struct fvm_mechanism_config {
     // For each instance index i, there are multiplicity[i] consecutive entries.
     std::vector<index_type> target;
 
+    // Gap junction peer CV index (gap junction mechanisms only)
+    std::vector<index_type> peer_cv;
+
+    // Gap junction weight, unit-less (gap junction mechanisms only)
+    std::vector<value_type> local_weight;
+
     // (Non-global) parameters and parameter values across the mechanism instance.
     std::vector<std::pair<std::string, std::vector<value_type>>> param_values;
 };
@@ -285,6 +291,19 @@ struct fvm_stimulus_config {
     std::vector<std::vector<double>> envelope_amplitude; // [A/m²]
 };
 
+// Maps gj {gid, lid} locations on a cell to their CV indices.
+std::unordered_map<cell_member_type, fvm_size_type> fvm_build_gap_junction_cv_map(
+    const std::vector<cable_cell>& cells,
+    const std::vector<cell_gid_type>& gids,
+    const fvm_cv_discretization& D);
+
+// Resolves gj_connections into {gid, lid} pairs, then to CV indices and a weight.
+std::unordered_map<cell_gid_type, std::vector<fvm_gap_junction>> fvm_resolve_gj_connections(
+    const std::vector<cell_gid_type>& gids,
+    const cell_label_range& gj_data,
+    const std::unordered_map<cell_member_type, fvm_size_type>& gj_cv,
+    const recipe& rec);
+
 struct fvm_mechanism_data {
     // Mechanism config, indexed by mechanism name.
     std::unordered_map<std::string, fvm_mechanism_config> mechanisms;
@@ -305,6 +324,12 @@ struct fvm_mechanism_data {
     bool post_events = false;
 };
 
-fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties& gprop, const std::vector<cable_cell>& cells, const fvm_cv_discretization& D, const arb::execution_context& ctx={});
+fvm_mechanism_data fvm_build_mechanism_data(
+    const cable_cell_global_properties& gprop,
+    const std::vector<cable_cell>& cells,
+    const std::vector<cell_gid_type>& gids,
+    const std::unordered_map<cell_gid_type, std::vector<fvm_gap_junction>>& gj_conns,
+    const fvm_cv_discretization& D,
+    const arb::execution_context& ctx={});
 
 } // namespace arb
diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp
index 9e4f5ed8..c3dbdfad 100644
--- a/arbor/fvm_lowered_cell_impl.hpp
+++ b/arbor/fvm_lowered_cell_impl.hpp
@@ -61,13 +61,6 @@ public:
         std::vector<deliverable_event> staged_events,
         std::vector<sample_event> staged_samples) override;
 
-    std::vector<fvm_gap_junction> fvm_gap_junctions(
-        const std::vector<cable_cell>& cells,
-        const std::vector<cell_gid_type>& gids,
-        const cell_label_range& gj_data,
-        const recipe& rec,
-        const fvm_cv_discretization& D);
-
     // Generates indom index for every gid, guarantees that gids belonging to the same supercell are in the same intdom
     // Fills cell_to_intdom map; returns number of intdoms
     fvm_size_type fvm_intdom(
@@ -252,9 +245,6 @@ fvm_integration_result fvm_lowered_cell_impl<Backend>::integrate(
             m->update_current();
         }
 
-        // Add current contribution from gap_junctions
-        state_->add_gj_current();
-
         PE(advance_integrate_events);
         state_->deliverable_events.drop_marked_events();
 
@@ -425,7 +415,7 @@ fvm_initialization_data fvm_lowered_cell_impl<Backend>::initialize(
         }
         fvm_info.num_targets[gid] = count;
 
-        for (const auto& [label, range]: c.gap_junction_ranges()) {
+        for (const auto& [label, range]: c.junction_ranges()) {
             fvm_info.gap_junction_data.add_label(label, range);
         }
     }
@@ -471,13 +461,14 @@ fvm_initialization_data fvm_lowered_cell_impl<Backend>::initialize(
                               D.cv_capacitance, D.face_conductance, D.cv_area, fvm_info.cell_to_intdom);
     sample_events_ = sample_event_stream(nintdom);
 
-    // Discretize mechanism data.
+    // Discretize and build gap junction info.
 
-    fvm_mechanism_data mech_data = fvm_build_mechanism_data(global_props, cells, D, context_);
+    auto gj_cvs = fvm_build_gap_junction_cv_map(cells, gids, D);
+    auto gj_conns = fvm_resolve_gj_connections(gids, fvm_info.gap_junction_data, gj_cvs, rec);
 
-    // Discretize and build gap junction info.
+    // Discretize mechanism data.
 
-    auto gj_vector = fvm_gap_junctions(cells, gids, fvm_info.gap_junction_data, rec, D);
+    fvm_mechanism_data mech_data = fvm_build_mechanism_data(global_props, cells, gids, gj_conns, D, context_);
 
     // Fill src_to_spike and cv_to_cell vectors only if mechanisms with post_events implemented are present.
     post_events_ = mech_data.post_events;
@@ -506,7 +497,7 @@ fvm_initialization_data fvm_lowered_cell_impl<Backend>::initialize(
             [&](const std::string& name) { return mech_instance(name).mech->data_alignment(); }));
 
     state_ = std::make_unique<shared_state>(
-                nintdom, ncell, max_detector, cv_to_intdom, std::move(cv_to_cell), gj_vector,
+                nintdom, ncell, max_detector, cv_to_intdom, std::move(cv_to_cell),
                 D.init_membrane_potential, D.temperature_K, D.diam_um, std::move(src_to_spike),
                 data_alignment? data_alignment: 1u);
 
@@ -540,6 +531,7 @@ fvm_initialization_data fvm_lowered_cell_impl<Backend>::initialize(
         mechanism_layout layout;
         layout.cv = config.cv;
         layout.multiplicity = config.multiplicity;
+        layout.peer_cv = config.peer_cv;
         layout.weight.resize(layout.cv.size());
 
         std::vector<fvm_index_type> multiplicity_divs;
@@ -571,6 +563,15 @@ fvm_initialization_data fvm_lowered_cell_impl<Backend>::initialize(
                 }
             }
             break;
+        case arb_mechanism_kind_gap_junction:
+            // Junction mechanism contributions are in [nA] (µS * mV); CV area A in [µm^2].
+            // F = 1/A * [nA/µm²] / [A/m²] = 1000/A.
+
+            for (auto i: count_along(layout.cv)) {
+                auto cv = layout.cv[i];
+                layout.weight[i] = config.local_weight[i] * 1000/D.cv_area[cv];
+            }
+            break;
         case arb_mechanism_kind_density:
             // Current density contributions from mechanism are already in [A/m²].
 
@@ -638,47 +639,6 @@ fvm_initialization_data fvm_lowered_cell_impl<Backend>::initialize(
     return fvm_info;
 }
 
-// Get vector of gap_junctions
-template <typename Backend>
-std::vector<fvm_gap_junction> fvm_lowered_cell_impl<Backend>::fvm_gap_junctions(
-        const std::vector<cable_cell>& cells,
-        const std::vector<cell_gid_type>& gids,
-        const cell_label_range& gap_junction_data,
-        const recipe& rec, const fvm_cv_discretization& D) {
-
-    std::vector<fvm_gap_junction> gj_vec;
-
-    std::unordered_map<cell_gid_type, std::vector<unsigned>> gid_to_cvs;
-    for (auto cell_idx: util::make_span(0, D.n_cell())) {
-        if (rec.gap_junctions_on(gids[cell_idx]).empty()) continue;
-
-        const auto& cell_gj = cells[cell_idx].gap_junction_sites();
-        gid_to_cvs[gids[cell_idx]].reserve(cell_gj.size());
-
-        for (auto gj : cell_gj) {
-            auto cv = D.geometry.location_cv(cell_idx, gj.loc, cv_prefer::cv_nonempty);
-            gid_to_cvs[gids[cell_idx]].push_back(cv);
-        }
-    }
-    label_resolution_map resolution_map({gap_junction_data, gids});
-    auto gj_resolver = resolver(&resolution_map);
-    for (auto gid: gids) {
-        auto gj_list = rec.gap_junctions_on(gid);
-        for (const auto& g: gj_list) {
-            if (g.local.policy != lid_selection_policy::assert_univalent) {
-                throw gj_unsupported_lid_selection_policy(gid, g.local.tag);
-            }
-            if (g.peer.label.policy != lid_selection_policy::assert_univalent) {
-                throw gj_unsupported_lid_selection_policy(g.peer.gid, g.peer.label.tag);
-            }
-            auto cv_local = gid_to_cvs[gid][gj_resolver.resolve({gid, g.local})];
-            auto cv_peer = gid_to_cvs[g.peer.gid][gj_resolver.resolve(g.peer)];
-            gj_vec.emplace_back(fvm_gap_junction(std::make_pair(cv_local, cv_peer), g.ggap * 1e3 / D.cv_area[cv_local]));
-        }
-    }
-    return gj_vec;
-}
-
 template <typename Backend>
 fvm_size_type fvm_lowered_cell_impl<Backend>::fvm_intdom(
         const recipe& rec,
@@ -756,7 +716,7 @@ struct probe_resolution_data {
 
     // Extent of density mechanism on cell.
     mextent mechanism_support(const std::string& name) const {
-        auto& mech_map = cell.region_assignments().template get<mechanism_desc>();
+        auto& mech_map = cell.region_assignments().template get<density>();
         auto opt_mm = util::value_by_key(mech_map, name);
 
         return opt_mm? opt_mm->support(): mextent{};
diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp
index 2eb12914..38331263 100644
--- a/arbor/include/arbor/cable_cell.hpp
+++ b/arbor/include/arbor/cable_cell.hpp
@@ -188,7 +188,7 @@ struct cable_cell_impl;
 template <typename T>
 using region_assignment =
     std::conditional_t<
-        std::is_same<T, mechanism_desc>::value || std::is_same<T, init_int_concentration>::value ||
+        std::is_same<T, density>::value || std::is_same<T, init_int_concentration>::value ||
         std::is_same<T, init_ext_concentration>::value || std::is_same<T, init_reversal_potential>::value,
         std::unordered_map<std::string, mcable_map<T>>,
         mcable_map<T>>;
@@ -207,17 +207,17 @@ using mlocation_map = std::vector<placed<T>>;
 template <typename T>
 using location_assignment =
     std::conditional_t<
-        std::is_same<T, mechanism_desc>::value,
+        std::is_same<T, synapse>::value || std::is_same<T, junction>::value,
         std::unordered_map<std::string, mlocation_map<T>>,
         mlocation_map<T>>;
 
 using cable_cell_region_map = static_typed_map<region_assignment,
-    mechanism_desc, init_membrane_potential, axial_resistivity,
+    density, init_membrane_potential, axial_resistivity,
     temperature_K, membrane_capacitance, init_int_concentration,
     init_ext_concentration, init_reversal_potential>;
 
 using cable_cell_location_map = static_typed_map<location_assignment,
-    mechanism_desc, i_clamp, gap_junction_site, threshold_detector>;
+    synapse, junction, i_clamp, threshold_detector>;
 
 // High-level abstract representation of a cell.
 class cable_cell {
@@ -226,8 +226,6 @@ public:
     using size_type = cell_local_size_type;
     using value_type = double;
 
-    using gap_junction_instance = mlocation;
-
     // Default constructor.
     cable_cell();
 
@@ -257,12 +255,12 @@ public:
 
     // Convenience access to placed items.
 
-    const std::unordered_map<std::string, mlocation_map<mechanism_desc>>& synapses() const {
-        return location_assignments().get<mechanism_desc>();
+    const std::unordered_map<std::string, mlocation_map<synapse>>& synapses() const {
+        return location_assignments().get<synapse>();
     }
 
-    const mlocation_map<gap_junction_site>& gap_junction_sites() const {
-        return location_assignments().get<gap_junction_site>();
+    const std::unordered_map<std::string, mlocation_map<junction>>& junctions() const {
+        return location_assignments().get<junction>();
     }
 
     const mlocation_map<threshold_detector>& detectors() const {
@@ -292,7 +290,7 @@ public:
     // The labeled lid_ranges of sources, targets and gap_junctions on the cell;
     const std::unordered_multimap<cell_tag_type, lid_range>& detector_ranges() const;
     const std::unordered_multimap<cell_tag_type, lid_range>& synapse_ranges() const;
-    const std::unordered_multimap<cell_tag_type, lid_range>& gap_junction_ranges() const;
+    const std::unordered_multimap<cell_tag_type, lid_range>& junction_ranges() const;
 
 private:
     std::unique_ptr<cable_cell_impl, void (*)(cable_cell_impl*)> impl_;
diff --git a/arbor/include/arbor/cable_cell_param.hpp b/arbor/include/arbor/cable_cell_param.hpp
index af7e6671..7d36e223 100644
--- a/arbor/include/arbor/cable_cell_param.hpp
+++ b/arbor/include/arbor/cable_cell_param.hpp
@@ -87,9 +87,6 @@ struct threshold_detector {
     double threshold;
 };
 
-// Tag type for dispatching cable_cell::place() calls that add gap junction sites.
-struct gap_junction_site {};
-
 // Setter types for painting physical and ion parameters or setting
 // cell-wide default:
 
@@ -187,6 +184,37 @@ private:
     std::unordered_map<std::string, double> param_;
 };
 
+// Tagged mechanism types for dispatching decor::place() and decor::paint() calls
+struct junction {
+    mechanism_desc mech;
+    explicit junction(mechanism_desc m): mech(std::move(m)) {}
+    junction(mechanism_desc m, const std::unordered_map<std::string, double>& params): mech(std::move(m)) {
+        for (const auto& [param, value]: params) {
+            mech.set(param, value);
+        }
+    }
+};
+
+struct synapse {
+    mechanism_desc mech;
+    explicit synapse(mechanism_desc m): mech(std::move(m)) {}
+    synapse(mechanism_desc m, const std::unordered_map<std::string, double>& params): mech(std::move(m)) {
+        for (const auto& [param, value]: params) {
+            mech.set(param, value);
+        }
+    }
+};
+
+struct density {
+    mechanism_desc mech;
+    explicit density(mechanism_desc m): mech(std::move(m)) {}
+    density(mechanism_desc m, const std::unordered_map<std::string, double>& params): mech(std::move(m)) {
+        for (const auto& [param, value]: params) {
+            mech.set(param, value);
+        }
+    }
+};
+
 struct ion_reversal_potential_method {
     std::string ion;
     mechanism_desc method;
@@ -200,13 +228,13 @@ using paintable =
                  init_int_concentration,
                  init_ext_concentration,
                  init_reversal_potential,
-                 mechanism_desc>;
+                 density>;
 
 using placeable =
-    std::variant<mechanism_desc,
-                 i_clamp,
+    std::variant<i_clamp,
                  threshold_detector,
-                 gap_junction_site>;
+                 synapse,
+                 junction>;
 
 using defaultable =
     std::variant<init_membrane_potential,
diff --git a/arbor/include/arbor/fvm_types.hpp b/arbor/include/arbor/fvm_types.hpp
index 0331b294..1b884965 100644
--- a/arbor/include/arbor/fvm_types.hpp
+++ b/arbor/include/arbor/fvm_types.hpp
@@ -10,14 +10,11 @@ using fvm_size_type  = arb_size_type;
 using fvm_index_type = arb_index_type;
 
 struct fvm_gap_junction {
-    using value_type = fvm_value_type;
-    using index_type = fvm_index_type;
-
-    std::pair<index_type, index_type> loc;
-    value_type weight = 0;
-
-    fvm_gap_junction() = default;
-    fvm_gap_junction(std::pair<index_type, index_type> l, value_type w): loc(l), weight(w) {}
+    cell_lid_type local_idx; // Index relative to other gap junction sites on the cell.
+    fvm_size_type local_cv;  // CV index of the local gap junction site.
+    fvm_size_type peer_cv;   // CV index of the peer gap junction site.
+    fvm_value_type weight;   // unit-less local weight of the connection.
 };
+ARB_DEFINE_LEXICOGRAPHIC_ORDERING(fvm_gap_junction, (a.local_cv, a.peer_cv, a.local_idx, a.weight), (b.local_cv, b.peer_cv, b.local_idx, b.weight))
 
 } // namespace arb
diff --git a/arbor/include/arbor/mechanism.hpp b/arbor/include/arbor/mechanism.hpp
index b6aa7919..62c90f99 100644
--- a/arbor/include/arbor/mechanism.hpp
+++ b/arbor/include/arbor/mechanism.hpp
@@ -74,6 +74,9 @@ struct mechanism_layout {
     // Maps in-instance index to CV index.
     std::vector<fvm_index_type> cv;
 
+    // Maps in-instance index to peer CV index (only for gap-junctions).
+    std::vector<fvm_index_type> peer_cv;
+
     // Maps in-instance index to compartment contribution.
     std::vector<fvm_value_type> weight;
 
diff --git a/arbor/include/arbor/mechanism_abi.h b/arbor/include/arbor/mechanism_abi.h
index 83770e2e..69e11c7c 100644
--- a/arbor/include/arbor/mechanism_abi.h
+++ b/arbor/include/arbor/mechanism_abi.h
@@ -21,6 +21,7 @@ typedef uint32_t arb_mechanism_kind;
 #define arb_mechanism_kind_point 1
 #define arb_mechanism_kind_density 2
 #define arb_mechanism_kind_reversal_potential 3
+#define arb_mechanism_kind_gap_junction 4
 
 typedef uint32_t arb_backend_kind;
 #define arb_backend_kind_nil 0
@@ -91,6 +92,7 @@ typedef struct arb_mechanism_ppack {
     arb_value_type* diam_um;
     arb_value_type* time_since_spike;
     arb_index_type* node_index;
+    arb_index_type* peer_index;
     arb_index_type* multiplicity;
     arb_value_type* weight;
     arb_size_type   mechanism_id;
diff --git a/arbor/include/arbor/mechinfo.hpp b/arbor/include/arbor/mechinfo.hpp
index 904fa90d..f0b55295 100644
--- a/arbor/include/arbor/mechinfo.hpp
+++ b/arbor/include/arbor/mechinfo.hpp
@@ -53,7 +53,6 @@ struct ion_dependency {
 using mechanism_fingerprint = std::string;
 
 struct mechanism_info {
-
     // mechanism_info is a convenient subset of the ABI mech description
     mechanism_info(const arb_mechanism_type&);
     mechanism_info() = default;
diff --git a/arbor/include/arbor/recipe.hpp b/arbor/include/arbor/recipe.hpp
index 76be9ddc..2853d4d0 100644
--- a/arbor/include/arbor/recipe.hpp
+++ b/arbor/include/arbor/recipe.hpp
@@ -54,11 +54,10 @@ struct cell_connection {
 struct gap_junction_connection {
     cell_global_label_type peer;
     cell_local_label_type local;
-
-    double ggap;
+    double weight; //unit-less
 
     gap_junction_connection(cell_global_label_type peer, cell_local_label_type local, double g):
-        peer(std::move(peer)), local(std::move(local)), ggap(g) {}
+        peer(std::move(peer)), local(std::move(local)), weight(g) {}
 };
 
 class recipe {
diff --git a/arbor/mechinfo.cpp b/arbor/mechinfo.cpp
index 341568a7..b2dbd7d0 100644
--- a/arbor/mechinfo.cpp
+++ b/arbor/mechinfo.cpp
@@ -4,7 +4,7 @@
 
 namespace arb {
 mechanism_info::mechanism_info(const arb_mechanism_type& m) {
-    kind = m.kind;
+    kind        = m.kind;
     post_events = m.has_post_events;
     linear      = m.is_linear;
     fingerprint = m.fingerprint;
diff --git a/arborio/cableio.cpp b/arborio/cableio.cpp
index 469fe96d..25f6d6e2 100644
--- a/arborio/cableio.cpp
+++ b/arborio/cableio.cpp
@@ -66,9 +66,6 @@ s_expr mksexp(const i_clamp& c) {
 s_expr mksexp(const threshold_detector& d) {
     return slist("threshold-detector"_symbol, d.threshold);
 }
-s_expr mksexp(const gap_junction_site& s) {
-    return slist("gap-junction-site"_symbol);
-}
 s_expr mksexp(const mechanism_desc& d) {
     std::vector<s_expr> mech;
     mech.push_back(s_expr(d.name()));
@@ -79,6 +76,15 @@ s_expr mksexp(const mechanism_desc& d) {
 s_expr mksexp(const ion_reversal_potential_method& e) {
     return slist("ion-reversal-potential-method"_symbol, s_expr(e.ion), mksexp(e.method));
 }
+s_expr mksexp(const junction& j) {
+    return slist("junction"_symbol, mksexp(j.mech));
+}
+s_expr mksexp(const synapse& j) {
+    return slist("synapse"_symbol, mksexp(j.mech));
+}
+s_expr mksexp(const density& j) {
+    return slist("density"_symbol, mksexp(j.mech));
+}
 s_expr mksexp(const mpoint& p) {
     return slist("point"_symbol, p.x, p.y, p.z, p.radius);
 }
@@ -215,15 +221,14 @@ pulse_tuple make_envelope_pulse(double delay, double duration, double amplitude)
 arb::i_clamp make_i_clamp_pulse(pulse_tuple p, double freq, double phase) {
     return arb::i_clamp::box(std::get<0>(p), std::get<1>(p), std::get<2>(p), freq, phase);
 }
-arb::gap_junction_site make_gap_junction_site() {
-    return arb::gap_junction_site{};
-}
 arb::cv_policy make_cv_policy(const cv_policy& p) {
     return p;
 }
 arb::ion_reversal_potential_method make_ion_reversal_potential_method(const std::string& ion, const arb::mechanism_desc& mech) {
     return ion_reversal_potential_method{ion, mech};
 }
+template <typename T>
+T make_wrapped_mechanism(const arb::mechanism_desc& mech) {return T(mech);}
 #undef ARBIO_DEFINE_SINGLE_ARG
 #undef ARBIO_DEFINE_DOUBLE_ARG
 
@@ -587,19 +592,19 @@ eval_map named_evals{
         "'current-clamp' with 3 arguments (env:envelope_pulse freq:real phase:real)")},
     {"threshold-detector", make_call<double>(make_threshold_detector,
         "'threshold-detector' with 1 argument (threshold:real)")},
-    {"gap-junction-site", make_call<>(make_gap_junction_site,
-        "'gap-junction-site' with 0 arguments")},
-    {"ion-reversal-potential-method", make_call<std::string, arb::mechanism_desc>(
-        make_ion_reversal_potential_method,
+    {"mechanism", make_mech_call("'mechanism' with a name argument, and 0 or more parameter settings"
+        "(name:string (param:string val:real))")},
+    {"ion-reversal-potential-method", make_call<std::string, arb::mechanism_desc>( make_ion_reversal_potential_method,
         "'ion-reversal-potential-method' with 2 arguments (ion:string mech:mechanism)")},
     {"cv-policy", make_call<cv_policy>(make_cv_policy,
         "'cv-policy' with 1 argument (p:policy)")},
-    {"mechanism", make_mech_call("'mechanism' with a name argument, and 0 or more parameter settings"
-                                 "(name:string (param:string val:real))")},
-    {"place", make_call<locset, gap_junction_site, std::string>(make_place, "'place' with 3 arguments (ls:locset gj:gap-junction-site name:string)")},
+    {"junction", make_call<arb::mechanism_desc>(make_wrapped_mechanism<junction>, "'junction' with 1 argumnet (m: mechanism)")},
+    {"synapse",  make_call<arb::mechanism_desc>(make_wrapped_mechanism<synapse>, "'synapse' with 1 argumnet (m: mechanism)")},
+    {"density",  make_call<arb::mechanism_desc>(make_wrapped_mechanism<density>, "'density' with 1 argumnet (m: mechanism)")},
     {"place", make_call<locset, i_clamp, std::string>(make_place, "'place' with 3 arguments (ls:locset c:current-clamp name:string)")},
     {"place", make_call<locset, threshold_detector, std::string>(make_place, "'place' with 3 arguments (ls:locset t:threshold-detector name:string)")},
-    {"place", make_call<locset, mechanism_desc, std::string>(make_place, "'place' with 3 arguments (ls:locset mech:mechanism name:string)")},
+    {"place", make_call<locset, junction, std::string>(make_place, "'place' with 3 arguments (ls:locset gj:junction name:string)")},
+    {"place", make_call<locset, synapse, std::string>(make_place, "'place' with 3 arguments (ls:locset mech:synapse name:string)")},
 
     {"paint", make_call<region, init_membrane_potential>(make_paint, "'paint' with 2 arguments (reg:region v:membrane-potential)")},
     {"paint", make_call<region, temperature_K>(make_paint, "'paint' with 2 arguments (reg:region v:temperature-kelvin)")},
@@ -608,7 +613,7 @@ eval_map named_evals{
     {"paint", make_call<region, init_int_concentration>(make_paint, "'paint' with 2 arguments (reg:region v:ion-internal-concentration)")},
     {"paint", make_call<region, init_ext_concentration>(make_paint, "'paint' with 2 arguments (reg:region v:ion-external-concentration)")},
     {"paint", make_call<region, init_reversal_potential>(make_paint, "'paint' with 2 arguments (reg:region v:ion-reversal-potential)")},
-    {"paint", make_call<region, mechanism_desc>(make_paint, "'paint' with 2 arguments (reg:region v:mechanism)")},
+    {"paint", make_call<region, density>(make_paint, "'paint' with 2 arguments (reg:region v:density)")},
 
     {"default", make_call<init_membrane_potential>(make_default, "'default' with 1 argument (v:membrane-potential)")},
     {"default", make_call<temperature_K>(make_default, "'default' with 1 argument (v:temperature-kelvin)")},
diff --git a/doc/concepts/cable_cell.rst b/doc/concepts/cable_cell.rst
index e81e218d..60ffdce9 100644
--- a/doc/concepts/cable_cell.rst
+++ b/doc/concepts/cable_cell.rst
@@ -5,7 +5,7 @@ Cable cells
 
 An Arbor *cable cell* is a full :ref:`description <modelcelldesc>` of a cell
 with morphology and cell dynamics like ion species and their properties, ion
-channels, synapses, gap junction sites, stimuli and spike detectors.
+channels, synapses, gap junction mechanisms, stimuli and spike detectors.
 
 Cable cells are constructed from three components:
 
@@ -17,7 +17,7 @@ When a cable cell is constructed the following steps are performed using the inp
 
 1. Concrete regions and locsets are generated for the morphology for each labelled region and locset in the dictionary
 2. The default values for parameters specified in the decor, such as ion species concentration, are instantiated.
-3. Dynamics (mechanisms, parameters, synapses, etc.) are instantiated on the regions and locsets as specified by the decor.
+3. Dynamics (mechanisms, parameters, synapses, gap junctions etc.) are instantiated on the regions and locsets as specified by the decor.
 
 Once constructed, the cable cell can be queried for specific information about the cell, but it can't be modified (it is *immutable*).
 
diff --git a/doc/concepts/decor.rst b/doc/concepts/decor.rst
index e1ec4bb2..3b15a827 100644
--- a/doc/concepts/decor.rst
+++ b/doc/concepts/decor.rst
@@ -20,8 +20,8 @@ The choice of region or locset is reflected in the two broad classes of dynamics
 * *Placed dynamics* are applied to locations on the cell, and are associated
   with entities that can be counted.
 
-  * :ref:`Synapses <cablecell-synapses>`.
-  * :ref:`Gap junction sites <cablecell-gj-sites>`.
+  * :ref:`Synapse mechanisms <cablecell-synapses>`.
+  * :ref:`Gap junction mechanisms <cablecell-gj-mechs>`.
   * :ref:`Threshold detectors <cablecell-threshold-detectors>` (spike detectors).
   * :ref:`Stimuli <cablecell-stimuli>`.
   * :ref:`Probes <cablecell-probes>`.
@@ -115,10 +115,10 @@ specialised on specific regions.
 ~~~~~~~~~~~~~~~~~~~~~
 
 Regions can have density mechanisms defined over their extents.
-Density mechanisms are :ref:`NMODL mechanisms <nmodl>`
-which describe biophysical processes. These are processes
-that are distributed in space, but whose behaviour is defined purely
-by the state of the cell and the process at any given point.
+:ref:`Density mechanisms <mechanisms-density>` are a kind of
+:ref:`NMODL mechanism <nmodl>` which describe biophysical processes.
+These are processes that are distributed in space, but whose behaviour is
+defined purely by the state of the cell and the process at any given point.
 
 The most common use for density mechanisms is to describe ion channel dynamics,
 for example the ``hh`` and ``pas`` mechanisms provided by NEURON and Arbor,
@@ -157,11 +157,12 @@ Take for example the built-in mechanism for passive leaky dynamics:
     # Create an instance of the same mechanism, that also sets conductance (range)
     m4 = arbor.mechanism('pas/e=-45', {'g': 0.1})
 
+    # And the mechanisms in `density` mechanism objects and add them to the decor.
     decor = arbor.decor()
-    decor.paint('"soma"', m1)
-    decor.paint('"soma"', m2) # error: can't place the same mechanism on overlapping regions
-    decor.paint('"soma"', m3) # error: can't have overlap between two instances of a mechanism
-                              #        with different values for a global parameter.
+    decor.paint('"soma"', arbor.density(m1))
+    decor.paint('"soma"', arbor.density(m2)) # error: can't place the same mechanism on overlapping regions
+    decor.paint('"soma"', arbor.density(m3)) # error: can't have overlap between two instances of a mechanism
+                                             #        with different values for a global parameter.
 
 .. _cablecell-ions:
 
@@ -271,20 +272,62 @@ locset.
 1. Connection sites
 ~~~~~~~~~~~~~~~~~~~
 
-Connections (synapses) are instances of NMODL POINT mechanisms. See also :term:`connection`.
+Similar to how regions can have density mechanisms defined over their extents,
+locsets can have point mechanisms placed on their individual locations.
+:ref:`Point mechanisms <mechanisms-point>` are a kind of :ref:`NMODL mechanism <nmodl>`
+which describe synaptic processes such as the ``expsyn`` mechanism provided by
+NEURON and Arbor, which models an exponential synapse.
 
-.. _cablecell-gj-sites:
+A point mechanism (synapse) can form the target of a :term:`connection` on a cell.
 
-2. Gap junction sites
-~~~~~~~~~~~~~~~~~~~~~
+.. code-block:: Python
+
+    decor = arbor.decor()
+
+    # Create an 'expsyn' mechanism with default parameter values (set in NMODL file).
+    expsyn = arbor.mechanism('expsyn')
+
+    # Wrap the 'expsyn' mechanism in a `synapse` object and add it to the decor.
+    decor.paint('"syn_loc_0"', arbor.synapse(expsyn))
+
+    # Create an 'expsyn' mechanism with default parameter values as a `synapse` object, and add it to the decor.
+    decor.paint('"syn_loc_1"', arbor.synapse("expsyn"))
+
+    # Create an 'expsyn' mechanism with modified 'tau' parameter as a `synapse` object, and add it to the decor.
+    decor.paint('"syn_loc_2"', arbor.synapse("expsyn", {"tau": 1.0}))
 
-See :term:`gap junction`.
 
 .. _cablecell-threshold-detectors:
 
-3. Threshold detectors (spike detectors).
+2. Threshold detectors (spike detectors).
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
+.. _cablecell-gj-mechs:
+
+3. Gap junction connection sites
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Locsets can also have junction mechanisms placed on their individual locations.
+:ref:`Junction mechanisms <mechanisms-junction>` are a kind of :ref:`NMODL mechanism <nmodl>`
+which describe gap-junction processes such as the ``gj`` mechanism provided by Arbor,
+which models a basic, linear, constant-conductance based gap-junction.
+
+A junction mechanism can form each of the endpoints of a :term:`gap junction connection`
+on two separate cells.
+
+.. code-block:: Python
+
+    decor = arbor.decor()
+
+    # Create a 'gj' mechanism with modified 'g' value.
+    gj = arbor.mechanism("gj", {"g": 2.0})
+
+    # Wrap the 'gj' mechanism in a `junction` object and add it to the decor.
+    decor.paint('"gj_loc_0"', arbor.junction(gj))
+
+    # Create a 'gj' mechanism with modified 'g' parameter as a `junction` object, and add it to the decor.
+    decor.paint('"gj_loc_1"', arbor.junction("gj", {"g": 1.5}))
+
 .. _cablecell-stimuli:
 
 4. Stimuli
diff --git a/doc/concepts/interconnectivity.rst b/doc/concepts/interconnectivity.rst
index cdc2c8ea..5ff303a2 100644
--- a/doc/concepts/interconnectivity.rst
+++ b/doc/concepts/interconnectivity.rst
@@ -5,8 +5,9 @@ Interconnectivity
 
 Networks can be regarded as graphs, where the nodes are locations on cells and the edges
 describe the communications between them. In Arbor, two sorts of edges are modelled: a
-:term:`connection` abstracts the propagation of action potentials (:term:`spikes <spike>`) through the network,
-while a :term:`gap junction` is used to describe a direct electrical connection between two locations on two cells.
+:term:`connection` abstracts the propagation of action potentials (:term:`spikes <spike>`)
+through the network, while a :term:`gap junction connection` is used to describe a direct
+electrical connection between two locations on two cells.
 Connections only capture the propagation delay and attenuation associated with spike
 connectivity: the biophysical modelling of the chemical synapses themselves is the
 responsibility of the target cell model.
@@ -51,7 +52,7 @@ A recipe lets you define which sites are connected to which.
 
 .. glossary::
 
-   gap junction
+   gap junction connection
       Gap junctions represent electrical synapses where transmission between cells is bidirectional and direct.
       They are modelled as a conductance between two **gap junction sites** on two cells.
 
@@ -63,7 +64,7 @@ A recipe lets you define which sites are connected to which.
          on the cell.
       2. Declare the Gap Junction connections in the recipe *on the local cell*: from a peer **gap junction site**
          identified using a :gen:`global_label`; to a local **gap junction site** identified using a
-         :gen:`local_label` (:gen:`gid` of the site is implicitly known); and a conductance in μS.
+         :gen:`local_label` (:gen:`gid` of the site is implicitly known); and a unit-less connection weight.
          Two of these connections are needed, on each of the peer and local cells.
 
    .. Note::
diff --git a/doc/concepts/mechanisms.rst b/doc/concepts/mechanisms.rst
index d98bc2d8..445347df 100644
--- a/doc/concepts/mechanisms.rst
+++ b/doc/concepts/mechanisms.rst
@@ -3,7 +3,7 @@
 Cable cell mechanisms
 =====================
 
-Mechanisms describe biophysical processes such as ion channels and synapses.
+Mechanisms describe biophysical processes such as ion channels, synapses and gap-junctions.
 Mechanisms are assigned to regions and locations on a cell morphology
 through the process of :ref:`decoration <cablecell-decoration>`.
 Mechanisms are described using a dialect of the :ref:`NMODL <nmodl>` domain
@@ -12,9 +12,9 @@ specific language that is similarly used in `NEURON <https://neuron.yale.edu/neu
 Arbor supports mechanism descriptions using the NMODL language through our ``modcc``
 compiler. ``modcc`` supports many of NMODL's features but there are a few
 additional :ref:`guidelines <formatnmodl>`.
-for users who wish to compile their own mechanisms for Arbor. Unfortunately, out-of-tree
-mechanism building is not yet supported in Arbor. However, we have many built-in mechanisms
-that can be used, which are oragnized in *mechanism catalogues*.
+for users who wish to compile their own mechanisms for Arbor. Out-of-tree mechanism
+building is available in Arbor (See: :ref:`mechanisms_dynamic`). We also have built-in
+mechanisms, which are organized in *mechanism catalogues*.
 
 Mechanism catalogues
 --------------------
@@ -68,6 +68,7 @@ Arbor provides the *default* catalogue with the following mechanisms:
   followed by an exponential decay (:ref:`point mechanism <mechanisms-point>`).
 * *exp2syn*: Bi-exponential conductance synapse described by two time constants:
   rise and decay (:ref:`point mechanism <mechanisms-point>`).
+* *gj*: Linear gap-junction mechanism with constant conductance (:ref:`junction mechanism <mechanisms-junction>`).
 
 With the exception of *nernst*, these mechanisms are the same as those available in NEURON.
 
@@ -178,9 +179,9 @@ shorthand ``("nernst/ca")`` can be used unambiguously.
 Mechanism types
 ---------------
 
-There are two broad categories of mechanism, density mechanisms and
-point mechanisms, and a third special density mechanism for
-computing ionic reversal potentials.
+There are three broad categories of mechanism: density mechanisms, point mechanisms,
+gap-junction mechanisms and a fourth special density mechanism for computing ionic
+reversal potential.
 
 .. _mechanisms-density:
 
@@ -195,6 +196,9 @@ Density mechanisms are commonly used to describe ion channel dynamics,
 for example the ``hh`` and ``pas`` mechanisms provided by NEURON and Arbor,
 which model classic Hodgkin-Huxley and passive leaky currents respectively.
 
+In NMODL, density mechanisms are identified using the ``SUFFIX`` keyword in the
+``NEURON`` block.
+
 .. _mechanisms-revpot:
 
 Ion reversal potential mechanisms
@@ -229,7 +233,7 @@ and ionic state.
 .. _mechanisms-point:
 
 Point mechanisms
-'''''''''''''''''''''''''''''''''
+''''''''''''''''
 
 *Point mechanisms*, which are associated with connection end points on a
 cable cell, are placed at discrete locations on the cell.
@@ -237,6 +241,27 @@ Unlike density mechanisms, whose behaviour is defined purely by the state of the
 and the process, their behaviour is additionally governed by the timing and weight of
 events delivered via incoming connections.
 
+In NMODL, point mechanisms are identified using the ``POINT_PROCESS`` keyword in the
+``NEURON`` block.
+
+.. _mechanisms-junction:
+
+Junction mechanisms
+'''''''''''''''''''
+
+*Junction mechanisms*, which are associated with gap-junction connection end points on a
+cable cell, are placed at discrete locations on the cell.
+A junction mechanism contributes a current at the discrete location of the cell on which it is placed.
+This current contribution depends on the state of the mechanism and the process, as well as the membrane
+potential at the discrete location which forms the other end of the gap-junction connection and the weight
+of that connection.
+
+In NMODL, junction mechanisms are identified using the ``JUNCTION_PROCESS`` keyword in the
+``NEURON`` block.
+
+.. note::
+    ``JUNCTION_PROCESS`` is an Arbor-specific extension to NMODL. The NMODL description of gap-junction
+    mechanisms in arbor is not identical to NEURON's though it is similar.
 
 API
 ---
diff --git a/doc/concepts/recipe.rst b/doc/concepts/recipe.rst
index ca1b4d45..7498471c 100644
--- a/doc/concepts/recipe.rst
+++ b/doc/concepts/recipe.rst
@@ -34,10 +34,10 @@ three cells:
   the **kind** of the cell.
 - ``Cell 1``: Is a soma and a single dendrite, with ``passive`` dynamics everywhere.
   It has a single synapse at the end of the dendrite labeled "syanpse_1" and a gap
-  junction site in the middle of the soma labeled "gap_junction_1".
+  junction mechanism in the middle of the soma labeled "gap_junction_1".
   This is the **description** of the cell. It's also a cable cell, which is its **cell kind**.
 - ``Cell 2``: Is a soma and a single dendrite, with ``passive`` dynamics everywhere.
-  It has a gap junction site in the middle of the soma labeled "gap_junction_2".
+  It has a gap junction mechanism in the middle of the soma labeled "gap_junction_2".
   This is the **description** of the cell. It's also a cable cell, which is its **cell kind**.
 
 The total **number of cells** in the model is 3. The **kind**, and **description** of each cell
diff --git a/doc/cpp/cable_cell.rst b/doc/cpp/cable_cell.rst
index aa2253a0..0a8b4366 100644
--- a/doc/cpp/cable_cell.rst
+++ b/doc/cpp/cable_cell.rst
@@ -47,7 +47,7 @@ given subsection of the morphology via the ``paint`` interface of the decor
 
 Ion channels and other distributed dynamical processes are also specified
 on the cell via the ``paint`` method; while synapses, current clamps,
-gap junction locations, and the site for testing the threshold potential
+gap junction mechanisms, and the site for testing the threshold potential
 are specified via the ``place`` method. See :ref:`cppcablecell-dynamics`, below.
 
 .. _cppcablecell-dynamics:
@@ -61,9 +61,10 @@ that are distributed in space, but whose behaviour is defined purely
 by the state of the cell and the process at any given point.
 
 Cells may also have *point* mechanisms, describing the dynamics
-at post-synaptic sites.
+at post-synaptic sites. And *junction* mechanisms, describing the
+dynamics at each site of the two sites of a gap-junction connection.
 
-A third type of mechanism, which describes ionic reversal potential
+A fourth type of mechanism, which describes ionic reversal potential
 behaviour, can be specified for cells or the whole model via cell parameter
 settings, described below.
 
@@ -87,21 +88,59 @@ mechanism name, and mechanism parameter values then set with the
    Returns a reference to the mechanism description, so that calls to
    :cpp:expr:`set` can be chained in a single expression.
 
+:cpp:type:`density`, :cpp:type:`synapse` and :cpp:type:`junction` objects are thin wrappers
+around a :cpp:type:`mechanism_desc`, needed for *painting* and *placing* mechanisms on a :cpp:type:`decor`.
+Relevant methods:
+
+.. cpp:function:: density::density(mechanism_desc mech)
+
+   Construct a density wrapper from the mechanism `mech`.
+
+.. cpp:function:: density::density(mechanism_desc mech, const std::unordered_map<std::string, double>& params)
+
+   For each ``{key, value}`` pair in `params`, set the parameter associated with ``key`` to ``value``
+   on mechanism ``mech``, then construct a density wrapper from the mechanism `mech`.
+
+.. cpp:function:: synapse::synapse(mechanism_desc mech)
+
+   Construct a synapse wrapper from the mechanism `mech`.
+
+.. cpp:function:: synapse::synapse(mechanism_desc mech, const std::unordered_map<std::string, double>& params)
+
+   For each ``{key, value}`` pair in `params`, set the parameter associated with ``key`` to ``value``
+   on mechanism ``mech``, then construct a synapse wrapper from the mechanism `mech`.
+
+.. cpp:function:: junction::junction(mechanism_desc mech)
+
+   Construct a junction wrapper from the mechanism `mech`.
+
+.. cpp:function:: junction::junction(mechanism_desc mech, const std::unordered_map<std::string, double>& params)
+
+   For each ``{key, value}`` pair in `params`, set the parameter associated with ``key`` to ``value``
+   on mechanism ``mech``, then construct a junction wrapper from the mechanism `mech`.
 
 Density mechanisms are associated with a cable cell object with:
 
-.. cpp:function:: void cable_cell::paint(const region&, mechanism_desc)
+.. cpp:function:: void cable_cell::paint(const region&, density)
 
 Point mechanisms, which are associated with connection end points on a
 cable cell, are placed on a set of locations given by a locset. The group
-of generated items requires a label. They are attached to a cell with:
+of generated items are given a label which can be used to create connections
+in the recipe. Point mechanisms are attached to a cell with:
+
+.. cpp:function:: void cable_cell::place(const locset&, synapse, cell_tag_type label)
+
+Gap-junction mechanisms, which are associated with gap-junction connection
+end points on a cable cell, are placed on a single location given by a locset
+(locsets with multiple locations will raise an exception). The generated item
+is given a label which can be used to create gap-junction connections in the
+recipe. Gap-junction mechanisms are attached to a cell with:
 
-.. cpp:function:: void cable_cell::place(const locset&, mechanism_desc, cell_tag_type label)
+.. cpp:function:: void cable_cell::place(const locset&, junction, cell_tag_type label)
 
 .. todo::
 
-   TODO: describe other ``place``-able things: current clamps, gap junction
-   sites, threshold potential measurement point.
+   TODO: describe other ``place``-able things: current clamps, threshold potential measurement point.
 
 .. _cppcablecell-electrical-properties:
 
@@ -228,7 +267,7 @@ A reversal potential mechanism described in NMODL:
 
 * May not maintain any STATE variables.
 * Can only write to the "eX" value associated with an ion.
-* Can not given as a POINT mechanism.
+* Can not be given as a POINT mechanism.
 
 Essentially, reversal potential mechanisms must be pure functions of cellular
 and ionic state.
diff --git a/doc/cpp/interconnectivity.rst b/doc/cpp/interconnectivity.rst
index 947498f4..e584902d 100644
--- a/doc/cpp/interconnectivity.rst
+++ b/doc/cpp/interconnectivity.rst
@@ -40,17 +40,16 @@ Interconnectivity
 
 .. cpp:class:: gap_junction_connection
 
-    Describes a gap junction between two gap junction sites. The :cpp:member:`local` site does not include
-    the gid of a cell, this is because a :cpp:class:`gap_junction_connection` is bound to the local
+    Describes a gap junction connection between two gap junction sites. The :cpp:member:`local` site does
+    not include the gid of a cell, this is because a :cpp:class:`gap_junction_connection` is bound to the local
     cell which means that the gid is implicitly known.
 
     .. note::
 
-       A bidirectional gap-junction between two cells ``c0`` and ``c1`` requires two
+       A bidirectional gap-junction connection between two cells ``c0`` and ``c1`` requires two
        :cpp:class:`gap_junction_connection` objects to be constructed: one where ``c0`` is the
        :cpp:member:`local` site, and ``c1`` is the :cpp:member:`peer` site; and another where ``c1`` is the
-       :cpp:member:`local` site, and ``c0`` is the :cpp:member:`peer` site. If :cpp:member:`ggap` is equal
-       in both connections, a symmetric gap-junction is formed, other wise the gap-junction is asymmetric.
+       :cpp:member:`local` site, and ``c0`` is the :cpp:member:`peer` site.
 
     .. cpp:member:: cell_global_label_type peer
 
@@ -63,6 +62,6 @@ Interconnectivity
         which packages a label of a group of gap junction sites on the cell and a selection policy.
         The gid of the local site's cell is implicitly known.
 
-    .. cpp:member:: float ggap
+    .. cpp:member:: float weight
 
-        gap junction conductance in μS.
+        unit-less gap junction connection weight.
diff --git a/doc/cpp/morphology.rst b/doc/cpp/morphology.rst
index ae8fa4bc..581290e2 100644
--- a/doc/cpp/morphology.rst
+++ b/doc/cpp/morphology.rst
@@ -115,7 +115,7 @@ by two stitches:
    stitched_morphology stitched(std::move(builder));
    cable_cell cell(stitched.morphology(), stitched.labels());
 
-   cell.paint("\"soma\"", "hh");
+   cell.paint("\"soma\"", density("hh"));
 
 
 .. _locsets-and-regions:
diff --git a/doc/fileformat/cable_cell.rst b/doc/fileformat/cable_cell.rst
index 95aebc30..745f35b6 100644
--- a/doc/fileformat/cable_cell.rst
+++ b/doc/fileformat/cable_cell.rst
@@ -95,9 +95,9 @@ of the cell):
    ion reversal potential method,      --,            --,             ✓
    density mechanism,                  --,             ✓,            --
    point mechanism,                    ✓,             --,            --
+   junction mechanism,                 ✓,             --,            --
    current clamp,                      ✓,             --,            --
    threshold detector,                 ✓,             --,            --
-   gap junction site,                  ✓,             --,            --
 
 The various properties and dynamics of the decor are described as follows:
 
@@ -153,6 +153,18 @@ The various properties and dynamics of the decor are described as follows:
 
       (ion-reversal-potential-method "ca" (mechanism "nernst/ca"))
 
+.. label:: (density method:mechanism)
+
+   This describes a *density* mechanism whose behavior is is defined by ``mechanism``.
+
+.. label:: (synapse method:mechanism)
+
+   This describes a *synapse* (point) mechanism whose behavior is is defined by ``mechanism``.
+
+.. label:: (junction method:mechanism)
+
+   This describes a *gap-junction* mechanism whose behavior is is defined by ``mechanism``.
+
 .. label:: (current-clamp (envelope-pulse delay:real duration:real amplitude:real) freq:real phase:real)
 
    This creates a *current clamp*. If the frequency ``freq`` (unit kHz) is zero, the current is a square
@@ -181,10 +193,6 @@ The various properties and dynamics of the decor are described as follows:
 
    This describes a *threshold-detector* object with value ``val`` (unit mV).
 
-.. label:: (gap-junction-site)
-
-   This describes a *gap-junction-site*.
-
 *Paintable* and *placeable* properties and dynamics are placed on regions (generated from :ref:`region expressions
 <labels-region-expr>`) and locsets (generated from :ref:`locset expressions <labels-locset-expr>`) respectively.
 *Defaultable* properties and dynamics apply to an entire cell.
@@ -243,10 +251,10 @@ Any number of paint, place and default expressions can be used to create a decor
         (default (membrane-potential -55.000000))
         (paint (region "custom") (temperature-kelvin 270))
         (paint (region "soma") (membrane-potential -50.000000))
-        (paint (all) (mechanism "pas"))
-        (paint (tag 4) (mechanism "Ih" ("gbar" 0.001)))
-        (place (locset "root") (mechanism "expsyn") "root_synapse")
-        (place (terminal) (gap-junction-site) "terminal_gj"))
+        (paint (all) (density (mechanism "pas")))
+        (paint (tag 4) (density (mechanism "Ih" ("gbar" 0.001))))
+        (place (locset "root") (synapse (mechanism "expsyn")) "root_synapse")
+        (place (terminal) (junction (mechanism "gj")) "terminal_gj"))
 
 Morphology
 ----------
@@ -337,9 +345,9 @@ expressions.
           (default (membrane-potential -55.000000))
           (paint (region "my_soma") (temperature-kelvin 270))
           (paint (region "my_region") (membrane-potential -50.000000))
-          (paint (tag 4) (mechanism "Ih" ("gbar" 0.001)))
-          (place (locset "root") (mechanism "expsyn") "root_synapse")
-          (place (location 1 0.2) (gap-junction-site) "terminal_gj"))
+          (paint (tag 4) (density (mechanism "Ih" ("gbar" 0.001))))
+          (place (locset "root") (synapse (mechanism "expsyn")) "root_synapse")
+          (place (location 1 0.2) (junction (mechanism "gj")) "terminal_gj"))
         (morphology
           (branch 0 -1
             (segment 0 (point 0 0 0 2) (point 4 0 0 2) 1)
@@ -405,7 +413,7 @@ Decoration
      (meta-data (version "0.1-dev"))
      (decor
        (default (membrane-potential -55.000000))
-       (place (locset "root") (mechanism "expsyn") "root_synapse")
+       (place (locset "root") (synapse (mechanism "expsyn")) "root_synapse")
        (paint (region "my_soma") (temperature-kelvin 270))))
 
 Morphology
@@ -434,7 +442,7 @@ Cable-cell
          (locset-def "root" (root)))
        (decor
          (default (membrane-potential -55.000000))
-         (place (locset "root") (mechanism "expsyn") "root_synapse")
+         (place (locset "root") (synapse (mechanism "expsyn")) "root_synapse")
          (paint (region "my_soma") (temperature-kelvin 270)))
        (morphology
           (branch 0 -1
diff --git a/doc/fileformat/nmodl.rst b/doc/fileformat/nmodl.rst
index fc40a000..961ce25c 100644
--- a/doc/fileformat/nmodl.rst
+++ b/doc/fileformat/nmodl.rst
@@ -109,3 +109,8 @@ Arbor-specific features
     POST_EVENT(t) {
        g = g + (0.1*t)
     }
+
+* Arbor allows a gap-junction mechanism to access the membrane potential at the peer site
+  of a gap-junction connection as well as the local site. The peer membrane potential is
+  made available through the ``v_peer`` variable while the local membrane potential
+  is available through ``v``, as usual.
\ No newline at end of file
diff --git a/doc/internals/mechanism_abi.rst b/doc/internals/mechanism_abi.rst
index 19ee56bf..df5d4250 100644
--- a/doc/internals/mechanism_abi.rst
+++ b/doc/internals/mechanism_abi.rst
@@ -59,6 +59,7 @@ This type collects all information independent of the backend.
     * point
     * density
     * reversal_potential
+    * junction
 
   .. c:member:: bool                      is_linear
 
@@ -222,6 +223,10 @@ fully formed to the interface. At this point:
 
     Indices of CVs covered by this mechanism, size is width
 
+  .. c:member:: arb_index_type* peer_index
+
+    Indices of peer CV of each CV in ``node_index``, needed for gap-junction connections, size is width.
+
   .. c:member:: arb_index_type* multiplicity
 
     [Unused]
diff --git a/doc/python/cable_cell.rst b/doc/python/cable_cell.rst
index 86fe945d..20018d83 100644
--- a/doc/python/cable_cell.rst
+++ b/doc/python/cable_cell.rst
@@ -47,9 +47,9 @@ Cable cells
 
         # Define decorations
         decor = arbor.decor()
-        decor.paint('"dend"', 'pas')
-        decor.paint('"axon"', 'hh')
-        decor.paint('"soma"', 'hh')
+        decor.paint('"dend"', arbor.density('pas'))
+        decor.paint('"axon"', arbor.density('hh'))
+        decor.paint('"soma"', arbor.density('hh'))
 
         # Construct a cable cell.
         cell = arbor.cable_cell(morph, labels, decor)
diff --git a/doc/python/decor.rst b/doc/python/decor.rst
index c9412966..33338141 100644
--- a/doc/python/decor.rst
+++ b/doc/python/decor.rst
@@ -67,7 +67,7 @@ Cable cell decoration
             decor.set_ion('na', int_con=5.0, rev_pot=70, method=None)
 
     Various specialisations of the ``paint`` method are available for setting properties
-    and mechanisms that are applied to regions.
+    and density mechanisms that are applied to regions.
 
     .. method:: paint(region, Vm=None, cm=None, rL=None, tempK=None)
 
@@ -104,59 +104,39 @@ Cable cell decoration
         :param float rev_pot: reversal potential [mV].
         :type rev_pot: float or None
 
-    .. method:: paint(region, mechanism)
+    .. method:: paint(region, density)
         :noindex:
 
-        Apply a mechanism with a region.
-        Returns a unique identifier that can be used to query the local indexes (see :gen:`index`) assigned to the placed items on the cable cell.
+        Apply a density mechanism on a region.
 
         :param str region: description of the region.
-        :param mechanism: the mechanism.
-        :type mechanism: :py:class:`mechanism`
+        :param density: the density mechanism.
+        :type density: :py:class:`density`
 
-    .. method:: paint(region, mech_name)
-        :noindex:
-
-        Apply a mechanism with a region using the name of the mechanism.
-        The mechanism will use the parameter values set in the mechanism catalogue.
-        Returns a unique identifier that can be used to query the local indexes (see :gen:`index`)
-        assigned to the placed items on the cable cell.
-
-        :param str region: description of the region.
-        :param str mechanism: the name of the mechanism.
-
-    .. method:: place(locations, mech_name, label)
-
-        Place one instance of the synapse named ``mech_name`` to each location in ``locations`` and label the
-        group of synapses with ``label``. The label can be used to form connections to one of the synapses
-        in the :py:class:`arbor.recipe` by creating a :py:class:`arbor.connection`.
-
-        :param str locations: description of the locset.
-        :param str mechanism: the name of the mechanism.
-        :param str label: the label of the group of synapses on the locset.
 
-    .. method:: place(locations, mechanism, label)
+    .. method:: place(locations, synapse, label)
         :noindex:
 
-        Place one instance of the synapse described by ``mechanism`` to each location in ``locations`` and label the
-        group of synapses with ``label``. The label can be used to form connections to one of the synapses
-        in the :py:class:`arbor.recipe` by creating a :py:class:`arbor.connection`.
+        Place one instance of the synapse mechanism described by ``synapse`` to each location in ``locations``
+        and label the group of synapses with ``label``. The label can be used to form connections to one of the
+        synapses in the :py:class:`arbor.recipe` by creating a :py:class:`arbor.connection`.
 
         :param str locations: description of the locset.
-        :param mechanism: the mechanism.
-        :type mechanism: :py:class:`mechanism`
+        :param synapse: the synapse.
+        :type synapse: :py:class:`synapse`
         :param str label: the label of the group of synapses on the locset.
 
-    .. method:: place(locations, site, label)
+    .. method:: place(locations, junction, label)
         :noindex:
 
-        Place one gap junction site at each location in ``locations`` and label the group of gap junction sites with
-        ``label``. The label can be used to form connections to/from one of the gap junction sites in the
-        :py:class:`arbor.recipe` by creating a :py:class:`arbor.gap_junction_connection`.
+        Place one instance of the gap junction mechanism described by ``junction`` at each location in ``locations``
+        and label the group of gap junction sites with ``label``. The label can be used to form gap junction
+        connections to/from one of labeled sites in the :py:class:`arbor.recipe` by creating a
+        :py:class:`arbor.gap_junction_connection`.
 
         :param str locations: description of the locset.
-        :param site: indicates a gap junction site..
-        :type site: :py:class:`gap_junction_site`
+        :param junction: the gap junction mechanism.
+        :type junction: :py:class:`junction`
         :param str label: the label of the group of gap junction sites on the locset.
 
     .. method:: place(locations, stim, label)
diff --git a/doc/python/interconnectivity.rst b/doc/python/interconnectivity.rst
index 7bb2e8b1..05d76e34 100644
--- a/doc/python/interconnectivity.rst
+++ b/doc/python/interconnectivity.rst
@@ -69,12 +69,11 @@ Interconnectivity
        A bidirectional gap-junction between two cells ``c0`` and ``c1`` requires two
        :class:`gap_junction_connection` objects to be constructed: one where ``c0`` is the
        :attr:`local` site, and ``c1`` is the :attr:`peer` site; and another where ``c1`` is the
-       :attr:`local` site, and ``c0`` is the :attr:`peer` site. If :attr:`ggap` is equal
-       in both connections, a symmetric gap-junction is formed, other wise the gap-junction is asymmetric.
+       :attr:`local` site, and ``c0`` is the :attr:`peer` site.
 
-    .. function::gap_junction_connection(peer, local, ggap)
+    .. function::gap_junction_connection(peer, local, weight)
 
-        Construct a gap junction connection between :attr:`peer` and :attr:`local` with conductance :attr:`ggap`.
+        Construct a gap junction connection between :attr:`peer` and :attr:`local` with weight :attr:`weight`.
 
     .. attribute:: peer
 
@@ -89,9 +88,9 @@ Interconnectivity
         the default :attr:`arbor.selection_policy.univalent` is used, or a (label, policy) tuple). The gid of the
         cell is implicitly known.
 
-    .. attribute:: ggap
+    .. attribute:: weight
 
-        The gap junction conductance [μS].
+        The unit-less weight of the gap junction connection.
 
 .. class:: spike_detector
 
diff --git a/doc/python/mechanisms.rst b/doc/python/mechanisms.rst
index 1a690de7..dbeb7903 100644
--- a/doc/python/mechanisms.rst
+++ b/doc/python/mechanisms.rst
@@ -3,9 +3,6 @@
 Cable cell mechanisms
 =====================
 
-When :ref:`decorating <cablecell-decoration>` a cable cell, we use a :py:class:`mechanism` type to describe a
-mechanism that is to be painted or placed on the cable cell.
-
 .. py:class:: mechanism
 
     Mechanisms describe physical processes, distributed over the membrane of the cell.
@@ -13,7 +10,9 @@ mechanism that is to be painted or placed on the cable cell.
     a function of the cell state and their own state where they are present.
     *Point mechanisms* are defined at discrete locations on the cell, which receive
     events from the network.
-    A third, specific type of density mechanism, which describes ionic reversal potential
+    *Junction mechanisms* are defined at discrete locations on the cell, which define the
+    behavior of a gap-junction mechanism.
+    A fourth, specific type of density mechanism, which describes ionic reversal potential
     behaviour, can be specified for cells or the whole model.
 
     The :class:`mechanism` type is a simple wrapper around a mechanism
@@ -37,44 +36,32 @@ mechanism that is to be painted or placed on the cable cell.
 
         import arbor
 
-        # hh dynamics with default parameters.
-        hh = arbor.mechanism('hh')
+        # A passive leaky channel with default parameter values (set in NOMDL file).
+        pas_0 = arbor.mechanism('pas')
+
+        # A passive leaky channel with custom conductance (range).
+        pas_1 = arbor.mechanism('pas', {'g': 0.02})
+
+        # A passive leaky channel with custom reversal potential (global).
+        pas_2 = arbor.mechanism('pas/e=-45')
+
+        # A passive leaky channel with custom reversal potential (global), and custom conductance (range).
+        pas_3 = arbor.mechanism('pas/e=-45', {'g', 0.1})
 
-        # A passive leaky channel with custom parameters
-        pas = arbor.mechanism('pas/e=-55.0', {'g': 0.02})
+        # This is an equivalent to pas_3, using set method to specify range parameters.
+        pas_4 = arbor.mechanism('pas/e=-45')
+        pas_4.set('g', 0.1)
 
         # Reversal potential using Nernst equation with GLOBAL parameter values
         # for Faraday's constant and the target ion species, set with a '/' followed
         # by comma-separated list of parameter after the base mechanism name.
         rev = arbor.mechanism('nernst/F=96485,x=ca')
 
-    Mechanisms can be painted to a region. Different mechanisms can be painted on top of each other.
-
-    .. code-block:: Python
-
-        import arbor
-
-        # Create pas mechanism with default parameter values (set in NOMDL file).
-        m1 = arbor.mechanism('pas')
-
-        # Create default mechainsm with custom conductance (range).
-        m2 = arbor.mechanism('pas', {'g', 0.1})
-
-        # Create a new pas mechanism with that changes reversal potential (global).
-        m3 = arbor.mechanism('pas/e=-45')
-
-        # Create an instance of the same mechanism, that also sets conductance (range).
-        m4 = arbor.mechanism('pas/e=-45', {'g', 0.1})
-
-        # This is an equivalent to m4, using set method to specify range parameters.
-        m5 = arbor.mechanism('pas/e=-45')
-        m5.set('g', 0.1)
+        # An exponential synapse with default parameter values (set in NOMDL file).
+        expsyn = arbor.mechanism("expsyn")
 
-        # Decorate the 'soma' with (multiple) mechanisms
-        decor.paint('"soma"', m1)
-        decor.paint('"soma"', m2) # Error: can't place the same mechanism on overlapping regions
-        decor.paint('"soma"', m3) # This is ok: m3 is a new, derived mechanism by virtue of
-                                  # having a different name, i.e. 'pas/e=-45' vs. 'pas'.
+        # A gap-junction mechanism with default parameter values (set in NOMDL file).
+        gj = arbor.mechanism("gj")
 
     .. method:: mechanism(name, params)
 
@@ -113,6 +100,156 @@ mechanism that is to be painted or placed on the cable cell.
 
        A dictionary of key-value pairs for the parameters.
 
+.. py:class:: density
+
+   When :ref:`decorating <cablecell-decoration>` a cable cell, we use a :py:class:`density` type to
+   wrap a density :py:class:`mechanism` that is to be painted on the cable cell.
+
+   Different :py:class:`density` mechanisms can be painted on top of each other.
+
+    .. code-block:: Python
+
+        import arbor
+
+        pas = arbor.mechanism('pas')
+        pas.set('g', 0.2)
+
+        # Decorate the 'soma' with (multiple) density mechanisms
+        decor.paint('"soma"', density(pas))
+        decor.paint('"soma"', density('pas', {'g': 0.1})) # Error: can't place the same mechanism on overlapping regions
+        decor.paint('"soma"', density('pas/e=-45'))       # This is ok: pas/e=-45 is a new, derived mechanism by virtue of
+                                                          # having a different name, i.e. 'pas/e=-45' vs. 'pas'.
+
+    .. py:attribute:: mech
+        :type: mechanism
+
+        The underlying mechanism.
+
+    .. method:: density(name)
+        :noindex:
+
+        constructs :attr:`mech` with *name* and default parameters.
+
+        :param name: name of mechanism.
+        :type name: str
+
+    .. method:: density(name, params)
+
+        constructs :attr:`mech` with *name* and range parameter overrides *params*.
+        for example: ``arbor.density('pas', {'g': 0.01})``.
+
+        :param name: name of mechanism.
+        :type name: str
+        :param params: A dictionary of parameter values, with parameter name as key.
+        :type params: dict[str, double]
+
+    .. method:: density(mech)
+        :noindex:
+
+        constructs :attr:`mech` from *mech*.
+
+        :param mech: mechanism description.
+        :type mech: :py:class:`mechanism`
+
+    .. method:: density(mech, params)
+
+        constructs :attr:`mech` from *mech* and sets the range parameter overrides *params*.
+
+        :param mech: mechanism description.
+        :type mech: :py:class:`mechanism`
+        :param params: A dictionary of parameter values, with parameter name as key.
+        :type params: dict[str, double]
+
+.. py:class:: synapse
+
+   When :ref:`decorating <cablecell-decoration>` a cable cell, we use a :py:class:`synapse` type to
+   wrap a point :py:class:`mechanism` that is to be placed on the cable cell.
+
+    .. py:attribute:: mech
+        :type: mechanism
+
+        The underlying mechanism.
+
+    .. method:: synapse(name)
+        :noindex:
+
+        constructs :attr:`mech` with *name* and default parameters.
+
+        :param name: name of mechanism.
+        :type name: str
+
+    .. method:: synapse(name, params)
+
+        constructs :attr:`mech` with *name* and range parameter overrides *params*.
+        for example: ``arbor.synapse('expsyn', {'tau': 0.01})``.
+
+        :param name: name of mechanism.
+        :type name: str
+        :param params: A dictionary of parameter values, with parameter name as key.
+        :type params: dict[str, double]
+
+    .. method:: synapse(mech)
+        :noindex:
+
+        constructs :attr:`mech` from *mech*.
+
+        :param mech: mechanism description.
+        :type mech: :py:class:`mechanism`
+
+    .. method:: synapse(mech, params)
+
+        constructs :attr:`mech` from *mech* and sets the range parameter overrides *params*.
+
+        :param mech: mechanism description.
+        :type mech: :py:class:`mechanism`
+        :param params: A dictionary of parameter values, with parameter name as key.
+        :type params: dict[str, double]
+
+
+.. py:class:: junction
+
+   When :ref:`decorating <cablecell-decoration>` a cable cell, we use a :py:class:`junction` type to
+   wrap a gap-junction :py:class:`mechanism` that is to be placed on the cable cell.
+
+    .. py:attribute:: mech
+        :type: mechanism
+
+        The underlying mechanism.
+
+    .. method:: junction(name)
+        :noindex:
+
+        constructs :attr:`mech` with *name* and default parameters.
+
+        :param name: name of mechanism.
+        :type name: str
+
+    .. method:: junction(name, params)
+
+        constructs :attr:`mech` with *name* and range parameter overrides *params*.
+        for example: ``arbor.junction('gj', {'g': 2})``.
+
+        :param name: name of mechanism.
+        :type name: str
+        :param params: A dictionary of parameter values, with parameter name as key.
+        :type params: dict[str, double]
+
+    .. method:: junction(mech)
+        :noindex:
+
+        constructs :attr:`mech` from *mech*.
+
+        :param mech: mechanism description.
+        :type mech: :py:class:`mechanism`
+
+    .. method:: junction(mech, params)
+
+        constructs :attr:`mech` from *mech* and sets the range parameter overrides *params*.
+
+        :param mech: mechanism description.
+        :type mech: :py:class:`mechanism`
+        :param params: A dictionary of parameter values, with parameter name as key.
+        :type params: dict[str, double]
 
 .. py:class:: mechanism_info
 
diff --git a/doc/tutorial/single_cell_detailed.rst b/doc/tutorial/single_cell_detailed.rst
index 7d24460b..a5db5411 100644
--- a/doc/tutorial/single_cell_detailed.rst
+++ b/doc/tutorial/single_cell_detailed.rst
@@ -249,9 +249,9 @@ We can override the default properties by *painting* new values on the relevant
    :lines: 59-61
 
 With the default and initial values taken care of, we now add some density mechanisms. Let's *paint*
-a *pas* mechanism everywhere on the cell using the previously defined "all" region; an *hh* mechanism
-on the "custom" region; and an *Ih* mechanism on the "dend" region. The *Ih* mechanism is explicitly
-constructed in order to change the default values of its 'gbar' parameter.
+a *pas* density mechanism everywhere on the cell using the previously defined "all" region; an *hh*
+density mechanism on the "custom" region; and an *Ih* density mechanism on the "dend" region. The *Ih*
+mechanism has a custom 'gbar' parameter.
 
 .. literalinclude:: ../../python/example/single_cell_detailed.py
    :language: python
diff --git a/doc/tutorial/single_cell_model.rst b/doc/tutorial/single_cell_model.rst
index 4ecd04f0..252c596a 100644
--- a/doc/tutorial/single_cell_model.rst
+++ b/doc/tutorial/single_cell_model.rst
@@ -68,8 +68,8 @@ following way:
   In the above example we set the initial membrane potential to -40 mV.
 * :meth:`arbor.decor.paint` is used to set properties or add dynamics to a region of the cell.
   We call this method 'painting' to convey that we are working on sections of a cell, as opposed to
-  precise locations: for example, we might want to *paint* an ion channel on all dendrites, and then
-  *place* a synapse at the tip of the axon. In the above example we paint
+  precise locations: for example, we might want to *paint* a density ion channel on all dendrites,
+  and then *place* a synapse at the tip of the axon. In the above example we paint
   HH dynamics on the region we previously named ``"soma"`` in our label dictionary.
 * :meth:`arbor.decor.place` is used to add objects on a precise
   :class:`arbor.location` on a cell. Examples of objects that are *placed* are synapses,
diff --git a/example/dryrun/branch_cell.hpp b/example/dryrun/branch_cell.hpp
index 0cae8cf0..9750c46c 100644
--- a/example/dryrun/branch_cell.hpp
+++ b/example/dryrun/branch_cell.hpp
@@ -112,8 +112,8 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param
 
     arb::decor decor;
 
-    decor.paint("soma"_lab, "hh");
-    decor.paint("dend"_lab, "pas");
+    decor.paint("soma"_lab, arb::density("hh"));
+    decor.paint("dend"_lab, arb::density("pas"));
 
     decor.set_default(arb::axial_resistivity{100}); // [Ω·cm]
 
@@ -121,11 +121,11 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param
     decor.place(arb::mlocation{0,0}, arb::threshold_detector{10}, "detector");
 
     // Add a synapse to the mid point of the first dendrite.
-    decor.place(arb::mlocation{0, 0.5}, "expsyn", "synapse");
+    decor.place(arb::mlocation{0, 0.5}, arb::synapse("expsyn"), "synapse");
 
     // Add additional synapses that will not be connected to anything.
     for (unsigned i=1u; i<params.synapses; ++i) {
-        decor.place(arb::mlocation{1, 0.5}, "expsyn", "dummy_synapses");
+        decor.place(arb::mlocation{1, 0.5}, arb::synapse("expsyn"), "dummy_synapses");
     }
 
     // Make a CV between every sample in the sample tree.
diff --git a/example/gap_junctions/gap_junctions.cpp b/example/gap_junctions/gap_junctions.cpp
index 28e4330b..c11802dc 100644
--- a/example/gap_junctions/gap_junctions.cpp
+++ b/example/gap_junctions/gap_junctions.cpp
@@ -303,17 +303,17 @@ arb::cable_cell gj_cell(cell_gid_type gid, unsigned ncell, double stim_duration)
     pas["g"] =  1.0/12000.0;
 
     // Paint density channels on all parts of the cell
-    decor.paint("(all)"_reg, nax);
-    decor.paint("(all)"_reg, kdrmt);
-    decor.paint("(all)"_reg, kamt);
-    decor.paint("(all)"_reg, pas);
+    decor.paint("(all)"_reg, arb::density{nax});
+    decor.paint("(all)"_reg, arb::density{kdrmt});
+    decor.paint("(all)"_reg, arb::density{kamt});
+    decor.paint("(all)"_reg, arb::density{pas});
 
     // Add a spike detector to the soma.
     decor.place(arb::mlocation{0,0}, arb::threshold_detector{10}, "detector");
 
     // Add two gap junction sites.
-    decor.place(arb::mlocation{0, 1}, arb::gap_junction_site{}, "local_0");
-    decor.place(arb::mlocation{0, 1}, arb::gap_junction_site{}, "local_1");
+    decor.place(arb::mlocation{0, 1}, arb::junction{"gj"}, "local_1");
+    decor.place(arb::mlocation{0, 0}, arb::junction{"gj"}, "local_0");
 
     // Attach a stimulus to the first cell of the first group
     if (!gid) {
@@ -322,7 +322,7 @@ arb::cable_cell gj_cell(cell_gid_type gid, unsigned ncell, double stim_duration)
     }
 
     // Add a synapse to the mid point of the first dendrite.
-    decor.place(arb::mlocation{0, 0.5}, "expsyn", "syn");
+    decor.place(arb::mlocation{0, 0.5}, arb::synapse{"expsyn"}, "syn");
 
     // Create the cell and set its electrical properties.
     return arb::cable_cell(tree, {}, decor);
diff --git a/example/generators/generators.cpp b/example/generators/generators.cpp
index 7fe3ea44..94bc75be 100644
--- a/example/generators/generators.cpp
+++ b/example/generators/generators.cpp
@@ -60,12 +60,12 @@ public:
         labels.set("soma", arb::reg::tagged(1));
 
         arb::decor decor;
-        decor.paint("soma"_lab, "pas");
+        decor.paint("soma"_lab, arb::density("pas"));
 
         // Add one synapse at the soma.
         // This synapse will be the target for all events, from both
         // event_generators.
-        decor.place(arb::mlocation{0, 0.5}, "expsyn", "syn");
+        decor.place(arb::mlocation{0, 0.5}, arb::synapse("expsyn"), "syn");
 
         return arb::cable_cell(tree, labels, decor);
     }
diff --git a/example/lfp/lfp.cpp b/example/lfp/lfp.cpp
index 61ddb27b..9669c26e 100644
--- a/example/lfp/lfp.cpp
+++ b/example/lfp/lfp.cpp
@@ -89,12 +89,12 @@ private:
         dec.set_default(cv_policy_fixed_per_branch(20, arb::reg::tagged(4)));
 
         // Add pas and hh mechanisms:
-        dec.paint(reg::tagged(1), "hh"); // (default parameters)
-        dec.paint(reg::tagged(4), mechanism_desc("pas/e=-70.0"));
+        dec.paint(reg::tagged(1), density("hh")); // (default parameters)
+        dec.paint(reg::tagged(4), density("pas/e=-70.0"));
 
         // Add exponential synapse at centre of soma.
         synapse_location_ = ls::on_components(0.5, reg::tagged(1));
-        dec.place(synapse_location_, mechanism_desc("expsyn").set("e", 0).set("tau", 2), "syn");
+        dec.place(synapse_location_, synapse("expsyn", {{"e", 0}, {"tau", 2}}), "syn");
 
         cell_ = cable_cell(tree, {}, dec);
     }
diff --git a/example/probe-demo/probe-demo.cpp b/example/probe-demo/probe-demo.cpp
index 6f67587a..585ea254 100644
--- a/example/probe-demo/probe-demo.cpp
+++ b/example/probe-demo/probe-demo.cpp
@@ -119,7 +119,7 @@ struct cable_recipe: public arb::recipe {
         tree.append(arb::mnpos, {0, 0, 0, 0.5*diam}, {length, 0, 0, 0.5*diam}, 1);
 
         arb::decor decor;
-        decor.paint(arb::reg::all(), "hh"); // HH mechanism over whole cell.
+        decor.paint(arb::reg::all(), arb::density("hh")); // HH mechanism over whole cell.
         decor.place(arb::mlocation{0, 0.}, arb::i_clamp{1.}, "iclamp"); // Inject a 1 nA current indefinitely.
 
         return arb::cable_cell(tree, {}, decor);
diff --git a/example/ring/branch_cell.hpp b/example/ring/branch_cell.hpp
index 789684ad..527f2601 100644
--- a/example/ring/branch_cell.hpp
+++ b/example/ring/branch_cell.hpp
@@ -112,8 +112,8 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param
 
     arb::decor decor;
 
-    decor.paint("soma"_lab, "hh");
-    decor.paint("dend"_lab, "pas");
+    decor.paint("soma"_lab, arb::density("hh"));
+    decor.paint("dend"_lab, arb::density("pas"));
 
     decor.set_default(arb::axial_resistivity{100}); // [Ω·cm]
 
@@ -121,12 +121,12 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param
     decor.place(arb::mlocation{0,0}, arb::threshold_detector{10}, "detector");
 
     // Add a synapse to the mid point of the first dendrite.
-    decor.place(arb::mlocation{0, 0.5}, "expsyn", "primary_syn");
+    decor.place(arb::mlocation{0, 0.5}, arb::synapse("expsyn"), "primary_syn");
 
     // Add additional synapses that will not be connected to anything.
 
     if (params.synapses > 1) {
-        decor.place(arb::ls::uniform("dend"_lab, 0, params.synapses - 2, gid), "expsyn", "extra_syns");
+        decor.place(arb::ls::uniform("dend"_lab, 0, params.synapses - 2, gid), arb::synapse("expsyn"), "extra_syns");
     }
 
     // Make a CV between every sample in the sample tree.
diff --git a/example/single/single.cpp b/example/single/single.cpp
index b7d3cff2..2721d2fd 100644
--- a/example/single/single.cpp
+++ b/example/single/single.cpp
@@ -68,14 +68,14 @@ struct single_recipe: public arb::recipe {
         arb::decor decor;
 
         // Add HH mechanism to soma, passive channels to dendrites.
-        decor.paint("soma"_lab, "hh");
-        decor.paint("dend"_lab, "pas");
+        decor.paint("soma"_lab, arb::density("hh"));
+        decor.paint("dend"_lab, arb::density("pas"));
 
         // Add synapse to last branch.
 
         arb::cell_lid_type last_branch = morpho.num_branches()-1;
         arb::mlocation end_last_branch = { last_branch, 1. };
-        decor.place(end_last_branch, "exp2syn", "synapse");
+        decor.place(end_last_branch, arb::synapse("exp2syn"), "synapse");
 
         return arb::cable_cell(morpho, dict, decor);
     }
diff --git a/mechanisms/CMakeLists.txt b/mechanisms/CMakeLists.txt
index f696245d..ac69ce94 100644
--- a/mechanisms/CMakeLists.txt
+++ b/mechanisms/CMakeLists.txt
@@ -25,7 +25,7 @@ make_catalogue(
   NAME default
   SOURCES "${CMAKE_CURRENT_SOURCE_DIR}/default"
   OUTPUT "CAT_DEFAULT_SOURCES"
-  MECHS exp2syn expsyn expsyn_stdp hh kamt kdrmt nax nernst pas
+  MECHS exp2syn expsyn expsyn_stdp hh kamt kdrmt nax nernst pas gj
   PREFIX "${PROJECT_SOURCE_DIR}/mechanisms"
   CXX_FLAGS_TARGET "${ARB_CXX_FLAGS_TARGET_FULL}"
   STANDALONE FALSE
diff --git a/mechanisms/default/gj.mod b/mechanisms/default/gj.mod
new file mode 100644
index 00000000..00d72cf8
--- /dev/null
+++ b/mechanisms/default/gj.mod
@@ -0,0 +1,14 @@
+NEURON {
+    JUNCTION_PROCESS gj
+    NONSPECIFIC_CURRENT i
+    RANGE g
+}
+INITIAL {}
+
+PARAMETER {
+    g = 1
+}
+
+BREAKPOINT {
+    i = g*(v - v_peer)
+}
diff --git a/modcc/identifier.hpp b/modcc/identifier.hpp
index 610a9b65..dd31222f 100644
--- a/modcc/identifier.hpp
+++ b/modcc/identifier.hpp
@@ -5,7 +5,7 @@
 #include <stdexcept>
 
 enum class moduleKind {
-    point, density, revpot
+    point, density, junction, revpot
 };
 
 /// indicate how a variable is accessed
@@ -41,6 +41,7 @@ enum class linkageKind {
 /// possible external data source for indexed variables
 enum class sourceKind {
     voltage,
+    peer_voltage,
     current_density,
     current,
     conductivity,
@@ -86,6 +87,7 @@ inline std::string to_string(linkageKind v) {
 inline std::string to_string(sourceKind v) {
     switch(v) {
     case sourceKind::voltage:             return "voltage";
+    case sourceKind::peer_voltage:        return "peer_voltage";
     case sourceKind::current_density:     return "current_density";
     case sourceKind::current:             return "current";
     case sourceKind::conductivity:        return "conductivity";
@@ -117,7 +119,7 @@ inline std::ostream& operator<< (std::ostream& os, linkageKind l) {
 
 inline sourceKind ion_source(const std::string& ion, const std::string& var, moduleKind mkind) {
     if (ion.empty()) return sourceKind::no_source;
-    else if (var=="i"+ion) return mkind==moduleKind::point? sourceKind::ion_current: sourceKind::ion_current_density;
+    else if (var=="i"+ion) return mkind==moduleKind::density? sourceKind::ion_current_density: sourceKind::ion_current;
     else if (var=="e"+ion) return sourceKind::ion_revpot;
     else if (var==ion+"i") return sourceKind::ion_iconc;
     else if (var==ion+"o") return sourceKind::ion_econc;
diff --git a/modcc/module.cpp b/modcc/module.cpp
index dc2c5b9d..5e46c307 100644
--- a/modcc/module.cpp
+++ b/modcc/module.cpp
@@ -586,13 +586,14 @@ void Module::add_variables_to_symbols() {
             make_symbol<IndexedVariable>(loc, name, data_source, acc, ch);
     };
 
-    sourceKind current_kind = kind_==moduleKind::point? sourceKind::current: sourceKind::current_density;
-    sourceKind conductance_kind = kind_==moduleKind::point? sourceKind::conductance: sourceKind::conductivity;
+    sourceKind current_kind = kind_==moduleKind::density? sourceKind::current_density: sourceKind::current;
+    sourceKind conductance_kind = kind_==moduleKind::density? sourceKind::conductivity: sourceKind::conductance;
 
     create_indexed_variable("current_", current_kind, accessKind::write, "", Location());
     create_indexed_variable("conductivity_", conductance_kind, accessKind::write, "", Location());
-    create_indexed_variable("v", sourceKind::voltage, accessKind::read,  "", Location());
-    create_indexed_variable("dt", sourceKind::dt, accessKind::read,  "", Location());
+    create_indexed_variable("v",      sourceKind::voltage, accessKind::read,  "", Location());
+    create_indexed_variable("v_peer", sourceKind::peer_voltage, accessKind::read,  "", Location());
+    create_indexed_variable("dt",     sourceKind::dt, accessKind::read,  "", Location());
 
     // If we put back support for accessing cell time again from NMODL code,
     // add indexed_variable also for "time" with appropriate cell-index based
@@ -604,9 +605,9 @@ void Module::add_variables_to_symbols() {
             accessKind::readwrite, visibilityKind::local, linkageKind::local, rangeKind::range, true);
     }
 
-    // Add parameters, ignoring built-in voltage variable "v".
+    // Add parameters, ignoring built-in voltage variables "v" and "v_peer".
     for (const Id& id: parameter_block_) {
-        if (id.name() == "v") {
+        if (id.name() == "v" || id.name() == "v_peer") {
             continue;
         }
 
@@ -653,9 +654,9 @@ void Module::add_variables_to_symbols() {
         parameter_block_.end()
     );
 
-    // Add 'assigned' variables, ignoring built-in voltage variable "v".
+    // Add 'assigned' variables, ignoring built-in voltage variables "v" and "v_peer".
     for (const Id& id: assigned_block_) {
-        if (id.name() == "v") {
+        if (id.name() == "v" || id.name() == "v_peer") {
             continue;
         }
 
diff --git a/modcc/parser.cpp b/modcc/parser.cpp
index 4b80992c..75bf1a57 100644
--- a/modcc/parser.cpp
+++ b/modcc/parser.cpp
@@ -233,8 +233,9 @@ void Parser::parse_neuron_block() {
 
         case tok::suffix:
         case tok::point_process:
-            neuron_block.kind = (token_.type == tok::suffix) ? moduleKind::density
-                                                             : moduleKind::point;
+        case tok::junction_process:
+            neuron_block.kind = (token_.type == tok::suffix) ? moduleKind::density :
+                                (token_.type == tok::point_process) ? moduleKind::point : moduleKind::junction;
 
             // set the modul kind
             module_->kind(neuron_block.kind);
diff --git a/modcc/printer/cprinter.cpp b/modcc/printer/cprinter.cpp
index 74d92372..23c4793f 100644
--- a/modcc/printer/cprinter.cpp
+++ b/modcc/printer/cprinter.cpp
@@ -280,6 +280,7 @@ std::string emit_cpp_source(const Module& module_, const printer_options& opt) {
                                    "[[maybe_unused]] auto* {0}diam_um           = pp->diam_um;\\\n"
                                    "[[maybe_unused]] auto* {0}time_since_spike  = pp->time_since_spike;\\\n"
                                    "[[maybe_unused]] auto* {0}node_index        = pp->node_index;\\\n"
+                                   "[[maybe_unused]] auto* {0}peer_index        = pp->peer_index;\\\n"
                                    "[[maybe_unused]] auto* {0}multiplicity      = pp->multiplicity;\\\n"
                                    "[[maybe_unused]] auto* {0}weight            = pp->weight;\\\n"
                                    "[[maybe_unused]] auto& {0}events            = pp->events;\\\n"
diff --git a/modcc/printer/gpuprinter.cpp b/modcc/printer/gpuprinter.cpp
index 73fda5a1..8c6bfc0c 100644
--- a/modcc/printer/gpuprinter.cpp
+++ b/modcc/printer/gpuprinter.cpp
@@ -99,7 +99,7 @@ std::string emit_gpu_cu_source(const Module& module_, const printer_options& opt
 
     auto ns_components = namespace_components(opt.cpp_namespace);
 
-    const bool is_point_proc = module_.kind() == moduleKind::point;
+    const bool is_point_proc = (module_.kind() == moduleKind::point) || (module_.kind() == moduleKind::junction);
 
     APIMethod* net_receive_api = find_api_method(module_, "net_rec_api");
     APIMethod* post_event_api  = find_api_method(module_, "post_event_api");
@@ -137,6 +137,7 @@ std::string emit_gpu_cu_source(const Module& module_, const printer_options& opt
                                    "auto* {0}diam_um           __attribute__((unused)) = params_.diam_um;\\\n"
                                    "auto* {0}time_since_spike  __attribute__((unused)) = params_.time_since_spike;\\\n"
                                    "auto* {0}node_index        __attribute__((unused)) = params_.node_index;\\\n"
+                                   "auto* {0}peer_index        __attribute__((unused)) = params_.peer_index;\\\n"
                                    "auto* {0}multiplicity      __attribute__((unused)) = params_.multiplicity;\\\n"
                                    "auto* {0}state_vars        __attribute__((unused)) = params_.state_vars;\\\n"
                                    "auto* {0}weight            __attribute__((unused)) = params_.weight;\\\n"
diff --git a/modcc/printer/printerutil.cpp b/modcc/printer/printerutil.cpp
index d36e9e6d..47d146a1 100644
--- a/modcc/printer/printerutil.cpp
+++ b/modcc/printer/printerutil.cpp
@@ -137,6 +137,11 @@ indexed_variable_info decode_indexed_variable(IndexedVariable* sym) {
         v.data_var="vec_v";
         v.readonly = true;
         break;
+    case sourceKind::peer_voltage:
+        v.data_var="vec_v";
+        v.node_index_var = "peer_index";
+        v.readonly = true;
+        break;
     case sourceKind::current_density:
         v.data_var = "vec_i";
         v.readonly = false;
diff --git a/modcc/printer/printerutil.hpp b/modcc/printer/printerutil.hpp
index 053c7fbc..cf872106 100644
--- a/modcc/printer/printerutil.hpp
+++ b/modcc/printer/printerutil.hpp
@@ -58,9 +58,10 @@ struct namespace_declaration_close {
 
 inline const char* module_kind_str(const Module& m) {
     switch (m.kind()) {
-    case moduleKind::density: return "arb_mechanism_kind_density";            break;
-    case moduleKind::point:   return "arb_mechanism_kind_point";              break;
-    case moduleKind::revpot:  return "arb_mechanism_kind_reversal_potential"; break;
+    case moduleKind::density:   return "arb_mechanism_kind_density";
+    case moduleKind::point:     return "arb_mechanism_kind_point";
+    case moduleKind::revpot:    return "arb_mechanism_kind_reversal_potential";
+    case moduleKind::junction:  return "arb_mechanism_kind_gap_junction";
     default: throw compiler_exception("Unknown module kind " + std::to_string((int)m.kind()));
     }
 }
diff --git a/modcc/token.cpp b/modcc/token.cpp
index e2226c7e..e884bc0d 100644
--- a/modcc/token.cpp
+++ b/modcc/token.cpp
@@ -51,7 +51,8 @@ static Keyword keywords[] = {
     {"SOLVE",       tok::solve},
     {"THREADSAFE",  tok::threadsafe},
     {"GLOBAL",      tok::global},
-    {"POINT_PROCESS", tok::point_process},
+    {"POINT_PROCESS",    tok::point_process},
+    {"JUNCTION_PROCESS", tok::junction_process},
     {"COMPARTMENT", tok::compartment},
     {"METHOD",      tok::method},
     {"STEADYSTATE", tok::steadystate},
@@ -132,7 +133,8 @@ static TokenString token_strings[] = {
     {"SOLVE",       tok::solve},
     {"THREADSAFE",  tok::threadsafe},
     {"GLOBAL",      tok::global},
-    {"POINT_PROCESS", tok::point_process},
+    {"POINT_PROCESS",    tok::point_process},
+    {"JUNCTION_PROCESS", tok::junction_process},
     {"COMPARTMENT", tok::compartment},
     {"METHOD",      tok::method},
     {"STEADYSTATE", tok::steadystate},
diff --git a/modcc/token.hpp b/modcc/token.hpp
index e3a7a12c..5e6466fe 100644
--- a/modcc/token.hpp
+++ b/modcc/token.hpp
@@ -63,7 +63,7 @@ enum class tok {
     range, local, conserve, compartment,
     solve, method, steadystate,
     threadsafe, global,
-    point_process,
+    point_process, junction_process,
     from, to,
 
     // prefix binary operators
diff --git a/python/cells.cpp b/python/cells.cpp
index 3bb1ea2b..1c26f783 100644
--- a/python/cells.cpp
+++ b/python/cells.cpp
@@ -345,14 +345,41 @@ void register_cells(pybind11::module& m) {
           "domain"_a="(all)", "the domain to which the policy is to be applied",
           "Policy to use the same number of CVs for each branch.");
 
-    // arb::gap_junction_site
-
-    pybind11::class_<arb::gap_junction_site> gjsite(m, "gap_junction",
-            "For marking a location on a cell morphology as a gap junction site.");
-    gjsite
-        .def(pybind11::init<>())
-        .def("__repr__", [](const arb::gap_junction_site&){return "<arbor.gap_junction>";})
-        .def("__str__", [](const arb::gap_junction_site&){return "<arbor.gap_junction>";});
+    // arb::density
+
+    pybind11::class_<arb::density> density(m, "density", "For painting a density mechanism on a region.");
+    density
+        .def(pybind11::init([](const std::string& name) {return arb::density(name);}))
+        .def(pybind11::init([](arb::mechanism_desc mech) {return arb::density(mech);}))
+        .def(pybind11::init([](const std::string& name, const std::unordered_map<std::string, double>& params) {return arb::density(name, params);}))
+        .def(pybind11::init([](arb::mechanism_desc mech, const std::unordered_map<std::string, double>& params) {return arb::density(mech, params);}))
+        .def_readonly("mech", &arb::density::mech, "The underlying mechanism.")
+        .def("__repr__", [](const arb::density& d){return "<arbor.density " + mechanism_desc_str(d.mech) + ">";})
+        .def("__str__", [](const arb::density& d){return "<arbor.density " + mechanism_desc_str(d.mech) + ">";});
+
+    // arb::synapse
+
+    pybind11::class_<arb::synapse> synapse(m, "synapse", "For placing a synaptic mechanism on a locset.");
+    synapse
+        .def(pybind11::init([](const std::string& name) {return arb::synapse(name);}))
+        .def(pybind11::init([](arb::mechanism_desc mech) {return arb::synapse(mech);}))
+        .def(pybind11::init([](const std::string& name, const std::unordered_map<std::string, double>& params) {return arb::synapse(name, params);}))
+        .def(pybind11::init([](arb::mechanism_desc mech, const std::unordered_map<std::string, double>& params) {return arb::synapse(mech, params);}))
+        .def_readonly("mech", &arb::synapse::mech, "The underlying mechanism.")
+        .def("__repr__", [](const arb::synapse& s){return "<arbor.synapse " + mechanism_desc_str(s.mech) + ">";})
+        .def("__str__", [](const arb::synapse& s){return "<arbor.synapse " + mechanism_desc_str(s.mech) + ">";});
+
+    // arb::junction
+
+    pybind11::class_<arb::junction> junction(m, "junction", "For placing a gap-junction mechanism on a locset.");
+    junction
+        .def(pybind11::init([](const std::string& name) {return arb::junction(name);}))
+        .def(pybind11::init([](arb::mechanism_desc mech) {return arb::junction(mech);}))
+        .def(pybind11::init([](const std::string& name, const std::unordered_map<std::string, double>& params) {return arb::junction(name, params);}))
+        .def(pybind11::init([](arb::mechanism_desc mech, const std::unordered_map<std::string, double>& params) {return arb::junction(mech, params);}))
+        .def_readonly("mech", &arb::junction::mech, "The underlying mechanism.")
+        .def("__repr__", [](const arb::junction& j){return "<arbor.junction " + mechanism_desc_str(j.mech) + ">";})
+        .def("__str__", [](const arb::junction& j){return "<arbor.junction " + mechanism_desc_str(j.mech) + ">";});
 
     // arb::i_clamp
 
@@ -531,17 +558,11 @@ void register_cells(pybind11::module& m) {
             "compartments in the cell, and can't be overriden locally.")
         // Paint mechanisms.
         .def("paint",
-            [](arb::decor& dec, const char* region, const arb::mechanism_desc& d) {
-                dec.paint(arborio::parse_region_expression(region).unwrap(), d);
+            [](arb::decor& dec, const char* region, const arb::density& mechanism) {
+                dec.paint(arborio::parse_region_expression(region).unwrap(), mechanism);
             },
             "region"_a, "mechanism"_a,
-            "Associate a mechanism with a region.")
-        .def("paint",
-            [](arb::decor& dec, const char* region, const char* mech_name) {
-                dec.paint(arborio::parse_region_expression(region).unwrap(), arb::mechanism_desc(mech_name));
-            },
-            "region"_a, "mechanism"_a,
-            "Associate a mechanism with a region.")
+            "Associate a density mechanism with a region.")
         // Paint membrane/static properties.
         .def("paint",
             [](arb::decor& dec,
@@ -570,33 +591,27 @@ void register_cells(pybind11::module& m) {
                 if (ext_con) dec.paint(r, arb::init_ext_concentration{name, *ext_con});
                 if (rev_pot) dec.paint(r, arb::init_reversal_potential{name, *rev_pot});
             },
-            "region"_a, "ion_name"_a,
+            "region"_a, pybind11::kw_only(), "ion_name"_a,
             pybind11::arg_v("int_con", pybind11::none(), "Initial internal concentration [mM]"),
             pybind11::arg_v("ext_con", pybind11::none(), "Initial external concentration [mM]"),
             pybind11::arg_v("rev_pot", pybind11::none(), "Initial reversal potential [mV]"),
             "Set ion species properties conditions on a region.")
         // Place synapses
         .def("place",
-            [](arb::decor& dec, const char* locset, const arb::mechanism_desc& d, const char* label_name) {
-                return dec.place(arborio::parse_locset_expression(locset).unwrap(), d, label_name); },
-            "locations"_a, "mechanism"_a, "label"_a,
-            "Place one instance of the synapse described by 'mechanism' on each location in 'locations'. "
-            "The group of synapses has the label 'label', used for forming connections between cells.")
-        .def("place",
-            [](arb::decor& dec, const char* locset, const char* mech_name, const char* label_name) {
-                return dec.place(arborio::parse_locset_expression(locset).unwrap(), mech_name, label_name);
+            [](arb::decor& dec, const char* locset, const arb::synapse& mechanism, const char* label_name) {
+                return dec.place(arborio::parse_locset_expression(locset).unwrap(), mechanism, label_name);
             },
-            "locations"_a, "mechanism"_a, "label"_a,
-            "Place one instance of the synapse described by 'mechanism' on each location in 'locations'."
+            "locations"_a, "synapse"_a, "label"_a,
+            "Place one instance of 'synapse' on each location in 'locations'."
             "The group of synapses has the label 'label', used for forming connections between cells.")
         // Place gap junctions.
         .def("place",
-            [](arb::decor& dec, const char* locset, const arb::gap_junction_site& site, const char* label_name) {
-                return dec.place(arborio::parse_locset_expression(locset).unwrap(), site, label_name);
+            [](arb::decor& dec, const char* locset, const arb::junction& mechanism, const char* label_name) {
+                return dec.place(arborio::parse_locset_expression(locset).unwrap(), mechanism, label_name);
             },
-            "locations"_a, "gapjunction"_a, "label"_a,
-            "Place one gap junction site labeled 'label' on each location in 'locations'."
-            "The group of gap junctions has the label 'label', used for forming connections between cells.")
+            "locations"_a, "junction"_a, "label"_a,
+            "Place one instance of 'junction' on each location in 'locations'."
+            "The group of junctions has the label 'label', used for forming gap-junction connections between cells.")
         // Place current clamp stimulus.
         .def("place",
             [](arb::decor& dec, const char* locset, const arb::i_clamp& stim, const char* label_name) {
diff --git a/python/example/network_ring.py b/python/example/network_ring.py
index 345238a8..691c05ad 100755
--- a/python/example/network_ring.py
+++ b/python/example/network_ring.py
@@ -45,11 +45,11 @@ def make_cable_cell(gid):
     decor = arbor.decor()
 
     # Put hh dynamics on soma, and passive properties on the dendrites.
-    decor.paint('"soma"', 'hh')
-    decor.paint('"dend"', 'pas')
+    decor.paint('"soma"', arbor.density('hh'))
+    decor.paint('"dend"', arbor.density('pas'))
 
     # (4) Attach a single synapse.
-    decor.place('"synapse_site"', 'expsyn', 'syn')
+    decor.place('"synapse_site"', arbor.synapse('expsyn'), 'syn')
 
     # Attach a spike detector with threshold of -10 mV.
     decor.place('"root"', arbor.spike_detector(-10), 'detector')
diff --git a/python/example/network_ring_mpi.py b/python/example/network_ring_mpi.py
index 6a38e2c6..c2deb198 100644
--- a/python/example/network_ring_mpi.py
+++ b/python/example/network_ring_mpi.py
@@ -47,11 +47,11 @@ def make_cable_cell(gid):
     decor = arbor.decor()
 
     # Put hh dynamics on soma, and passive properties on the dendrites.
-    decor.paint('"soma"', 'hh')
-    decor.paint('"dend"', 'pas')
+    decor.paint('"soma"', arbor.density('hh'))
+    decor.paint('"dend"', arbor.density('pas'))
 
     # (4) Attach a single synapse.
-    decor.place('"synapse_site"', 'expsyn', 'syn')
+    decor.place('"synapse_site"', arbor.synapse('expsyn'), 'syn')
 
     # Attach a spike detector with threshold of -10 mV.
     decor.place('"root"', arbor.spike_detector(-10), 'detector')
diff --git a/python/example/single_cell_cable.py b/python/example/single_cell_cable.py
index b44c5662..f22d26f0 100755
--- a/python/example/single_cell_cable.py
+++ b/python/example/single_cell_cable.py
@@ -88,7 +88,7 @@ class Cable(arbor.recipe):
         decor.set_property(rL=self.rL)
 
         decor.paint('"cable"',
-                    arbor.mechanism(f'pas/e={self.Vm}', {'g': self.g}))
+                    arbor.density(f'pas/e={self.Vm}', {'g': self.g}))
 
         decor.place('"start"', arbor.iclamp(self.stimulus_start, self.stimulus_duration, self.stimulus_amplitude), "iclamp")
 
diff --git a/python/example/single_cell_detailed.py b/python/example/single_cell_detailed.py
index 0ff3af50..0bd1d3f2 100755
--- a/python/example/single_cell_detailed.py
+++ b/python/example/single_cell_detailed.py
@@ -5,7 +5,7 @@ import arbor
 import pandas
 import seaborn
 import sys
-from arbor import mechanism as mech
+from arbor import density
 
 # (1) Read the morphology from an SWC file.
 
@@ -61,9 +61,9 @@ decor.paint('"custom"', tempK=270)
 decor.paint('"soma"',   Vm=-50)
 
 # Paint density mechanisms.
-decor.paint('"all"', 'pas')
-decor.paint('"custom"', 'hh')
-decor.paint('"dend"',  mech('Ih', {'gbar': 0.001}))
+decor.paint('"all"', density('pas'))
+decor.paint('"custom"', density('hh'))
+decor.paint('"dend"',  density('Ih', {'gbar': 0.001}))
 
 # Place stimuli and spike detectors.
 decor.place('"root"', arbor.iclamp(10, 1, current=2), 'iclamp0')
diff --git a/python/example/single_cell_detailed_recipe.py b/python/example/single_cell_detailed_recipe.py
index d9275e8d..ee0f2857 100644
--- a/python/example/single_cell_detailed_recipe.py
+++ b/python/example/single_cell_detailed_recipe.py
@@ -5,7 +5,7 @@ import arbor
 import pandas
 import seaborn
 import sys
-from arbor import mechanism as mech
+from arbor import density
 
 # (1) Read the morphology from an SWC file.
 
@@ -61,9 +61,9 @@ decor.paint('"custom"', tempK=270)
 decor.paint('"soma"',   Vm=-50)
 
 # Paint density mechanisms.
-decor.paint('"all"', 'pas')
-decor.paint('"custom"', 'hh')
-decor.paint('"dend"',  mech('Ih', {'gbar': 0.001}))
+decor.paint('"all"', density('pas'))
+decor.paint('"custom"', density('hh'))
+decor.paint('"dend"', density('Ih', {'gbar': 0.001}))
 
 # Place stimuli and spike detectors.
 decor.place('"root"', arbor.iclamp(10, 1, current=2), 'iclamp0')
diff --git a/python/example/single_cell_extracellular_potentials.py b/python/example/single_cell_extracellular_potentials.py
index 21508c41..adf3d629 100644
--- a/python/example/single_cell_extracellular_potentials.py
+++ b/python/example/single_cell_extracellular_potentials.py
@@ -92,7 +92,7 @@ decor.set_property(
 # passive mech w. leak reversal potential (mV)
 pas = arbor.mechanism('pas/e=-65')
 pas.set('g', 0.0001)  # leak conductivity (S/cm2)
-decor.paint('(all)', pas)
+decor.paint('(all)', arbor.density(pas))
 
 # set sinusoid input current at mid point of terminating CV (segment)
 iclamp = arbor.iclamp(5,  # stimulation onset (ms)
diff --git a/python/example/single_cell_model.py b/python/example/single_cell_model.py
index a21445d3..2410971b 100755
--- a/python/example/single_cell_model.py
+++ b/python/example/single_cell_model.py
@@ -15,7 +15,7 @@ labels = arbor.label_dict({'soma':   '(tag 1)',
 # (3) Create and set up a decor object
 decor = arbor.decor()
 decor.set_property(Vm=-40)
-decor.paint('"soma"', 'hh')
+decor.paint('"soma"', arbor.density('hh'))
 decor.place('"midpoint"', arbor.iclamp( 10, 2, 0.8), "iclamp")
 decor.place('"midpoint"', arbor.spike_detector(-10), "detector")
 
diff --git a/python/example/single_cell_nml.py b/python/example/single_cell_nml.py
index d582b36b..008a0b4f 100755
--- a/python/example/single_cell_nml.py
+++ b/python/example/single_cell_nml.py
@@ -50,10 +50,10 @@ decor.set_property(Vm=-55)
 decor.set_ion('ca', method=mech('nernst/x=ca'))
 #decor.set_ion('ca', method='nernst/x=ca')
 # hh mechanism on the soma and axon.
-decor.paint('"soma"', 'hh')
-decor.paint('"axon"', 'hh')
+decor.paint('"soma"', arbor.density('hh'))
+decor.paint('"axon"', arbor.density('hh'))
 # pas mechanism the dendrites.
-decor.paint('"dend"', 'pas')
+decor.paint('"dend"', arbor.density('pas'))
 # Increase resistivity on dendrites.
 decor.paint('"dend"', rL=500)
 # Attach stimuli that inject 4 nA current for 1 ms, starting at 3 and 8 ms.
@@ -74,7 +74,7 @@ decor.discretization(policy)
 # Combine morphology with region and locset definitions to make a cable cell.
 cell = arbor.cable_cell(morpho, labels, decor)
 
-print(cell.locations('"axon_end"'))
+print(cell.locations('axon_end'))
 
 # Make single cell model.
 m = arbor.single_cell_model(cell)
diff --git a/python/example/single_cell_recipe.py b/python/example/single_cell_recipe.py
index 193e258a..8ef071c8 100644
--- a/python/example/single_cell_recipe.py
+++ b/python/example/single_cell_recipe.py
@@ -17,7 +17,7 @@ labels = arbor.label_dict({'soma':   '(tag 1)',
 # (3) Create cell and set properties
 decor = arbor.decor()
 decor.set_property(Vm=-40)
-decor.paint('"soma"', 'hh')
+decor.paint('"soma"', arbor.density('hh'))
 decor.place('"midpoint"', arbor.iclamp( 10, 2, 0.8), "iclamp")
 decor.place('"midpoint"', arbor.spike_detector(-10), "detector")
 cell = arbor.cable_cell(tree, labels, decor)
diff --git a/python/example/single_cell_stdp.py b/python/example/single_cell_stdp.py
index b2db196e..761a331e 100755
--- a/python/example/single_cell_stdp.py
+++ b/python/example/single_cell_stdp.py
@@ -32,15 +32,16 @@ class single_recipe(arbor.recipe):
 
         decor = arbor.decor()
         decor.set_property(Vm=-40)
-        decor.paint('(all)', 'hh')
+        decor.paint('(all)', arbor.density('hh'))
 
         decor.place('"center"', arbor.spike_detector(-10), "detector")
-        decor.place('"center"', 'expsyn', "synapse")
+        decor.place('"center"', arbor.synapse('expsyn'), "synapse")
 
-        mech_syn = arbor.mechanism('expsyn_stdp')
-        mech_syn.set("max_weight", 1.)
+        mech = arbor.mechanism('expsyn_stdp')
+        mech.set("max_weight", 1.)
+        syn = arbor.synapse(mech)
 
-        decor.place('"center"', mech_syn, "stpd_synapse")
+        decor.place('"center"', syn, "stpd_synapse")
 
         cell = arbor.cable_cell(tree, labels, decor)
 
diff --git a/python/example/single_cell_swc.py b/python/example/single_cell_swc.py
index b6e91ce7..dbeacbfa 100755
--- a/python/example/single_cell_swc.py
+++ b/python/example/single_cell_swc.py
@@ -42,10 +42,10 @@ decor.set_property(Vm=-55)
 decor.set_ion('ca', method=mech('nernst/x=ca'))
 #decor.set_ion('ca', method='nernst/x=ca')
 # hh mechanism on the soma and axon.
-decor.paint('"soma"', 'hh')
-decor.paint('"axon"', 'hh')
+decor.paint('"soma"', arbor.density('hh'))
+decor.paint('"axon"', arbor.density('hh'))
 # pas mechanism the dendrites.
-decor.paint('"dend"', 'pas')
+decor.paint('"dend"', arbor.density('pas'))
 # Increase resistivity on dendrites.
 decor.paint('"dend"', rL=500)
 # Attach stimuli that inject 4 nA current for 1 ms, starting at 3 and 8 ms.
diff --git a/python/recipe.cpp b/python/recipe.cpp
index b5d23897..de2cb132 100644
--- a/python/recipe.cpp
+++ b/python/recipe.cpp
@@ -123,8 +123,8 @@ std::string con_to_string(const arb::cell_connection& c) {
 }
 
 std::string gj_to_string(const arb::gap_junction_connection& gc) {
-    return util::pprintf("<arbor.gap_junction_connection: peer ({}, \"{}\", {}), local (\"{}\", {}), ggap {}>",
-         gc.peer.gid, gc.peer.label.tag, gc.peer.label.policy, gc.local.tag, gc.local.policy, gc.ggap);
+    return util::pprintf("<arbor.gap_junction_connection: peer ({}, \"{}\", {}), local (\"{}\", {}), weight {}>",
+         gc.peer.gid, gc.peer.label.tag, gc.peer.label.policy, gc.local.tag, gc.local.policy, gc.weight);
 }
 
 void register_recipe(pybind11::module& m) {
@@ -158,17 +158,17 @@ void register_recipe(pybind11::module& m) {
         "Describes a gap junction between two gap junction sites.");
     gap_junction_connection
         .def(pybind11::init<arb::cell_global_label_type, arb::cell_local_label_type, double>(),
-            "peer"_a, "local"_a, "ggap"_a,
+            "peer"_a, "local"_a, "weight"_a,
             "Construct a gap junction connection with arguments:\n"
             "  peer:  remote half of the gap junction connection.\n"
             "  local: local half of the gap junction connection.\n"
-            "  ggap:  Gap junction conductance [μS].")
+            "  weight:  Gap junction connection weight [unit-less].")
         .def_readwrite("peer", &arb::gap_junction_connection::peer,
             "Remote gid and label of the gap junction connection.")
         .def_readwrite("local", &arb::gap_junction_connection::local,
             "Local label of the gap junction connection.")
-        .def_readwrite("ggap", &arb::gap_junction_connection::ggap,
-            "Gap junction conductance [μS].")
+        .def_readwrite("weight", &arb::gap_junction_connection::weight,
+            "Gap junction connection weight [unit-less].")
         .def("__str__",  &gj_to_string)
         .def("__repr__", &gj_to_string);
 
diff --git a/python/test/fixtures.py b/python/test/fixtures.py
index 64145597..38142fb3 100644
--- a/python/test/fixtures.py
+++ b/python/test/fixtures.py
@@ -159,7 +159,7 @@ def cable_cell():
     # (3) Create cell and set properties
     decor = arbor.decor()
     decor.set_property(Vm=-40)
-    decor.paint('"soma"', 'hh')
+    decor.paint('"soma"', arbor.density('hh'))
     decor.place('"midpoint"', arbor.iclamp( 10, 2, 0.8), "iclamp")
     decor.place('"midpoint"', arbor.spike_detector(-10), "detector")
     return arbor.cable_cell(tree, labels, decor)
diff --git a/python/test/unit/test_cable_probes.py b/python/test/unit/test_cable_probes.py
index b992b4a6..6050e661 100644
--- a/python/test/unit/test_cable_probes.py
+++ b/python/test/unit/test_cable_probes.py
@@ -19,10 +19,10 @@ class cc_recipe(A.recipe):
 
         dec = A.decor()
 
-        dec.place('(location 0 0.08)', "expsyn", "syn0")
-        dec.place('(location 0 0.09)', "exp2syn", "syn1")
+        dec.place('(location 0 0.08)', A.synapse("expsyn"), "syn0")
+        dec.place('(location 0 0.09)', A.synapse("exp2syn"), "syn1")
         dec.place('(location 0 0.1)', A.iclamp(20.), "iclamp")
-        dec.paint('(all)', "hh")
+        dec.paint('(all)', A.density("hh"))
 
         self.cell = A.cable_cell(st, A.label_dict(), dec)
 
diff --git a/python/test/unit/test_catalogues.py b/python/test/unit/test_catalogues.py
index f94dc737..550dda8a 100644
--- a/python/test/unit/test_catalogues.py
+++ b/python/test/unit/test_catalogues.py
@@ -20,7 +20,7 @@ class recipe(arb.recipe):
             raise
 
         d = arb.decor()
-        d.paint('(all)', 'pas')
+        d.paint('(all)', arb.density('pas'))
         d.set_property(Vm=0.0)
         self.cell = arb.cable_cell(self.tree, arb.label_dict(), d)
 
diff --git a/test/common_cells.cpp b/test/common_cells.cpp
index 49db2965..5d1f3aa2 100644
--- a/test/common_cells.cpp
+++ b/test/common_cells.cpp
@@ -179,7 +179,7 @@ cable_cell_description make_cell_soma_only(bool with_stim) {
     soma_cell_builder builder(18.8/2.0);
 
     auto c = builder.make_cell();
-    c.decorations.paint("soma"_lab, "hh");
+    c.decorations.paint("soma"_lab, density("hh"));
     if (with_stim) {
         c.decorations.place(builder.location({0,0.5}), i_clamp{10., 100., 0.1}, "cc");
     }
@@ -213,8 +213,8 @@ cable_cell_description make_cell_ball_and_stick(bool with_stim) {
     builder.add_branch(0, 200, 1.0/2, 1.0/2, 4, "dend");
 
     auto c = builder.make_cell();
-    c.decorations.paint("soma"_lab, "hh");
-    c.decorations.paint("dend"_lab, "pas");
+    c.decorations.paint("soma"_lab, density("hh"));
+    c.decorations.paint("dend"_lab, density("pas"));
     if (with_stim) {
         c.decorations.place(builder.location({1,1}), i_clamp{5, 80, 0.3}, "cc");
     }
@@ -251,8 +251,8 @@ cable_cell_description make_cell_ball_and_3stick(bool with_stim) {
     builder.add_branch(1, 100, 0.5, 0.5, 4, "dend");
 
     auto c = builder.make_cell();
-    c.decorations.paint("soma"_lab, "hh");
-    c.decorations.paint("dend"_lab, "pas");
+    c.decorations.paint("soma"_lab, density("hh"));
+    c.decorations.paint("dend"_lab, density("pas"));
     if (with_stim) {
         c.decorations.place(builder.location({2,1}), i_clamp{5.,  80., 0.45}, "cc0");
         c.decorations.place(builder.location({3,1}), i_clamp{40., 10.,-0.2}, "cc1");
diff --git a/test/unit-distributed/test_communicator.cpp b/test/unit-distributed/test_communicator.cpp
index e7e081d5..e71590ba 100644
--- a/test/unit-distributed/test_communicator.cpp
+++ b/test/unit-distributed/test_communicator.cpp
@@ -200,7 +200,7 @@ namespace {
                 arb::decor decor;
                 decor.set_default(arb::cv_policy_fixed_per_branch(10));
                 decor.place(arb::mlocation{0, 0.5}, arb::threshold_detector{10}, "src");
-                decor.place(arb::mlocation{0, 0.5}, "expsyn", "tgt");
+                decor.place(arb::mlocation{0, 0.5}, arb::synapse("expsyn"), "tgt");
                 return arb::cable_cell(arb::morphology(tree), {}, decor);
             }
             return arb::lif_cell("src", "tgt");
@@ -273,7 +273,7 @@ namespace {
             arb::decor decor;
             decor.set_default(arb::cv_policy_fixed_per_branch(10));
             decor.place(arb::mlocation{0, 0.5}, arb::threshold_detector{10}, "src");
-            decor.place(arb::ls::uniform(arb::reg::all(), 0, size_, gid), "expsyn", "tgt");
+            decor.place(arb::ls::uniform(arb::reg::all(), 0, size_, gid), arb::synapse("expsyn"), "tgt");
             return arb::cable_cell(arb::morphology(tree), {}, decor);
         }
         cell_kind get_cell_kind(cell_gid_type gid) const override {
@@ -345,8 +345,8 @@ namespace {
             tree.append(arb::mnpos, {0, 0, 0.0, 1.0}, {0, 0, 200, 1.0}, 1);
             arb::decor decor;
             if (gid%3 != 1) {
-                decor.place(arb::ls::uniform(arb::reg::all(), 0, 1, gid), "expsyn", "synapses_0");
-                decor.place(arb::ls::uniform(arb::reg::all(), 2, 2, gid), "expsyn", "synapses_1");
+                decor.place(arb::ls::uniform(arb::reg::all(), 0, 1, gid), arb::synapse("expsyn"), "synapses_0");
+                decor.place(arb::ls::uniform(arb::reg::all(), 2, 2, gid), arb::synapse("expsyn"), "synapses_1");
             }
             else {
                 decor.place(arb::ls::uniform(arb::reg::all(), 0, 2, gid), arb::threshold_detector{10}, "detectors_0");
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index c52c27e6..dbb6ddfa 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -5,6 +5,8 @@ set(test_mechanisms
     celsius_test
     diam_test
     fixed_ica_current
+    gj0
+    gj1
     linear_ca_conc
     non_linear
     param_as_state
diff --git a/test/unit/mod/gj0.mod b/test/unit/mod/gj0.mod
new file mode 100644
index 00000000..b4eae217
--- /dev/null
+++ b/test/unit/mod/gj0.mod
@@ -0,0 +1,14 @@
+NEURON {
+    JUNCTION_PROCESS gj0
+    NONSPECIFIC_CURRENT i
+    RANGE g
+}
+INITIAL {}
+
+PARAMETER {
+    g = 1
+}
+
+BREAKPOINT {
+    i = g*(v - v_peer)
+}
diff --git a/test/unit/mod/gj1.mod b/test/unit/mod/gj1.mod
new file mode 100644
index 00000000..1083fe13
--- /dev/null
+++ b/test/unit/mod/gj1.mod
@@ -0,0 +1,15 @@
+NEURON {
+    JUNCTION_PROCESS gj1
+    NONSPECIFIC_CURRENT i
+    RANGE g, e
+}
+INITIAL {}
+
+PARAMETER {
+    g = 1
+    e = 0
+}
+
+BREAKPOINT {
+    i = g*(v - v_peer - e)
+}
diff --git a/test/unit/test_abi.cpp b/test/unit/test_abi.cpp
index 5cc736e1..29e533a8 100644
--- a/test/unit/test_abi.cpp
+++ b/test/unit/test_abi.cpp
@@ -47,12 +47,11 @@ TEST(abi, multicore_initialisation) {
     std::vector<arb_value_type> temp(ncv, 23);
     std::vector<arb_value_type> diam(ncv, 1.);
     std::vector<arb_value_type> vinit(ncv, -65);
-    std::vector<arb::fvm_gap_junction> gj = {};
     std::vector<arb_index_type> src_to_spike = {};
 
     arb::multicore::shared_state shared_state(ncell, ncell, 0,
                                               cv_to_intdom, cv_to_intdom,
-                                              gj, vinit, temp, diam, src_to_spike,
+                                              vinit, temp, diam, src_to_spike,
                                               mech.data_alignment());
 
     arb::mechanism_layout layout;
@@ -125,12 +124,11 @@ TEST(abi, multicore_null) {
     std::vector<arb_value_type> temp(ncv, 23);
     std::vector<arb_value_type> diam(ncv, 1.);
     std::vector<arb_value_type> vinit(ncv, -65);
-    std::vector<arb::fvm_gap_junction> gj = {};
     std::vector<arb_index_type> src_to_spike = {};
 
     arb::multicore::shared_state shared_state(ncell, ncell, 0,
                                               cv_to_intdom, cv_to_intdom,
-                                              gj, vinit, temp, diam, src_to_spike,
+                                              vinit, temp, diam, src_to_spike,
                                               mech.data_alignment());
 
     arb::mechanism_layout layout;
@@ -190,12 +188,11 @@ TEST(abi, gpu_initialisation) {
     std::vector<arb_value_type> temp(ncv, 23);
     std::vector<arb_value_type> diam(ncv, 1.);
     std::vector<arb_value_type> vinit(ncv, -65);
-    std::vector<arb::fvm_gap_junction> gj = {};
     std::vector<arb_index_type> src_to_spike = {};
 
     arb::gpu::shared_state shared_state(ncell, ncell, 0,
                                         cv_to_intdom, cv_to_intdom,
-                                        gj, vinit, temp, diam, src_to_spike,
+                                        vinit, temp, diam, src_to_spike,
                                         1);
 
     arb::mechanism_layout layout;
@@ -267,12 +264,11 @@ TEST(abi, gpu_null) {
     std::vector<arb_value_type> temp(ncv, 23);
     std::vector<arb_value_type> diam(ncv, 1.);
     std::vector<arb_value_type> vinit(ncv, -65);
-    std::vector<arb::fvm_gap_junction> gj = {};
     std::vector<arb_index_type> src_to_spike = {};
 
     arb::gpu::shared_state shared_state(ncell, ncell, 0,
                                         cv_to_intdom, cv_to_intdom,
-                                        gj, vinit, temp, diam, src_to_spike,
+                                        vinit, temp, diam, src_to_spike,
                                         1);
 
     arb::mechanism_layout layout;
diff --git a/test/unit/test_cable_cell.cpp b/test/unit/test_cable_cell.cpp
index 9019330c..4dee0af7 100644
--- a/test/unit/test_cable_cell.cpp
+++ b/test/unit/test_cable_cell.cpp
@@ -31,13 +31,13 @@ TEST(cable_cell, lid_ranges) {
 
     // Place synapses and threshold detectors in interleaved order.
     // Note: there are 2 terminal points.
-    decorations.place("term"_lab, "expsyn", "t0");
-    decorations.place("term"_lab, "expsyn", "t1");
+    decorations.place("term"_lab, synapse("expsyn"), "t0");
+    decorations.place("term"_lab, synapse("expsyn"), "t1");
     decorations.place("term"_lab, threshold_detector{-10}, "s0");
-    decorations.place(empty_sites, "expsyn", "t2");
+    decorations.place(empty_sites, synapse("expsyn"), "t2");
     decorations.place("term"_lab, threshold_detector{-20}, "s1");
-    decorations.place(three_sites, "expsyn", "t3");
-    decorations.place("term"_lab, "exp2syn", "t3");
+    decorations.place(three_sites, synapse("expsyn"), "t3");
+    decorations.place("term"_lab, synapse("exp2syn"), "t3");
 
     cable_cell cell(morph, dict, decorations);
 
diff --git a/test/unit/test_domain_decomposition.cpp b/test/unit/test_domain_decomposition.cpp
index 223f894f..b2f0db19 100644
--- a/test/unit/test_domain_decomposition.cpp
+++ b/test/unit/test_domain_decomposition.cpp
@@ -68,7 +68,7 @@ namespace {
 
         arb::util::unique_any get_cell_description(cell_gid_type) const override {
             auto c = arb::make_cell_soma_only(false);
-            c.decorations.place(mlocation{0,1}, gap_junction_site{}, "gj");
+            c.decorations.place(mlocation{0,1}, junction("gj"), "gj");
             return {arb::cable_cell(c)};
         }
 
diff --git a/test/unit/test_event_delivery.cpp b/test/unit/test_event_delivery.cpp
index 9ffcb086..caca91e3 100644
--- a/test/unit/test_event_delivery.cpp
+++ b/test/unit/test_event_delivery.cpp
@@ -35,9 +35,9 @@ struct test_recipe: public n_cable_cell_recipe {
         labels.set("soma", arb::reg::tagged(1));
 
         decor decorations;
-        decorations.place(mlocation{0, 0.5}, "expsyn", "synapse");
+        decorations.place(mlocation{0, 0.5}, synapse("expsyn"), "synapse");
         decorations.place(mlocation{0, 0.5}, threshold_detector{-64}, "detector");
-        decorations.place(mlocation{0, 0.5}, gap_junction_site{}, "gapjunction");
+        decorations.place(mlocation{0, 0.5}, junction("gj"), "gapjunction");
         cable_cell c(st, labels, decorations);
 
         return c;
@@ -110,7 +110,7 @@ struct test_recipe_gj: public test_recipe {
             if (p.first == i) gjs.push_back({{p.second, "gapjunction", lid_selection_policy::assert_univalent},
                                              {"gapjunction", lid_selection_policy::assert_univalent}, 0.});
             if (p.second == i) gjs.push_back({{p.first, "gapjunction", lid_selection_policy::assert_univalent},
-                                                    {"gapjunction", lid_selection_policy::assert_univalent}, 0.});
+                                             {"gapjunction", lid_selection_policy::assert_univalent}, 0.});
         }
         return gjs;
     }
diff --git a/test/unit/test_fvm_layout.cpp b/test/unit/test_fvm_layout.cpp
index 981b291a..0d227ee2 100644
--- a/test/unit/test_fvm_layout.cpp
+++ b/test/unit/test_fvm_layout.cpp
@@ -1,15 +1,22 @@
 #include <limits>
 #include <string>
+#include <tuple>
 #include <vector>
 
-#include <arborio/label_parse.hpp>
-
 #include <arbor/cable_cell.hpp>
 #include <arbor/math.hpp>
 #include <arbor/mechcat.hpp>
-#include "arbor/cable_cell_param.hpp"
-#include "arbor/morph/morphology.hpp"
-#include "arbor/morph/segment_tree.hpp"
+#include <arbor/cable_cell_param.hpp>
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/segment_tree.hpp>
+
+#include <arborio/label_parse.hpp>
+
+#include <arborenv/concurrency.hpp>
+
+#include "backends/multicore/fvm.hpp"
+#include "fvm_lowered_cell.hpp"
+#include "fvm_lowered_cell_impl.hpp"
 #include "fvm_layout.hpp"
 #include "util/maputil.hpp"
 #include "util/rangeutil.hpp"
@@ -30,6 +37,9 @@ using util::count_along;
 using util::ptr_by_key;
 using util::value_by_key;
 
+using backend = arb::multicore::backend;
+using fvm_cell = arb::fvm_lowered_cell_impl<backend>;
+
 namespace {
     struct system {
         std::vector<soma_cell_builder> builders;
@@ -56,8 +66,8 @@ namespace {
             builder.add_branch(0, 200, 1.0/2, 1.0/2, 4, "dend");
 
             auto description = builder.make_cell();
-            description.decorations.paint("soma"_lab, "hh");
-            description.decorations.paint("dend"_lab, "pas");
+            description.decorations.paint("soma"_lab, density("hh"));
+            description.decorations.paint("dend"_lab, density("pas"));
             description.decorations.place(builder.location({1,1}), i_clamp{5, 80, 0.3}, "clamp");
 
             s.builders.push_back(std::move(builder));
@@ -99,8 +109,8 @@ namespace {
             auto b3 = b.add_branch(1, 180, 0.35, 0.35, 4, "dend");
             auto desc = b.make_cell();
 
-            desc.decorations.paint("soma"_lab, "hh");
-            desc.decorations.paint("dend"_lab, "pas");
+            desc.decorations.paint("soma"_lab, density("hh"));
+            desc.decorations.paint("dend"_lab, density("pas"));
 
             using ::arb::reg::branch;
             auto c1 = reg::cable(b1-1, b.location({b1, 0}).pos, 1);
@@ -127,6 +137,77 @@ namespace {
         ASSERT_EQ(1u, cells[0].morphology().num_branches());
         ASSERT_EQ(3u, cells[1].morphology().num_branches());
     }
+
+    system six_cell_gj_system() {
+        system s;
+        auto& descriptions = s.descriptions;
+
+        // Cell 0: simple soma cell, 1 CV.
+        {
+            soma_cell_builder builder(12.6157/2.0);
+
+            auto desc = builder.make_cell();
+            desc.decorations.paint("soma"_lab, density("hh"));
+
+            s.builders.push_back(std::move(builder));
+            descriptions.push_back(desc);
+        }
+        // Cell 1: ball and 3-stick, 2 CVs per dendrite, 1 CV for the soma.
+        {
+            soma_cell_builder b(7.);
+            b.add_branch(0, 200, 0.5,  0.5, 2,  "dend");
+            b.add_branch(1, 300, 0.4,  0.4, 2,  "dend");
+            b.add_branch(1, 180, 0.35, 0.35, 2, "dend");
+            auto desc = b.make_cell();
+
+            desc.decorations.paint("soma"_lab, density("hh"));
+            desc.decorations.paint("dend"_lab, density("pas"));
+
+            s.builders.push_back(std::move(b));
+            descriptions.push_back(desc);
+        }
+        // Cell 2: ball and stick, 1 CV for the soma, 2 CVs for the dendrite.
+        {
+            soma_cell_builder b(7.);
+            b.add_branch(0, 200, 0.5,  0.5, 2,  "dend");
+            auto desc = b.make_cell();
+
+            s.builders.push_back(std::move(b));
+            descriptions.push_back(desc);
+        }
+        // Cell 3: simple soma cell, 1 CV. 1 gap junction.
+        {
+            soma_cell_builder b(7.);
+            auto desc = b.make_cell();
+
+            desc.decorations.paint("soma"_lab, density("hh"));
+
+            s.builders.push_back(std::move(b));
+            descriptions.push_back(desc);
+        }
+        // Cell 4: ball and stick, 1 CV for the soma, 3 CVs for the dendrite.
+        {
+            soma_cell_builder b(7.);
+            b.add_branch(0, 200, 0.5,  0.5, 3,  "dend");
+            auto desc = b.make_cell();
+
+            desc.decorations.paint("soma"_lab, density("pas"));
+
+            s.builders.push_back(std::move(b));
+            descriptions.push_back(desc);
+        }
+        // Cell 5: ball and stick, 1 CV for the soma, 1 CV for the dendrite.
+        {
+            soma_cell_builder b(7.);
+            b.add_branch(0, 200, 0.5,  0.5, 1,  "dend");
+            auto desc = b.make_cell();
+
+            s.builders.push_back(std::move(b));
+            descriptions.push_back(desc);
+        }
+        return s;
+    }
+
 } // namespace
 
 template<typename P>
@@ -137,17 +218,19 @@ void check_compatible_mechanism_failure(cable_cell_global_properties gprop, P pa
 
     auto cells = system.cells();
     check_two_cell_system(cells);
+    std::vector<cell_gid_type> gids = {0,1};
+    std::unordered_map<cell_gid_type, std::vector<fvm_gap_junction>> gj_conns = {{0, {}}, {1, {}}};
     fvm_cv_discretization D = fvm_cv_discretize(cells, gprop.default_parameters);
 
-    EXPECT_THROW(fvm_build_mechanism_data(gprop, cells, D), arb::cable_cell_error);
+    EXPECT_THROW(fvm_build_mechanism_data(gprop, cells, gids, gj_conns, D), arb::cable_cell_error);
 }
 
 TEST(fvm_layout, compatible_mechanisms) {
     cable_cell_global_properties gprop;
     gprop.default_parameters = neuron_parameter_defaults;
 
-    check_compatible_mechanism_failure(gprop, [](auto& sys) {sys.descriptions[0].decorations.place(sys.builders[0].location({1, 0.4}), "hh", "syn0"); });
-    check_compatible_mechanism_failure(gprop, [](auto& sys) {sys.descriptions[1].decorations.paint(sys.builders[1].cable(mcable{0}), "expsyn"); });
+    check_compatible_mechanism_failure(gprop, [](auto& sys) {sys.descriptions[0].decorations.place(sys.builders[0].location({1, 0.4}), synapse("hh"), "syn0"); });
+    check_compatible_mechanism_failure(gprop, [](auto& sys) {sys.descriptions[1].decorations.paint(sys.builders[1].cable(mcable{0}), density("expsyn")); });
     check_compatible_mechanism_failure(gprop, [](auto& sys) {sys.descriptions[0].decorations.set_default(ion_reversal_potential_method{"na", "expsyn"}); });
 
     gprop.default_parameters.reversal_potential_method["na"] = "pas";
@@ -160,18 +243,20 @@ TEST(fvm_layout, mech_index) {
     auto& builders = system.builders;
 
     // Add four synapses of two varieties across the cells.
-    descriptions[0].decorations.place(builders[0].location({1, 0.4}), "expsyn", "syn0");
-    descriptions[0].decorations.place(builders[0].location({1, 0.4}), "expsyn", "syn1");
-    descriptions[1].decorations.place(builders[1].location({2, 0.4}), "exp2syn", "syn3");
-    descriptions[1].decorations.place(builders[1].location({3, 0.4}), "expsyn", "syn4");
+    descriptions[0].decorations.place(builders[0].location({1, 0.4}), synapse("expsyn"), "syn0");
+    descriptions[0].decorations.place(builders[0].location({1, 0.4}), synapse("expsyn"), "syn1");
+    descriptions[1].decorations.place(builders[1].location({2, 0.4}), synapse("exp2syn"), "syn3");
+    descriptions[1].decorations.place(builders[1].location({3, 0.4}), synapse("expsyn"), "syn4");
 
     cable_cell_global_properties gprop;
     gprop.default_parameters = neuron_parameter_defaults;
 
     auto cells = system.cells();
     check_two_cell_system(cells);
+    std::vector<cell_gid_type> gids = {0,1};
+    std::unordered_map<cell_gid_type, std::vector<fvm_gap_junction>> gj_conns = {{0, {}}, {1, {}}};
     fvm_cv_discretization D = fvm_cv_discretize(cells, gprop.default_parameters);
-    fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, D);
+    fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, gids, gj_conns, D);
 
     auto& hh_config = M.mechanisms.at("hh");
     auto& expsyn_config = M.mechanisms.at("expsyn");
@@ -280,14 +365,14 @@ TEST(fvm_layout, coalescing_synapses) {
     {
         auto desc = builder.make_cell();
 
-        desc.decorations.place(builder.location({1, 0.3}), "expsyn", "syn0");
-        desc.decorations.place(builder.location({1, 0.5}), "expsyn", "syn1");
-        desc.decorations.place(builder.location({1, 0.7}), "expsyn", "syn2");
-        desc.decorations.place(builder.location({1, 0.9}), "expsyn", "syn3");
+        desc.decorations.place(builder.location({1, 0.3}), synapse("expsyn"), "syn0");
+        desc.decorations.place(builder.location({1, 0.5}), synapse("expsyn"), "syn1");
+        desc.decorations.place(builder.location({1, 0.7}), synapse("expsyn"), "syn2");
+        desc.decorations.place(builder.location({1, 0.9}), synapse("expsyn"), "syn3");
 
         cable_cell cell(desc);
         fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
-        fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
+        fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, {0}, {{0, {}}}, D);
 
         auto &expsyn_config = M.mechanisms.at("expsyn");
         EXPECT_EQ(ivec({2, 3, 4, 5}), expsyn_config.cv);
@@ -297,14 +382,14 @@ TEST(fvm_layout, coalescing_synapses) {
         auto desc = builder.make_cell();
 
         // Add synapses of two varieties.
-        desc.decorations.place(builder.location({1, 0.3}), "expsyn", "syn0");
-        desc.decorations.place(builder.location({1, 0.5}), "exp2syn", "syn1");
-        desc.decorations.place(builder.location({1, 0.7}), "expsyn", "syn2");
-        desc.decorations.place(builder.location({1, 0.9}), "exp2syn", "syn3");
+        desc.decorations.place(builder.location({1, 0.3}), synapse("expsyn"), "syn0");
+        desc.decorations.place(builder.location({1, 0.5}), synapse("exp2syn"), "syn1");
+        desc.decorations.place(builder.location({1, 0.7}), synapse("expsyn"), "syn2");
+        desc.decorations.place(builder.location({1, 0.9}), synapse("exp2syn"), "syn3");
 
         cable_cell cell(desc);
         fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
-        fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
+        fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, {0}, {{0, {}}}, D);
 
         auto &expsyn_config = M.mechanisms.at("expsyn");
         EXPECT_EQ(ivec({2, 4}), expsyn_config.cv);
@@ -317,14 +402,14 @@ TEST(fvm_layout, coalescing_synapses) {
     {
         auto desc = builder.make_cell();
 
-        desc.decorations.place(builder.location({1, 0.3}), "expsyn", "syn0");
-        desc.decorations.place(builder.location({1, 0.5}), "expsyn", "syn1");
-        desc.decorations.place(builder.location({1, 0.7}), "expsyn", "syn2");
-        desc.decorations.place(builder.location({1, 0.9}), "expsyn", "syn3");
+        desc.decorations.place(builder.location({1, 0.3}), synapse("expsyn"), "syn0");
+        desc.decorations.place(builder.location({1, 0.5}), synapse("expsyn"), "syn1");
+        desc.decorations.place(builder.location({1, 0.7}), synapse("expsyn"), "syn2");
+        desc.decorations.place(builder.location({1, 0.9}), synapse("expsyn"), "syn3");
 
         cable_cell cell(desc);
         fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
-        fvm_mechanism_data M = fvm_build_mechanism_data(gprop_no_coalesce, {cell}, D);
+        fvm_mechanism_data M = fvm_build_mechanism_data(gprop_no_coalesce, {cell}, {0}, {{0, {}}}, D);
 
         auto &expsyn_config = M.mechanisms.at("expsyn");
         EXPECT_EQ(ivec({2, 3, 4, 5}), expsyn_config.cv);
@@ -334,14 +419,14 @@ TEST(fvm_layout, coalescing_synapses) {
         auto desc = builder.make_cell();
 
         // Add synapses of two varieties.
-        desc.decorations.place(builder.location({1, 0.3}), "expsyn", "syn0");
-        desc.decorations.place(builder.location({1, 0.5}), "exp2syn", "syn1");
-        desc.decorations.place(builder.location({1, 0.7}), "expsyn", "syn2");
-        desc.decorations.place(builder.location({1, 0.9}), "exp2syn", "syn3");
+        desc.decorations.place(builder.location({1, 0.3}), synapse("expsyn"), "syn0");
+        desc.decorations.place(builder.location({1, 0.5}), synapse("exp2syn"), "syn1");
+        desc.decorations.place(builder.location({1, 0.7}), synapse("expsyn"), "syn2");
+        desc.decorations.place(builder.location({1, 0.9}), synapse("exp2syn"), "syn3");
 
         cable_cell cell(desc);
         fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
-        fvm_mechanism_data M = fvm_build_mechanism_data(gprop_no_coalesce, {cell}, D);
+        fvm_mechanism_data M = fvm_build_mechanism_data(gprop_no_coalesce, {cell}, {0}, {{0, {}}}, D);
 
         auto &expsyn_config = M.mechanisms.at("expsyn");
         EXPECT_EQ(ivec({2, 4}), expsyn_config.cv);
@@ -355,14 +440,14 @@ TEST(fvm_layout, coalescing_synapses) {
         auto desc = builder.make_cell();
 
         // Add synapses of two varieties.
-        desc.decorations.place(builder.location({1, 0.3}), "expsyn", "syn0");
-        desc.decorations.place(builder.location({1, 0.3}), "expsyn", "syn1");
-        desc.decorations.place(builder.location({1, 0.7}), "expsyn", "syn2");
-        desc.decorations.place(builder.location({1, 0.7}), "expsyn", "syn3");
+        desc.decorations.place(builder.location({1, 0.3}), synapse("expsyn"), "syn0");
+        desc.decorations.place(builder.location({1, 0.3}), synapse("expsyn"), "syn1");
+        desc.decorations.place(builder.location({1, 0.7}), synapse("expsyn"), "syn2");
+        desc.decorations.place(builder.location({1, 0.7}), synapse("expsyn"), "syn3");
 
         cable_cell cell(desc);
         fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
-        fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
+        fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, {0}, {{0, {}}}, D);
 
         auto &expsyn_config = M.mechanisms.at("expsyn");
         EXPECT_EQ(ivec({2, 4}), expsyn_config.cv);
@@ -372,14 +457,14 @@ TEST(fvm_layout, coalescing_synapses) {
         auto desc = builder.make_cell();
 
         // Add synapses of two varieties.
-        desc.decorations.place(builder.location({1, 0.3}), syn_desc("expsyn", 0, 0.2), "syn0");
-        desc.decorations.place(builder.location({1, 0.3}), syn_desc("expsyn", 0, 0.2), "syn1");
-        desc.decorations.place(builder.location({1, 0.3}), syn_desc("expsyn", 0.1, 0.2), "syn2");
-        desc.decorations.place(builder.location({1, 0.7}), syn_desc("expsyn", 0.1, 0.2), "syn3");
+        desc.decorations.place(builder.location({1, 0.3}), synapse(syn_desc("expsyn", 0, 0.2)), "syn0");
+        desc.decorations.place(builder.location({1, 0.3}), synapse(syn_desc("expsyn", 0, 0.2)), "syn1");
+        desc.decorations.place(builder.location({1, 0.3}), synapse(syn_desc("expsyn", 0.1, 0.2)), "syn2");
+        desc.decorations.place(builder.location({1, 0.7}), synapse(syn_desc("expsyn", 0.1, 0.2)), "syn3");
 
         cable_cell cell(desc);
         fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
-        fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
+        fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, {0}, {{0, {}}}, D);
 
         std::vector<exp_instance> instances{
             exp_instance(2, L{0, 1}, 0., 0.2),
@@ -395,18 +480,18 @@ TEST(fvm_layout, coalescing_synapses) {
         auto desc = builder.make_cell();
 
         // Add synapses of two varieties.
-        desc.decorations.place(builder.location({1, 0.7}), syn_desc("expsyn", 0, 3), "syn0");
-        desc.decorations.place(builder.location({1, 0.7}), syn_desc("expsyn", 1, 3), "syn1");
-        desc.decorations.place(builder.location({1, 0.7}), syn_desc("expsyn", 0, 3), "syn2");
-        desc.decorations.place(builder.location({1, 0.7}), syn_desc("expsyn", 1, 3), "syn3");
-        desc.decorations.place(builder.location({1, 0.3}), syn_desc("expsyn", 0, 2), "syn4");
-        desc.decorations.place(builder.location({1, 0.3}), syn_desc("expsyn", 1, 2), "syn5");
-        desc.decorations.place(builder.location({1, 0.3}), syn_desc("expsyn", 0, 2), "syn6");
-        desc.decorations.place(builder.location({1, 0.3}), syn_desc("expsyn", 1, 2), "syn7");
+        desc.decorations.place(builder.location({1, 0.7}), synapse(syn_desc("expsyn", 0, 3)), "syn0");
+        desc.decorations.place(builder.location({1, 0.7}), synapse(syn_desc("expsyn", 1, 3)), "syn1");
+        desc.decorations.place(builder.location({1, 0.7}), synapse(syn_desc("expsyn", 0, 3)), "syn2");
+        desc.decorations.place(builder.location({1, 0.7}), synapse(syn_desc("expsyn", 1, 3)), "syn3");
+        desc.decorations.place(builder.location({1, 0.3}), synapse(syn_desc("expsyn", 0, 2)), "syn4");
+        desc.decorations.place(builder.location({1, 0.3}), synapse(syn_desc("expsyn", 1, 2)), "syn5");
+        desc.decorations.place(builder.location({1, 0.3}), synapse(syn_desc("expsyn", 0, 2)), "syn6");
+        desc.decorations.place(builder.location({1, 0.3}), synapse(syn_desc("expsyn", 1, 2)), "syn7");
 
         cable_cell cell(desc);
         fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
-        fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
+        fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, {0}, {{0, {}}}, D);
 
         std::vector<exp_instance> instances{
             exp_instance(2, L{4, 6}, 0.0, 2.0),
@@ -423,20 +508,20 @@ TEST(fvm_layout, coalescing_synapses) {
         auto desc = builder.make_cell();
 
         // Add synapses of two varieties.
-        desc.decorations.place(builder.location({1, 0.3}), syn_desc("expsyn",  1, 2), "syn0");
-        desc.decorations.place(builder.location({1, 0.3}), syn_desc_2("exp2syn", 4, 1), "syn1");
-        desc.decorations.place(builder.location({1, 0.3}), syn_desc("expsyn",  1, 2), "syn2");
-        desc.decorations.place(builder.location({1, 0.3}), syn_desc("expsyn",  5, 1), "syn3");
-        desc.decorations.place(builder.location({1, 0.3}), syn_desc_2("exp2syn", 1, 3), "syn4");
-        desc.decorations.place(builder.location({1, 0.3}), syn_desc("expsyn",  1, 2), "syn5");
-        desc.decorations.place(builder.location({1, 0.7}), syn_desc_2("exp2syn", 2, 2), "syn6");
-        desc.decorations.place(builder.location({1, 0.7}), syn_desc_2("exp2syn", 2, 1), "syn7");
-        desc.decorations.place(builder.location({1, 0.7}), syn_desc_2("exp2syn", 2, 1), "syn8");
-        desc.decorations.place(builder.location({1, 0.7}), syn_desc_2("exp2syn", 2, 2), "syn9");
+        desc.decorations.place(builder.location({1, 0.3}), synapse(syn_desc("expsyn",  1, 2)), "syn0");
+        desc.decorations.place(builder.location({1, 0.3}), synapse(syn_desc_2("exp2syn", 4, 1)), "syn1");
+        desc.decorations.place(builder.location({1, 0.3}), synapse(syn_desc("expsyn",  1, 2)), "syn2");
+        desc.decorations.place(builder.location({1, 0.3}), synapse(syn_desc("expsyn",  5, 1)), "syn3");
+        desc.decorations.place(builder.location({1, 0.3}), synapse(syn_desc_2("exp2syn", 1, 3)), "syn4");
+        desc.decorations.place(builder.location({1, 0.3}), synapse(syn_desc("expsyn",  1, 2)), "syn5");
+        desc.decorations.place(builder.location({1, 0.7}), synapse(syn_desc_2("exp2syn", 2, 2)), "syn6");
+        desc.decorations.place(builder.location({1, 0.7}), synapse(syn_desc_2("exp2syn", 2, 1)), "syn7");
+        desc.decorations.place(builder.location({1, 0.7}), synapse(syn_desc_2("exp2syn", 2, 1)), "syn8");
+        desc.decorations.place(builder.location({1, 0.7}), synapse(syn_desc_2("exp2syn", 2, 2)), "syn9");
 
         cable_cell cell(desc);
         fvm_cv_discretization D = fvm_cv_discretize({cell}, neuron_parameter_defaults);
-        fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, D);
+        fvm_mechanism_data M = fvm_build_mechanism_data(gprop_coalesce, {cell}, {0}, {{0, {}}}, D);
 
         for (auto &instance: {exp_instance(2, L{0,2,5}, 1, 2),
                               exp_instance(2, L{3},     5, 1)}) {
@@ -472,21 +557,23 @@ TEST(fvm_layout, synapse_targets) {
         return mechanism_desc(name).set("e", syn_e.at(idx));
     };
 
-    descriptions[0].decorations.place(builders[0].location({1, 0.9}), syn_desc("expsyn", 0), "syn0");
-    descriptions[0].decorations.place(builders[0].location({0, 0.5}), syn_desc("expsyn", 1), "syn1");
-    descriptions[0].decorations.place(builders[0].location({1, 0.4}), syn_desc("expsyn", 2), "syn2");
+    descriptions[0].decorations.place(builders[0].location({1, 0.9}), synapse(syn_desc("expsyn", 0)), "syn0");
+    descriptions[0].decorations.place(builders[0].location({0, 0.5}), synapse(syn_desc("expsyn", 1)), "syn1");
+    descriptions[0].decorations.place(builders[0].location({1, 0.4}), synapse(syn_desc("expsyn", 2)), "syn2");
 
-    descriptions[1].decorations.place(builders[1].location({2, 0.4}), syn_desc("exp2syn", 3), "syn3");
-    descriptions[1].decorations.place(builders[1].location({1, 0.4}), syn_desc("exp2syn", 4), "syn4");
-    descriptions[1].decorations.place(builders[1].location({3, 0.4}), syn_desc("expsyn", 5), "syn5");
-    descriptions[1].decorations.place(builders[1].location({3, 0.7}), syn_desc("exp2syn", 6), "syn6");
+    descriptions[1].decorations.place(builders[1].location({2, 0.4}), synapse(syn_desc("exp2syn", 3)), "syn3");
+    descriptions[1].decorations.place(builders[1].location({1, 0.4}), synapse(syn_desc("exp2syn", 4)), "syn4");
+    descriptions[1].decorations.place(builders[1].location({3, 0.4}), synapse(syn_desc("expsyn", 5)), "syn5");
+    descriptions[1].decorations.place(builders[1].location({3, 0.7}), synapse(syn_desc("exp2syn", 6)), "syn6");
 
     cable_cell_global_properties gprop;
     gprop.default_parameters = neuron_parameter_defaults;
 
     auto cells = system.cells();
+    std::vector<cell_gid_type> gids = {0,1};
+    std::unordered_map<cell_gid_type, std::vector<fvm_gap_junction>> gj_conns = {{0, {}}, {1, {}}};
     fvm_cv_discretization D = fvm_cv_discretize(cells, gprop.default_parameters);
-    fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, D);
+    fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, gids, gj_conns, D);
 
     ASSERT_EQ(1u, M.mechanisms.count("expsyn"));
     ASSERT_EQ(1u, M.mechanisms.count("exp2syn"));
@@ -521,6 +608,657 @@ TEST(fvm_layout, synapse_targets) {
     }
 }
 
+TEST(fvm_lowered, gj_example_0) {
+    arb::proc_allocation resources;
+    if (auto nt = arbenv::get_env_num_threads()) {
+        resources.num_threads = nt;
+    } else {
+        resources.num_threads = arbenv::thread_concurrency();
+    }
+    arb::execution_context context(resources);
+
+    class gap_recipe: public recipe {
+    public:
+        gap_recipe(std::vector<cable_cell> cells, arb::cable_cell_global_properties gprop) : gprop_(gprop), cells_(cells) {}
+
+        cell_size_type num_cells() const override { return n_; }
+        cell_kind get_cell_kind(cell_gid_type) const override { return cell_kind::cable; }
+        util::unique_any get_cell_description(cell_gid_type gid) const override {
+            return cells_[gid];
+        }
+        std::vector<arb::gap_junction_connection> gap_junctions_on(cell_gid_type gid) const override{
+            std::vector<gap_junction_connection> conns;
+            conns.push_back(gap_junction_connection({(gid+1)%2, "gj", lid_selection_policy::assert_univalent}, {"gj", lid_selection_policy::assert_univalent}, 0.5));
+            return conns;
+        }
+        std::any get_global_properties(cell_kind) const override {
+            return gprop_;
+        }
+    protected:
+        arb::cable_cell_global_properties gprop_;
+        std::vector<cable_cell> cells_;
+        cell_size_type n_ = 2;
+    };
+
+    std::vector<cable_cell> cells;
+    soma_cell_builder b0(2.1);
+    b0.add_branch(0, 10, 0.3, 0.2, 5, "dend");
+    auto c0 = b0.make_cell();
+    auto loc_0 = b0.location({1, 0.8});
+    c0.decorations.place(loc_0, junction("gj"), "gj");
+    cells.push_back(c0);
+
+    soma_cell_builder b1(2.4);
+    b1.add_branch(0, 10, 0.3, 0.2, 2, "dend");
+    auto c1 = b1.make_cell();
+    auto loc_1 = b1.location({1, 1});
+    c1.decorations.place(loc_1, junction("gj", {{"g", 0.5}}), "gj");
+    cells.push_back(c1);
+
+    // Check the GJ CV map
+    cable_cell_global_properties gprop;
+    gprop.default_parameters = neuron_parameter_defaults;
+
+    std::vector<cell_gid_type> gids = {0, 1};
+
+    auto D = fvm_cv_discretize(cells, gprop.default_parameters, context);
+    auto gj_cvs = fvm_build_gap_junction_cv_map(cells, gids, D);
+
+    auto cv_0 = D.geometry.location_cv(0, loc_0, cv_prefer::cv_nonempty);
+    auto cv_1 = D.geometry.location_cv(1, loc_1, cv_prefer::cv_nonempty);
+
+    EXPECT_EQ(2u, gj_cvs.size());
+
+    EXPECT_EQ(cv_0, gj_cvs.at(cell_member_type{0, 0}));
+    EXPECT_EQ(cv_1, gj_cvs.at(cell_member_type{1, 0}));
+
+    // Check the resolved GJ connections
+    fvm_cell fvcell(context);
+    gap_recipe rec(cells, gprop);
+
+    auto fvm_info = fvcell.initialize(gids, rec);
+    auto gj_conns = fvm_resolve_gj_connections(gids, fvm_info.gap_junction_data, gj_cvs, rec);
+
+    EXPECT_EQ(1u, gj_conns.at(0).size());
+    EXPECT_EQ(1u, gj_conns.at(1).size());
+
+    auto gj0 = fvm_gap_junction{0, cv_0, cv_1, 0.5};
+    auto gj1 = fvm_gap_junction{0, cv_1, cv_0, 0.5};
+
+    EXPECT_EQ(gj0, gj_conns.at(0).front());
+    EXPECT_EQ(gj1, gj_conns.at(1).front());
+
+    // Check the GJ mechanism data
+    auto M = fvm_build_mechanism_data(gprop, cells, gids, gj_conns, D, context);
+
+    EXPECT_EQ(1u, M.mechanisms.size());
+    ASSERT_EQ(1u, M.mechanisms.count("gj"));
+
+    const auto& gj_data = M.mechanisms["gj"];
+    const auto& gj_g_param = *ptr_by_key(gj_data.param_values, "g"s);
+
+    std::vector<std::tuple<int, int, double, double>> expected_gj_values = {
+        {cv_0, cv_1, 0.5, 1.},
+        {cv_1, cv_0, 0.5, 0.5},
+    };
+    util::sort(expected_gj_values);
+
+    std::vector<int> expected_gj_cv, expected_gj_peer_cv;
+    std::vector<double> expected_gj_weight, expected_gj_g;
+    util::assign(expected_gj_cv,      util::transform_view(expected_gj_values, [](const auto& x){return std::get<0>(x);}));
+    util::assign(expected_gj_peer_cv, util::transform_view(expected_gj_values, [](const auto& x){return std::get<1>(x);}));
+    util::assign(expected_gj_weight,  util::transform_view(expected_gj_values, [](const auto& x){return std::get<2>(x);}));
+    util::assign(expected_gj_g,       util::transform_view(expected_gj_values, [](const auto& x){return std::get<3>(x);}));
+
+    EXPECT_EQ(2u, gj_data.cv.size());
+    EXPECT_EQ(expected_gj_cv,      gj_data.cv);
+    EXPECT_EQ(expected_gj_peer_cv, gj_data.peer_cv);
+    EXPECT_EQ(expected_gj_weight,  gj_data.local_weight);
+    EXPECT_EQ(expected_gj_g,       gj_g_param);
+}
+
+TEST(fvm_lowered, gj_example_1) {
+    arb::proc_allocation resources;
+    if (auto nt = arbenv::get_env_num_threads()) {
+        resources.num_threads = nt;
+    } else {
+        resources.num_threads = arbenv::thread_concurrency();
+    }
+    arb::execution_context context(resources);
+
+    class gap_recipe: public recipe {
+    public:
+        gap_recipe(std::vector<cable_cell> cells, arb::cable_cell_global_properties gprop) : gprop_(gprop), cells_(cells) {}
+
+        cell_size_type num_cells() const override { return n_; }
+        cell_kind get_cell_kind(cell_gid_type) const override { return cell_kind::cable; }
+        util::unique_any get_cell_description(cell_gid_type gid) const override {
+            return cells_[gid];
+        }
+        std::vector<arb::gap_junction_connection> gap_junctions_on(cell_gid_type gid) const override{
+            std::vector<gap_junction_connection> conns;
+            switch (gid) {
+            case 0:
+                return {
+                    gap_junction_connection({2, "gj0"}, {"gj1"}, 0.01),
+                    gap_junction_connection({1, "gj0"}, {"gj0"}, 0.03),
+                    gap_junction_connection({1, "gj1"}, {"gj0"}, 0.04)
+                };
+            case 1:
+                return {
+                    gap_junction_connection({0, "gj0"}, {"gj0"}, 0.03),
+                    gap_junction_connection({0, "gj0"}, {"gj1"}, 0.04),
+                    gap_junction_connection({2, "gj1"}, {"gj1"}, 0.02),
+                    gap_junction_connection({2, "gj2"}, {"gj3"}, 0.01)
+                };
+            case 2:
+                return {
+                    gap_junction_connection({0, "gj1"}, {"gj0"}, 0.01),
+                    gap_junction_connection({1, "gj1"}, {"gj1"}, 0.02),
+                    gap_junction_connection({1, "gj3"}, {"gj2"}, 0.01)
+                };
+            default : return {};
+            }
+            return conns;
+        }
+        std::any get_global_properties(cell_kind) const override {
+            return gprop_;
+        }
+    protected:
+        arb::cable_cell_global_properties gprop_;
+        std::vector<cable_cell> cells_;
+        cell_size_type n_ = 3;
+    };
+
+    soma_cell_builder b0(2.1);
+    b0.add_branch(0, 8, 0.3, 0.2, 4, "dend");
+
+    auto c0 = b0.make_cell();
+    mlocation c0_gj[2] = {b0.location({1, 1}), b0.location({1, 0.5})};
+
+    c0.decorations.place(c0_gj[0], junction{"gj"}, "gj0");
+    c0.decorations.place(c0_gj[1], junction{"gj"}, "gj1");
+
+    soma_cell_builder b1(1.4);
+    b1.add_branch(0, 12, 0.3, 0.5, 6, "dend");
+    b1.add_branch(1,  9, 0.3, 0.2, 3, "dend");
+    b1.add_branch(1,  5, 0.2, 0.2, 5, "dend");
+
+    auto c1 = b1.make_cell();
+    mlocation c1_gj[4] = {b1.location({2, 1}), b1.location({1, 1}), b1.location({1, 0.45}), b1.location({1, 0.1})};
+
+    c1.decorations.place(c1_gj[0], junction{"gj"}, "gj0");
+    c1.decorations.place(c1_gj[1], junction{"gj"}, "gj1");
+    c1.decorations.place(c1_gj[2], junction{"gj"}, "gj2");
+    c1.decorations.place(c1_gj[3], junction{"gj"}, "gj3");
+
+    soma_cell_builder b2(2.9);
+    b2.add_branch(0, 4, 0.3, 0.5, 2, "dend");
+    b2.add_branch(1, 6, 0.4, 0.2, 2, "dend");
+    b2.add_branch(1, 8, 0.1, 0.2, 2, "dend");
+    b2.add_branch(2, 4, 0.2, 0.2, 2, "dend");
+    b2.add_branch(2, 4, 0.2, 0.2, 2, "dend");
+
+    auto c2 = b2.make_cell();
+    mlocation c2_gj[3] = {b2.location({1, 0.5}), b2.location({4, 1}), b2.location({2, 1})};
+
+    c2.decorations.place(c2_gj[0], junction{"gj"}, "gj0");
+    c2.decorations.place(c2_gj[1], junction{"gj"}, "gj1");
+    c2.decorations.place(c2_gj[2], junction{"gj"}, "gj2");
+
+    // Check the GJ CV map
+    cable_cell_global_properties gprop;
+    gprop.default_parameters = neuron_parameter_defaults;
+
+    std::vector<cable_cell> cells{c0, c1, c2};
+    std::vector<cell_gid_type> gids = {0, 1, 2};
+
+    auto D = fvm_cv_discretize(cells, neuron_parameter_defaults, context);
+    unsigned c0_gj_cv[2], c1_gj_cv[4], c2_gj_cv[3];
+    for (int i = 0; i<2; ++i) c0_gj_cv[i] = D.geometry.location_cv(0, c0_gj[i], cv_prefer::cv_nonempty);
+    for (int i = 0; i<4; ++i) c1_gj_cv[i] = D.geometry.location_cv(1, c1_gj[i], cv_prefer::cv_nonempty);
+    for (int i = 0; i<3; ++i) c2_gj_cv[i] = D.geometry.location_cv(2, c2_gj[i], cv_prefer::cv_nonempty);
+
+    auto gj_cvs = fvm_build_gap_junction_cv_map(cells, gids, D);
+
+    EXPECT_EQ(9u, gj_cvs.size());
+    EXPECT_EQ(c0_gj_cv[0], gj_cvs.at(cell_member_type{0, 0}));
+    EXPECT_EQ(c0_gj_cv[1], gj_cvs.at(cell_member_type{0, 1}));
+    EXPECT_EQ(c1_gj_cv[0], gj_cvs.at(cell_member_type{1, 0}));
+    EXPECT_EQ(c1_gj_cv[1], gj_cvs.at(cell_member_type{1, 1}));
+    EXPECT_EQ(c1_gj_cv[2], gj_cvs.at(cell_member_type{1, 2}));
+    EXPECT_EQ(c1_gj_cv[3], gj_cvs.at(cell_member_type{1, 3}));
+    EXPECT_EQ(c2_gj_cv[0], gj_cvs.at(cell_member_type{2, 0}));
+    EXPECT_EQ(c2_gj_cv[1], gj_cvs.at(cell_member_type{2, 1}));
+    EXPECT_EQ(c2_gj_cv[2], gj_cvs.at(cell_member_type{2, 2}));
+
+    // Check the resolved GJ connections
+    fvm_cell fvcell(context);
+    gap_recipe rec(cells, gprop);
+
+    auto fvm_info = fvcell.initialize(gids, rec);
+    auto gj_conns = fvm_resolve_gj_connections(gids, fvm_info.gap_junction_data, gj_cvs, rec);
+
+    std::unordered_map<cell_gid_type, std::vector<fvm_gap_junction>> expected;
+    expected[0] = {
+        {0, c0_gj_cv[0], c1_gj_cv[0], 0.03},
+        {0, c0_gj_cv[0], c1_gj_cv[1], 0.04},
+        {1, c0_gj_cv[1], c2_gj_cv[0], 0.01},
+    };
+    expected[1] = {
+        {0, c1_gj_cv[0], c0_gj_cv[0], 0.03},
+        {1, c1_gj_cv[1], c0_gj_cv[0], 0.04},
+        {1, c1_gj_cv[1], c2_gj_cv[1], 0.02},
+        {3, c1_gj_cv[3], c2_gj_cv[2], 0.01}
+    };
+    expected[2] = {
+        {0, c2_gj_cv[0], c0_gj_cv[1], 0.01},
+        {1, c2_gj_cv[1], c1_gj_cv[1], 0.02},
+        {2, c2_gj_cv[2], c1_gj_cv[3], 0.01}
+    };
+
+    util::sort(expected.at(0));
+    util::sort(expected.at(1));
+    util::sort(expected.at(2));
+
+    EXPECT_EQ(expected.at(0), gj_conns.at(0));
+    EXPECT_EQ(expected.at(1), gj_conns.at(1));
+    EXPECT_EQ(expected.at(2), gj_conns.at(2));
+
+    // Check the GJ mechanism data
+    auto M = fvm_build_mechanism_data(gprop, cells, gids, gj_conns, D, context);
+
+    EXPECT_EQ(1u, M.mechanisms.size());
+    ASSERT_EQ(1u, M.mechanisms.count("gj"));
+
+    const auto& gj_data = M.mechanisms["gj"];
+    const auto& gj_g_param = *ptr_by_key(gj_data.param_values, "g"s);
+
+    std::vector<std::tuple<int, int, double, double>> expected_gj_values = {
+        {c0_gj_cv[0], c1_gj_cv[0], 0.03, 1.},
+        {c0_gj_cv[0], c1_gj_cv[1], 0.04, 1.},
+        {c0_gj_cv[1], c2_gj_cv[0], 0.01, 1.},
+        {c1_gj_cv[0], c0_gj_cv[0], 0.03, 1.},
+        {c1_gj_cv[1], c0_gj_cv[0], 0.04, 1.},
+        {c1_gj_cv[1], c2_gj_cv[1], 0.02, 1.},
+        {c1_gj_cv[3], c2_gj_cv[2], 0.01, 1.},
+        {c2_gj_cv[0], c0_gj_cv[1], 0.01, 1.},
+        {c2_gj_cv[1], c1_gj_cv[1], 0.02, 1.},
+        {c2_gj_cv[2], c1_gj_cv[3], 0.01, 1.}
+    };
+    util::sort(expected_gj_values);
+
+    std::vector<int> expected_gj_cv, expected_gj_peer_cv;
+    std::vector<double> expected_gj_weight, expected_gj_g;
+    util::assign(expected_gj_cv,      util::transform_view(expected_gj_values, [](const auto& x){return std::get<0>(x);}));
+    util::assign(expected_gj_peer_cv, util::transform_view(expected_gj_values, [](const auto& x){return std::get<1>(x);}));
+    util::assign(expected_gj_weight,  util::transform_view(expected_gj_values, [](const auto& x){return std::get<2>(x);}));
+    util::assign(expected_gj_g,       util::transform_view(expected_gj_values, [](const auto& x){return std::get<3>(x);}));
+
+    EXPECT_EQ(10u, gj_data.cv.size());
+    EXPECT_EQ(expected_gj_cv,      gj_data.cv);
+    EXPECT_EQ(expected_gj_peer_cv, gj_data.peer_cv);
+    EXPECT_EQ(expected_gj_weight,  gj_data.local_weight);
+    EXPECT_EQ(expected_gj_g,       gj_g_param);
+}
+
+TEST(fvm_layout, gj_example_2) {
+    arb::proc_allocation resources;
+    if (auto nt = arbenv::get_env_num_threads()) {
+        resources.num_threads = nt;
+    } else {
+        resources.num_threads = arbenv::thread_concurrency();
+    }
+    arb::execution_context context(resources);
+
+    class gap_recipe: public recipe {
+    public:
+        gap_recipe(std::vector<cable_cell> cells, arb::cable_cell_global_properties gprop) : gprop_(gprop), cells_(cells) {}
+
+        cell_size_type num_cells() const override { return n_; }
+        cell_kind get_cell_kind(cell_gid_type) const override { return cell_kind::cable; }
+        util::unique_any get_cell_description(cell_gid_type gid) const override {
+            return cells_[gid];
+        }
+        std::vector<arb::gap_junction_connection> gap_junctions_on(cell_gid_type gid) const override{
+            std::vector<gap_junction_connection> conns;
+            switch (gid) {
+            case 0:
+                return {};
+            case 1:
+                return {
+                    gap_junction_connection({2, "j5"}, {"j2"}, 0.03),
+                    gap_junction_connection({2, "j6"}, {"j3"}, 0.04),
+                    gap_junction_connection({3, "j7"}, {"j0"}, 0.02),
+                    gap_junction_connection({5, "j8"}, {"j1"}, 0.01),
+                    gap_junction_connection({5, "j9"}, {"j4"}, 0.01)
+                };
+            case 2:
+                return {
+                    gap_junction_connection({1, "j2"}, {"j5"}, 0.03),
+                    gap_junction_connection({1, "j3"}, {"j6"}, 0.04),
+                    gap_junction_connection({3, "j7"}, {"j5"}, 0.02),
+                    gap_junction_connection({5, "j9"}, {"j6"}, 0.01),
+                };
+            case 3:
+                return {
+                    gap_junction_connection({1, "j0"}, {"j7"}, 0.03),
+                    gap_junction_connection({2, "j5"}, {"j7"}, 0.04),
+                    gap_junction_connection({5, "j8"}, {"j7"}, 0.02),
+                };
+            case 4:
+                return {};
+            case 5:
+                return {
+                    gap_junction_connection({1, "j1"}, {"j8"}, 0.03),
+                    gap_junction_connection({1, "j4"}, {"j9"}, 0.04),
+                    gap_junction_connection({2, "j6"}, {"j9"}, 0.02),
+                    gap_junction_connection({3, "j7"}, {"j8"}, 0.01),
+                };
+            default : return {};
+            }
+            return conns;
+        }
+        std::any get_global_properties(cell_kind) const override {
+            return gprop_;
+        }
+    protected:
+        arb::cable_cell_global_properties gprop_;
+        std::vector<cable_cell> cells_;
+        cell_size_type n_ = 6;
+    };
+
+    auto system = six_cell_gj_system();
+    auto& desc = system.descriptions;
+    auto& builders = system.builders;
+
+    std::vector<mlocation> locs_1 = {
+        builders[1].location({2, 0}),
+        builders[1].location({3, 0.3}),
+        builders[1].location({0, 0.5}),
+        builders[1].location({0, 0.6}),
+        builders[1].location({3, 0.8})};
+    desc[1].decorations.place(locs_1[0], junction("gj1", {{"g", 0.3}, {"e", 1e-3}}), "j0");
+    desc[1].decorations.place(locs_1[1], junction("gj0", {{"g", 0.2}}), "j1");
+    desc[1].decorations.place(locs_1[2], junction("gj0", {{"g", 1.5}}), "j2");
+    desc[1].decorations.place(locs_1[3], junction("gj0"), "j3");
+    desc[1].decorations.place(locs_1[4], junction("gj0"), "j4");
+
+    std::vector<mlocation> locs_2 = {
+        builders[2].location({0, 0.5}),
+        builders[2].location({1, 0})};
+    desc[2].decorations.place(locs_2[0], junction("gj0", {{"g", 2.3}}), "j5");
+    desc[2].decorations.place(locs_2[1], junction("gj1", {{"e", 5e-4}}), "j6");
+
+    std::vector<mlocation> locs_3 = {
+        builders[3].location({0, 0.5})};
+    desc[3].decorations.place(locs_3[0], junction("gj1", {{"g", 0.4}}), "j7");
+
+    std::vector<mlocation> locs_5 = {
+        builders[5].location({0, 0.1}),
+        builders[5].location({0, 0.2})};
+    desc[5].decorations.place(locs_5[0], junction("gj1"), "j8");
+    desc[5].decorations.place(locs_5[1], junction("gj0", {{"g", 1.6}}), "j9");
+
+    // Check the GJ CV map
+    cable_cell_global_properties gprop;
+    auto cat = make_unit_test_catalogue();
+    cat.import(arb::global_default_catalogue(), "");
+    gprop.catalogue = &cat;
+    gprop.default_parameters = neuron_parameter_defaults;
+
+    auto cells = system.cells();
+    std::vector<cell_gid_type> gids = {0, 1, 2, 3, 4, 5};
+
+    fvm_cv_discretization D = fvm_cv_discretize(cells, gprop.default_parameters);
+
+    unsigned cvs_1[5], cvs_2[2], cvs_3[1], cvs_5[2];
+    for (int i = 0; i<5; ++i) cvs_1[i] = D.geometry.location_cv(1, locs_1[i], cv_prefer::cv_nonempty);
+    for (int i = 0; i<2; ++i) cvs_2[i] = D.geometry.location_cv(2, locs_2[i], cv_prefer::cv_nonempty);
+    for (int i = 0; i<1; ++i) cvs_3[i] = D.geometry.location_cv(3, locs_3[i], cv_prefer::cv_nonempty);
+    for (int i = 0; i<2; ++i) cvs_5[i] = D.geometry.location_cv(5, locs_5[i], cv_prefer::cv_nonempty);
+
+    auto gj_cvs = fvm_build_gap_junction_cv_map(cells, gids, D);
+
+    EXPECT_EQ(10u, gj_cvs.size());
+    EXPECT_EQ(cvs_1[0], gj_cvs.at(cell_member_type{1, 0}));
+    EXPECT_EQ(cvs_1[1], gj_cvs.at(cell_member_type{1, 1}));
+    EXPECT_EQ(cvs_1[2], gj_cvs.at(cell_member_type{1, 2}));
+    EXPECT_EQ(cvs_1[3], gj_cvs.at(cell_member_type{1, 3}));
+    EXPECT_EQ(cvs_1[4], gj_cvs.at(cell_member_type{1, 4}));
+    EXPECT_EQ(cvs_2[0], gj_cvs.at(cell_member_type{2, 0}));
+    EXPECT_EQ(cvs_2[1], gj_cvs.at(cell_member_type{2, 1}));
+    EXPECT_EQ(cvs_3[0], gj_cvs.at(cell_member_type{3, 0}));
+    EXPECT_EQ(cvs_5[0], gj_cvs.at(cell_member_type{5, 0}));
+    EXPECT_EQ(cvs_5[1], gj_cvs.at(cell_member_type{5, 1}));
+
+    // Check the resolved GJ connections
+    fvm_cell fvcell(context);
+    gap_recipe rec(cells, gprop);
+
+    auto fvm_info = fvcell.initialize(gids, rec);
+    auto gj_conns = fvm_resolve_gj_connections(gids, fvm_info.gap_junction_data, gj_cvs, rec);
+
+    std::unordered_map<cell_gid_type, std::vector<fvm_gap_junction>> expected;
+    expected[0] = {};
+    expected[1] = {
+        {2, cvs_1[2], cvs_2[0], 0.03},
+        {3, cvs_1[3], cvs_2[1], 0.04},
+        {0, cvs_1[0], cvs_3[0], 0.02},
+        {1, cvs_1[1], cvs_5[0], 0.01},
+        {4, cvs_1[4], cvs_5[1], 0.01},
+    };
+    expected[2] = {
+        {0, cvs_2[0], cvs_1[2], 0.03},
+        {1, cvs_2[1], cvs_1[3], 0.04},
+        {0, cvs_2[0], cvs_3[0], 0.02},
+        {1, cvs_2[1], cvs_5[1], 0.01}
+    };
+    expected[3] = {
+        {0, cvs_3[0], cvs_1[0], 0.03},
+        {0, cvs_3[0], cvs_2[0], 0.04},
+        {0, cvs_3[0], cvs_5[0], 0.02}
+    };
+    expected[4] = {};
+    expected[5] = {
+        {0, cvs_5[0], cvs_1[1], 0.03},
+        {1, cvs_5[1], cvs_1[4], 0.04},
+        {1, cvs_5[1], cvs_2[1], 0.02},
+        {0, cvs_5[0], cvs_3[0], 0.01}
+    };
+
+    util::sort(expected.at(1));
+    util::sort(expected.at(2));
+    util::sort(expected.at(3));
+    util::sort(expected.at(5));
+
+    EXPECT_EQ(expected.at(0), gj_conns.at(0));
+    EXPECT_EQ(expected.at(1), gj_conns.at(1));
+    EXPECT_EQ(expected.at(2), gj_conns.at(2));
+    EXPECT_EQ(expected.at(3), gj_conns.at(3));
+    EXPECT_EQ(expected.at(4), gj_conns.at(4));
+    EXPECT_EQ(expected.at(5), gj_conns.at(5));
+
+    // Check the GJ mechanism data
+    auto M = fvm_build_mechanism_data(gprop, cells, gids, gj_conns, D, context);
+
+    EXPECT_EQ(4u, M.mechanisms.size());
+    ASSERT_EQ(1u, M.mechanisms.count("gj0"));
+    ASSERT_EQ(1u, M.mechanisms.count("gj1"));
+
+    const auto& gj0_data = M.mechanisms["gj0"];
+    const auto& gj1_data = M.mechanisms["gj1"];
+
+    const auto& gj0_params = gj0_data.param_values;
+    const auto& gj1_params = gj1_data.param_values;
+
+    const auto& gj0_g_param = *ptr_by_key(gj0_params, "g"s);
+    const auto& gj1_g_param = *ptr_by_key(gj1_params, "g"s);
+    const auto& gj1_e_param = *ptr_by_key(gj1_params, "e"s);
+
+    std::vector<std::tuple<int, int, double, double>> expected_gj0_values = {
+        {cvs_1[2], cvs_2[0], 0.03, 1.5},
+        {cvs_1[3], cvs_2[1], 0.04, 1. },
+        {cvs_1[1], cvs_5[0], 0.01, 0.2},
+        {cvs_1[4], cvs_5[1], 0.01, 1.0},
+        {cvs_2[0], cvs_1[2], 0.03, 2.3},
+        {cvs_2[0], cvs_3[0], 0.02, 2.3},
+        {cvs_5[1], cvs_1[4], 0.04, 1.6},
+        {cvs_5[1], cvs_2[1], 0.02, 1.6}
+    };
+    util::sort(expected_gj0_values);
+
+    std::vector<int> expected_gj0_cv, expected_gj0_peer_cv;
+    std::vector<double> expected_gj0_weight, expected_gj0_g;
+    util::assign(expected_gj0_cv,      util::transform_view(expected_gj0_values, [](const auto& x){return std::get<0>(x);}));
+    util::assign(expected_gj0_peer_cv, util::transform_view(expected_gj0_values, [](const auto& x){return std::get<1>(x);}));
+    util::assign(expected_gj0_weight,  util::transform_view(expected_gj0_values, [](const auto& x){return std::get<2>(x);}));
+    util::assign(expected_gj0_g,       util::transform_view(expected_gj0_values, [](const auto& x){return std::get<3>(x);}));
+
+    EXPECT_EQ(8u, gj0_data.cv.size());
+    EXPECT_EQ(expected_gj0_cv,      gj0_data.cv);
+    EXPECT_EQ(expected_gj0_peer_cv, gj0_data.peer_cv);
+    EXPECT_EQ(expected_gj0_weight,  gj0_data.local_weight);
+    EXPECT_EQ(expected_gj0_g,       gj0_g_param);
+
+    std::vector<std::tuple<int, int, double, double, double>> expected_gj1_values = {
+        {cvs_1[0], cvs_3[0], 0.02, 0.3, 1e-3},
+        {cvs_2[1], cvs_1[3], 0.04, 1.0, 5e-4},
+        {cvs_2[1], cvs_5[1], 0.01, 1.0, 5e-4},
+        {cvs_3[0], cvs_1[0], 0.03, 0.4, 0.},
+        {cvs_3[0], cvs_2[0], 0.04, 0.4, 0.},
+        {cvs_3[0], cvs_5[0], 0.02, 0.4, 0.},
+        {cvs_5[0], cvs_1[1], 0.03, 1.0, 0.},
+        {cvs_5[0], cvs_3[0], 0.01, 1.0, 0.}
+    };
+    util::sort(expected_gj1_values);
+
+    std::vector<int> expected_gj1_cv, expected_gj1_peer_cv;
+    std::vector<double> expected_gj1_weight, expected_gj1_g, expected_gj1_e;
+    util::assign(expected_gj1_cv,      util::transform_view(expected_gj1_values, [](const auto& x){return std::get<0>(x);}));
+    util::assign(expected_gj1_peer_cv, util::transform_view(expected_gj1_values, [](const auto& x){return std::get<1>(x);}));
+    util::assign(expected_gj1_weight,  util::transform_view(expected_gj1_values, [](const auto& x){return std::get<2>(x);}));
+    util::assign(expected_gj1_g,       util::transform_view(expected_gj1_values, [](const auto& x){return std::get<3>(x);}));
+    util::assign(expected_gj1_e,       util::transform_view(expected_gj1_values, [](const auto& x){return std::get<4>(x);}));
+
+    EXPECT_EQ(8u, gj1_data.cv.size());
+    EXPECT_EQ(expected_gj1_cv,      gj1_data.cv);
+    EXPECT_EQ(expected_gj1_peer_cv, gj1_data.peer_cv);
+    EXPECT_EQ(expected_gj1_weight,  gj1_data.local_weight);
+    EXPECT_EQ(expected_gj1_g,       gj1_g_param);
+    EXPECT_EQ(expected_gj1_e,       gj1_e_param);
+}
+
+TEST(fvm_lowered, cell_group_gj) {
+    arb::proc_allocation resources;
+    if (auto nt = arbenv::get_env_num_threads()) {
+        resources.num_threads = nt;
+    }
+    else {
+        resources.num_threads = arbenv::thread_concurrency();
+    }
+    arb::execution_context context(resources);
+
+    class gap_recipe: public recipe {
+    public:
+        gap_recipe(const std::vector<cable_cell>& cg0, const std::vector<cable_cell>& cg1) {
+            cells_ = cg0;
+            cells_.insert(cells_.end(), cg1.begin(), cg1.end());
+            gprop_.default_parameters = neuron_parameter_defaults;
+        }
+
+        cell_size_type num_cells() const override { return n_; }
+        cell_kind get_cell_kind(cell_gid_type) const override { return cell_kind::cable; }
+        util::unique_any get_cell_description(cell_gid_type gid) const override {
+            return cells_[gid];
+        }
+        std::vector<arb::gap_junction_connection> gap_junctions_on(cell_gid_type gid) const override{
+            std::vector<gap_junction_connection> conns;
+            if (gid % 2 == 0) {
+                // connect 5 of the first 10 cells in a ring; connect 5 of the second 10 cells in a ring
+                auto next_cell = gid == 8 ? 0 : (gid == 18 ? 10 : gid + 2);
+                auto prev_cell = gid == 0 ? 8 : (gid == 10 ? 18 : gid - 2);
+                conns.push_back(gap_junction_connection({next_cell, "gj", lid_selection_policy::assert_univalent},
+                                                        {"gj", lid_selection_policy::assert_univalent}, 0.03));
+                conns.push_back(gap_junction_connection({prev_cell, "gj", lid_selection_policy::assert_univalent},
+                                                        {"gj", lid_selection_policy::assert_univalent}, 0.03));
+            }
+            return conns;
+        }
+        std::any get_global_properties(cell_kind) const override {
+            return gprop_;
+        }
+
+    protected:
+        arb::cable_cell_global_properties gprop_;
+        std::vector<cable_cell> cells_;
+        cell_size_type n_ = 20;
+    };
+
+    std::vector<cable_cell> cell_group0;
+    std::vector<cable_cell> cell_group1;
+
+    // Make 20 cells
+    for (unsigned i = 0; i < 20; i++) {
+        cable_cell_description c = soma_cell_builder(2.1).make_cell();
+        if (i % 2 == 0) {
+            c.decorations.place(mlocation{0, 1}, junction{"gj"}, "gj");
+        }
+        if (i < 10) {
+            cell_group0.push_back(c);
+        }
+        else {
+            cell_group1.push_back(c);
+        }
+    }
+
+    std::vector<cell_gid_type> gids_cg0 = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
+    std::vector<cell_gid_type> gids_cg1 = {10, 11, 12, 13, 14, 15, 16, 17, 18, 19};
+
+    gap_recipe rec(cell_group0, cell_group1);
+
+    fvm_cell fvcell0(context);
+    fvm_cell fvcell1(context);
+
+    auto fvm_info_0 = fvcell0.initialize(gids_cg0, rec);
+    auto fvm_info_1 = fvcell1.initialize(gids_cg1, rec);
+
+    auto num_dom0 = fvcell0.fvm_intdom(rec, gids_cg0, fvm_info_0.cell_to_intdom);
+    auto num_dom1 = fvcell1.fvm_intdom(rec, gids_cg1, fvm_info_1.cell_to_intdom);
+
+    fvm_cv_discretization D0 = fvm_cv_discretize(cell_group0, neuron_parameter_defaults, context);
+    fvm_cv_discretization D1 = fvm_cv_discretize(cell_group1, neuron_parameter_defaults, context);
+
+    auto gj_cvs_0 = fvm_build_gap_junction_cv_map(cell_group0, gids_cg0, D0);
+    auto gj_cvs_1 = fvm_build_gap_junction_cv_map(cell_group1, gids_cg1, D1);
+
+    auto GJ0 = fvm_resolve_gj_connections(gids_cg0, fvm_info_0.gap_junction_data, gj_cvs_0, rec);
+    auto GJ1 = fvm_resolve_gj_connections(gids_cg1, fvm_info_1.gap_junction_data, gj_cvs_1, rec);
+
+    EXPECT_EQ(gids_cg0.size(), GJ0.size());
+    EXPECT_EQ(gids_cg1.size(), GJ1.size());
+
+    std::vector<std::vector<fvm_gap_junction>> expected = {
+        {{0, 0, 2, 0.03}, {0, 0, 8, 0.03}}, {},
+        {{0, 2, 0, 0.03}, {0, 2, 4, 0.03}}, {},
+        {{0, 4, 2, 0.03} ,{0, 4, 6, 0.03}}, {},
+        {{0, 6, 4, 0.03}, {0, 6, 8, 0.03}}, {},
+        {{0, 8, 0, 0.03}, {0, 8, 6, 0.03}}, {}
+    };
+
+    for (unsigned i = 0; i < GJ0.size(); i++) {
+        EXPECT_EQ(expected[i], GJ0[i]);
+        EXPECT_EQ(expected[i], GJ1[i+10]);
+    }
+
+    std::vector<fvm_index_type> expected_doms= {0u, 1u, 0u, 2u, 0u, 3u, 0u, 4u, 0u, 5u};
+    EXPECT_EQ(6u, num_dom0);
+    EXPECT_EQ(6u, num_dom1);
+
+    EXPECT_EQ(expected_doms, fvm_info_0.cell_to_intdom);
+    EXPECT_EQ(expected_doms, fvm_info_1.cell_to_intdom);
+}
+
 namespace {
     double wm_impl(double wa, double xa) {
         return wa? xa/wa: 0;
@@ -590,10 +1328,10 @@ TEST(fvm_layout, density_norm_area) {
     hh_3["gl"] = seg3_gl;
 
     auto desc = builder.make_cell();
-    desc.decorations.paint("soma"_lab, std::move(hh_0));
-    desc.decorations.paint("reg1"_lab, std::move(hh_1));
-    desc.decorations.paint("reg2"_lab, std::move(hh_2));
-    desc.decorations.paint("reg3"_lab, std::move(hh_3));
+    desc.decorations.paint("soma"_lab, density(std::move(hh_0)));
+    desc.decorations.paint("reg1"_lab, density(std::move(hh_1)));
+    desc.decorations.paint("reg2"_lab, density(std::move(hh_2)));
+    desc.decorations.paint("reg3"_lab, density(std::move(hh_3)));
 
     std::vector<cable_cell> cells{desc};
 
@@ -636,8 +1374,10 @@ TEST(fvm_layout, density_norm_area) {
     cable_cell_global_properties gprop;
     gprop.default_parameters = neuron_parameter_defaults;
 
+    std::vector<cell_gid_type> gids = {0};
+    std::unordered_map<cell_gid_type, std::vector<fvm_gap_junction>> gj_conns = {{0, {}}};
     fvm_cv_discretization D = fvm_cv_discretize(cells, gprop.default_parameters);
-    fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, D);
+    fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, gids, gj_conns, D);
 
     // Grab the HH parameters from the mechanism.
 
@@ -690,8 +1430,8 @@ TEST(fvm_layout, density_norm_area_partial) {
     auto desc = builder.make_cell();
     desc.decorations.set_default(cv_policy_fixed_per_branch(1));
 
-    desc.decorations.paint(builder.cable({1, 0., 0.3}), hh_begin);
-    desc.decorations.paint(builder.cable({1, 0.4, 1.}), hh_end);
+    desc.decorations.paint(builder.cable({1, 0., 0.3}), density(hh_begin));
+    desc.decorations.paint(builder.cable({1, 0.4, 1.}), density(hh_end));
 
     std::vector<cable_cell> cells{desc};
 
@@ -710,8 +1450,10 @@ TEST(fvm_layout, density_norm_area_partial) {
     cable_cell_global_properties gprop;
     gprop.default_parameters = neuron_parameter_defaults;
 
+    std::vector<cell_gid_type> gids = {0};
+    std::unordered_map<cell_gid_type, std::vector<fvm_gap_junction>> gj_conns = {{0, {}}};
     fvm_cv_discretization D = fvm_cv_discretize(cells, gprop.default_parameters);
-    fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, D);
+    fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, gids, gj_conns, D);
 
     // Grab the HH parameters from the mechanism.
 
@@ -739,8 +1481,10 @@ TEST(fvm_layout, density_norm_area_partial) {
 
 TEST(fvm_layout, valence_verify) {
     auto desc = soma_cell_builder(6).make_cell();
-    desc.decorations.paint("soma"_lab, "test_cl_valence");
+    desc.decorations.paint("soma"_lab, density("test_cl_valence"));
     std::vector<cable_cell> cells{desc};
+    std::vector<cell_gid_type> gids = {0};
+    std::unordered_map<cell_gid_type, std::vector<fvm_gap_junction>> gj_conns = {{0, {}}};
 
     cable_cell_global_properties gprop;
     gprop.default_parameters = neuron_parameter_defaults;
@@ -751,16 +1495,16 @@ TEST(fvm_layout, valence_verify) {
     gprop.catalogue = &testcat;
 
     // Missing the 'cl' ion:
-    EXPECT_THROW(fvm_build_mechanism_data(gprop, cells, D), cable_cell_error);
+    EXPECT_THROW(fvm_build_mechanism_data(gprop, cells, gids, gj_conns, D), cable_cell_error);
 
     // Adding ion, should be fine now:
     gprop.default_parameters.ion_data["cl"] = { 1., 1., 0. };
     gprop.ion_species["cl"] = -1;
-    EXPECT_NO_THROW(fvm_build_mechanism_data(gprop, cells, D));
+    EXPECT_NO_THROW(fvm_build_mechanism_data(gprop, cells, gids, gj_conns, D));
 
     // 'cl' ion has wrong charge:
     gprop.ion_species["cl"] = -2;
-    EXPECT_THROW(fvm_build_mechanism_data(gprop, cells, D), cable_cell_error);
+    EXPECT_THROW(fvm_build_mechanism_data(gprop, cells, gids, gj_conns, D), cable_cell_error);
 }
 
 TEST(fvm_layout, ion_weights) {
@@ -832,13 +1576,14 @@ TEST(fvm_layout, ion_weights) {
 
         for (auto i: mech_branches[run]) {
             auto cab = builder.cable({i, 0, 1});
-            desc.decorations.paint(reg::cable(cab.branch, cab.prox_pos, cab.dist_pos), "test_ca");
+            desc.decorations.paint(reg::cable(cab.branch, cab.prox_pos, cab.dist_pos), density("test_ca"));
         }
 
         std::vector<cable_cell> cells{desc};
-
+        std::vector<cell_gid_type> gids = {0};
+        std::unordered_map<cell_gid_type, std::vector<fvm_gap_junction>> gj_conns = {{0, {}}};
         fvm_cv_discretization D = fvm_cv_discretize(cells, gprop.default_parameters);
-        fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, D);
+        fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, gids, gj_conns, D);
 
         ASSERT_EQ(1u, M.ions.count("ca"s));
         auto& ca = M.ions.at("ca"s);
@@ -868,9 +1613,9 @@ TEST(fvm_layout, revpot) {
     builder.add_branch(1, 200, 0.5, 0.5, 1, "dend");
     builder.add_branch(1, 100, 0.5, 0.5, 1, "dend");
     auto desc = builder.make_cell();
-    desc.decorations.paint("soma"_lab, "read_eX/c");
-    desc.decorations.paint("soma"_lab, "read_eX/a");
-    desc.decorations.paint("dend"_lab, "read_eX/a");
+    desc.decorations.paint("soma"_lab, density("read_eX/c"));
+    desc.decorations.paint("soma"_lab, density("read_eX/a"));
+    desc.decorations.paint("dend"_lab, density("read_eX/a"));
 
     std::vector<cable_cell_description> descriptions{desc, desc};
 
@@ -892,8 +1637,10 @@ TEST(fvm_layout, revpot) {
         test_gprop.default_parameters.reversal_potential_method["b"] = write_eb_ec;
 
         std::vector<cable_cell> cells{descriptions[0], descriptions[1]};
+        std::vector<cell_gid_type> gids = {0,1};
+        std::unordered_map<cell_gid_type, std::vector<fvm_gap_junction>> gj_conns = {{0, {}}, {1, {}}};
         fvm_cv_discretization D = fvm_cv_discretize(cells, test_gprop.default_parameters);
-        EXPECT_THROW(fvm_build_mechanism_data(test_gprop, cells, D), cable_cell_error);
+        EXPECT_THROW(fvm_build_mechanism_data(test_gprop, cells, gids, gj_conns, D), cable_cell_error);
     }
 
     {
@@ -902,10 +1649,12 @@ TEST(fvm_layout, revpot) {
         test_gprop.default_parameters.reversal_potential_method["b"] = write_eb_ec;
         test_gprop.default_parameters.reversal_potential_method["c"] = write_eb_ec;
         descriptions[1].decorations.set_default(ion_reversal_potential_method{"c", "write_eX/c"});
-        std::vector<cable_cell> cells{descriptions[0], descriptions[1]};
 
+        std::vector<cable_cell> cells{descriptions[0], descriptions[1]};
+        std::vector<cell_gid_type> gids = {0,1};
+        std::unordered_map<cell_gid_type, std::vector<fvm_gap_junction>> gj_conns = {{0, {}}, {1, {}}};
         fvm_cv_discretization D = fvm_cv_discretize(cells, test_gprop.default_parameters);
-        EXPECT_THROW(fvm_build_mechanism_data(test_gprop, cells, D), cable_cell_error);
+        EXPECT_THROW(fvm_build_mechanism_data(test_gprop, cells, gids, gj_conns, D), cable_cell_error);
     }
 
     {
@@ -915,8 +1664,10 @@ TEST(fvm_layout, revpot) {
         descriptions[1].decorations.set_default(ion_reversal_potential_method{"c", write_eb_ec});
 
         std::vector<cable_cell> cells{descriptions[0], descriptions[1]};
+        std::vector<cell_gid_type> gids = {0,1};
+        std::unordered_map<cell_gid_type, std::vector<fvm_gap_junction>> gj_conns = {{0, {}}, {1, {}}};
         fvm_cv_discretization D = fvm_cv_discretize(cells, gprop.default_parameters);
-        fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, D);
+        fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, gids, gj_conns, D);
 
         // Only CV which needs write_multiple_eX/x=b,y=c is the soma (first CV)
         // of the second cell.
diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp
index d8295388..cee09538 100644
--- a/test/unit/test_fvm_lowered.cpp
+++ b/test/unit/test_fvm_lowered.cpp
@@ -80,7 +80,7 @@ public:
 
     arb::util::unique_any get_cell_description(cell_gid_type gid) const override {
         auto c = soma_cell_builder(20).make_cell();
-        c.decorations.place(mlocation{0, 1}, gap_junction_site{}, "gj");
+        c.decorations.place(mlocation{0, 1}, junction("gj"), "gj");
         return {cable_cell{c}};
     }
 
@@ -145,7 +145,7 @@ public:
 
     arb::util::unique_any get_cell_description(cell_gid_type) const override {
         auto c = soma_cell_builder(20).make_cell();
-        c.decorations.place(mlocation{0,1}, gap_junction_site{}, "gj");
+        c.decorations.place(mlocation{0,1}, junction("gj"), "gj");
         return {cable_cell{c}};
     }
 
@@ -247,10 +247,10 @@ TEST(fvm_lowered, target_handles) {
     };
 
     // (in increasing target order)
-    descriptions[0].decorations.place(mlocation{0, 0.7}, "expsyn", "syn0");
-    descriptions[0].decorations.place(mlocation{0, 0.3}, "expsyn", "syn1");
-    descriptions[1].decorations.place(mlocation{2, 0.2}, "exp2syn", "syn2");
-    descriptions[1].decorations.place(mlocation{2, 0.8}, "expsyn", "syn3");
+    descriptions[0].decorations.place(mlocation{0, 0.7}, synapse("expsyn"), "syn0");
+    descriptions[0].decorations.place(mlocation{0, 0.3}, synapse("expsyn"), "syn1");
+    descriptions[1].decorations.place(mlocation{2, 0.2}, synapse("exp2syn"), "syn2");
+    descriptions[1].decorations.place(mlocation{2, 0.8}, synapse("expsyn"), "syn3");
 
     descriptions[1].decorations.place(mlocation{0, 0}, threshold_detector{3.3}, "detector");
 
@@ -457,14 +457,14 @@ TEST(fvm_lowered, derived_mechs) {
 
         switch (i) {
             case 0:
-                cell.decorations.paint(reg::all(), "test_kin1");
+                cell.decorations.paint(reg::all(), density("test_kin1"));
                 break;
             case 1:
-                cell.decorations.paint(reg::all(), "custom_kin1");
+                cell.decorations.paint(reg::all(), density("custom_kin1"));
                 break;
             case 2:
-                cell.decorations.paint(reg::all(), "test_kin1");
-                cell.decorations.paint(reg::all(), "custom_kin1");
+                cell.decorations.paint(reg::all(), density("test_kin1"));
+                cell.decorations.paint(reg::all(), density("custom_kin1"));
                 break;
         }
         cells.push_back(cell);
@@ -548,8 +548,8 @@ TEST(fvm_lowered, null_region) {
     builder.add_branch(0, 100, 0.5, 0.5, 4, "dend");
     auto cell = builder.make_cell();
 
-    cell.decorations.paint(reg::nil(), "test_kin1");
-    cell.decorations.paint(reg::nil(), "custom_kin1");
+    cell.decorations.paint(reg::nil(), density("test_kin1"));
+    cell.decorations.paint(reg::nil(), density("custom_kin1"));
 
     cable1d_recipe rec(cable_cell{cell});
     rec.catalogue() = make_unit_test_catalogue();
@@ -578,7 +578,7 @@ TEST(fvm_lowered, read_valence) {
 
         soma_cell_builder builder(6);
         auto cell = builder.make_cell();
-        cell.decorations.paint("soma"_lab, "test_ca_read_valence");
+        cell.decorations.paint("soma"_lab, density("test_ca_read_valence"));
         cable1d_recipe rec(cable_cell{cell});
         rec.catalogue() = make_unit_test_catalogue();
 
@@ -598,7 +598,7 @@ TEST(fvm_lowered, read_valence) {
         // Check ion renaming.
         soma_cell_builder builder(6);
         auto cell = builder.make_cell();
-        cell.decorations.paint("soma"_lab, "cr_read_valence");
+        cell.decorations.paint("soma"_lab, density("cr_read_valence"));
         cable1d_recipe rec(cable_cell{cell});
         rec.catalogue() = make_unit_test_catalogue();
         rec.catalogue() = make_unit_test_catalogue();
@@ -628,7 +628,6 @@ TEST(fvm_lowered, ionic_concentrations) {
     std::vector<fvm_value_type> temp(ncv, 23);
     std::vector<fvm_value_type> diam(ncv, 1.);
     std::vector<fvm_value_type> vinit(ncv, -65);
-    std::vector<fvm_gap_junction> gj = {};
     std::vector<fvm_index_type> src_to_spike = {};
 
     fvm_ion_config ion_config;
@@ -653,7 +652,7 @@ TEST(fvm_lowered, ionic_concentrations) {
     auto& write_cai_mech = write_cai.mech;
 
     auto shared_state = std::make_unique<typename backend::shared_state>(
-            ncell, ncell, 0, cv_to_intdom, cv_to_intdom, gj, vinit, temp, diam, src_to_spike, read_cai_mech->data_alignment());
+            ncell, ncell, 0, cv_to_intdom, cv_to_intdom, vinit, temp, diam, src_to_spike, read_cai_mech->data_alignment());
     shared_state->add_ion("ca", 2, ion_config);
 
     shared_state->instantiate(*read_cai_mech, 0, overrides, layout);
@@ -712,8 +711,8 @@ TEST(fvm_lowered, ionic_currents) {
     m2["coeff"] = coeff;
 
     auto c = b.make_cell();
-    c.decorations.paint("soma"_lab, m1);
-    c.decorations.paint("soma"_lab, m2);
+    c.decorations.paint("soma"_lab, density(m1));
+    c.decorations.paint("soma"_lab, density(m2));
 
     cable1d_recipe rec({cable_cell{c}});
     rec.catalogue() = make_unit_test_catalogue();
@@ -754,7 +753,7 @@ TEST(fvm_lowered, point_ionic_current) {
     double soma_area_m2 = 4*math::pi<double>*r*r*1e-12; // [m²]
 
     // Event weight is translated by point_ica_current into a current contribution in nA.
-    c.decorations.place(mlocation{0u, 0.5}, "point_ica_current", "syn");
+    c.decorations.place(mlocation{0u, 0.5}, synapse("point_ica_current"), "syn");
 
     cable1d_recipe rec({cable_cell{c}});
     rec.catalogue() = make_unit_test_catalogue();
@@ -827,10 +826,10 @@ TEST(fvm_lowered, weighted_write_ion) {
     const double con_ext = 120;
 
     // Ca ion reader test_kinlva on CV 2 and 3 via branch 2:
-    c.decorations.paint(reg::branch(1), "test_kinlva");
+    c.decorations.paint(reg::branch(1), density("test_kinlva"));
 
     // Ca ion writer test_ca on CV 2 and 4 via branch 3:
-    c.decorations.paint(reg::branch(2), "test_ca");
+    c.decorations.paint(reg::branch(2), density("test_ca"));
 
     cable1d_recipe rec({cable_cell{c}});
     rec.catalogue() = make_unit_test_catalogue();
@@ -874,332 +873,6 @@ TEST(fvm_lowered, weighted_write_ion) {
     EXPECT_TRUE(testing::seq_almost_eq<double>(expected_iconc, ion_iconc));
 }
 
-TEST(fvm_lowered, gj_coords_simple) {
-    arb::proc_allocation resources;
-    if (auto nt = arbenv::get_env_num_threads()) {
-        resources.num_threads = nt;
-    }
-    else {
-        resources.num_threads = arbenv::thread_concurrency();
-    }
-    arb::execution_context context(resources);
-
-    using pair = std::pair<int, int>;
-
-    class gap_recipe: public recipe {
-    public:
-        gap_recipe(std::vector<cable_cell> cells) : cells_(cells) {
-            gprop_.default_parameters = neuron_parameter_defaults;
-        }
-
-        cell_size_type num_cells() const override { return n_; }
-        cell_kind get_cell_kind(cell_gid_type) const override { return cell_kind::cable; }
-        util::unique_any get_cell_description(cell_gid_type gid) const override {
-            return cells_[gid];
-        }
-        std::vector<arb::gap_junction_connection> gap_junctions_on(cell_gid_type gid) const override{
-            std::vector<gap_junction_connection> conns;
-            conns.push_back(gap_junction_connection({(gid+1)%2, "gj", lid_selection_policy::assert_univalent}, {"gj", lid_selection_policy::assert_univalent}, 0.5));
-            return conns;
-        }
-        std::any get_global_properties(cell_kind) const override {
-            return gprop_;
-        }
-    protected:
-        arb::cable_cell_global_properties gprop_;
-        std::vector<cable_cell> cells_;
-        cell_size_type n_ = 2;
-    };
-
-    fvm_cell fvcell(context);
-
-    std::vector<cable_cell> cells;
-    {
-        soma_cell_builder b(2.1);
-        b.add_branch(0, 10, 0.3, 0.2, 5, "dend");
-        auto c = b.make_cell();
-        c.decorations.place(b.location({1, 0.8}), gap_junction_site{}, "gj");
-        cells.push_back(c);
-    }
-
-    {
-        soma_cell_builder b(2.4);
-        b.add_branch(0, 10, 0.3, 0.2, 2, "dend");
-        auto c = b.make_cell();
-        c.decorations.place(b.location({1, 1}), gap_junction_site{}, "gj");
-        cells.push_back(c);
-    }
-    gap_recipe rec(cells);
-
-    fvm_cv_discretization D = fvm_cv_discretize(cells, neuron_parameter_defaults, context);
-
-    std::vector<cell_gid_type> gids = {0, 1};
-    auto fvm_info = fvcell.initialize(gids, rec);
-
-    auto GJ = fvcell.fvm_gap_junctions(cells, gids, fvm_info.gap_junction_data, rec, D);
-
-    auto weight = [&](fvm_value_type g, fvm_index_type i){
-        return g * 1e3 / D.cv_area[i];
-    };
-
-    EXPECT_EQ(pair({5,10}), GJ[0].loc);
-    EXPECT_EQ(weight(0.5, 5), GJ[0].weight);
-
-    EXPECT_EQ(pair({10,5}), GJ[1].loc);
-    EXPECT_EQ(weight(0.5, 10), GJ[1].weight);
-}
-
-TEST(fvm_lowered, gj_coords_complex) {
-    arb::proc_allocation resources;
-    if (auto nt = arbenv::get_env_num_threads()) {
-        resources.num_threads = nt;
-    }
-    else {
-        resources.num_threads = arbenv::thread_concurrency();
-    }
-    arb::execution_context context(resources);
-
-    class gap_recipe: public recipe {
-    public:
-        gap_recipe(std::vector<cable_cell> cells) : cells_(cells) {
-            gprop_.default_parameters = neuron_parameter_defaults;
-        }
-
-        cell_size_type num_cells() const override { return n_; }
-        cell_kind get_cell_kind(cell_gid_type) const override { return cell_kind::cable; }
-        util::unique_any get_cell_description(cell_gid_type gid) const override {
-            return cells_[gid];
-        }
-        std::vector<arb::gap_junction_connection> gap_junctions_on(cell_gid_type gid) const override{
-            std::vector<gap_junction_connection> conns;
-            switch (gid) {
-            case 0:
-                return {
-                    gap_junction_connection({2, "gj0", lid_selection_policy::assert_univalent}, {"gj1", lid_selection_policy::assert_univalent}, 0.01),
-                    gap_junction_connection({1, "gj0", lid_selection_policy::assert_univalent}, {"gj0", lid_selection_policy::assert_univalent}, 0.03),
-                    gap_junction_connection({1, "gj1", lid_selection_policy::assert_univalent}, {"gj0", lid_selection_policy::assert_univalent}, 0.04)
-                };
-            case 1:
-                return {
-                    gap_junction_connection({0, "gj0", lid_selection_policy::assert_univalent}, {"gj0", lid_selection_policy::assert_univalent}, 0.03),
-                    gap_junction_connection({0, "gj0", lid_selection_policy::assert_univalent}, {"gj1", lid_selection_policy::assert_univalent}, 0.04),
-                    gap_junction_connection({2, "gj1", lid_selection_policy::assert_univalent}, {"gj2", lid_selection_policy::assert_univalent}, 0.02),
-                    gap_junction_connection({2, "gj2", lid_selection_policy::assert_univalent}, {"gj3", lid_selection_policy::assert_univalent}, 0.01)
-                };
-            case 2:
-                return {
-                    gap_junction_connection({0, "gj1", lid_selection_policy::assert_univalent}, {"gj0", lid_selection_policy::assert_univalent}, 0.01),
-                    gap_junction_connection({1, "gj2", lid_selection_policy::assert_univalent}, {"gj1", lid_selection_policy::assert_univalent}, 0.02),
-                    gap_junction_connection({1, "gj3", lid_selection_policy::assert_univalent}, {"gj2", lid_selection_policy::assert_univalent}, 0.01)
-                };
-            default : return {};
-            }
-            return conns;
-        }
-        std::any get_global_properties(cell_kind) const override {
-            return gprop_;
-        }
-    protected:
-        arb::cable_cell_global_properties gprop_;
-        std::vector<cable_cell> cells_;
-        cell_size_type n_ = 3;
-    };
-
-    // Add 5 gap junctions
-    soma_cell_builder b0(2.1);
-    b0.add_branch(0, 8, 0.3, 0.2, 4, "dend");
-
-    auto c0 = b0.make_cell();
-    mlocation c0_gj[2] = {b0.location({1, 1}), b0.location({1, 0.5})};
-
-    c0.decorations.place(c0_gj[0], gap_junction_site{}, "gj0");
-    c0.decorations.place(c0_gj[1], gap_junction_site{}, "gj1");
-
-    soma_cell_builder b1(1.4);
-    b1.add_branch(0, 12, 0.3, 0.5, 6, "dend");
-    b1.add_branch(1,  9, 0.3, 0.2, 3, "dend");
-    b1.add_branch(1,  5, 0.2, 0.2, 5, "dend");
-
-    auto c1 = b1.make_cell();
-    mlocation c1_gj[4] = {b1.location({2, 1}), b1.location({1, 1}), b1.location({1, 0.45}), b1.location({1, 0.1})};
-
-    c1.decorations.place(c1_gj[0], gap_junction_site{}, "gj0");
-    c1.decorations.place(c1_gj[1], gap_junction_site{}, "gj1");
-    c1.decorations.place(c1_gj[2], gap_junction_site{}, "gj2");
-    c1.decorations.place(c1_gj[3], gap_junction_site{}, "gj3");
-
-
-    soma_cell_builder b2(2.9);
-    b2.add_branch(0, 4, 0.3, 0.5, 2, "dend");
-    b2.add_branch(1, 6, 0.4, 0.2, 2, "dend");
-    b2.add_branch(1, 8, 0.1, 0.2, 2, "dend");
-    b2.add_branch(2, 4, 0.2, 0.2, 2, "dend");
-    b2.add_branch(2, 4, 0.2, 0.2, 2, "dend");
-
-    auto c2 = b2.make_cell();
-    mlocation c2_gj[3] = {b2.location({1, 0.5}), b2.location({4, 1}), b2.location({2, 1})};
-
-    c2.decorations.place(c2_gj[0], gap_junction_site{}, "gj0");
-    c2.decorations.place(c2_gj[1], gap_junction_site{}, "gj1");
-    c2.decorations.place(c2_gj[2], gap_junction_site{}, "gj2");
-
-    std::vector<cable_cell> cells{c0, c1, c2};
-    std::vector<cell_gid_type> gids = {0, 1, 2};
-
-    gap_recipe rec(cells);
-    fvm_cell fvcell(context);
-
-    auto fvm_info = fvcell.initialize(gids, rec);
-    fvcell.fvm_intdom(rec, gids, fvm_info.cell_to_intdom);
-    fvm_cv_discretization D = fvm_cv_discretize(cells, neuron_parameter_defaults, context);
-
-    using namespace cv_prefer;
-    int c0_gj_cv[2];
-    for (int i = 0; i<2; ++i) c0_gj_cv[i] = D.geometry.location_cv(0, c0_gj[i], cv_nonempty);
-
-    int c1_gj_cv[4];
-    for (int i = 0; i<4; ++i) c1_gj_cv[i] = D.geometry.location_cv(1, c1_gj[i], cv_nonempty);
-
-    int c2_gj_cv[3];
-    for (int i = 0; i<3; ++i) c2_gj_cv[i] = D.geometry.location_cv(2, c2_gj[i], cv_nonempty);
-
-    std::vector<fvm_gap_junction> GJ = fvcell.fvm_gap_junctions(cells, gids, fvm_info.gap_junction_data, rec, D);
-    EXPECT_EQ(10u, GJ.size());
-
-    auto weight = [&](fvm_value_type g, fvm_index_type i){
-        return g * 1e3 / D.cv_area[i];
-    };
-
-    std::vector<fvm_gap_junction> expected = {
-        {{c0_gj_cv[0], c1_gj_cv[0]}, weight(0.03, c0_gj_cv[0])},
-        {{c0_gj_cv[0], c1_gj_cv[1]}, weight(0.04, c0_gj_cv[0])},
-        {{c0_gj_cv[1], c2_gj_cv[0]}, weight(0.01, c0_gj_cv[1])},
-        {{c1_gj_cv[0], c0_gj_cv[0]}, weight(0.03, c1_gj_cv[0])},
-        {{c1_gj_cv[1], c0_gj_cv[0]}, weight(0.04, c1_gj_cv[1])},
-        {{c1_gj_cv[2], c2_gj_cv[1]}, weight(0.02, c1_gj_cv[2])},
-        {{c1_gj_cv[3], c2_gj_cv[2]}, weight(0.01, c1_gj_cv[3])},
-        {{c2_gj_cv[0], c0_gj_cv[1]}, weight(0.01, c2_gj_cv[0])},
-        {{c2_gj_cv[1], c1_gj_cv[2]}, weight(0.02, c2_gj_cv[1])},
-        {{c2_gj_cv[2], c1_gj_cv[3]}, weight(0.01, c2_gj_cv[2])}
-    };
-
-    using util::sort_by;
-    using util::transform_view;
-
-    auto gj_loc = [](const fvm_gap_junction gj) { return gj.loc; };
-    auto gj_weight = [](const fvm_gap_junction gj) { return gj.weight; };
-
-    sort_by(GJ, [](fvm_gap_junction gj) { return gj.loc; });
-    sort_by(expected, [](fvm_gap_junction gj) { return gj.loc; });
-
-    EXPECT_TRUE(testing::seq_eq(transform_view(expected, gj_loc), transform_view(GJ, gj_loc)));
-    EXPECT_TRUE(testing::seq_almost_eq<double>(transform_view(expected, gj_weight), transform_view(GJ, gj_weight)));
-}
-
-TEST(fvm_lowered, cell_group_gj) {
-    arb::proc_allocation resources;
-    if (auto nt = arbenv::get_env_num_threads()) {
-        resources.num_threads = nt;
-    }
-    else {
-        resources.num_threads = arbenv::thread_concurrency();
-    }
-    arb::execution_context context(resources);
-
-    using pair = std::pair<int, int>;
-
-    class gap_recipe: public recipe {
-    public:
-        gap_recipe(const std::vector<cable_cell>& cg0, const std::vector<cable_cell>& cg1) {
-            cells_ = cg0;
-            cells_.insert(cells_.end(), cg1.begin(), cg1.end());
-            gprop_.default_parameters = neuron_parameter_defaults;
-        }
-
-        cell_size_type num_cells() const override { return n_; }
-        cell_kind get_cell_kind(cell_gid_type) const override { return cell_kind::cable; }
-        util::unique_any get_cell_description(cell_gid_type gid) const override {
-            return cells_[gid];
-        }
-        std::vector<arb::gap_junction_connection> gap_junctions_on(cell_gid_type gid) const override{
-            std::vector<gap_junction_connection> conns;
-            if (gid % 2 == 0) {
-                // connect 5 of the first 10 cells in a ring; connect 5 of the second 10 cells in a ring
-                auto next_cell = gid == 8 ? 0 : (gid == 18 ? 10 : gid + 2);
-                auto prev_cell = gid == 0 ? 8 : (gid == 10 ? 18 : gid - 2);
-                conns.push_back(gap_junction_connection({next_cell, "gj", lid_selection_policy::assert_univalent},
-                                                             {"gj", lid_selection_policy::assert_univalent}, 0.03));
-                conns.push_back(gap_junction_connection({prev_cell, "gj", lid_selection_policy::assert_univalent},
-                                                             {"gj", lid_selection_policy::assert_univalent}, 0.03));
-            }
-            return conns;
-        }
-        std::any get_global_properties(cell_kind) const override {
-            return gprop_;
-        }
-
-    protected:
-        arb::cable_cell_global_properties gprop_;
-        std::vector<cable_cell> cells_;
-        cell_size_type n_ = 20;
-    };
-
-    std::vector<cable_cell> cell_group0;
-    std::vector<cable_cell> cell_group1;
-
-    // Make 20 cells
-    for (unsigned i = 0; i < 20; i++) {
-        cable_cell_description c = soma_cell_builder(2.1).make_cell();
-        if (i % 2 == 0) {
-            c.decorations.place(mlocation{0, 1}, gap_junction_site{}, "gj");
-        }
-        if (i < 10) {
-            cell_group0.push_back(c);
-        }
-        else {
-            cell_group1.push_back(c);
-        }
-    }
-
-    std::vector<cell_gid_type> gids_cg0 = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
-    std::vector<cell_gid_type> gids_cg1 = {10, 11, 12, 13, 14, 15, 16, 17, 18, 19};
-
-    gap_recipe rec(cell_group0, cell_group1);
-
-    fvm_cell fvcell0(context);
-    fvm_cell fvcell1(context);
-
-    auto fvm_info_0 = fvcell0.initialize(gids_cg0, rec);
-    auto fvm_info_1 = fvcell1.initialize(gids_cg1, rec);
-
-    auto num_dom0 = fvcell0.fvm_intdom(rec, gids_cg0, fvm_info_0.cell_to_intdom);
-    auto num_dom1 = fvcell1.fvm_intdom(rec, gids_cg1, fvm_info_1.cell_to_intdom);
-
-    fvm_cv_discretization D0 = fvm_cv_discretize(cell_group0, neuron_parameter_defaults, context);
-    fvm_cv_discretization D1 = fvm_cv_discretize(cell_group1, neuron_parameter_defaults, context);
-
-    auto GJ0 = fvcell0.fvm_gap_junctions(cell_group0, gids_cg0, fvm_info_0.gap_junction_data, rec, D0);
-    auto GJ1 = fvcell1.fvm_gap_junctions(cell_group1, gids_cg1, fvm_info_1.gap_junction_data, rec, D1);
-
-    EXPECT_EQ(10u, GJ0.size());
-    EXPECT_EQ(10u, GJ1.size());
-
-    std::vector<pair> expected_loc = {{0, 2}, {0, 8}, {2, 4}, {2, 0}, {4, 6} ,{4, 2}, {6, 8}, {6, 4}, {8, 0}, {8, 6}};
-
-    for (unsigned i = 0; i < GJ0.size(); i++) {
-        EXPECT_EQ(expected_loc[i], GJ0[i].loc);
-        EXPECT_EQ(expected_loc[i], GJ1[i].loc);
-    }
-
-    std::vector<fvm_index_type> expected_doms= {0u, 1u, 0u, 2u, 0u, 3u, 0u, 4u, 0u, 5u};
-    EXPECT_EQ(6u, num_dom0);
-    EXPECT_EQ(6u, num_dom1);
-
-    EXPECT_EQ(expected_doms, fvm_info_0.cell_to_intdom);
-    EXPECT_EQ(expected_doms, fvm_info_1.cell_to_intdom);
-
-}
-
 TEST(fvm_lowered, integration_domains) {
     {
         execution_context context;
@@ -1299,7 +972,7 @@ TEST(fvm_lowered, post_events_shared_state) {
         unsigned ncell_;
         unsigned ncv_;
         std::vector<unsigned> detectors_per_cell_;
-        arb::mechanism_desc synapse_;
+        arb::synapse synapse_;
         mechanism_catalogue cat_;
     };
 
@@ -1375,8 +1048,8 @@ TEST(fvm_lowered, label_data) {
             {
                 arb::decor decor;
                 decor.set_default(arb::cv_policy_fixed_per_branch(10));
-                decor.place(uniform(all(), 0, 3, 42), "expsyn", "4_synapses");
-                decor.place(uniform(all(), 4, 4, 42), "expsyn", "1_synapse");
+                decor.place(uniform(all(), 0, 3, 42), arb::synapse("expsyn"), "4_synapses");
+                decor.place(uniform(all(), 4, 4, 42), arb::synapse("expsyn"), "1_synapse");
                 decor.place(uniform(all(), 5, 5, 42), arb::threshold_detector{10}, "1_detector");
 
                 cells_.push_back(arb::cable_cell(arb::morphology(tree), {}, decor));
@@ -1386,8 +1059,8 @@ TEST(fvm_lowered, label_data) {
                 decor.set_default(arb::cv_policy_fixed_per_branch(10));
                 decor.place(uniform(all(), 0, 2, 24), arb::threshold_detector{10}, "3_detectors");
                 decor.place(uniform(all(), 3, 4, 24), arb::threshold_detector{10}, "2_detectors");
-                decor.place(uniform(all(), 5, 6, 24), arb::gap_junction_site(), "2_gap_junctions");
-                decor.place(uniform(all(), 7, 7, 24), arb::gap_junction_site(), "1_gap_junction");
+                decor.place(uniform(all(), 5, 6, 24), arb::junction("gj"), "2_gap_junctions");
+                decor.place(uniform(all(), 7, 7, 24), arb::junction("gj"), "1_gap_junction");
 
                 cells_.push_back(arb::cable_cell(arb::morphology(tree), {}, decor));
             }
diff --git a/test/unit/test_kinetic_linear.cpp b/test/unit/test_kinetic_linear.cpp
index 07b2f130..45a1330f 100644
--- a/test/unit/test_kinetic_linear.cpp
+++ b/test/unit/test_kinetic_linear.cpp
@@ -37,7 +37,6 @@ void run_test(std::string mech_name,
     fvm_size_type ncv = 1;
     std::vector<fvm_index_type> cv_to_intdom(ncv, 0);
 
-    std::vector<fvm_gap_junction> gj = {};
     auto instance = cat.instance(backend::kind, mech_name);
     auto& test = instance.mech;
 
@@ -47,7 +46,7 @@ void run_test(std::string mech_name,
     std::vector<fvm_index_type> src_to_spike = {};
 
     auto shared_state = std::make_unique<typename backend::shared_state>(
-            ncell, ncell, 0, cv_to_intdom, cv_to_intdom, gj, vinit, temp, diam, src_to_spike, test->data_alignment());
+            ncell, ncell, 0, cv_to_intdom, cv_to_intdom, vinit, temp, diam, src_to_spike, test->data_alignment());
 
     mechanism_layout layout;
     mechanism_overrides overrides;
diff --git a/test/unit/test_mc_cell_group.cpp b/test/unit/test_mc_cell_group.cpp
index 035c7d67..37449682 100644
--- a/test/unit/test_mc_cell_group.cpp
+++ b/test/unit/test_mc_cell_group.cpp
@@ -27,8 +27,8 @@ namespace {
         soma_cell_builder builder(12.6157/2.0);
         builder.add_branch(0, 200, 0.5, 0.5, 101, "dend");
         auto d = builder.make_cell();
-        d.decorations.paint("soma"_lab, "hh");
-        d.decorations.paint("dend"_lab, "pas");
+        d.decorations.paint("soma"_lab, density("hh"));
+        d.decorations.paint("dend"_lab, density("pas"));
         d.decorations.place(builder.location({1,1}), i_clamp::box(5, 80, 0.3), "clamp0");
         d.decorations.place(builder.location({0, 0}), threshold_detector{0}, "detector0");
         return d;
diff --git a/test/unit/test_mc_cell_group_gpu.cpp b/test/unit/test_mc_cell_group_gpu.cpp
index 1170c80b..d2bbcd18 100644
--- a/test/unit/test_mc_cell_group_gpu.cpp
+++ b/test/unit/test_mc_cell_group_gpu.cpp
@@ -27,8 +27,8 @@ namespace {
         soma_cell_builder builder(12.6157/2.0);
         builder.add_branch(0, 200, 0.5, 0.5, 101, "dend");
         auto d = builder.make_cell();
-        d.decorations.paint("soma"_lab, "hh");
-        d.decorations.paint("dend"_lab, "pas");
+        d.decorations.paint("soma"_lab, density("hh"));
+        d.decorations.paint("dend"_lab, density("pas"));
         d.decorations.place(builder.location({1,1}), i_clamp::box(5, 80, 0.3), "clamp0");
         d.decorations.place(builder.location({0, 0}), threshold_detector{0}, "detector0");
         return d;
diff --git a/test/unit/test_mech_temp_diam.cpp b/test/unit/test_mech_temp_diam.cpp
index cb872d84..49a720aa 100644
--- a/test/unit/test_mech_temp_diam.cpp
+++ b/test/unit/test_mech_temp_diam.cpp
@@ -24,7 +24,6 @@ void run_celsius_test() {
     fvm_size_type ncv = 3;
     std::vector<fvm_index_type> cv_to_intdom(ncv, 0);
 
-    std::vector<fvm_gap_junction> gj = {};
     auto instance = cat.instance(backend::kind, "celsius_test");
     auto& celsius_test = instance.mech;
 
@@ -37,7 +36,7 @@ void run_celsius_test() {
     std::vector<fvm_index_type> src_to_spike = {};
 
     auto shared_state = std::make_unique<typename backend::shared_state>(
-        ncell, ncell, 0, cv_to_intdom, cv_to_intdom, gj, vinit, temp, diam, src_to_spike, celsius_test->data_alignment());
+        ncell, ncell, 0, cv_to_intdom, cv_to_intdom, vinit, temp, diam, src_to_spike, celsius_test->data_alignment());
 
     mechanism_layout layout;
     mechanism_overrides overrides;
@@ -75,7 +74,6 @@ void run_diam_test() {
     fvm_size_type ncv = 3;
     std::vector<fvm_index_type> cv_to_intdom(ncv, 0);
 
-    std::vector<fvm_gap_junction> gj = {};
     auto instance = cat.instance(backend::kind, "diam_test");
     auto& celsius_test = instance.mech;
 
@@ -95,7 +93,7 @@ void run_diam_test() {
     }
 
     auto shared_state = std::make_unique<typename backend::shared_state>(
-            ncell, ncell, 0, cv_to_intdom, cv_to_intdom, gj, vinit, temp, diam, src_to_spike, celsius_test->data_alignment());
+            ncell, ncell, 0, cv_to_intdom, cv_to_intdom, vinit, temp, diam, src_to_spike, celsius_test->data_alignment());
 
 
     shared_state->instantiate(*celsius_test, 0, overrides, layout);
diff --git a/test/unit/test_mechcat.cpp b/test/unit/test_mechcat.cpp
index 4d3d0cd3..05077ca2 100644
--- a/test/unit/test_mechcat.cpp
+++ b/test/unit/test_mechcat.cpp
@@ -402,7 +402,7 @@ TEST(mechcat, instantiate) {
     // write its specialized global variables to shared state, but we do in
     // these tests for testing purposes.
 
-    mechanism_layout layout = {{0u, 1u, 2u}, {1., 2., 1.}, {1u, 1u, 1u}};
+    mechanism_layout layout = {{0u, 1u, 2u}, {}, {1., 2., 1.}, {1u, 1u, 1u}};
     bar_backend bar;
 
     auto cat = build_fake_catalogue();
diff --git a/test/unit/test_probe.cpp b/test/unit/test_probe.cpp
index 13fc9889..dca9870b 100644
--- a/test/unit/test_probe.cpp
+++ b/test/unit/test_probe.cpp
@@ -272,8 +272,8 @@ void run_expsyn_g_probe_test(const context& ctx) {
     builder.add_branch(0, 200, 1.0/2, 1.0/2, 1, "dend");
     builder.add_branch(0, 200, 1.0/2, 1.0/2, 1, "dend");
     auto bs = builder.make_cell();
-    bs.decorations.place(loc0, "expsyn", "syn0");
-    bs.decorations.place(loc1, "expsyn", "syn1");
+    bs.decorations.place(loc0, synapse("expsyn"), "syn0");
+    bs.decorations.place(loc1, synapse("expsyn"), "syn1");
     bs.decorations.set_default(cv_policy_fixed_per_branch(2));
 
     auto run_test = [&](bool coalesce_synapses) {
@@ -376,9 +376,9 @@ void run_expsyn_g_cell_probe_test(const context& ctx) {
         for (unsigned j = 0; j<10; ++j) {
             auto idx = (bid*10+j)*2;
             mlocation expsyn_loc{bid, 0.1*j};
-            d.place(expsyn_loc, "expsyn", "syn"+std::to_string(idx));
+            d.place(expsyn_loc, synapse("expsyn"), "syn"+std::to_string(idx));
             expsyn_target_loc_map[2*n_expsyn] = expsyn_loc;
-            d.place(mlocation{bid, 0.1*j+0.05}, "exp2syn", "syn"+std::to_string(idx+1));
+            d.place(mlocation{bid, 0.1*j+0.05}, synapse("exp2syn"), "syn"+std::to_string(idx+1));
             ++n_expsyn;
         }
     }
@@ -503,9 +503,9 @@ void run_ion_density_probe_test(const context& ctx) {
     // Calcium ions everywhere, half written by write_ca1, half by write_ca2.
     // Sodium ions only on distal half.
 
-    d.paint(mcable{0, 0., 0.5}, "write_ca1");
-    d.paint(mcable{0, 0.5, 1.}, "write_ca2");
-    d.paint(mcable{0, 0.5, 1.}, "write_na3");
+    d.paint(mcable{0, 0., 0.5}, density("write_ca1"));
+    d.paint(mcable{0, 0.5, 1.}, density("write_ca2"));
+    d.paint(mcable{0, 0.5, 1.}, density("write_na3"));
 
     // Place probes in each CV.
 
@@ -680,7 +680,7 @@ void run_partial_density_probe_test(const context& ctx) {
     //    CV 4:   8.6
     //    CV 5:  10.5
 
-    auto mk_mech = [](double param) { return mechanism_desc("param_as_state").set("p", param); };
+    auto mk_mech = [](double param) { return density(mechanism_desc("param_as_state").set("p", param)); };
 
     d0.paint(mcable{0, 0.0, 0.1}, mk_mech(2));
     d0.paint(mcable{0, 0.2, 0.3}, mk_mech(3));
@@ -781,7 +781,7 @@ void run_axial_and_ion_current_sampled_probe_test(const context& ctx) {
     // For τ = 0.1 ms, set conductance to 0.01 S/cm² and membrance capacitance
     // to 0.01 F/m².
 
-    d.paint(reg::all(), mechanism_desc("ca_linear").set("g", 0.01)); // [S/cm²]
+    d.paint(reg::all(), density("ca_linear", {{"g", 0.01}})); // [S/cm²]
     d.set_default(membrane_capacitance{0.01}); // [F/m²]
     const double tau = 0.1; // [ms]
 
@@ -933,9 +933,9 @@ void run_multi_probe_test(const context& ctx) {
     decor d;
 
     // Paint mechanism on branches 1, 2, and 5, omitting branch 4.
-    d.paint(reg::branch(1), mechanism_desc("param_as_state").set("p", 10.));
-    d.paint(reg::branch(2), mechanism_desc("param_as_state").set("p", 20.));
-    d.paint(reg::branch(5), mechanism_desc("param_as_state").set("p", 50.));
+    d.paint(reg::branch(1), density("param_as_state", {{"p", 10.}}));
+    d.paint(reg::branch(2), density("param_as_state", {{"p", 20.}}));
+    d.paint(reg::branch(5), density("param_as_state", {{"p", 50.}}));
 
     auto tracev = run_simple_sampler<double, mlocation>(ctx, 0.1, {cable_cell{m, {}, d}}, 0, cable_probe_density_state{ls::terminal(), "param_as_state", "s"}, {0.});
 
@@ -1025,7 +1025,7 @@ void run_total_current_probe_test(const context& ctx) {
     const double tau = 0.1;     // [ms]
     d0.place(mlocation{0, 0}, i_clamp(0.3), "clamp0");
 
-    d0.paint(reg::all(), mechanism_desc("ca_linear").set("g", 0.01)); // [S/cm²]
+    d0.paint(reg::all(), density("ca_linear", {{"g", 0.01}})); // [S/cm²]
     d0.set_default(membrane_capacitance{0.01}); // [F/m²]
     // Tweak membrane capacitance on cells[1] so as to change dynamics a bit.
     auto d1 = d0;
@@ -1198,13 +1198,13 @@ void run_exact_sampling_probe_test(const context& ctx) {
             std::vector<cable_cell_description> cd;
             cd.assign(4, builder.make_cell());
 
-            cd[0].decorations.place(mlocation{1, 0.1}, "expsyn", "syn");
-            cd[1].decorations.place(mlocation{1, 0.1}, "exp2syn", "syn");
-            cd[2].decorations.place(mlocation{1, 0.9}, "expsyn", "syn");
-            cd[3].decorations.place(mlocation{1, 0.9}, "exp2syn", "syn");
+            cd[0].decorations.place(mlocation{1, 0.1}, synapse("expsyn"), "syn");
+            cd[1].decorations.place(mlocation{1, 0.1}, synapse("exp2syn"), "syn");
+            cd[2].decorations.place(mlocation{1, 0.9}, synapse("expsyn"), "syn");
+            cd[3].decorations.place(mlocation{1, 0.9}, synapse("exp2syn"), "syn");
 
-            cd[1].decorations.place(mlocation{1, 0.2}, gap_junction_site{}, "gj");
-            cd[3].decorations.place(mlocation{1, 0.2}, gap_junction_site{}, "gj");
+            cd[1].decorations.place(mlocation{1, 0.2}, junction("gj"), "gj");
+            cd[3].decorations.place(mlocation{1, 0.2}, junction("gj"), "gj");
 
             for (auto& d: cd) cells_.push_back(d);
         }
@@ -1341,9 +1341,9 @@ TEST(probe, get_probe_metadata) {
     decor d;
 
     // Paint mechanism on branches 1, 2, and 5, omitting branch 4.
-    d.paint(reg::branch(1), mechanism_desc("param_as_state").set("p", 10.));
-    d.paint(reg::branch(2), mechanism_desc("param_as_state").set("p", 20.));
-    d.paint(reg::branch(5), mechanism_desc("param_as_state").set("p", 50.));
+    d.paint(reg::branch(1), density("param_as_state", {{"p", 10.}}));
+    d.paint(reg::branch(2), density("param_as_state", {{"p", 20.}}));
+    d.paint(reg::branch(5), density("param_as_state", {{"p", 50.}}));
 
     cable1d_recipe rec(cable_cell{m, {}, d}, false);
     rec.catalogue() = make_unit_test_catalogue(global_default_catalogue());
diff --git a/test/unit/test_recipe.cpp b/test/unit/test_recipe.cpp
index 0ecd2e33..37d1a7c7 100644
--- a/test/unit/test_recipe.cpp
+++ b/test/unit/test_recipe.cpp
@@ -80,12 +80,12 @@ namespace {
 
         // Add a num_synapses synapses to the cell.
         for (auto i: util::make_span(num_synapses)) {
-            decorations.place(arb::mlocation{0,(double)i/num_synapses}, "expsyn", "synapse"+std::to_string(i));
+            decorations.place(arb::mlocation{0,(double)i/num_synapses}, arb::synapse("expsyn"), "synapse"+std::to_string(i));
         }
 
         // Add a num_gj gap_junctions to the cell.
         for (auto i: util::make_span(num_gj)) {
-            decorations.place(arb::mlocation{0,(double)i/num_gj}, arb::gap_junction_site{}, "gapjunction"+std::to_string(i));
+            decorations.place(arb::mlocation{0,(double)i/num_gj}, arb::junction("gj"), "gapjunction"+std::to_string(i));
         }
 
         return arb::cable_cell(tree, {}, decorations);
diff --git a/test/unit/test_s_expr.cpp b/test/unit/test_s_expr.cpp
index 7309e247..f8aa256b 100644
--- a/test/unit/test_s_expr.cpp
+++ b/test/unit/test_s_expr.cpp
@@ -464,9 +464,6 @@ std::ostream& operator<<(std::ostream& o, const i_clamp& c) {
 std::ostream& operator<<(std::ostream& o, const threshold_detector& p) {
     return o << "(threshold-detector " << p.threshold << ')';
 }
-std::ostream& operator<<(std::ostream& o, const gap_junction_site& p) {
-    return o << "(gap-junction-site)";
-}
 std::ostream& operator<<(std::ostream& o, const init_membrane_potential& p) {
     return o << "(membrane-potential " << p.value << ')';
 }
@@ -495,6 +492,15 @@ std::ostream& operator<<(std::ostream& o, const mechanism_desc& m) {
     }
     return o << ')';
 }
+std::ostream& operator<<(std::ostream& o, const junction& p) {
+    return o << "(junction " << p.mech << ')';
+}
+std::ostream& operator<<(std::ostream& o, const synapse& p) {
+    return o << "(synapse " << p.mech << ')';
+}
+std::ostream& operator<<(std::ostream& o, const density& p) {
+    return o << "(density " << p.mech << ')';
+}
 std::ostream& operator<<(std::ostream& o, const ion_reversal_potential_method& p) {
     return o << "(ion-reversal-potential-method \"" << p.ion << "\" " << p.method << ')';
 }
@@ -604,8 +610,8 @@ TEST(decor_literals, round_tripping) {
         "(ion-external-concentration \"h\" -50.1)",
         "(ion-reversal-potential \"na\" 30)"};
     auto paint_literals = {
-        "(mechanism \"hh\")",
-        "(mechanism \"pas\" (\"g\" 0.02))",
+        "(density (mechanism \"hh\"))",
+        "(density (mechanism \"pas\" (\"g\" 0.02)))",
     };
     auto default_literals = {
         "(ion-reversal-potential-method \"ca\" (mechanism \"nernst/ca\"))",
@@ -614,8 +620,8 @@ TEST(decor_literals, round_tripping) {
     auto place_literals = {
         "(current-clamp (envelope (10 0.5) (110 0.5) (110 0)) 10 0.25)",
         "(threshold-detector -10)",
-        "(gap-junction-site)",
-        "(mechanism \"expsyn\")"};
+        "(junction (mechanism \"gj\"))",
+        "(synapse (mechanism \"expsyn\"))"};
     for (auto l: paint_default_literals) {
         EXPECT_EQ(l, round_trip_variant<defaultable>(l));
         EXPECT_EQ(l, round_trip_variant<paintable>(l));
@@ -659,8 +665,8 @@ TEST(decor_expressions, round_tripping) {
         "(paint (radius-gt (tag 3) 1) (ion-internal-concentration \"ca\" 75.1))",
         "(paint (intersect (cable 2 0 0.5) (region \"axon\")) (ion-external-concentration \"h\" -50.1))",
         "(paint (region \"my_region\") (ion-reversal-potential \"na\" 30))",
-        "(paint (cable 2 0.1 0.4) (mechanism \"hh\"))",
-        "(paint (all) (mechanism \"pas\" (\"g\" 0.02)))"
+        "(paint (cable 2 0.1 0.4) (density (mechanism \"hh\")))",
+        "(paint (all) (density (mechanism \"pas\" (\"g\" 0.02))))"
     };
     auto decorate_default_literals = {
         "(default (membrane-potential -65.1))",
@@ -676,8 +682,8 @@ TEST(decor_expressions, round_tripping) {
     auto decorate_place_literals = {
         "(place (location 3 0.2) (current-clamp (envelope (10 0.5) (110 0.5) (110 0)) 0.5 0.25) \"clamp\")",
         "(place (terminal) (threshold-detector -10) \"detector\")",
-        "(place (root) (gap-junction-site) \"gap_junction\")",
-        "(place (locset \"my!ls\") (mechanism \"expsyn\") \"synapse\")"};
+        "(place (root) (junction (mechanism \"gj\")) \"gap_junction\")",
+        "(place (locset \"my!ls\") (synapse (mechanism \"expsyn\")) \"synapse\")"};
 
     for (auto l: decorate_paint_literals) {
         EXPECT_EQ(l, round_trip<paint_pair>(l));
@@ -744,10 +750,12 @@ TEST(decor, round_tripping) {
                                 "          1)))\n"
                                 "    (paint \n"
                                 "      (region \"dend\")\n"
-                                "      (mechanism \"pas\"))\n"
+                                "      (density \n"
+                                "        (mechanism \"pas\")))\n"
                                 "    (paint \n"
                                 "      (region \"soma\")\n"
-                                "      (mechanism \"hh\"))\n"
+                                "      (density \n"
+                                "        (mechanism \"hh\")))\n"
                                 "    (paint \n"
                                 "      (join \n"
                                 "        (tag 1)\n"
@@ -755,7 +763,8 @@ TEST(decor, round_tripping) {
                                 "      (ion-internal-concentration \"ca\" 0.500000))\n"
                                 "    (place \n"
                                 "      (location 0 0)\n"
-                                "      (gap-junction-site)\n"
+                                "      (junction \n"
+                                "        (mechanism \"gj\"))\n"
                                 "      \"gap-junction\")\n"
                                 "    (place \n"
                                 "      (location 0 0)\n"
@@ -763,8 +772,9 @@ TEST(decor, round_tripping) {
                                 "      \"detector\")\n"
                                 "    (place \n"
                                 "      (location 0 0.5)\n"
-                                "      (mechanism \"expsyn\" \n"
-                                "        (\"tau\" 1.500000))\n"
+                                "      (synapse \n"
+                                "        (mechanism \"expsyn\" \n"
+                                "          (\"tau\" 1.500000)))\n"
                                 "      \"synapse\")))";
 
     EXPECT_EQ(component_str, round_trip_component(component_str.c_str()));
@@ -961,11 +971,13 @@ TEST(cable_cell, round_tripping) {
                                 "    (decor \n"
                                 "      (paint \n"
                                 "        (region \"dend\")\n"
-                                "        (mechanism \"pas\"))\n"
+                                "        (density \n"
+                                "          (mechanism \"pas\")))\n"
                                 "      (paint \n"
                                 "        (region \"soma\")\n"
-                                "        (mechanism \"hh\" \n"
-                                "          (\"el\" 0.500000)))\n"
+                                "        (density \n"
+                                "          (mechanism \"hh\" \n"
+                                "            (\"el\" 0.500000))))\n"
                                 "      (place \n"
                                 "        (location 0 1)\n"
                                 "        (current-clamp \n"
@@ -988,19 +1000,19 @@ TEST(cable_cell_literals, errors) {
                      "(membrane-capacitance ",       // syntax error
                      "(mem-capacitance 3.5)",        // invalid function
                      "(ion-initial-concentration ca 0.1)",   // unquoted ion
-                     "(mechanism \"hh\" (gl 3.5))",          // unqouted parameter
-                     "(mechanism \"pas\" ((\"g\" 0.5) (\"e\" 0.2)))",   // paranthesis around params
-                     "(mechanism \"pas\" (\"g\" 0.5 0.1) (\"e\" 0.2))", // too many values
-                     "(gap-junction-site 0)",                // too many arguments
-                     "(current-clamp (envelope (10 0.5) (110 0.5) (110 0)))",  // too few arguments
+                     "(density (mechanism \"hh\" (gl 3.5)))",// unqouted parameter
+                     "(density (mechanism \"pas\" ((\"g\" 0.5) (\"e\" 0.2))))",   // paranthesis around params
+                     "(density (mechanism \"pas\" (\"g\" 0.5 0.1) (\"e\" 0.2)))", // too many values
+                     "(current-clamp (envelope (10 0.5) (110 0.5) (110 0)))",     // too few arguments
                      "(current-clamp (envelope (10 0.5) (110 0.5) (110 0)) 10)",  // too few arguments
                      "(paint (region) (mechanism \"hh\"))",  // invalid region
                      "(paint (tag 1) (mechanims hh))",       // invalid painting
                      "(paint (terminal) (membrance-capacitance 0.2))", // can't paint a locset
                      "(paint (tag 3))",                      // too few arguments
-                     "(place (locset) (gap-junction-site) \"gj\")",        // invalid locset
-                     "(place (gap-junction-site) (location 0 1), \"gj\")", // swapped argument order
-                     "(place (location 0 1) (mechanism \"expsyn\"))",      // missing label
+                     "(place (locset) (junction (mechanism \"gj\")) \"gj\")",        // invalid locset
+                     "(place (junction (mechanism \"gj\")) (location 0 1), \"gj\")", // swapped argument order
+                     "(place (location 0 1) (mechanism \"expsyn\"))",                // missing synapse
+                     "(place (location 0 1) (synapse (mechanism \"expsyn\")))",      // missing label
                      "(region-def my_region (tag 3))",       // unquoted region name
                      "(locset-def \"my_ls\" (tag 3))",       // invalid locset
                      "(locset-def \"my_ls\")",               // too few arguments
@@ -1062,10 +1074,10 @@ TEST(doc_expressions, parse) {
                      "  (default (membrane-potential -55.000000))\n"
                      "  (paint (region \"custom\") (temperature-kelvin 270))\n"
                      "  (paint (region \"soma\") (membrane-potential -50.000000))\n"
-                     "  (paint (all) (mechanism \"pas\"))\n"
-                     "  (paint (tag 4) (mechanism \"Ih\" (\"gbar\" 0.001)))\n"
-                     "  (place (locset \"root\") (mechanism \"expsyn\") \"root_synapse\")\n"
-                     "  (place (terminal) (gap-junction-site) \"terminal_gj\"))",
+                     "  (paint (all) (density (mechanism \"pas\")))\n"
+                     "  (paint (tag 4) (density (mechanism \"Ih\" (\"gbar\" 0.001))))\n"
+                     "  (place (locset \"root\") (synapse (mechanism \"expsyn\")) \"root_synapse\")\n"
+                     "  (place (terminal) (junction (mechanism \"gj\")) \"terminal_gj\"))",
                      "(morphology\n"
                      "  (branch 0 -1\n"
                      "    (segment 0 (point 0 0 0 2) (point 4 0 0 2) 1)\n"
@@ -1095,9 +1107,9 @@ TEST(doc_expressions, parse) {
                      "    (default (membrane-potential -55.000000))\n"
                      "    (paint (region \"my_soma\") (temperature-kelvin 270))\n"
                      "    (paint (region \"my_region\") (membrane-potential -50.000000))\n"
-                     "    (paint (tag 4) (mechanism \"Ih\" (\"gbar\" 0.001)))\n"
-                     "    (place (locset \"root\") (mechanism \"expsyn\") \"root_synapse\")\n"
-                     "    (place (location 1 0.2) (gap-junction-site) \"terminal_gj\"))\n"
+                     "    (paint (tag 4) (density (mechanism \"Ih\" (\"gbar\" 0.001))))\n"
+                     "    (place (locset \"root\") (synapse (mechanism \"expsyn\")) \"root_synapse\")\n"
+                     "    (place (location 1 0.2) (junction (mechanism \"gj\")) \"terminal_gj\"))\n"
                      "  (morphology\n"
                      "    (branch 0 -1\n"
                      "      (segment 0 (point 0 0 0 2) (point 4 0 0 2) 1)\n"
@@ -1134,7 +1146,7 @@ TEST(doc_expressions, parse) {
                             "  (meta-data (version \"" + arborio::acc_version() +"\"))\n"
                             "  (decor\n"
                             "    (default (membrane-potential -55.000000))\n"
-                            "    (place (locset \"root\") (mechanism \"expsyn\") \"root_synapse\")\n"
+                            "    (place (locset \"root\") (synapse (mechanism \"expsyn\")) \"root_synapse\")\n"
                             "    (paint (region \"my_soma\") (temperature-kelvin 270))))",
                             "(arbor-component\n"
                             "  (meta-data (version \"" + arborio::acc_version() +"\"))\n"
@@ -1151,7 +1163,7 @@ TEST(doc_expressions, parse) {
                             "      (locset-def \"root\" (root)))\n"
                             "    (decor\n"
                             "      (default (membrane-potential -55.000000))\n"
-                            "      (place (locset \"root\") (mechanism \"expsyn\") \"root_synapse\")\n"
+                            "      (place (locset \"root\") (synapse (mechanism \"expsyn\")) \"root_synapse\")\n"
                             "      (paint (region \"my_soma\") (temperature-kelvin 270)))\n"
                             "    (morphology\n"
                             "       (branch 0 -1\n"
diff --git a/test/unit/test_spikes.cpp b/test/unit/test_spikes.cpp
index 292c1fa3..882de042 100644
--- a/test/unit/test_spikes.cpp
+++ b/test/unit/test_spikes.cpp
@@ -226,7 +226,7 @@ TEST(SPIKES_TEST_CLASS, threshold_watcher_interpolation) {
         decor.set_default(arb::cv_policy_every_segment());
         decor.place("mid"_lab, arb::threshold_detector{10}, "detector");
         decor.place("mid"_lab, arb::i_clamp::box(0.01+i*dt, duration, 0.5), "clamp");
-        decor.place("mid"_lab, arb::mechanism_desc("expsyn"), "synapse");
+        decor.place("mid"_lab, arb::synapse("expsyn"), "synapse");
 
         arb::cable_cell cell(morpho, dict, decor);
         cable1d_recipe rec({cell});
diff --git a/test/unit/test_synapses.cpp b/test/unit/test_synapses.cpp
index 2fc897ae..cd82579f 100644
--- a/test/unit/test_synapses.cpp
+++ b/test/unit/test_synapses.cpp
@@ -31,9 +31,9 @@ TEST(synapses, add_to_cell) {
 
     auto description = make_cell_soma_only(false);
 
-    description.decorations.place(mlocation{0, 0.1}, "expsyn", "synapse0");
-    description.decorations.place(mlocation{0, 0.2}, "exp2syn", "synapse1");
-    description.decorations.place(mlocation{0, 0.3}, "expsyn", "synapse2");
+    description.decorations.place(mlocation{0, 0.1}, synapse("expsyn"), "synapse0");
+    description.decorations.place(mlocation{0, 0.2}, synapse("exp2syn"), "synapse1");
+    description.decorations.place(mlocation{0, 0.3}, synapse("expsyn"), "synapse2");
 
     cable_cell cell(description);
 
@@ -43,16 +43,16 @@ TEST(synapses, add_to_cell) {
     ASSERT_EQ(1u, syns["exp2syn"].size());
 
     EXPECT_EQ((mlocation{0, 0.1}), syns["expsyn"][0].loc);
-    EXPECT_EQ("expsyn", syns["expsyn"][0].item.name());
+    EXPECT_EQ("expsyn", syns["expsyn"][0].item.mech.name());
 
     EXPECT_EQ((mlocation{0, 0.3}), syns["expsyn"][1].loc);
-    EXPECT_EQ("expsyn", syns["expsyn"][1].item.name());
+    EXPECT_EQ("expsyn", syns["expsyn"][1].item.mech.name());
 
     EXPECT_EQ((mlocation{0, 0.2}), syns["exp2syn"][0].loc);
-    EXPECT_EQ("exp2syn", syns["exp2syn"][0].item.name());
+    EXPECT_EQ("exp2syn", syns["exp2syn"][0].item.mech.name());
 
     // adding a synapse to an invalid branch location should throw.
-    description.decorations.place(mlocation{1, 0.3}, "expsyn", "synapse3");
+    description.decorations.place(mlocation{1, 0.3}, synapse("expsyn"), "synapse3");
     EXPECT_THROW((cell=description), std::runtime_error);
 }
 
@@ -85,7 +85,6 @@ TEST(synapses, syn_basic_state) {
     auto exp2syn = unique_cast<mechanism>(global_default_catalogue().instance(backend::kind, "exp2syn").mech);
     ASSERT_TRUE(exp2syn);
 
-    std::vector<fvm_gap_junction> gj = {};
     auto align = std::max(expsyn->data_alignment(), exp2syn->data_alignment());
 
     shared_state state(num_intdom,
@@ -93,7 +92,6 @@ TEST(synapses, syn_basic_state) {
         0,
         std::vector<index_type>(num_comp, 0),
         std::vector<index_type>(num_comp, 0),
-        {},
         std::vector<value_type>(num_comp, -65),
         std::vector<value_type>(num_comp, temp_K),
         std::vector<value_type>(num_comp, 1.),
@@ -109,8 +107,8 @@ TEST(synapses, syn_basic_state) {
     std::vector<index_type> syn_mult(num_syn, 1);
     std::vector<value_type> syn_weight(num_syn, 1.0);
 
-    state.instantiate(*expsyn,  0, {}, {syn_cv, syn_weight, syn_mult});
-    state.instantiate(*exp2syn, 1, {}, {syn_cv, syn_weight, syn_mult});
+    state.instantiate(*expsyn,  0, {}, {syn_cv, {}, syn_weight, syn_mult});
+    state.instantiate(*exp2syn, 1, {}, {syn_cv, {}, syn_weight, syn_mult});
 
     // Parameters initialized to default values?
 
diff --git a/test/unit/unit_test_catalogue.cpp b/test/unit/unit_test_catalogue.cpp
index aef54594..d9e3adb3 100644
--- a/test/unit/unit_test_catalogue.cpp
+++ b/test/unit/unit_test_catalogue.cpp
@@ -10,6 +10,8 @@
 #include "mechanisms/ca_linear.hpp"
 #include "mechanisms/celsius_test.hpp"
 #include "mechanisms/diam_test.hpp"
+#include "mechanisms/gj0.hpp"
+#include "mechanisms/gj1.hpp"
 #include "mechanisms/non_linear.hpp"
 #include "mechanisms/param_as_state.hpp"
 #include "mechanisms/post_events_syn.hpp"
@@ -60,6 +62,8 @@ using namespace arb;
 mechanism_catalogue make_unit_test_catalogue(const mechanism_catalogue& from) {
     mechanism_catalogue cat(from);
 
+    ADD_MECH(cat, gj0)
+    ADD_MECH(cat, gj1)
     ADD_MECH(cat, test_ca)
     ADD_MECH(cat, test_kin1)
     ADD_MECH(cat, test_kinlva)
-- 
GitLab