From 07e0bb624f44076bdcaa563169ce846efa08b601 Mon Sep 17 00:00:00 2001
From: Thorsten Hater <24411438+thorstenhater@users.noreply.github.com>
Date: Mon, 5 Sep 2022 15:02:02 +0200
Subject: [PATCH] Mutable connection table (#1919)

Add functionality, docs, and examples on editing the connection table.  This is a first, small PR
on the topic, further functionality will come as requested.

Teaser example
```cxx
    // This recipe -- omitted -- implements a mutable connection table
    // and returns it via the `connections_on` callback. New entries are
    // added using `add_connection` which makes one link from a central
    // `spike_source_cell` to the given `gid`
    auto rec = recipe();
    rec.add_connection(1);

    auto sim = arb::simulation(rec, ctx);

    // Record spikes
    std::vector<arb::spike> spikes;
    sim.set_local_spike_callback([&spikes](const auto& s) { spikes.insert(spikes.end(), s.begin(), s.end()); });

    // Run and print all spikes per gid
    sim.run(0.25, 0.025);
    show_spikes(spikes);

    // Add a new connection
    rec.add_connection(2);
    sim.update_connections(rec);

    // Run and print all spikes per gid
    sim.run(0.5, 0.025);
    show_spikes(spikes);
```

# Update July '22

- added re-wiring of Event Generators
- implemented a split between a recipe and a recipe-update (cheekily dubbed `topping` until further notice)
- `recipes` are `toppings` via inheritance
---
 arbor/communication/communicator.cpp          | 114 ++++---
 arbor/communication/communicator.hpp          |  13 +-
 arbor/connection.hpp                          |  55 ++--
 arbor/domain_decomposition.cpp                |   2 +-
 arbor/execution_context.cpp                   |  22 +-
 arbor/include/arbor/context.hpp               |  22 +-
 arbor/include/arbor/domain_decomposition.hpp  |   5 +-
 arbor/include/arbor/load_balance.hpp          |   2 +-
 arbor/include/arbor/profile/meter_manager.hpp |   8 +-
 arbor/include/arbor/profile/profiler.hpp      |   2 +-
 arbor/include/arbor/recipe.hpp                |  37 ++-
 arbor/include/arbor/simulation.hpp            |   8 +-
 arbor/include/arbor/spike.hpp                 |   2 +-
 arbor/label_resolution.hpp                    |   2 +-
 arbor/partition_load_balance.cpp              |   2 +-
 arbor/profile/meter_manager.cpp               |   8 +-
 arbor/profile/profiler.cpp                    |   2 +-
 arbor/simulation.cpp                          |  67 ++--
 doc/concepts/domdec.rst                       |  56 +++-
 doc/concepts/interconnectivity.rst            | 108 ++++++-
 doc/dev/communication.rst                     | 300 ++++++++++++++++++
 doc/dev/index.rst                             |   1 +
 doc/python/simulation.rst                     |  14 +-
 example/CMakeLists.txt                        |   1 +
 example/brunel/brunel.cpp                     |   6 +-
 example/plasticity/CMakeLists.txt             |   4 +
 example/plasticity/branch_cell.hpp            | 138 ++++++++
 example/plasticity/plasticity.cpp             | 104 ++++++
 example/plasticity/readme.md                  |  22 ++
 python/context.cpp                            |   2 +-
 python/context.hpp                            |   2 +-
 python/example/plasticity.py                  | 111 +++++++
 python/simulation.cpp                         |  10 +-
 scripts/run_cpp_examples.sh                   |   2 +-
 scripts/run_python_examples.sh                |   1 +
 test/unit-distributed/test_communicator.cpp   |  16 +-
 test/unit/test_probe.cpp                      |  24 +-
 37 files changed, 1058 insertions(+), 237 deletions(-)
 create mode 100644 doc/dev/communication.rst
 create mode 100644 example/plasticity/CMakeLists.txt
 create mode 100644 example/plasticity/branch_cell.hpp
 create mode 100644 example/plasticity/plasticity.cpp
 create mode 100644 example/plasticity/readme.md
 create mode 100644 python/example/plasticity.py

diff --git a/arbor/communication/communicator.cpp b/arbor/communication/communicator.cpp
index 9fa4a0d7..6f816bf7 100644
--- a/arbor/communication/communicator.cpp
+++ b/arbor/communication/communicator.cpp
@@ -26,15 +26,23 @@ communicator::communicator(const recipe& rec,
                            const domain_decomposition& dom_dec,
                            const label_resolution_map& source_resolution_map,
                            const label_resolution_map& target_resolution_map,
-                           execution_context& ctx)
-{
-    distributed_ = ctx.distributed;
-    thread_pool_ = ctx.thread_pool;
+                           execution_context& ctx):  num_total_cells_{rec.num_cells()},
+                                                     num_local_cells_{dom_dec.num_local_cells()},
+                                                     num_local_groups_{dom_dec.num_groups()},
+                                                     num_domains_{(cell_size_type) ctx.distributed->size()},
+                                                     distributed_{ctx.distributed},
+                                                     thread_pool_{ctx.thread_pool} {
+    update_connections(rec, dom_dec, source_resolution_map, target_resolution_map);
+}
 
-    num_domains_ = distributed_->size();
-    num_local_groups_ = dom_dec.num_groups();
-    num_local_cells_ = dom_dec.num_local_cells();
-    auto num_total_cells = rec.num_cells();
+void communicator::update_connections(const connectivity& rec,
+                                      const domain_decomposition& dom_dec,
+                                      const label_resolution_map& source_resolution_map,
+                                      const label_resolution_map& target_resolution_map) {
+    // Forget all lingering information
+    connections_.clear();
+    connection_part_.clear();
+    index_divisions_.clear();
 
     // For caching information about each cell
     struct gid_info {
@@ -81,10 +89,11 @@ communicator::communicator(const recipe& rec,
 
     for (const auto& cell: gid_infos) {
         for (auto c: cell.conns) {
-            if (c.source.gid >= num_total_cells) {
-                throw arb::bad_connection_source_gid(cell.gid, c.source.gid, num_total_cells);
+            auto sgid = c.source.gid;
+            if (sgid >= num_total_cells_) {
+                throw arb::bad_connection_source_gid(cell.gid, sgid, num_total_cells_);
             }
-            const auto src = dom_dec.gid_domain(c.source.gid);
+            const auto src = dom_dec.gid_domain(sgid);
             src_domains.push_back(src);
             src_counts[src]++;
         }
@@ -95,17 +104,18 @@ communicator::communicator(const recipe& rec,
     // the connections as partitioned by the domain of their source gid.
     connections_.resize(n_cons);
     util::make_partition(connection_part_, src_counts);
-    auto offsets = connection_part_;
-    std::size_t pos = 0;
+    auto offsets = connection_part_; // Copy, as we use this as the list of current target indices to write into
+    auto src_domain = src_domains.begin();
     auto target_resolver = resolver(&target_resolution_map);
     for (const auto& cell: gid_infos) {
+        auto index = cell.index_on_domain;
         auto source_resolver = resolver(&source_resolution_map);
         for (const auto& c: cell.conns) {
-            const auto i = offsets[src_domains[pos]]++;
             auto src_lid = source_resolver.resolve(c.source);
             auto tgt_lid = target_resolver.resolve({cell.gid, c.dest});
-            connections_[i] = {{c.source.gid, src_lid}, tgt_lid, c.weight, c.delay, cell.index_on_domain};
-            ++pos;
+            auto offset  = offsets[*src_domain]++;
+            ++src_domain;
+            connections_[offset] = {{c.source.gid, src_lid}, tgt_lid, c.weight, c.delay, index};
         }
     }
 
@@ -119,9 +129,7 @@ communicator::communicator(const recipe& rec,
     // This is num_domains_ independent sorts, so it can be parallelized trivially.
     const auto& cp = connection_part_;
     threading::parallel_for::apply(0, num_domains_, thread_pool_.get(),
-        [&](cell_size_type i) {
-            util::sort(util::subrange_view(connections_, cp[i], cp[i+1]));
-        });
+                                   [&](cell_size_type i) { util::sort(util::subrange_view(connections_, cp[i], cp[i+1])); });
 }
 
 std::pair<cell_size_type, cell_size_type> communicator::group_queue_range(cell_size_type i) {
@@ -130,11 +138,9 @@ std::pair<cell_size_type, cell_size_type> communicator::group_queue_range(cell_s
 }
 
 time_type communicator::min_delay() {
-    auto local_min = std::numeric_limits<time_type>::max();
-    for (auto& con : connections_) {
-        local_min = std::min(local_min, con.delay());
-    }
-
+    auto local_min = std::accumulate(connections_.begin(), connections_.end(),
+                                     std::numeric_limits<time_type>::max(),
+                                     [](auto&& acc, auto&& el) { return std::min(acc, time_type(el.delay)); });
     return distributed_->min(local_min);
 }
 
@@ -153,29 +159,21 @@ gathered_vector<spike> communicator::exchange(std::vector<spike> local_spikes) {
     return global_spikes;
 }
 
-void communicator::make_event_queues(
-        const gathered_vector<spike>& global_spikes,
-        std::vector<pse_vector>& queues)
-{
+void communicator::make_event_queues(const gathered_vector<spike>& global_spikes,
+                                     std::vector<pse_vector>& queues) {
     arb_assert(queues.size()==num_local_cells_);
-
-    using util::subrange_view;
-    using util::make_span;
-    using util::make_range;
-
+    // Predicate for partitioning
+    struct spike_pred {
+        bool operator()(const spike& spk, const cell_member_type& src) { return spk.source < src; }
+        bool operator()(const cell_member_type& src, const spike& spk) { return src < spk.source; }
+    };
     const auto& sp = global_spikes.partition();
     const auto& cp = connection_part_;
-    for (auto dom: make_span(num_domains_)) {
-        auto cons = subrange_view(connections_, cp[dom], cp[dom+1]);
-        auto spks = subrange_view(global_spikes.values(), sp[dom], sp[dom+1]);
-
-        struct spike_pred {
-            bool operator()(const spike& spk, const cell_member_type& src)
-                {return spk.source<src;}
-            bool operator()(const cell_member_type& src, const spike& spk)
-                {return src<spk.source;}
-        };
-
+    for (auto dom: util::make_span(num_domains_)) {
+        auto cons = util::subrange_view(connections_,           cp[dom], cp[dom+1]);
+        auto spks = util::subrange_view(global_spikes.values(), sp[dom], sp[dom+1]);
+        auto sp = spks.begin(), se = spks.end();
+        auto cn = cons.begin(), ce = cons.end();
         // We have a choice of whether to walk spikes or connections:
         // i.e., we can iterate over the spikes, and for each spike search
         // the for connections that have the same source; or alternatively
@@ -186,27 +184,19 @@ void communicator::make_event_queues(
         // complexity of order max(S log(C), C log(S)), where S is the
         // number of spikes, and C is the number of connections.
         if (cons.size()<spks.size()) {
-            auto sp = spks.begin();
-            auto cn = cons.begin();
-            while (cn!=cons.end() && sp!=spks.end()) {
-                auto sources = std::equal_range(sp, spks.end(), cn->source(), spike_pred());
-                for (auto s: make_range(sources)) {
-                    queues[cn->index_on_domain()].push_back(cn->make_event(s));
-                }
-
+            while (cn!=ce && sp!=se) {
+                auto& queue = queues[cn->index_on_domain];
+                auto src = cn->source;
+                auto sources = std::equal_range(sp, se, src, spike_pred());
+                for (auto s: util::make_range(sources)) queue.push_back(cn->make_event(s));
                 sp = sources.first;
                 ++cn;
             }
         }
         else {
-            auto cn = cons.begin();
-            auto sp = spks.begin();
-            while (cn!=cons.end() && sp!=spks.end()) {
-                auto targets = std::equal_range(cn, cons.end(), sp->source);
-                for (auto c: make_range(targets)) {
-                    queues[c.index_on_domain()].push_back(c.make_event(*sp));
-                }
-
+            while (cn!=ce && sp!=se) {
+                auto targets = std::equal_range(cn, ce, sp->source);
+                for (auto c: util::make_range(targets)) queues[c.index_on_domain].push_back(c.make_event(*sp));
                 cn = targets.first;
                 ++sp;
             }
@@ -218,6 +208,10 @@ std::uint64_t communicator::num_spikes() const {
     return num_spikes_;
 }
 
+void communicator::set_num_spikes(std::uint64_t n) {
+    num_spikes_ = n;
+}
+
 cell_size_type communicator::num_local_cells() const {
     return num_local_cells_;
 }
diff --git a/arbor/communication/communicator.hpp b/arbor/communication/communicator.hpp
index 61872c86..7a4becd3 100644
--- a/arbor/communication/communicator.hpp
+++ b/arbor/communication/communicator.hpp
@@ -62,6 +62,7 @@ public:
 
     /// Returns the total number of global spikes over the duration of the simulation
     std::uint64_t num_spikes() const;
+    void set_num_spikes(std::uint64_t n);
 
     cell_size_type num_local_cells() const;
 
@@ -69,10 +70,16 @@ public:
 
     void reset();
 
+    void update_connections(const connectivity& rec,
+                            const domain_decomposition& dom_dec,
+                            const label_resolution_map& source_resolution_map,
+                            const label_resolution_map& target_resolution_map);
+
 private:
-    cell_size_type num_local_cells_;
-    cell_size_type num_local_groups_;
-    cell_size_type num_domains_;
+    cell_size_type num_total_cells_ = 0;
+    cell_size_type num_local_cells_ = 0;
+    cell_size_type num_local_groups_ = 0;
+    cell_size_type num_domains_ = 0;
     std::vector<connection> connections_;
     std::vector<cell_size_type> connection_part_;
     std::vector<cell_size_type> index_divisions_;
diff --git a/arbor/connection.hpp b/arbor/connection.hpp
index 0779fda6..7fab0162 100644
--- a/arbor/connection.hpp
+++ b/arbor/connection.hpp
@@ -15,52 +15,33 @@ public:
                float w,
                float d,
                cell_gid_type didx=cell_gid_type(-1)):
-        source_(src),
-        destination_(dest),
-        weight_(w),
-        delay_(d),
-        index_on_domain_(didx)
+        source(src),
+        destination(dest),
+        weight(w),
+        delay(d),
+        index_on_domain(didx)
     {}
 
-    float weight() const { return weight_; }
-    time_type delay() const { return delay_; }
+    spike_event make_event(const spike& s) { return { destination, s.time + delay, weight}; }
 
-    cell_member_type source() const { return source_; }
-    cell_lid_type destination() const { return destination_; }
-    cell_size_type index_on_domain() const { return index_on_domain_; }
-
-    spike_event make_event(const spike& s) {
-        return {destination_, s.time + delay_, weight_};
-    }
-
-private:
-    cell_member_type source_;
-    cell_lid_type destination_;
-    float weight_;
-    float delay_;
-    cell_size_type index_on_domain_;
+    cell_member_type source;
+    cell_lid_type destination;
+    float weight;
+    float delay;
+    cell_size_type index_on_domain;
 };
 
 // connections are sorted by source id
 // these operators make for easy interopability with STL algorithms
-
-static inline bool operator<(const connection& lhs, const connection& rhs) {
-    return lhs.source() < rhs.source();
-}
-
-static inline bool operator<(const connection& lhs, cell_member_type rhs) {
-    return lhs.source() < rhs;
-}
-
-static inline bool operator<(cell_member_type lhs, const connection& rhs) {
-    return lhs < rhs.source();
-}
+static inline bool operator<(const connection& lhs, const connection& rhs) { return lhs.source < rhs.source; }
+static inline bool operator<(const connection& lhs, cell_member_type rhs)  { return lhs.source < rhs; }
+static inline bool operator<(cell_member_type lhs, const connection& rhs)  { return lhs < rhs.source; }
 
 } // namespace arb
 
 static inline std::ostream& operator<<(std::ostream& o, arb::connection const& con) {
-    return o << "con [" << con.source() << " -> " << con.destination()
-             << " : weight " << con.weight()
-             << ", delay " << con.delay()
-             << ", index " << con.index_on_domain() << "]";
+    return o << "con [" << con.source << " -> " << con.destination
+             << " : weight " << con.weight
+             << ", delay " << con.delay
+             << ", index " << con.index_on_domain << "]";
 }
diff --git a/arbor/domain_decomposition.cpp b/arbor/domain_decomposition.cpp
index de443ab3..aa220821 100644
--- a/arbor/domain_decomposition.cpp
+++ b/arbor/domain_decomposition.cpp
@@ -15,7 +15,7 @@
 namespace arb {
 domain_decomposition::domain_decomposition(
     const recipe& rec,
-    const context& ctx,
+    context ctx,
     const std::vector<group_description>& groups)
 {
     struct partition_gid_domain {
diff --git a/arbor/execution_context.cpp b/arbor/execution_context.cpp
index 0437024b..4b555f1c 100644
--- a/arbor/execution_context.cpp
+++ b/arbor/execution_context.cpp
@@ -14,10 +14,6 @@
 
 namespace arb {
 
-void execution_context_deleter::operator()(execution_context* p) const {
-    delete p;
-}
-
 execution_context::execution_context(const proc_allocation& resources):
     distributed(make_local_context()),
     thread_pool(std::make_shared<threading::task_system>(resources.num_threads)),
@@ -26,7 +22,7 @@ execution_context::execution_context(const proc_allocation& resources):
 {}
 
 ARB_ARBOR_API context make_context(const proc_allocation& p) {
-    return context(new execution_context(p));
+    return std::make_shared<execution_context>(p);
 }
 
 #ifdef ARB_HAVE_MPI
@@ -40,7 +36,7 @@ execution_context::execution_context(const proc_allocation& resources, MPI_Comm
 
 template <>
 ARB_ARBOR_API context make_context<MPI_Comm>(const proc_allocation& p, MPI_Comm comm) {
-    return context(new execution_context(p, comm));
+    return std::make_shared<execution_context>(p, comm);
 }
 #endif
 template <>
@@ -55,30 +51,30 @@ execution_context::execution_context(
 
 template <>
 ARB_ARBOR_API context make_context(const proc_allocation& p, dry_run_info d) {
-    return context(new execution_context(p, d));
+    return std::make_shared<execution_context>(p, d);
 }
 
-ARB_ARBOR_API std::string distribution_type(const context& ctx) {
+ARB_ARBOR_API std::string distribution_type(context ctx) {
     return ctx->distributed->name();
 }
 
-ARB_ARBOR_API bool has_gpu(const context& ctx) {
+ARB_ARBOR_API bool has_gpu(context ctx) {
     return ctx->gpu->has_gpu();
 }
 
-ARB_ARBOR_API unsigned num_threads(const context& ctx) {
+ARB_ARBOR_API unsigned num_threads(context ctx) {
     return ctx->thread_pool->get_num_threads();
 }
 
-ARB_ARBOR_API unsigned num_ranks(const context& ctx) {
+ARB_ARBOR_API unsigned num_ranks(context ctx) {
     return ctx->distributed->size();
 }
 
-ARB_ARBOR_API unsigned rank(const context& ctx) {
+ARB_ARBOR_API unsigned rank(context ctx) {
     return ctx->distributed->id();
 }
 
-ARB_ARBOR_API bool has_mpi(const context& ctx) {
+ARB_ARBOR_API bool has_mpi(context ctx) {
     return ctx->distributed->name() == "MPI";
 }
 
diff --git a/arbor/include/arbor/context.hpp b/arbor/include/arbor/context.hpp
index 75ba004b..f4fadc38 100644
--- a/arbor/include/arbor/context.hpp
+++ b/arbor/include/arbor/context.hpp
@@ -47,14 +47,8 @@ struct proc_allocation {
 struct execution_context;
 
 // arb::context is an opaque handle for the execution context for use
-// in the public API, implemented as a unique pointer.
-//
-// As execution_context is an incomplete type, an explicit deleter must be
-// provided.
-struct ARB_ARBOR_API execution_context_deleter {
-    void operator()(execution_context*) const;
-};
-using context = std::unique_ptr<execution_context, execution_context_deleter>;
+// in the public API, implemented as a shared pointer.
+using context = std::shared_ptr<execution_context>;
 
 // Helpers for creating contexts. These are implemented in the back end.
 
@@ -68,11 +62,11 @@ ARB_ARBOR_API context make_context(const proc_allocation& resources, Comm comm);
 
 // Queries for properties of execution resources in a context.
 
-ARB_ARBOR_API std::string distribution_type(const context&);
-ARB_ARBOR_API bool has_gpu(const context&);
-ARB_ARBOR_API unsigned num_threads(const context&);
-ARB_ARBOR_API bool has_mpi(const context&);
-ARB_ARBOR_API unsigned num_ranks(const context&);
-ARB_ARBOR_API unsigned rank(const context&);
+ARB_ARBOR_API std::string distribution_type(context);
+ARB_ARBOR_API bool has_gpu(context);
+ARB_ARBOR_API unsigned num_threads(context);
+ARB_ARBOR_API bool has_mpi(context);
+ARB_ARBOR_API unsigned num_ranks(context);
+ARB_ARBOR_API unsigned rank(context);
 
 }
diff --git a/arbor/include/arbor/domain_decomposition.hpp b/arbor/include/arbor/domain_decomposition.hpp
index 0de6383a..21accd76 100644
--- a/arbor/include/arbor/domain_decomposition.hpp
+++ b/arbor/include/arbor/domain_decomposition.hpp
@@ -36,7 +36,10 @@ struct group_description {
 class ARB_ARBOR_API domain_decomposition {
 public:
     domain_decomposition() = delete;
-    domain_decomposition(const recipe& rec, const context& ctx, const std::vector<group_description>& groups);
+    domain_decomposition(const recipe& rec, context ctx, const std::vector<group_description>& groups);
+
+    domain_decomposition(const domain_decomposition&) = default;
+    domain_decomposition& operator=(const domain_decomposition&) = default;
 
     int gid_domain(cell_gid_type gid) const;
     int num_domains() const;
diff --git a/arbor/include/arbor/load_balance.hpp b/arbor/include/arbor/load_balance.hpp
index 1c073980..31614403 100644
--- a/arbor/include/arbor/load_balance.hpp
+++ b/arbor/include/arbor/load_balance.hpp
@@ -22,6 +22,6 @@ using partition_hint_map = std::unordered_map<cell_kind, partition_hint>;
 
 ARB_ARBOR_API domain_decomposition partition_load_balance(
     const recipe& rec,
-    const context& ctx,
+    context ctx,
     const partition_hint_map& hint_map = {});
 } // namespace arb
diff --git a/arbor/include/arbor/profile/meter_manager.hpp b/arbor/include/arbor/profile/meter_manager.hpp
index 985e43c6..1b57c44f 100644
--- a/arbor/include/arbor/profile/meter_manager.hpp
+++ b/arbor/include/arbor/profile/meter_manager.hpp
@@ -26,7 +26,7 @@ struct measurement {
     std::string name;
     std::string units;
     std::vector<std::vector<double>> measurements;
-    measurement(std::string, std::string, const std::vector<double>&, const context&);
+    measurement(std::string, std::string, const std::vector<double>&, context);
 };
 
 class ARB_ARBOR_API meter_manager {
@@ -41,8 +41,8 @@ private:
 
 public:
     meter_manager();
-    void start(const context& ctx);
-    void checkpoint(std::string name, const context& ctx);
+    void start(context ctx);
+    void checkpoint(std::string name, context ctx);
 
     const std::vector<std::unique_ptr<meter>>& meters() const;
     const std::vector<std::string>& checkpoint_names() const;
@@ -59,7 +59,7 @@ struct meter_report {
     std::vector<std::string> hosts;
 };
 
-ARB_ARBOR_API meter_report make_meter_report(const meter_manager& manager, const context& ctx);
+ARB_ARBOR_API meter_report make_meter_report(const meter_manager& manager, context ctx);
 ARB_ARBOR_API std::ostream& operator<<(std::ostream& o, const meter_report& report);
 
 } // namespace profile
diff --git a/arbor/include/arbor/profile/profiler.hpp b/arbor/include/arbor/profile/profiler.hpp
index 8d7d02c7..fd8dd980 100644
--- a/arbor/include/arbor/profile/profiler.hpp
+++ b/arbor/include/arbor/profile/profiler.hpp
@@ -36,7 +36,7 @@ struct profile {
 
 // TODO: remove declaration and update the docs
 void profiler_clear();
-ARB_ARBOR_API void profiler_initialize(context& ctx);
+ARB_ARBOR_API void profiler_initialize(context ctx);
 ARB_ARBOR_API void profiler_enter(std::size_t region_id);
 ARB_ARBOR_API void profiler_leave();
 
diff --git a/arbor/include/arbor/recipe.hpp b/arbor/include/arbor/recipe.hpp
index ed5e6218..bf17317e 100644
--- a/arbor/include/arbor/recipe.hpp
+++ b/arbor/include/arbor/recipe.hpp
@@ -61,28 +61,43 @@ struct gap_junction_connection {
         peer(std::move(peer)), local(std::move(local)), weight(g) {}
 };
 
-class ARB_ARBOR_API recipe {
-public:
-    virtual cell_size_type num_cells() const = 0;
-
-    // Cell description type will be specific to cell kind of cell with given gid.
-    virtual util::unique_any get_cell_description(cell_gid_type gid) const = 0;
-    virtual cell_kind get_cell_kind(cell_gid_type) const = 0;
-
-    virtual std::vector<event_generator> event_generators(cell_gid_type) const {
+struct ARB_ARBOR_API has_gap_junctions {
+    virtual std::vector<gap_junction_connection> gap_junctions_on(cell_gid_type) const {
         return {};
     }
+};
+
+struct ARB_ARBOR_API has_synapses {
     virtual std::vector<cell_connection> connections_on(cell_gid_type) const {
         return {};
     }
-    virtual std::vector<gap_junction_connection> gap_junctions_on(cell_gid_type) const {
+};
+
+struct ARB_ARBOR_API has_probes {
+    virtual std::vector<probe_info> get_probes(cell_gid_type gid) const {
         return {};
     }
+};
 
-    virtual std::vector<probe_info> get_probes(cell_gid_type gid) const {
+struct ARB_ARBOR_API has_generators {
+    virtual std::vector<event_generator> event_generators(cell_gid_type) const {
         return {};
     }
+};
 
+// Toppings allow updating a simulation
+struct ARB_ARBOR_API connectivity: public has_synapses, has_generators {
+    virtual ~connectivity() {}
+};
+
+// Recipes allow building a simulation by lazy queries
+struct ARB_ARBOR_API recipe: public has_gap_junctions, has_probes, connectivity {
+    // number of cells to build
+    virtual cell_size_type num_cells() const = 0;
+    // Cell description type will be specific to cell kind of cell with given gid.
+    virtual util::unique_any get_cell_description(cell_gid_type gid) const = 0;
+    // Query cell kind per gid
+    virtual cell_kind get_cell_kind(cell_gid_type) const = 0;
     // Global property type will be specific to given cell kind.
     virtual std::any get_global_properties(cell_kind) const { return std::any{}; };
 
diff --git a/arbor/include/arbor/simulation.hpp b/arbor/include/arbor/simulation.hpp
index 205c8f94..a96e0cd7 100644
--- a/arbor/include/arbor/simulation.hpp
+++ b/arbor/include/arbor/simulation.hpp
@@ -28,11 +28,13 @@ class simulation_state;
 class ARB_ARBOR_API simulation {
 public:
 
-    simulation(const recipe& rec, const context& ctx, const domain_decomposition& decomp);
+    simulation(const recipe& rec, context ctx, const domain_decomposition& decomp);
 
     simulation(const recipe& rec,
-               const context& ctx=make_context(),
-               std::function<domain_decomposition(const recipe&, const context&)> balancer=[](auto& r, auto& c) { return partition_load_balance(r, c); }): simulation(rec, ctx, balancer(rec, ctx)) {}
+               context ctx=make_context(),
+               std::function<domain_decomposition(const recipe&, context)> balancer=[](auto& r, auto c) { return partition_load_balance(r, c); }): simulation(rec, ctx, balancer(rec, ctx)) {}
+
+    void update(const connectivity& rec);
 
     void reset();
 
diff --git a/arbor/include/arbor/spike.hpp b/arbor/include/arbor/spike.hpp
index 69b1558e..f804dfb5 100644
--- a/arbor/include/arbor/spike.hpp
+++ b/arbor/include/arbor/spike.hpp
@@ -17,7 +17,7 @@ struct basic_spike {
     basic_spike() = default;
 
     basic_spike(id_type s, time_type t):
-        source(s), time(t)
+        source(std::move(s)), time(t)
     {}
 
     friend bool operator==(const basic_spike& l, const basic_spike& r) {
diff --git a/arbor/label_resolution.hpp b/arbor/label_resolution.hpp
index c015524e..dbfe2014 100644
--- a/arbor/label_resolution.hpp
+++ b/arbor/label_resolution.hpp
@@ -77,7 +77,7 @@ public:
         lid_hopefully at(unsigned idx) const;
     };
 
-    label_resolution_map() = delete;
+    label_resolution_map() = default;
     explicit label_resolution_map(const cell_labels_and_gids&);
 
     const range_set& at(const cell_gid_type& gid, const cell_tag_type& tag) const;
diff --git a/arbor/partition_load_balance.cpp b/arbor/partition_load_balance.cpp
index 71467aab..8718c743 100644
--- a/arbor/partition_load_balance.cpp
+++ b/arbor/partition_load_balance.cpp
@@ -22,7 +22,7 @@ namespace arb {
 
 ARB_ARBOR_API domain_decomposition partition_load_balance(
     const recipe& rec,
-    const context& ctx,
+    context ctx,
     const partition_hint_map& hint_map)
 {
     using util::make_span;
diff --git a/arbor/profile/meter_manager.cpp b/arbor/profile/meter_manager.cpp
index c2540c84..2886c1c2 100644
--- a/arbor/profile/meter_manager.cpp
+++ b/arbor/profile/meter_manager.cpp
@@ -24,7 +24,7 @@ double mean(const C& c) {
 
 measurement::measurement(std::string n, std::string u,
                          const std::vector<double>& readings,
-                         const context& ctx):
+                         context ctx):
     name(std::move(n)), units(std::move(u))
 {
     auto dist = ctx->distributed;
@@ -54,7 +54,7 @@ meter_manager::meter_manager() {
     }
 };
 
-void meter_manager::start(const context& ctx) {
+void meter_manager::start(context ctx) {
     arb_assert(!started_);
 
     started_ = true;
@@ -70,7 +70,7 @@ void meter_manager::start(const context& ctx) {
     start_time_ = timer_type::tic();
 };
 
-void meter_manager::checkpoint(std::string name, const context& ctx) {
+void meter_manager::checkpoint(std::string name, context ctx) {
     arb_assert(started_);
 
     // Record the time taken on this domain since the last checkpoint
@@ -101,7 +101,7 @@ const std::vector<double>& meter_manager::times() const {
 
 // Build a report of meters, for use at the end of a simulation
 // for output to file or analysis.
-ARB_ARBOR_API meter_report make_meter_report(const meter_manager& manager, const context& ctx) {
+ARB_ARBOR_API meter_report make_meter_report(const meter_manager& manager, context ctx) {
     meter_report report;
 
     // Add the times to the meter outputs
diff --git a/arbor/profile/profiler.cpp b/arbor/profile/profiler.cpp
index c5e0b641..be5b2aff 100644
--- a/arbor/profile/profiler.cpp
+++ b/arbor/profile/profiler.cpp
@@ -357,7 +357,7 @@ ARB_ARBOR_API void profiler_enter(region_id_type region_id) {
     profiler::get_global_profiler().enter(region_id);
 }
 
-ARB_ARBOR_API void profiler_initialize(context& ctx) {
+ARB_ARBOR_API void profiler_initialize(context ctx) {
     profiler::get_global_profiler().initialize(ctx->thread_pool);
 }
 
diff --git a/arbor/simulation.cpp b/arbor/simulation.cpp
index 562db162..a822ae61 100644
--- a/arbor/simulation.cpp
+++ b/arbor/simulation.cpp
@@ -90,7 +90,9 @@ ARB_ARBOR_API void merge_cell_events(
 
 class simulation_state {
 public:
-    simulation_state(const recipe& rec, const domain_decomposition& decomp, execution_context ctx);
+    simulation_state(const recipe& rec, const domain_decomposition& decomp, context ctx);
+
+    void update(const connectivity& rec);
 
     void reset();
 
@@ -116,6 +118,8 @@ public:
     spike_export_function global_export_callback_;
     spike_export_function local_export_callback_;
     epoch_function epoch_callback_;
+    label_resolution_map source_resolution_map_;
+    label_resolution_map target_resolution_map_;
 
 private:
     // Record last computed epoch (integration interval).
@@ -137,6 +141,8 @@ private:
     std::unordered_map<cell_gid_type, gid_local_info> gid_to_local_;
 
     communicator communicator_;
+    context ctx_;
+    domain_decomposition ddc_;
 
     task_system_handle task_system_;
 
@@ -183,11 +189,13 @@ private:
 simulation_state::simulation_state(
         const recipe& rec,
         const domain_decomposition& decomp,
-        execution_context ctx
+        context ctx
     ):
-    task_system_(ctx.thread_pool),
-    local_spikes_({thread_private_spike_store(ctx.thread_pool), thread_private_spike_store(ctx.thread_pool)})
-{
+    ctx_{ctx},
+    ddc_{decomp},
+    task_system_(ctx->thread_pool),
+    local_spikes_({thread_private_spike_store(ctx->thread_pool),
+                  thread_private_spike_store(ctx->thread_pool)}) {
     // Generate the cell groups in parallel, with one task per cell group.
     auto num_groups = decomp.num_groups();
     cell_groups_.resize(num_groups);
@@ -197,7 +205,7 @@ simulation_state::simulation_state(
         [&](cell_group_ptr& group, int i) {
           const auto& group_info = decomp.group(i);
           cell_label_range sources, targets;
-          auto factory = cell_kind_implementation(group_info.kind, group_info.backend, ctx);
+          auto factory = cell_kind_implementation(group_info.kind, group_info.backend, *ctx_);
           group = factory(group_info.gids, rec, sources, targets);
 
           cg_sources[i] = cell_labels_and_gids(std::move(sources), group_info.gids);
@@ -209,44 +217,44 @@ simulation_state::simulation_state(
         local_sources.append(cg_sources.at(i));
         local_targets.append(cg_targets.at(i));
     }
-    auto global_sources = ctx.distributed->gather_cell_labels_and_gids(local_sources);
-
-    auto source_resolution_map = label_resolution_map(std::move(global_sources));
-    auto target_resolution_map = label_resolution_map(std::move(local_targets));
+    auto global_sources = ctx->distributed->gather_cell_labels_and_gids(local_sources);
 
-    communicator_ = arb::communicator(rec, decomp, source_resolution_map, target_resolution_map, ctx);
-
-    const auto num_local_cells = communicator_.num_local_cells();
+    source_resolution_map_ = label_resolution_map(std::move(global_sources));
+    target_resolution_map_ = label_resolution_map(std::move(local_targets));
+    communicator_ = communicator(rec, ddc_, source_resolution_map_, target_resolution_map_, *ctx_);
+    update(rec);
+    epoch_.reset();
+}
 
+void simulation_state::update(const connectivity& rec) {
+    communicator_.update_connections(rec, ddc_, source_resolution_map_, target_resolution_map_);
     // Use half minimum delay of the network for max integration interval.
     t_interval_ = communicator_.min_delay()/2;
 
+    const auto num_local_cells = communicator_.num_local_cells();
     // Initialize empty buffers for pending events for each local cell
     pending_events_.resize(num_local_cells);
-
+    // Forget old generators, if present
+    event_generators_.clear();
     event_generators_.resize(num_local_cells);
     cell_size_type lidx = 0;
     cell_size_type grpidx = 0;
-
-    auto target_resolution_map_ptr = std::make_shared<label_resolution_map>(std::move(target_resolution_map));
-    for (const auto& group_info: decomp.groups()) {
+    auto target_resolution_map_ptr = std::make_shared<label_resolution_map>(target_resolution_map_);
+    for (const auto& group_info: ddc_.groups()) {
         for (auto gid: group_info.gids) {
             // Store mapping of gid to local cell index.
-            gid_to_local_[gid] = gid_local_info{lidx, grpidx};
-
-            // Resolve event_generator targets.
-            // Each event generator gets their own resolver state.
+            gid_to_local_[gid] = {lidx, grpidx};
+            // Resolve event_generator targets; each event generator gets their own resolver state.
             auto event_gens = rec.event_generators(gid);
             for (auto& g: event_gens) {
-                g.resolve_label([target_resolution_map_ptr, event_resolver=resolver(target_resolution_map_ptr.get()), gid]
-                    (const cell_local_label_type& label) mutable {
+                g.resolve_label([target_resolution_map_ptr,
+                                 event_resolver=resolver(target_resolution_map_ptr.get()),
+                                 gid] (const cell_local_label_type& label) mutable {
                         return event_resolver.resolve({gid, label});
                     });
             }
-
             // Set up the event generators for cell gid.
             event_generators_[lidx] = event_gens;
-
             ++lidx;
         }
         ++grpidx;
@@ -257,10 +265,9 @@ simulation_state::simulation_state(
     // the following epoch. In each buffer there is one lane for each local cell.
     event_lanes_[0].resize(num_local_cells);
     event_lanes_[1].resize(num_local_cells);
-
-    epoch_.reset();
 }
 
+
 void simulation_state::reset() {
     epoch_ = epoch();
 
@@ -510,16 +517,18 @@ void simulation_state::inject_events(const cse_vector& events) {
 
 simulation::simulation(
     const recipe& rec,
-    const context& ctx,
+    context ctx,
     const domain_decomposition& decomp)
 {
-    impl_.reset(new simulation_state(rec, decomp, *ctx));
+    impl_.reset(new simulation_state(rec, decomp, ctx));
 }
 
 void simulation::reset() {
     impl_->reset();
 }
 
+void simulation::update(const connectivity& rec) { impl_->update(rec); }
+
 time_type simulation::run(time_type tfinal, time_type dt) {
     if (dt <= 0.0) {
         throw domain_error("Finite time-step must be supplied.");
diff --git a/doc/concepts/domdec.rst b/doc/concepts/domdec.rst
index 835152f2..f4dac781 100644
--- a/doc/concepts/domdec.rst
+++ b/doc/concepts/domdec.rst
@@ -3,30 +3,64 @@
 Domain decomposition
 ====================
 
-An Arbor simulation requires a :ref:`modelrecipe`, a :ref:`(hardware) context <modelhardware>`, and a domain decomposition. The Recipe contains the neuroscientific model, the hardware context describes the computational resources you are going to execute the simulation on, and the domain decomposition describes how Arbor will use the hardware. Since the context and domain decomposition may seem closely related at first, it might be instructive to see how recipes are used by Arbor: 
+An Arbor simulation requires a :ref:`modelrecipe`, a :ref:`(hardware) context
+<modelhardware>`, and a domain decomposition. The Recipe contains the
+neuroscientific model, the hardware context describes the computational
+resources you are going to execute the simulation on, and the domain
+decomposition describes how Arbor will use the hardware. Since the context and
+domain decomposition may seem closely related at first, it might be instructive
+to see how recipes are used by Arbor:
 
 .. raw:: html
    :file: domdec-diag-1.html
 
-A *domain decomposition* describes the distribution of the model over the available computational resources.
-The description partitions the cells in the model as follows:
+A *domain decomposition* describes the distribution of the model over the
+available computational resources. The description partitions the cells in the
+model as follows:
 
-* group the cells into cell groups of the same :ref:`kind <modelcellkind>` of cell;
+* group the cells into cell groups of the same :ref:`kind <modelcellkind>` of
+  cell;
 * assign each cell group to either a CPU core or GPU on a specific MPI rank.
 
-The number of cells in each cell group depends on different factors, including the type of the cell, and whether the
-cell group will run on a CPU core or the GPU. The domain decomposition is solely responsible for describing the distribution
-of cells across cell groups and domains.
+The number of cells in each cell group depends on different factors, including
+the type of the cell, and whether the cell group will run on a CPU core or the
+GPU. The domain decomposition is solely responsible for describing the
+distribution of cells across cell groups and domains.
 
-The domain decomposition can be built manually by the modeler, or an automatic load balancer can be used.
+The domain decomposition can be built manually by the modeler, or an automatic
+load balancer can be used.
 
+We define some terms as used in the context of connectivity
+
+.. glossary::
+
+   connection
+      Tuple of ``(source, target, weight, delay)`` describing an
+      axon/synapse connection as travelling time (`delay`) and attenuation
+      (`weight`) between two sites `source = (gid, spike_detector)` and `target
+      = (gid, synapse)` where `spike_detector` and `synapse` are string labels.
+
+   cell_group
+      List of same-kinded cells that share some information. Must not be split
+      across domains.
+
+   domain
+      Produced by a `load_balancer`, a list of all `cell_groups`
+      located on the same hardware. A `communicator` deals with the full set of
+      cells of one `domain`.
+
+   domain_decomposition
+      List of domains; distributed across MPI processes.
 
 Load balancers
 --------------
 
-A *load balancer* generates the domain decomposition using the model :ref:`recipe <modelrecipe>` and a description of
-the available computational resources on which the model will run described by an :ref:`execution context <modelcontext>`.
-Currently Arbor provides one load balancer and more will be added over time.
+A *load balancer* generates the domain decomposition using the model
+:ref:`recipe <modelrecipe>` and a description of the available computational
+resources on which the model will run described by an :ref:`execution context
+<modelcontext>`. Currently Arbor provides one automatic load balancer and a
+method for manually listing out the cells per MPI task. More approaches might be
+added over time.
 
 API
 ---
diff --git a/doc/concepts/interconnectivity.rst b/doc/concepts/interconnectivity.rst
index 4fee5f69..8c560648 100644
--- a/doc/concepts/interconnectivity.rst
+++ b/doc/concepts/interconnectivity.rst
@@ -12,9 +12,67 @@ Connections only capture the propagation delay and attenuation associated with s
 connectivity: the biophysical modelling of the chemical synapses themselves is the
 responsibility of the target cell model.
 
-Connection sites and gap junction sites are defined on locations on cells as part of the
-:ref:`cell description <modelcelldesc>`.
-A recipe lets you define which sites are connected to which.
+Connection sites and gap junction sites are defined on locations on cells as
+part of the :ref:`cell description <modelcelldesc>`.
+These sites as such are not connected yet, however the :ref:`recipe <modelrecipe>`
+exposes a number of callbacks to form connections and gap junctions between sites.
+The recipe callbacks are interrogated during simulation creation.
+
+In addition, simulations may update their connectivity by building a new
+connection table outside calls to `run`, for example
+
+.. code-block:: python
+
+    rec = recipe()
+    dec = arb.domain_decomposition(rec, ctx)
+    sim = arb.simulation(rec, ctx, dec)
+
+    # run simulation for 0.25ms with the basic connectivity
+    sim.run(0.25, 0.025)
+
+    # extend the recipe to more connections
+    rec.add_connections()
+    #  use `connections_on` to build a new connection table
+    sim.update_connections(rec)
+
+    # run simulation for 0.25ms with the extended connectivity
+    sim.run(0.5, 0.025)
+
+This will completely replace the old table, previous connections to be retained
+must be explicitly included in the updated callback. This can also be used to
+update connection weights and delays. Note, however, that there is currently no
+way to introduce new sites to the simulation, nor any changes to gap junctions.
+
+The ``update_connections`` method accepts either a full ``recipe`` (but will
+**only** use the ``connections_on`` and ``events_generators`` callbacks) or a
+``connectivity``, which is a reduced recipe exposing only the relevant callbacks.
+Currently ``connectivity`` is only available in C++; Python users have to pass a
+full recipe.
+
+.. warn::
+
+   The semantics of connection updates are subtle and might produce surprising
+   results if handled carelessly. In particular, spikes in-flight over a
+   connection will *always* be delivered, even if the connection has been
+   deleted before the time of delivery has passed (`= t_emitted +
+   connection_delay`). As Arbor's connection model joins processes on the axon,
+   the synaptic cleft, and the receiving synapse into a simple pair `(weight,
+   delay)` it is unclear 'where' the action potential is located at the time of
+   deletion relative to the locus of disconnection. Thus, it was decided to
+   deliver spike events regardless. This is will not cause issues when the
+   transition is slow and smooth, ie weights decays over time towards a small
+   value and then the connection is removed. However, drastic and/or frequent
+   changes across busy synapses might cause unexpected behaviour.
+
+.. note::
+
+   Arbor uses a lazily constructed network (from the ``recipe`` callbacks) for
+   good reason; storing the full connectivity (for all ``gids``) in the
+   ``recipe`` can lead to probitivively large memory footprints. Keep this in
+   mind when designing your connectivity and heed the consequences of doing I/O
+   in these callbacks. This is doubly important when using models with dynamic
+   connectivity where the temptation to store all connections is even larger and
+   each call to ``update`` will re-evaluate the corresponding callbacks.
 
 .. _modelconnections:
 
@@ -35,6 +93,14 @@ A recipe lets you define which sites are connected to which.
          a :gen:`global_label`; a target identified using a :gen:`local_label` (:gen:`gid` of target is
          the argument of the recipe method); a connection delay and a connection weight.
 
+         .. code-block:: python
+
+             def connections_on(self, gid):
+                 if gid + 1 < self.num_cells():
+                     return [arbor.connection((gid + 1, "spike-source"), "synapse", weight, delay)]
+                 else:
+                     return []
+
    spike
    action potential
       Spikes travel over :term:`connections <connection>`. In a synapse, they generate an event.
@@ -71,14 +137,34 @@ A recipe lets you define which sites are connected to which.
 
       Similarly to `Connections`, Gap Junctions in Arbor are defined in two steps:
 
-      1. Create labeled **gap junction sites** on two separate cells as part of their
-         :ref:`cell descriptions <modelcelldesc>` in the :ref:`recipe <modelrecipe>`.
-         Each labeled group of gap junctions may contain multiple items on possibly multiple locations
-         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 unit-less connection weight.
-         Two of these connections are needed, on each of the peer and local cells.
+      1. Create labeled **gap junction sites** on two separate cells as part of
+         their :ref:`cell descriptions <modelcelldesc>` in the :ref:`recipe
+         <modelrecipe>`. Each labeled group of gap junctions may contain multiple
+         items on possibly multiple locations 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 unit-less connection weight. Two of these connections are needed, on
+         each of the peer and local cells. The callback `gap_junctions_on`
+         returns a list of these items, eg
+
+         .. code-block:: python
+
+             def gap_junctions_on(self, gid):
+                 n = self.num_cells
+                 if gid + 1 < n and gid > 0:
+                     return [arbor.gap_junction_connection((gid + 1, "gj"), "gj", weight),
+                             arbor.gap_junction_connection((gid - 1, "gj"), "gj", weight),]
+                 elif gid + 1 < n:
+                     return [arbor.gap_junction_connection((gid + 1, "gj"), "gj", weight),]
+                 if gid > 0:
+                     return [arbor.gap_junction_connection((gid - 1, "gj"), "gj", weight),]
+                 else:
+                     return []
+
+         Note that gap junction connections are symmetrical and thus the above
+         example generates two connections, one incoming and one outgoing.
 
    .. Note::
       Only cable cells support gap junctions as of now.
diff --git a/doc/dev/communication.rst b/doc/dev/communication.rst
new file mode 100644
index 00000000..81661345
--- /dev/null
+++ b/doc/dev/communication.rst
@@ -0,0 +1,300 @@
+.. _communication:
+
+Communication in Arbor
+======================
+
+Communication between cells (and thus cell groups) is facilitated using discrete
+events which we call `spikes` in analogy to the underlying neurobiological
+process. Spikes are the only way to communicate between different cell kinds, eg
+cable cells and integrate-and-fire. In accordance to current theory, spikes are
+reduced to a simple time value. The connection is abstracted into a weight and a
+delay; modelling all axonal processes. While this may seem crude, it is a well
+supported model and commonly used in neuronal simulations.
+
+Connections are formed between sources (cable cell: spike detectors) and targets
+(cable cell: synapses). During runtime, events from all sources are concatenated
+on all MPI ranks using ``Allgatherv`` and targets are responsible for selecting
+events they have subscribed to. This is optimised for by sorting events locally
+(by source) and relying on the process layout to convert this into a *globally*
+sorted array.
+
+Exchange of Spikes
+==================
+
+We will start by discussing the exchange of spikes, an extended summary of the
+full communicator class can be found in the next section. Communication of
+spikes is facilitated through the ``communicator::exchange`` method. which is
+
+.. code-block:: c++
+
+   // aka cell_member_type, lexicographic order provided
+   struct id_type {
+       gid_type gid;
+       lid_type index;
+   };
+
+   struct spike {
+       id_type source;
+       double time;
+   };
+
+   using spike_vec = std::vector<spike>;
+
+   // From gathered_vector
+   template <typename T>
+   struct gathered_vector {
+       std::vector<T> values_;
+       std::vector<int> partition_;
+   };
+
+   using g_spike_vec = gathered_vector<spike>;
+
+   g_spike_vec communicator::exchange(spike_vec local_spikes) { // Take by value, since we modify (sort) anyhow.
+       // sort the spikes in ascending order of source gid
+       util::sort_by(local_spikes, [](spike s){ return s.source; });
+       // global all-to-all to gather a local copy of the global spike list on each node.
+       auto global_spikes = distributed_->gather_spikes(local_spikes);
+       num_spikes_ += global_spikes.size();
+       return global_spikes;
+   }
+
+   // From mpi_context
+   g_spike_vec gather_spikes(const spike_ve& local_spikes) const {
+       return mpi::gather_all_with_partition(local_spikes, comm_);
+   }
+
+   // From mpi
+   g_spike_vec gather_all_with_partition(const spike_vec& values, MPI_Comm comm) {
+       using traits = mpi_traits<spike>;
+       // Collect all sender's counts, scale to bytes
+       std::vector<int> counts = gather_all(int(values.size()), comm);
+       for (auto& c: counts) c *= traits::count();
+       // Scan-left to get offsets
+       std::vector<int> displs;
+       util::make_partition(displs, counts);
+       // Allocate recv buffer and exchange
+       spike_vec buffer(displs.back()/traits::count());
+       MPI_Allgatherv(// send buffer  count to send       MPI datatype of spike
+                      values.data(),  counts[rank(comm)], traits::mpi_type(),
+                      // recv buffer  count of each sender offset of senders  MPI datatype of spike
+                      buffer.data(),  counts.data(),       displs.data(),     traits::mpi_type(),
+                      comm);
+       // Scale back to sizeof(spike)
+       for (auto& d: displs) d /= traits::count();
+       return {std::move(buffer), std::move(displs)};
+   }
+
+Note that these snippets have been simplified and shortened from the
+actual code, they are intended for education only.
+
+After ``exchange`` is done, each process has received an object like
+this
+
+.. code-block:: c++
+
+   {
+     spikes:  { {s0, t0}, ...},
+     offsets: { 0, ... }
+   }
+
+Now ``spikes`` is the array of all spikes during the last *epoch* where
+each sub-array of spikes is sorted, ie between ``offsets[ix]`` and
+``offest[ix+1]``. The ``offsets`` array has a length of MPI task count
+and its ``i``'th element gives the position of the first spike sent by
+task ``i``.
+
+Distribution of Events to Targets
+=================================
+
+Having received the generated spikes, the concatenated data is converted
+into events on each local cell group. This is done asynchronously with
+computation of the next cell state. In ``simulation.cpp`` we find
+
+.. code-block:: c++
+
+   auto exchange = [this](epoch prev) {
+       // Collate locally generated spikes.
+       auto all_local_spikes = local_spikes(prev.id).gather();
+       // Gather generated spikes across all ranks.
+       auto global_spikes = communicator_.exchange(all_local_spikes);
+       // Append events formed from global spikes to per-cell pending event queues.
+       communicator_.make_event_queues(global_spikes, pending_events_);
+   };
+
+which uses this
+
+.. code-block:: c++
+
+   // Check each global spike in turn to see it generates local events.
+   // If so, make the events and insert them into the appropriate event list.
+   //
+   // Takes reference to a vector of event lists as an argument, with one list
+   // for each local cell group. On completion, the events in each list are
+   // all events that must be delivered to targets in that cell group as a
+   // result of the global spike exchange, plus any events that were already
+   // in the list.
+   void make_event_queues(const gathered_vector<spike>& global_spikes,
+                          std::vector<pse_vector>& queues) {
+       // Predicate for partitioning
+       struct spike_pred {
+           bool operator()(const spike& spk, const cell_member_type& src) { return spk.source < src; }
+           bool operator()(const cell_member_type& src, const spike& spk) { return src < spk.source; }
+       };
+
+       const auto& sp = global_spikes.partition();
+       for (auto dom: util::make_span(num_domains_)) {
+           // Fetch connections and spikes per integration domain
+           auto cons = util::subrange_view(connections_,           connection_part_[dom], connection_part_[dom + 1]);
+           auto spks = util::subrange_view(global_spikes.values(), sp[dom],               sp[dom + 1]);
+           auto sp = spks.begin(), se = spks.end();
+           auto cn = cons.begin(), ce = cons.end();
+           // We have a choice of whether to walk spikes or connections:
+           // i.e., we can iterate over the spikes, and for each spike search
+           // the for connections that have the same source; or alternatively
+           // for each connection, we can search the list of spikes for spikes
+           // with the same source.
+           //
+           // We iterate over whichever set is the smallest, which has
+           // complexity of order max(S log(C), C log(S)), where S is the
+           // number of spikes, and C is the number of connections.
+           if (cons.size() < spks.size()) {
+               while (cn != ce && sp != ce) {
+                   auto src = cn->source();           // Source for connection
+                   auto cix = cn->index_on_domain();  // Queue index for connection
+                   // Given a source src split the range [sp, spks.end) into a pair sources=[l, h]
+                   // st  *l is the last element not smaller than src
+                   // and *h is the first element greater than src.
+                   // 'Greater' and 'smaller' are defined via the predicate above.
+                   // The range [sp, spks.end) must be (partially) ordered wrt the predicate.
+                   auto sources = std::equal_range(sp, se, src, spike_pred());
+                   // Consequently, the range returned is the range of equal spike sources,
+                   // we pick out ours and add all of them to the appropriate queue.
+                   for (auto s: util::make_range(sources)) queues[cix].push_back(cn->make_event(s));
+                   // now, move to next
+                   sp = sources.first;
+                   ++cn;
+               }
+           }
+           else {
+               while (cn != ce && sp != se) {
+                   auto targets = std::equal_range(cn, ce, sp->source);
+                   for (auto c: util::make_range(targets)) queues[c.index_on_domain()].push_back(c.make_event(*sp));
+                   cn = targets.first;
+                   ++sp;
+               }
+           }
+       }
+   }
+
+After ``make_event_queues`` there is one queue per cell and each queue
+is filled with a time ordered list of events for that cell. We now need
+to understand the actual connection table stored in
+
+.. code-block:: c++
+
+   struct connection {
+       spike_event make_event(const spike& s) {
+           return { destination_, s.time + delay_, weight_};
+       }
+
+       cell_member_type source;
+       cell_lid_type destination;
+       float weight;
+       float delay;
+       cell_size_type index_on_domain;
+   };
+
+   struct communicator {
+       // [...]
+       cell_size_type num_domains_;
+       std::vector<connection> connections_;
+       std::vector<cell_size_type> connection_part_;
+       // [...]
+   };
+
+The ``connections`` vector is a list of connections partitioned by the
+domain (as in domain decomposition) of their source's ``gid``, while
+``connection_part`` stores the partioning indices.
+
+Building the Connection Table
+=============================
+
+The table of connections on the local rank is built during the construction of
+the ``communicator`` object
+
+.. code-block:: c++
+
+   communicator::communicator(const recipe& rec,
+                              const domain_decomposition& dom_dec,
+                              const label_resolution_map& source_resolution_map,
+                              const label_resolution_map& target_resolution_map,
+                              execution_context& ctx);
+
+After that process,
+
+.. code-block:: c++
+
+   struct communicator {
+       // ...
+       std::vector<connection> connections_;
+       std::vector<cell_size_type> connection_part_;
+   };
+
+will contain all connections in ``connections_`` partitioned by the
+domain of the source's ``gid`` in ``dom_dec``. Beginnings of the
+respective partitions are pointed to by the indices in
+``connection_part_``.
+
+The algorithm for building is slightly obscured by caching and the use
+of labels and resolving them via ``target_/source_resolution_map`` to
+local ids on the respective source and target cells.
+
+.. note::
+
+   The ``label_resolution_map`` class is used to translate from labels at the
+   user facing API layers to Arbor's internal mappings in the vein of
+   ``(cell_gid, item_offset)``, where ``item_offset`` is an automatically
+   assigned integer ID. Textual labels are created by calls to ``place``
+   as in this example
+
+   .. code-block:: c++
+
+      auto d = arb::decor{};
+      d.place("..."_ls, arb::synapse{"..."}, "synapse-label");
+
+The construction is performed in-place
+
+.. code-block:: c++
+
+    // Allocate space for our connections
+    connections_.resize(n_cons);
+    // We have pre-computed `src_counts`, connection_part_ will now hold the starting indices
+    // of each `domain`.
+    util::make_partition(connection_part_, src_counts);
+    // Copy, as we use this as the list of the currently available next free target slots in
+    // `connections_`
+    auto offsets = connection_part_;
+    auto target_resolver = resolver(&target_resolution_map);
+    auto src_domain = src_domains.begin();
+    for (const auto& cell: gid_infos) {
+        auto index = cell.index_on_domain;
+        auto source_resolver = resolver(&source_resolution_map);
+        for (const auto& c: cell.conns) {
+            // Compute index representation of labels
+            auto src_lid = source_resolver.resolve(c.source);
+            auto tgt_lid = target_resolver.resolve({cell.gid, c.dest});
+            // Get offset of current source and bump to next free slot
+            auto offset  = offsets[*src_domain]++;
+            // Write connection info into slot
+            connections_[offset] = {{c.source.gid, src_lid}, tgt_lid, c.weight, c.delay, index};
+            // Next source domain
+            ++src_domain;
+        }
+    }
+    // Now
+    // * all slots in `connections_` are filled.
+    // * `offsets` points at the ends of each partition.
+
+
+Next, each *partition* is sorted independently according to their
+source's ``gid``.
diff --git a/doc/dev/index.rst b/doc/dev/index.rst
index 0174b15f..9de6651c 100644
--- a/doc/dev/index.rst
+++ b/doc/dev/index.rst
@@ -16,6 +16,7 @@ Here we document internal components of Arbor. These pages can be useful if you'
 
    cable_cell
    cell_groups
+   communication
    debug
    matrix_solver
    simd_api
diff --git a/doc/python/simulation.rst b/doc/python/simulation.rst
index 533041d7..818c7363 100644
--- a/doc/python/simulation.rst
+++ b/doc/python/simulation.rst
@@ -63,12 +63,22 @@ over the local and distributed hardware resources (see :ref:`pydomdec`). Then, t
 
     **Constructor:**
 
-    .. function:: simulation(recipe, domain_decomposition, context)
+    .. function:: simulation(recipe, context, domain_decomposition)
 
-        Initialize the model described by an :py:class:`arbor.recipe`, with cells and network distributed according to :py:class:`arbor.domain_decomposition`, and computational resources described by :py:class:`arbor.context`.
+        Initialize the model described by an :py:class:`arbor.recipe`, with
+        cells and network distributed according to
+        :py:class:`arbor.domain_decomposition`, and computational resources
+        described by :py:class:`arbor.context`.
 
     **Updating Model State:**
 
+    .. function:: update_connections(recipe)
+
+        Rebuild the connection table as described by
+        :py:class:`arbor.recipe::connections_on` The recipe must differ **only**
+        in the return value of its :py:func:`connections_on` when compared to
+        the original recipe used to construct the simulation object.
+
     .. function:: reset()
 
         Reset the state of the simulation to its initial state.
diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt
index 752737ec..f1393268 100644
--- a/example/CMakeLists.txt
+++ b/example/CMakeLists.txt
@@ -11,5 +11,6 @@ add_subdirectory(ring)
 add_subdirectory(gap_junctions)
 add_subdirectory(single)
 add_subdirectory(probe-demo)
+add_subdirectory(plasticity)
 add_subdirectory(lfp)
 add_subdirectory(diffusion)
diff --git a/example/brunel/brunel.cpp b/example/brunel/brunel.cpp
index 351e656f..e1e19b3b 100644
--- a/example/brunel/brunel.cpp
+++ b/example/brunel/brunel.cpp
@@ -70,7 +70,7 @@ std::ostream& operator<<(std::ostream& o, const cl_options& opt);
 
 std::optional<cl_options> read_options(int argc, char** argv);
 
-void banner(const context& ctx);
+void banner(context ctx);
 
 // Samples m unique values in interval [start, end) - gid.
 // We exclude gid because we don't want self-loops.
@@ -249,7 +249,7 @@ int main(int argc, char** argv) {
 
         simulation sim(recipe,
                        context,
-                       [&hints](auto& r, auto& c) { return partition_load_balance(r, c, hints); });
+                       [&hints](auto& r, auto c) { return partition_load_balance(r, c, hints); });
 
         // Set up spike recording.
         std::vector<arb::spike> recorded_spikes;
@@ -296,7 +296,7 @@ int main(int argc, char** argv) {
     return 0;
 }
 
-void banner(const context& ctx) {
+void banner(context ctx) {
     std::cout << "==========================================\n";
     std::cout << "  Brunel model miniapp\n";
     std::cout << "  - distributed : " << arb::num_ranks(ctx)
diff --git a/example/plasticity/CMakeLists.txt b/example/plasticity/CMakeLists.txt
new file mode 100644
index 00000000..904ea094
--- /dev/null
+++ b/example/plasticity/CMakeLists.txt
@@ -0,0 +1,4 @@
+add_executable(plasticity EXCLUDE_FROM_ALL plasticity.cpp)
+add_dependencies(examples plasticity)
+
+target_link_libraries(plasticity PRIVATE arbor arborio arborenv)
diff --git a/example/plasticity/branch_cell.hpp b/example/plasticity/branch_cell.hpp
new file mode 100644
index 00000000..527f2601
--- /dev/null
+++ b/example/plasticity/branch_cell.hpp
@@ -0,0 +1,138 @@
+#pragma once
+
+#include <array>
+#include <random>
+
+#include <nlohmann/json.hpp>
+
+#include <arborio/label_parse.hpp>
+
+#include <arbor/cable_cell.hpp>
+#include <arbor/cable_cell_param.hpp>
+#include <arbor/common_types.hpp>
+#include <arbor/morph/segment_tree.hpp>
+
+#include <string>
+#include <sup/json_params.hpp>
+
+using namespace arborio::literals;
+
+// Parameters used to generate the random cell morphologies.
+struct cell_parameters {
+    cell_parameters() = default;
+
+    //  Maximum number of levels in the cell (not including the soma)
+    unsigned max_depth = 5;
+
+    // The following parameters are described as ranges.
+    // The first value is at the soma, and the last value is used on the last level.
+    // Values at levels in between are found by linear interpolation.
+    std::array<double,2> branch_probs = {1.0, 0.5}; //  Probability of a branch occuring.
+    std::array<unsigned,2> compartments = {20, 2};  //  Compartment count on a branch.
+    std::array<double,2> lengths = {200, 20};       //  Length of branch in μm.
+
+    // The number of synapses per cell.
+    unsigned synapses = 1;
+};
+
+cell_parameters parse_cell_parameters(nlohmann::json& json) {
+    cell_parameters params;
+    sup::param_from_json(params.max_depth, "depth", json);
+    sup::param_from_json(params.branch_probs, "branch-probs", json);
+    sup::param_from_json(params.compartments, "compartments", json);
+    sup::param_from_json(params.lengths, "lengths", json);
+    sup::param_from_json(params.synapses, "synapses", json);
+
+    return params;
+}
+
+// Helper used to interpolate in branch_cell.
+template <typename T>
+double interp(const std::array<T,2>& r, unsigned i, unsigned n) {
+    double p = i * 1./(n-1);
+    double r0 = r[0];
+    double r1 = r[1];
+    return r[0] + p*(r1-r0);
+}
+
+arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params) {
+    arb::segment_tree tree;
+
+    // Add soma.
+    double srad = 12.6157/2.0; // soma radius
+    int stag = 1; // soma tag
+    tree.append(arb::mnpos, {0, 0,-srad, srad}, {0, 0, srad, srad}, stag); // For area of 500 μm².
+
+    std::vector<std::vector<unsigned>> levels;
+    levels.push_back({0});
+
+    // Standard mersenne_twister_engine seeded with gid.
+    std::mt19937 gen(gid);
+    std::uniform_real_distribution<double> dis(0, 1);
+
+    double drad = 0.5; // Diameter of 1 μm for each dendrite cable.
+    int dtag = 3;      // Dendrite tag.
+
+    double dist_from_soma = srad; // Start dendrite at the edge of the soma.
+    for (unsigned i=0; i<params.max_depth; ++i) {
+        // Branch prob at this level.
+        double bp = interp(params.branch_probs, i, params.max_depth);
+        // Length at this level.
+        double l = interp(params.lengths, i, params.max_depth);
+        // Number of compartments at this level.
+        unsigned nc = std::round(interp(params.compartments, i, params.max_depth));
+
+        std::vector<unsigned> sec_ids;
+        for (unsigned sec: levels[i]) {
+            for (unsigned j=0; j<2; ++j) {
+                if (dis(gen)<bp) {
+                    auto z = dist_from_soma;
+                    auto dz = l/nc;
+                    auto p = sec;
+                    for (unsigned k=0; k<nc; ++k) {
+                        p = tree.append(p, {0,0,z+(k+1)*dz, drad}, dtag);
+                    }
+                    sec_ids.push_back(p);
+                }
+            }
+        }
+        if (sec_ids.empty()) {
+            break;
+        }
+        levels.push_back(sec_ids);
+
+        dist_from_soma += l;
+    }
+
+    arb::label_dict labels;
+
+    using arb::reg::tagged;
+    labels.set("soma", tagged(stag));
+    labels.set("dend", tagged(dtag));
+
+    arb::decor decor;
+
+    decor.paint("soma"_lab, arb::density("hh"));
+    decor.paint("dend"_lab, arb::density("pas"));
+
+    decor.set_default(arb::axial_resistivity{100}); // [Ω·cm]
+
+    // Add spike threshold detector at the soma.
+    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}, 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), arb::synapse("expsyn"), "extra_syns");
+    }
+
+    // Make a CV between every sample in the sample tree.
+    decor.set_default(arb::cv_policy_every_segment());
+
+    arb::cable_cell cell(arb::morphology(tree), labels, decor);
+
+    return cell;
+}
diff --git a/example/plasticity/plasticity.cpp b/example/plasticity/plasticity.cpp
new file mode 100644
index 00000000..4b97d458
--- /dev/null
+++ b/example/plasticity/plasticity.cpp
@@ -0,0 +1,104 @@
+#include <any>
+#include <iostream>
+#include <iomanip>
+#include <unordered_map>
+
+#include <arborio/label_parse.hpp>
+
+#include <arbor/spike_source_cell.hpp>
+#include <arbor/cable_cell.hpp>
+#include <arbor/simulation.hpp>
+#include <arbor/util/any_cast.hpp>
+#include <arbor/util/any_ptr.hpp>
+
+using namespace arborio::literals;
+
+// Fan out network with n members
+// - one spike source at gid 0
+// - n-1 passive, soma-only cable cells
+// Starts out disconnected, but new connections from the source to the cable cells may be added.
+struct recipe: public arb::recipe {
+    const std::string                           // Textuals labels for
+        syn = "synapse",                        //     exponential synapse on cable cells
+        det = "detector",                       //     spike detector on cable cells
+        src = "src";                            //     spike source
+    const double                                // Parameters
+        r_soma  = 12.6157/2.0,                  //     soma radius; A=500 μm²
+        f_spike = 0.25,                         //     spike source interval
+        weight  = 5.5,                          //     connection weight
+        delay   = 0.05;                         //                delay
+    arb::locset center = "(location 0 0.5)"_ls; // Soma center
+    arb::region all    = "(all)"_reg;           // Whole cell
+    arb::cell_size_type n_ = 0;                 // Cell count
+
+    mutable std::unordered_map<arb::cell_gid_type, std::vector<arb::cell_connection>> connected; // lookup table for connections
+    // Required but uninteresting methods
+    recipe(arb::cell_size_type n): n_{n} {}
+    arb::cell_size_type num_cells() const override { return n_; }
+    arb::cell_kind get_cell_kind(arb::cell_gid_type gid) const override { return 0 == gid ? arb::cell_kind::spike_source : arb::cell_kind::cable; }
+    std::any get_global_properties(arb::cell_kind kind) const override {
+        if (kind == arb::cell_kind::spike_source) return {};
+        arb::cable_cell_global_properties ccp;
+        ccp.default_parameters = arb::neuron_parameter_defaults;
+        return ccp;
+    }
+    // For printing Um
+    std::vector<arb::probe_info> get_probes(arb::cell_gid_type gid) const override {
+        if (gid == 0) return {};
+        return {arb::cable_probe_membrane_voltage{center}};
+    }
+    // Look up the (potential) connection to this cell
+    std::vector<arb::cell_connection> connections_on(arb::cell_gid_type gid) const override { return connected[gid]; }
+    // Connect cell `to` to the spike source
+    void add_connection(arb::cell_gid_type to) { assert(to > 0); connected[to] = {arb::cell_connection({0, src}, {syn}, weight, delay)}; }
+    // Return the cell at gid
+    arb::util::unique_any get_cell_description(arb::cell_gid_type gid) const override {
+        // source at gid 0
+        if (gid == 0) return arb::spike_source_cell{src, arb::regular_schedule(f_spike)};
+        // all others are receiving cable cells; single CV w/ HH
+        arb::segment_tree tree; tree.append(arb::mnpos, {-r_soma, 0, 0, r_soma}, {r_soma, 0, 0, r_soma}, 1);
+        auto decor = arb::decor{};
+        decor.paint(all, arb::density("hh", {{"gl", 5}}));
+        decor.place(center, arb::synapse("expsyn"), syn);
+        decor.place(center, arb::threshold_detector{-10.0}, det);
+        decor.set_default(arb::cv_policy_every_segment());
+        return arb::cable_cell({tree}, {}, decor);
+    }
+};
+
+void sampler(arb::probe_metadata pm, std::size_t n, const arb::sample_record* samples) {
+    auto* loc = arb::util::any_cast<const arb::mlocation*>(pm.meta);
+    std::cout << std::fixed << std::setprecision(4);
+    for (std::size_t i = 0; i<n; ++i) {
+        auto* value = arb::util::any_cast<const double*>(samples[i].data);
+        std::cout << "|  " << samples[i].time << " |      " << loc->pos << " | " << *value << " |\n";
+    }
+}
+
+void spike_cb(const std::vector<arb::spike>& spikes) {
+    for(const auto& spike: spikes) std::cout << " * " << spike.source << "@" << spike.time << '\n';
+}
+
+void print_header(double from, double to) {
+    std::cout << "\n"
+              << "Running simulation from " << from << "ms to " << to << "ms\n"
+              << "Spikes are marked: *\n"
+              << "\n"
+              << "| Time/ms | Position/um | Um/mV    |\n"
+              << "|---------+-------------+----------|\n";
+}
+
+const double dt = 0.05;
+
+int main(int argc, char** argv) {
+    auto rec = recipe(3);
+    rec.add_connection(1);
+    auto sim = arb::simulation(rec);
+    sim.add_sampler(arb::all_probes, arb::regular_schedule(dt), sampler, arb::sampling_policy::exact);
+    sim.set_local_spike_callback(spike_cb);
+    print_header(0, 1);
+    sim.run(1.0, dt);
+    rec.add_connection(2);
+    print_header(1, 2);
+    sim.run(2.0, dt);
+}
diff --git a/example/plasticity/readme.md b/example/plasticity/readme.md
new file mode 100644
index 00000000..e2cac0b5
--- /dev/null
+++ b/example/plasticity/readme.md
@@ -0,0 +1,22 @@
+# Plasticity Example
+
+A miniapp that demonstrates how to use mutable connectivity. The simulation is
+run in two parts with different connectivity.
+
+All cells will print out their _generated_ spikes and (cable cells only!) their
+membrane potential at the soma. The topology is this
+```
+         +--------------+
+         | 0 Spike      |
+         |   Source     |
+         | f=1/0.0125ms |
+         +--------------+
+        /                \ 
+       /                  \    <--- This connection will not be present for the first part
+      /                    \        and enabled for the second part of the simulation.
+     /                      \
+     +---------+             +---------+
+     | 1 Cable |             | 2 Cable |
+     |   Cell  |             |   Cell  |
+     +---------+             +---------+
+```
diff --git a/python/context.cpp b/python/context.cpp
index 1de01e62..2f0e6ffd 100644
--- a/python/context.cpp
+++ b/python/context.cpp
@@ -92,7 +92,7 @@ context_shim create_context(unsigned threads, pybind11::object gpu, pybind11::ob
         return context_shim(arb::make_context(alloc, c->comm));
     }
 #endif
-    return context_shim(arb::make_context(alloc));
+    return context_shim{arb::make_context(alloc)};
 }
 
 std::ostream& operator<<(std::ostream& o, const proc_allocation_shim& alloc) {
diff --git a/python/context.hpp b/python/context.hpp
index 8fff4bc3..3ba8bfe4 100644
--- a/python/context.hpp
+++ b/python/context.hpp
@@ -6,7 +6,7 @@ namespace pyarb {
 
 struct context_shim {
     arb::context context;
-    context_shim(arb::context&& c): context(std::move(c)) {}
+    context_shim(arb::context c): context{c} {}
 };
 
 } // namespace pyarb
diff --git a/python/example/plasticity.py b/python/example/plasticity.py
new file mode 100644
index 00000000..ee20edcc
--- /dev/null
+++ b/python/example/plasticity.py
@@ -0,0 +1,111 @@
+import arbor as A
+
+
+class recipe(A.recipe):
+    def __init__(self, n):
+        # Call our parent; needed for proper initialization.
+        A.recipe.__init__(self)
+        # Cell count; should be at least 3 for the example to work.
+        assert n >= 3
+        self.cells = n
+        # Uniform weights and delays for the connectivity.
+        self.weight = 0.75
+        self.delay = 0.1
+        # Track the connections from the source cell to a given gid.
+        self.connected = set()
+
+    def num_cells(self):
+        return self.cells
+
+    def cell_kind(self, gid):
+        # Cell 0 is the spike source, all others cable cells.
+        if gid == 0:
+            return A.cell_kind.spike_source
+        else:
+            return A.cell_kind.cable
+
+    def global_properties(self, kind):
+        # Cable cells return the NRN defaults.
+        if kind == A.cell_kind.cable:
+            return A.neuron_cable_properties()
+        # Spike source cells have nothing to report.
+        return None
+
+    def cell_description(self, gid):
+        # Cell 0 is reserved for the spike source, spiking every 0.0125ms.
+        if gid == 0:
+            return A.spike_source_cell("source", A.regular_schedule(0.0125))
+        # All cells >= 1 are cable cells w/ a simple, soma-only morphology
+        # comprising two segments of radius r=3um and length l=3um *each*.
+        #
+        #   +-- --+  ^
+        #   |  *  |  3um
+        #   +-- --+  v
+        #   < 6um >
+        #
+        #  * Detectors and Synapses are here, at the midpoint.
+        r = 3
+        tree = A.segment_tree()
+        tree.append(A.mnpos, A.mpoint(-r, 0, 0, r), A.mpoint(r, 0, 0, r), tag=1)
+        # ... and a synapse/detector pair
+        decor = A.decor()
+        #   - just have a leaky membrane here.
+        decor.paint("(all)", A.density("pas"))
+        #   - synapse to receive incoming spikes from the source cell.
+        decor.place("(location 0 0.5)", A.synapse("expsyn"), "synapse")
+        #   - detector for reporting spikes on the cable cells.
+        decor.place("(location 0 0.5)", A.spike_detector(-10.0), "detector")
+        # return the cable cell description
+        return A.cable_cell(tree, A.label_dict(), decor)
+
+    def connections_on(self, gid):
+        # If we have added a connection to this cell, return it, else nothing.
+        if gid in self.connected:
+            # Connect spike source to synapse(s)
+            return [A.connection((0, "source"), "synapse", self.weight, self.delay)]
+        return []
+
+    def add_connection_to_spike_source(self, to):
+        """Add a connection from the spike source at gid=0 to the cable cell
+        gid=<to>. Note that we try to minimize the information stored here,
+        which is important with large networks. Also we cannot connect the
+        source back to itself.
+        """
+        assert to != 0
+        self.connected.add(to)
+
+
+# Make an unconnected network with 2 cable cells and one spike source,
+rec = recipe(3)
+
+# but before setting up anything, connect cable cell gid=1 to spike source gid=0
+# and make the simulation of the simple network
+#
+#    spike_source <gid=0> ----> cable_cell <gid=1>
+#
+#                               cable_cell <gid=2>
+#
+# Note that the connection is just _recorded_ in the recipe, the actual connectivity
+# is set up in the simulation construction.
+rec.add_connection_to_spike_source(1)
+sim = A.simulation(rec)
+sim.record(A.spike_recording.all)
+
+# then run the simulation for a bit
+sim.run(0.25, 0.025)
+
+# update the simulation to
+#
+#    spike_source <gid=0> ----> cable_cell <gid=1>
+#                        \
+#                         ----> cable_cell <gid=2>
+rec.add_connection_to_spike_source(2)
+sim.update(rec)
+
+# and run the simulation for another bit.
+sim.run(0.5, 0.025)
+
+# When finished, print spike times and locations.
+print("spikes:")
+for sp in sim.spikes():
+    print(" ", sp)
diff --git a/python/simulation.cpp b/python/simulation.cpp
index 647522d0..03b80bd0 100644
--- a/python/simulation.cpp
+++ b/python/simulation.cpp
@@ -94,6 +94,10 @@ public:
         sim_->set_binning_policy(policy, bin_interval);
     }
 
+    void update(std::shared_ptr<py_recipe>& rec) {
+        sim_->update(py_recipe_shim(rec));
+    }
+
     void record(spike_recording policy) {
         auto spike_recorder = [this](const std::vector<arb::spike>& spikes) {
             auto old_size = spike_record_.size();
@@ -215,6 +219,10 @@ void register_simulation(pybind11::module& m, pyarb_global_ptr global_ptr) {
              "recipe"_a,
              pybind11::arg_v("context", pybind11::none(), "Execution context"),
              pybind11::arg_v("domains", pybind11::none(), "Domain decomposition"))
+        .def("update", &simulation_shim::update,
+             "Rebuild the connection table from recipe::connections_on and the event"
+             "generators based on recipe::event_generators.",
+             "recipe"_a)
         .def("reset", &simulation_shim::reset,
             pybind11::call_guard<pybind11::gil_scoped_release>(),
             "Reset the state of the simulation to its initial state.")
@@ -227,7 +235,7 @@ void register_simulation(pybind11::module& m, pyarb_global_ptr global_ptr) {
             "tfinal"_a, "dt"_a=0.025)
         .def("set_binning_policy", &simulation_shim::set_binning_policy,
             "Set the binning policy for event delivery, and the binning time interval if applicable [ms].",
-            "policy"_a, "bin_interval"_a)
+             "policy"_a, "bin_interval"_a)
         .def("record", &simulation_shim::record,
             "Disable or enable local or global spike recording.")
         .def("spikes", &simulation_shim::spikes,
diff --git a/scripts/run_cpp_examples.sh b/scripts/run_cpp_examples.sh
index 36c70669..efbcf815 100755
--- a/scripts/run_cpp_examples.sh
+++ b/scripts/run_cpp_examples.sh
@@ -26,7 +26,7 @@ check () {
     fi
 }
 
-for ex in bench brunel gap_junctions generators lfp ring single-cell "probe-demo v"
+for ex in bench brunel gap_junctions generators lfp ring single-cell "probe-demo v" plasticity
 do
     echo "   - $ex"
     dir=`echo $ex | tr ' ' '_'`
diff --git a/scripts/run_python_examples.sh b/scripts/run_python_examples.sh
index ff269083..4277d7b7 100755
--- a/scripts/run_python_examples.sh
+++ b/scripts/run_python_examples.sh
@@ -30,3 +30,4 @@ $PREFIX python python/example/network_ring.py
 $PREFIX python python/example/network_ring_gpu.py # by default, gpu_id=None
 $PREFIX python python/example/network_two_cells_gap_junctions.py
 $PREFIX python python/example/diffusion.py
+$PREFIX python python/example/plasticity.py
diff --git a/test/unit-distributed/test_communicator.cpp b/test/unit-distributed/test_communicator.cpp
index eeefa667..e3ac5250 100644
--- a/test/unit-distributed/test_communicator.cpp
+++ b/test/unit-distributed/test_communicator.cpp
@@ -644,10 +644,10 @@ TEST(communicator, all2all)
     for (auto i: util::make_span(0, n_global)) {
         for (unsigned j = 0; j < n_local; ++j) {
             auto c = connections[i*n_local+j];
-            EXPECT_EQ(i, c.source().gid);
-            EXPECT_EQ(0u, c.source().index);
-            EXPECT_EQ(i, c.destination());
-            EXPECT_LT(c.index_on_domain(), n_local);
+            EXPECT_EQ(i, c.source.gid);
+            EXPECT_EQ(0u, c.source.index);
+            EXPECT_EQ(i, c.destination);
+            EXPECT_LT(c.index_on_domain, n_local);
         }
     }
 
@@ -689,7 +689,7 @@ TEST(communicator, mini_network)
     // sort connections by source then target
     auto connections = C.connections();
     util::sort(connections, [](const connection& lhs, const connection& rhs) {
-      return std::forward_as_tuple(lhs.source(), lhs.index_on_domain(), lhs.destination()) < std::forward_as_tuple(rhs.source(), rhs.index_on_domain(), rhs.destination());
+      return std::forward_as_tuple(lhs.source, lhs.index_on_domain, lhs.destination) < std::forward_as_tuple(rhs.source, rhs.index_on_domain, rhs.destination);
     });
 
     // Expect one set of 22 connections from every rank: these have been sorted.
@@ -701,9 +701,9 @@ TEST(communicator, mini_network)
         std::vector<cell_gid_type> ex_source_gids(22u, i*3 + 1);
         for (unsigned j = 0; j < 22u; ++j) {
             auto c = connections[i*22 + j];
-            EXPECT_EQ(ex_source_gids[j], c.source().gid);
-            EXPECT_EQ(ex_source_lids[j], c.source().index);
-            EXPECT_EQ(ex_target_lids[i%2][j], c.destination());
+            EXPECT_EQ(ex_source_gids[j], c.source.gid);
+            EXPECT_EQ(ex_source_lids[j], c.source.index);
+            EXPECT_EQ(ex_target_lids[i%2][j], c.destination);
         }
     }
 }
diff --git a/test/unit/test_probe.cpp b/test/unit/test_probe.cpp
index a122552a..61f5fef2 100644
--- a/test/unit/test_probe.cpp
+++ b/test/unit/test_probe.cpp
@@ -99,7 +99,7 @@ static morphology make_stick_morphology() {
 }
 
 template <typename Backend>
-void run_v_i_probe_test(const context& ctx) {
+void run_v_i_probe_test(context ctx) {
     using fvm_cell = typename backend_access<Backend>::fvm_cell;
     auto deref = [](const arb_value_type* p) { return backend_access<Backend>::deref(p); };
 
@@ -198,7 +198,7 @@ void run_v_i_probe_test(const context& ctx) {
 }
 
 template <typename Backend>
-void run_v_cell_probe_test(const context& ctx) {
+void run_v_cell_probe_test(context ctx) {
     using fvm_cell = typename backend_access<Backend>::fvm_cell;
 
     // Take the per-cable voltage over a Y-shaped cell with and without
@@ -258,7 +258,7 @@ void run_v_cell_probe_test(const context& ctx) {
 }
 
 template <typename Backend>
-void run_expsyn_g_probe_test(const context& ctx) {
+void run_expsyn_g_probe_test(context ctx) {
     using fvm_cell = typename backend_access<Backend>::fvm_cell;
     auto deref = [](const arb_value_type* p) { return backend_access<Backend>::deref(p); };
 
@@ -353,7 +353,7 @@ void run_expsyn_g_probe_test(const context& ctx) {
 }
 
 template <typename Backend>
-void run_expsyn_g_cell_probe_test(const context& ctx) {
+void run_expsyn_g_cell_probe_test(context ctx) {
     using fvm_cell = typename backend_access<Backend>::fvm_cell;
     auto deref = [](const auto* p) { return backend_access<Backend>::deref(p); };
 
@@ -483,7 +483,7 @@ void run_expsyn_g_cell_probe_test(const context& ctx) {
 }
 
 template <typename Backend>
-void run_ion_density_probe_test(const context& ctx) {
+void run_ion_density_probe_test(context ctx) {
     using fvm_cell = typename backend_access<Backend>::fvm_cell;
     auto deref = [](const arb_value_type* p) { return backend_access<Backend>::deref(p); };
 
@@ -656,7 +656,7 @@ void run_ion_density_probe_test(const context& ctx) {
 }
 
 template <typename Backend>
-void run_partial_density_probe_test(const context& ctx) {
+void run_partial_density_probe_test(context ctx) {
     using fvm_cell = typename backend_access<Backend>::fvm_cell;
     auto deref = [](const arb_value_type* p) { return backend_access<Backend>::deref(p); };
 
@@ -769,7 +769,7 @@ void run_partial_density_probe_test(const context& ctx) {
 }
 
 template <typename Backend>
-void run_axial_and_ion_current_sampled_probe_test(const context& ctx) {
+void run_axial_and_ion_current_sampled_probe_test(context ctx) {
     // On a passive cable in steady-state, the capacitive membrane current will be zero,
     // and the axial currents should balance the stimulus and ionic membrane currents in any CV.
     //
@@ -937,7 +937,7 @@ auto run_simple_sampler(
 }
 
 template <typename Backend>
-void run_multi_probe_test(const context& ctx) {
+void run_multi_probe_test(context ctx) {
     // Construct and run thorugh simple sampler a probe defined over
     // cell terminal points; check metadata and values.
 
@@ -971,7 +971,7 @@ void run_multi_probe_test(const context& ctx) {
 }
 
 template <typename Backend>
-void run_v_sampled_probe_test(const context& ctx) {
+void run_v_sampled_probe_test(context ctx) {
     soma_cell_builder builder(12.6157/2.0);
     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");
@@ -1015,7 +1015,7 @@ void run_v_sampled_probe_test(const context& ctx) {
 
 
 template <typename Backend>
-void run_total_current_probe_test(const context& ctx) {
+void run_total_current_probe_test(context ctx) {
     // Model two passive Y-shaped cells with a similar but not identical
     // time constant Ï„.
     //
@@ -1147,7 +1147,7 @@ void run_total_current_probe_test(const context& ctx) {
 
 
 template <typename Backend>
-void run_stimulus_probe_test(const context& ctx) {
+void run_stimulus_probe_test(context ctx) {
     // Model two simple stick cable cells, 3 CVs each, and stimuli on cell 0, cv 1
     // and cell 1, cv 2. Run both cells in the same cell group.
 
@@ -1192,7 +1192,7 @@ void run_stimulus_probe_test(const context& ctx) {
 }
 
 template <typename Backend>
-void run_exact_sampling_probe_test(const context& ctx) {
+void run_exact_sampling_probe_test(context ctx) {
     // As the exact sampling implementation interacts with the event delivery
     // implementation within in cable cell groups, construct a somewhat
     // elaborate model with 4 cells and a gap junction between cell 1 and 3.
-- 
GitLab