diff --git a/CMakeLists.txt b/CMakeLists.txt
index bed72ed2c6e51a33d2579e2ab23e667eda06a651..25b81d6d6b94d4f0396c34a5c1c2879b71c8f590 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -7,6 +7,10 @@ string(REGEX MATCH "^[0-9]+(\\.[0-9]+)?(\\.[0-9]+)?(\\.[0-9]+)?" numeric_version
 project(arbor VERSION ${numeric_version})
 enable_language(CXX)
 
+# Turn on this option to force the compilers to produce color output when output is
+# redirected from the terminal (e.g. when using ninja or a pager).
+
+option(ARBDEV_COLOR "Always produce ANSI-colored output (GNU/Clang only)." OFF)
 
 #----------------------------------------------------------
 # Configure-time build options for Arbor:
@@ -48,7 +52,6 @@ option(ARB_WITH_PROFILING "use built-in profiling" OFF)
 
 option(ARB_WITH_ASSERTIONS "enable arb_assert() assertions in code" OFF)
 
-
 #----------------------------------------------------------
 # Python front end for Arbor:
 #----------------------------------------------------------
diff --git a/arbor/CMakeLists.txt b/arbor/CMakeLists.txt
index 0174af4b445c94d0116bf1c6fbfd36da64e42d0b..0368fc88613d5213bdc9582a473ad4ae6dbcf39f 100644
--- a/arbor/CMakeLists.txt
+++ b/arbor/CMakeLists.txt
@@ -33,6 +33,7 @@ set(arbor_sources
     morph/embed_pwlin.cpp
     morph/stitch.cpp
     morph/label_dict.cpp
+    morph/label_parse.cpp
     morph/locset.cpp
     morph/morphexcept.cpp
     morph/morphology.cpp
@@ -53,6 +54,7 @@ set(arbor_sources
     spike_event_io.cpp
     spike_source_cell_group.cpp
     swcio.cpp
+    s_expr.cpp
     threading/threading.cpp
     thread_private_spike_store.cpp
     tree.cpp
diff --git a/arbor/include/arbor/morph/label_parse.hpp b/arbor/include/arbor/morph/label_parse.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ecc82d533d0b33333cd2d1e01a723529b3ae4e6a
--- /dev/null
+++ b/arbor/include/arbor/morph/label_parse.hpp
@@ -0,0 +1,23 @@
+#pragma once
+
+#include <any>
+
+#include <arbor/morph/region.hpp>
+#include <arbor/arbexcept.hpp>
+#include <arbor/util/expected.hpp>
+
+namespace arb {
+
+struct label_parse_error: arb::arbor_exception {
+    label_parse_error(const std::string& msg);
+};
+
+template <typename T>
+using parse_hopefully = arb::util::expected<T, label_parse_error>;
+
+parse_hopefully<std::any> parse_label_expression(const std::string&);
+
+parse_hopefully<arb::region> parse_region_expression(const std::string& s);
+parse_hopefully<arb::locset> parse_locset_expression(const std::string& s);
+
+} // namespace arb
diff --git a/arbor/include/arbor/morph/locset.hpp b/arbor/include/arbor/morph/locset.hpp
index 170e1fbb15191c5c5c0594eca11556fe303768fd..b66f268cdfa2f2ffba309f1ead5237134f49d206 100644
--- a/arbor/include/arbor/morph/locset.hpp
+++ b/arbor/include/arbor/morph/locset.hpp
@@ -49,7 +49,7 @@ public:
     locset(mlocation_list other);
 
     // Implicitly convert string to named locset expression.
-    locset(std::string label);
+    locset(const std::string& label);
     locset(const char* label);
 
     template <typename Impl,
diff --git a/arbor/include/arbor/morph/region.hpp b/arbor/include/arbor/morph/region.hpp
index f875dcdf888b983d6779fac626e2fab6044e4049..62afa9e0541013cf2b0ed2520c4672ac8779ca99 100644
--- a/arbor/include/arbor/morph/region.hpp
+++ b/arbor/include/arbor/morph/region.hpp
@@ -59,7 +59,7 @@ public:
     region(mcable_list);
 
     // Implicitly convert string to named region expression.
-    region(std::string label);
+    region(const std::string& label);
     region(const char* label);
 
     friend mextent thingify(const region& r, const mprovider& m) {
diff --git a/arbor/include/arbor/string_literals.hpp b/arbor/include/arbor/string_literals.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..056c25231abd0078387c15b12d0cbc2207885c41
--- /dev/null
+++ b/arbor/include/arbor/string_literals.hpp
@@ -0,0 +1,12 @@
+#pragma once
+
+#include <string>
+
+namespace arb::literals {
+
+inline
+std::string operator "" _lab(const char* s, std::size_t) {
+    return std::string("\"") + s + "\"";
+}
+
+} // namespace arb
diff --git a/python/morph_parse.cpp b/arbor/morph/label_parse.cpp
similarity index 68%
rename from python/morph_parse.cpp
rename to arbor/morph/label_parse.cpp
index 33caa33eac9a331055b27e229b4c9091194f9e91..a7832d5c36b2e6579335093cd1c1a919a08995a0 100644
--- a/python/morph_parse.cpp
+++ b/arbor/morph/label_parse.cpp
@@ -1,16 +1,22 @@
 #include <any>
 #include <limits>
 
+#include <arbor/arbexcept.hpp>
+#include <arbor/morph/region.hpp>
+#include <arbor/morph/locset.hpp>
+#include <arbor/morph/label_parse.hpp>
 #include <arbor/morph/region.hpp>
 #include <arbor/morph/locset.hpp>
 
-#include "error.hpp"
+#include "arbor/util/expected.hpp"
 #include "s_expr.hpp"
-#include "morph_parse.hpp"
+#include "util/strprintf.hpp"
 
-using std::any_cast;
+namespace arb {
 
-namespace pyarb {
+label_parse_error::label_parse_error(const std::string& msg):
+    arb::arbor_exception(msg)
+{}
 
 struct nil_tag {};
 
@@ -36,24 +42,24 @@ bool match<arb::locset>(const std::type_info& info) {
 
 template <typename T>
 T eval_cast(std::any arg) {
-    return std::move(any_cast<T&>(arg));
+    return std::move(std::any_cast<T&>(arg));
 }
 
 template <>
 double eval_cast<double>(std::any arg) {
-    if (arg.type()==typeid(int)) return any_cast<int>(arg);
-    return any_cast<double>(arg);
+    if (arg.type()==typeid(int)) return std::any_cast<int>(arg);
+    return std::any_cast<double>(arg);
 }
 
 template <>
 arb::region eval_cast<arb::region>(std::any arg) {
-    if (arg.type()==typeid(arb::region)) return any_cast<arb::region>(arg);
+    if (arg.type()==typeid(arb::region)) return std::any_cast<arb::region>(arg);
     return arb::reg::nil();
 }
 
 template <>
 arb::locset eval_cast<arb::locset>(std::any arg) {
-    if (arg.type()==typeid(arb::locset)) return any_cast<arb::locset>(arg);
+    if (arg.type()==typeid(arb::locset)) return std::any_cast<arb::locset>(arg);
     return arb::ls::nil();
 }
 
@@ -189,34 +195,34 @@ std::unordered_multimap<std::string, evaluator> eval_map {
                             "'cable' with 3 arguments: (branch_id:integer prox:real dist:real)")},
     {"region",  make_call<std::string>(arb::reg::named,
                             "'region' with 1 argument: (name:string)")},
-    {"distal_interval",  make_call<arb::locset, double>(arb::reg::distal_interval,
-                            "'distal_interval' with 2 arguments: (start:locset extent:real)")},
-    {"distal_interval", make_call<arb::locset>(
+    {"distal-interval",  make_call<arb::locset, double>(arb::reg::distal_interval,
+                            "'distal-interval' with 2 arguments: (start:locset extent:real)")},
+    {"distal-interval", make_call<arb::locset>(
                             [](arb::locset ls){return arb::reg::distal_interval(std::move(ls), std::numeric_limits<double>::max());},
-                            "'distal_interval' with 1 argument: (start:locset)")},
-    {"proximal_interval",make_call<arb::locset, double>(arb::reg::proximal_interval,
-                            "'proximal_interval' with 2 arguments: (start:locset extent:real)")},
-    {"proximal_interval", make_call<arb::locset>(
+                            "'distal-interval' with 1 argument: (start:locset)")},
+    {"proximal-interval",make_call<arb::locset, double>(arb::reg::proximal_interval,
+                            "'proximal-interval' with 2 arguments: (start:locset extent:real)")},
+    {"proximal-interval", make_call<arb::locset>(
                             [](arb::locset ls){return arb::reg::proximal_interval(std::move(ls), std::numeric_limits<double>::max());},
                             "'proximal_interval' with 1 argument: (start:locset)")},
     {"complete", make_call<arb::region>(arb::reg::complete,
-                            "'super' with 1 argment: (reg:region)")},
-    {"radius_lt",make_call<arb::region, double>(arb::reg::radius_lt,
-                            "'radius_lt' with 2 arguments: (reg:region radius:real)")},
-    {"radius_le",make_call<arb::region, double>(arb::reg::radius_le,
-                            "'radius_le' with 2 arguments: (reg:region radius:real)")},
-    {"radius_gt",make_call<arb::region, double>(arb::reg::radius_gt,
-                            "'radius_gt' with 2 arguments: (reg:region radius:real)")},
-    {"radius_ge",make_call<arb::region, double>(arb::reg::radius_ge,
-                            "'radius_ge' with 2 arguments: (reg:region radius:real)")},
-    {"z_dist_from_root_lt",make_call<double>(arb::reg::z_dist_from_root_lt,
-                            "'z_dist_from_root_lt' with 1 arguments: (distance:real)")},
-    {"z_dist_from_root_le",make_call<double>(arb::reg::z_dist_from_root_le,
-                            "'z_dist_from_root_le' with 1 arguments: (distance:real)")},
-    {"z_dist_from_root_gt",make_call<double>(arb::reg::z_dist_from_root_gt,
-                            "'z_dist_from_root_gt' with 1 arguments: (distance:real)")},
-    {"z_dist_from_root_ge",make_call<double>(arb::reg::z_dist_from_root_ge,
-                            "'z_dist_from_root_ge' with 1 arguments: (distance:real)")},
+                            "'complete' with 1 argment: (reg:region)")},
+    {"radius-lt",make_call<arb::region, double>(arb::reg::radius_lt,
+                            "'radius-lt' with 2 arguments: (reg:region radius:real)")},
+    {"radius-le",make_call<arb::region, double>(arb::reg::radius_le,
+                            "'radius-le' with 2 arguments: (reg:region radius:real)")},
+    {"radius-gt",make_call<arb::region, double>(arb::reg::radius_gt,
+                            "'radius-gt' with 2 arguments: (reg:region radius:real)")},
+    {"radius-ge",make_call<arb::region, double>(arb::reg::radius_ge,
+                            "'radius-ge' with 2 arguments: (reg:region radius:real)")},
+    {"z-dist-from-root-lt",make_call<double>(arb::reg::z_dist_from_root_lt,
+                            "'z-dist-from-root-lt' with 1 arguments: (distance:real)")},
+    {"z-dist-from-root-le",make_call<double>(arb::reg::z_dist_from_root_le,
+                            "'z-dist-from-root-le' with 1 arguments: (distance:real)")},
+    {"z-dist-from-root-gt",make_call<double>(arb::reg::z_dist_from_root_gt,
+                            "'z-dist-from-root-gt' with 1 arguments: (distance:real)")},
+    {"z-dist-from-root-ge",make_call<double>(arb::reg::z_dist_from_root_ge,
+                            "'z-dist-from-root-ge' with 1 arguments: (distance:real)")},
     {"join",    make_fold<arb::region>(static_cast<arb::region(*)(arb::region, arb::region)>(arb::join),
                             "'join' with at least 2 arguments: (region region [...region])")},
     {"intersect",make_fold<arb::region>(static_cast<arb::region(*)(arb::region, arb::region)>(arb::intersect),
@@ -236,8 +242,8 @@ std::unordered_multimap<std::string, evaluator> eval_map {
                             "'proximal' with 1 argument: (reg:region)")},
     {"uniform",make_call<arb::region, int, int, int>(arb::ls::uniform,
                             "'uniform' with 4 arguments: (reg:region, first:int, last:int, seed:int)")},
-    {"on_branches",make_call<double>(arb::ls::on_branches,
-                            "'on_branches' with 1 argument: (pos:double)")},
+    {"on-branches",make_call<double>(arb::ls::on_branches,
+                            "'on-branches' with 1 argument: (pos:double)")},
     {"locset",  make_call<std::string>(arb::ls::named,
                             "'locset' with 1 argument: (name:string)")},
     {"restrict",  make_call<arb::locset, arb::region>(arb::ls::restrict,
@@ -252,13 +258,14 @@ parse_hopefully<std::any> eval(const s_expr& e);
 
 parse_hopefully<std::vector<std::any>> eval_args(const s_expr& e) {
     if (!e) return {std::vector<std::any>{}}; // empty argument list
-    const s_expr* h = &e;
     std::vector<std::any> args;
-    while (*h) {
-        auto arg = eval(h->head());
-        if (!arg) return std::move(arg.error());
-        args.push_back(std::move(*arg));
-        h = &h->tail();
+    for (auto& h: e) {
+        if (auto arg=eval(h)) {
+            args.push_back(std::move(*arg));
+        }
+        else {
+            return util::unexpected(std::move(arg.error()));
+        }
     }
     return args;
 }
@@ -297,15 +304,19 @@ std::string eval_description(const char* name, const std::vector<std::any>& args
     return msg;
 }
 
+label_parse_error parse_error(std::string const& msg, src_location loc) {
+    return {util::pprintf("error in label description at {}: {}.", loc, msg)};
+}
 // Evaluate an s expression.
-// On success the result is wrapped in util::any, where the result is one of:
+// On success the result is wrapped in std::any, where the result is one of:
 //      int         : an integer atom
 //      double      : a real atom
+//      std::string : a string atom: to be treated as a label
 //      arb::region : a region
 //      arb::locset : a locset
 //
 // If there invalid input is detected, hopefully return value contains
-// a parse_error_state with an error string and location.
+// a label_error_state with an error string and location.
 //
 // If there was an unexpected/fatal error, an exception will be thrown.
 parse_hopefully<std::any> eval(const s_expr& e) {
@@ -319,12 +330,17 @@ parse_hopefully<std::any> eval(const s_expr& e) {
             case tok::nil:
                 return {nil_tag()};
             case tok::string:
-                return {std::string(t.spelling)};
+                return std::any{std::string(t.spelling)};
+            // An arbitrary symbol in a region/locset expression is an error, and is
+            // often a result of not quoting a label correctly.
+            case tok::symbol:
+                return util::unexpected(parse_error(
+                        util::pprintf("Unexpected symbol '{}' in a region or locset definition. If '{}' is a label, it must be quoted {}{}{}", e, e, '"', e, '"'),
+                        location(e)));
             case tok::error:
-                return parse_error_state{e.atom().spelling, location(e)};
+                return util::unexpected(parse_error(e.atom().spelling, location(e)));
             default:
-                return parse_error_state{
-                    util::pprintf("Unexpected term: {}", e), location(e)};
+                return util::unexpected(parse_error(util::pprintf("Unexpected term '{}' in a region or locset definition", e), location(e)));
         }
     }
     if (e.head().is_atom()) {
@@ -356,14 +372,51 @@ parse_hopefully<std::any> eval(const s_expr& e) {
         for (auto i=matches.first; i!=matches.second; ++i) {
             msg += util::pprintf("\n  Candidate {}  {}", ++count, i->second.message);
         }
-        return parse_error_state{std::move(msg), location(e)};
+        return util::unexpected(parse_error(msg, location(e)));
     }
 
-    return parse_error_state{
-        util::pprintf("Unable to evaluate '{}': expression must be either integer, real expression of the form (op <args>)", e),
-            location(e)};
+    return util::unexpected(parse_error(
+            util::pprintf("'{}' is not either integer, real expression of the form (op <args>)", e),
+            location(e)));
+}
+
+parse_hopefully<std::any> parse_label_expression(const std::string& e) {
+    return eval(parse_s_expr(e));
+}
+
+parse_hopefully<arb::region> parse_region_expression(const std::string& s) {
+    if (auto e = eval(parse_s_expr(s))) {
+        if (e->type() == typeid(region)) {
+            return {std::move(std::any_cast<region&>(*e))};
+        }
+        if (e->type() == typeid(std::string)) {
+            return {reg::named(std::move(std::any_cast<std::string&>(*e)))};
+        }
+        return util::unexpected(
+                label_parse_error(
+                util::pprintf("Invalid region description: '{}' is neither a valid region expression or region label string.", s)));
+    }
+    else {
+        return util::unexpected(label_parse_error(std::string()+e.error().what()));
+    }
 }
 
-} // namespace pyarb
+parse_hopefully<arb::locset> parse_locset_expression(const std::string& s) {
+    if (auto e = eval(parse_s_expr(s))) {
+        if (e->type() == typeid(locset)) {
+            return {std::move(std::any_cast<locset&>(*e))};
+        }
+        if (e->type() == typeid(std::string)) {
+            return {ls::named(std::move(std::any_cast<std::string&>(*e)))};
+        }
+        return util::unexpected(
+                label_parse_error(
+                util::pprintf("Invalid locset description: '{}' is neither a valid locset expression or locset label string.", s)));
+    }
+    else {
+        return util::unexpected(label_parse_error(std::string()+e.error().what()));
+    }
+}
 
+} // namespace arb
 
diff --git a/arbor/morph/locset.cpp b/arbor/morph/locset.cpp
index 67fce43b61914bb4dee5431bf163d1641d51c565..17a27addb585c5a3a1f54becaaddbd806490bf27 100644
--- a/arbor/morph/locset.cpp
+++ b/arbor/morph/locset.cpp
@@ -3,6 +3,7 @@
 #include <numeric>
 
 #include <arbor/math.hpp>
+#include <arbor/morph/label_parse.hpp>
 #include <arbor/morph/locset.hpp>
 #include <arbor/morph/morphexcept.hpp>
 #include <arbor/morph/morphology.hpp>
@@ -212,7 +213,7 @@ mlocation_list thingify_(const most_distal_& n, const mprovider& p) {
 }
 
 std::ostream& operator<<(std::ostream& o, const most_distal_& x) {
-    return o << "(distal \"" << x.reg << "\")";
+    return o << "(distal " << x.reg << ")";
 }
 
 // Most proximal points of a region
@@ -237,7 +238,7 @@ mlocation_list thingify_(const most_proximal_& n, const mprovider& p) {
 }
 
 std::ostream& operator<<(std::ostream& o, const most_proximal_& x) {
-    return o << "(proximal \"" << x.reg << "\")";
+    return o << "(proximal " << x.reg << ")";
 }
 
 // Boundary points of a region.
@@ -584,12 +585,15 @@ locset::locset(mlocation_list ll) {
     *this = ls::location_list(std::move(ll));
 }
 
-locset::locset(std::string name) {
-    *this = ls::named(std::move(name));
+locset::locset(const std::string& desc) {
+    if (auto r=parse_locset_expression(desc)) {
+        *this = *r;
+    }
+    else {
+        throw r.error();
+    }
 }
 
-locset::locset(const char* name) {
-    *this = ls::named(name);
-}
+locset::locset(const char* label): locset(std::string(label)) {}
 
 } // namespace arb
diff --git a/arbor/morph/region.cpp b/arbor/morph/region.cpp
index e78a59b4c8da8ce974a5e0323876d1927bf7c3b8..895666dcfacf4d367796af2b6699c35b6a4416b1 100644
--- a/arbor/morph/region.cpp
+++ b/arbor/morph/region.cpp
@@ -1,7 +1,9 @@
+#include <any>
+#include <stack>
 #include <string>
 #include <vector>
-#include <stack>
 
+#include <arbor/morph/label_parse.hpp>
 #include <arbor/morph/locset.hpp>
 #include <arbor/morph/primitives.hpp>
 #include <arbor/morph/morphexcept.hpp>
@@ -9,6 +11,7 @@
 #include <arbor/morph/region.hpp>
 #include <arbor/util/optional.hpp>
 
+#include "s_expr.hpp"
 #include "util/mergeview.hpp"
 #include "util/span.hpp"
 #include "util/strprintf.hpp"
@@ -279,7 +282,9 @@ mextent thingify_(const distal_interval_& reg, const mprovider& p) {
 }
 
 std::ostream& operator<<(std::ostream& o, const distal_interval_& d) {
-    return o << "(distal_interval " << d.start << " " << d.distance << ")";
+    return d.distance==std::numeric_limits<double>::max()?
+        o << "(distal-interval " << d.start << ")":
+        o << "(distal-interval " << d.start << " " << d.distance << ")";
 }
 
 // Region comprising points up to `distance` proximal to a point in `end`.
@@ -333,7 +338,9 @@ mextent thingify_(const proximal_interval_& reg, const mprovider& p) {
 }
 
 std::ostream& operator<<(std::ostream& o, const proximal_interval_& d) {
-    return o << "(proximal_interval " << d.end << " " << d.distance << ")";
+    return d.distance==std::numeric_limits<double>::max()?
+        o << "(proximal-interval " << d.end << ")":
+        o << "(proximal-interval " << d.end << " " << d.distance << ")";
 }
 
 mextent radius_cmp(const mprovider& p, region r, double val, comp_op op) {
@@ -368,7 +375,7 @@ mextent thingify_(const radius_lt_& r, const mprovider& p) {
 }
 
 std::ostream& operator<<(std::ostream& o, const radius_lt_& r) {
-    return o << "(radius_lt " << r.reg << " " << r.val << ")";
+    return o << "(radius-lt " << r.reg << " " << r.val << ")";
 }
 
 // Region with all segments with radius less than r
@@ -386,7 +393,7 @@ mextent thingify_(const radius_le_& r, const mprovider& p) {
 }
 
 std::ostream& operator<<(std::ostream& o, const radius_le_& r) {
-    return o << "(radius_le " << r.reg << " " << r.val << ")";
+    return o << "(radius-le " << r.reg << " " << r.val << ")";
 }
 
 // Region with all segments with radius greater than r
@@ -404,7 +411,7 @@ mextent thingify_(const radius_gt_& r, const mprovider& p) {
 }
 
 std::ostream& operator<<(std::ostream& o, const radius_gt_& r) {
-    return o << "(radius_gt " << r.reg << " " << r.val << ")";
+    return o << "(radius-gt " << r.reg << " " << r.val << ")";
 }
 
 // Region with all segments with radius greater than or equal to r
@@ -422,7 +429,7 @@ mextent thingify_(const radius_ge_& r, const mprovider& p) {
 }
 
 std::ostream& operator<<(std::ostream& o, const radius_ge_& r) {
-    return o << "(radius_ge " << r.reg << " " << r.val << ")";
+    return o << "(radius-ge " << r.reg << " " << r.val << ")";
 }
 
 mextent projection_cmp(const mprovider& p, double v, comp_op op) {
@@ -451,7 +458,7 @@ mextent thingify_(const projection_lt_& r, const mprovider& p) {
 }
 
 std::ostream& operator<<(std::ostream& o, const projection_lt_& r) {
-    return o << "(projection_lt " << r.val << ")";
+    return o << "(projection-lt " << r.val << ")";
 }
 
 // Region with all segments with projection less than or equal to val
@@ -468,7 +475,7 @@ mextent thingify_(const projection_le_& r, const mprovider& p) {
 }
 
 std::ostream& operator<<(std::ostream& o, const projection_le_& r) {
-    return o << "(projection_le " << r.val << ")";
+    return o << "(projection-le " << r.val << ")";
 }
 
 // Region with all segments with projection greater than val
@@ -485,7 +492,7 @@ mextent thingify_(const projection_gt_& r, const mprovider& p) {
 }
 
 std::ostream& operator<<(std::ostream& o, const projection_gt_& r) {
-    return o << "(projection_gt " << r.val << ")";
+    return o << "(projection-gt " << r.val << ")";
 }
 
 // Region with all segments with projection greater than val
@@ -502,7 +509,7 @@ mextent thingify_(const projection_ge_& r, const mprovider& p) {
 }
 
 std::ostream& operator<<(std::ostream& o, const projection_ge_& r) {
-    return o << "(projection_ge " << r.val << ")";
+    return o << "(projection-ge " << r.val << ")";
 }
 
 region z_dist_from_root_lt(double r0) {
@@ -547,7 +554,7 @@ mextent thingify_(const named_& n, const mprovider& p) {
 }
 
 std::ostream& operator<<(std::ostream& o, const named_& x) {
-    return o << "(region " << x.name << ")";
+    return o << "(region \"" << x.name << "\")";
 }
 
 // Adds all cover points to a region.
@@ -738,13 +745,16 @@ region::region() {
 
 // Implicit constructors/converters.
 
-region::region(std::string label) {
-    *this = reg::named(std::move(label));
+region::region(const std::string& desc) {
+    if (auto r=parse_region_expression(desc)) {
+        *this = *r;
+    }
+    else {
+        throw r.error();
+    }
 }
 
-region::region(const char* label) {
-    *this = reg::named(label);
-}
+region::region(const char* label): region(std::string(label)) {}
 
 region::region(mcable c) {
     *this = reg::cable(c.branch, c.prox_pos, c.dist_pos);
diff --git a/python/s_expr.cpp b/arbor/s_expr.cpp
similarity index 53%
rename from python/s_expr.cpp
rename to arbor/s_expr.cpp
index e9db04781fb1d179c9b4d48f1fb4e7e9967af3cb..b5c9073e795dac77af8a5ba9afe319c718c29936 100644
--- a/python/s_expr.cpp
+++ b/arbor/s_expr.cpp
@@ -2,16 +2,16 @@
 #include <cstring>
 #include <string>
 #include <memory>
+#include <unordered_map>
 #include <ostream>
 #include <variant>
 #include <vector>
 
 #include <arbor/arbexcept.hpp>
 
+#include "util/strprintf.hpp"
 #include "s_expr.hpp"
-#include "strprintf.hpp"
-
-namespace pyarb {
+namespace arb {
 
 inline bool is_alphanumeric(char c) {
     return std::isdigit(c) || std::isalpha(c);
@@ -19,6 +19,32 @@ inline bool is_alphanumeric(char c) {
 inline bool is_plusminus(char c) {
     return (c=='-' || c=='+');
 }
+inline bool is_valid_symbol_char(char c) {
+    switch (c) {
+        case '+':
+        case '-':
+        case '*':
+        case '/':
+        case '@':
+        case '$':
+        case '%':
+        case '^':
+        case '&':
+        case '_':
+        case '=':
+        case '<':
+        case '>':
+        case '~':
+        case '.':
+            return true;
+        default:
+            return is_alphanumeric(c);
+    }
+}
+
+std::ostream& operator<<(std::ostream& o, const src_location& l) {
+    return o << l.line << ":" << l.column;
+}
 
 std::ostream& operator<<(std::ostream& o, const tok& t) {
     switch (t) {
@@ -27,7 +53,7 @@ std::ostream& operator<<(std::ostream& o, const tok& t) {
         case tok::rparen: return o << "rparen";
         case tok::real:   return o << "real";
         case tok::integer:return o << "integer";
-        case tok::name:   return o << "name";
+        case tok::symbol: return o << "symbol";
         case tok::string: return o << "string";
         case tok::eof:    return o << "eof";
         case tok::error:  return o << "error";
@@ -46,10 +72,9 @@ std::ostream& operator<<(std::ostream& o, const token& t) {
 // lexer
 //
 
-struct parser_error: public arb::arbor_exception {
-    int loc;
-    parser_error(const std::string& msg, int l):
-        arbor_exception(msg), loc(l)
+struct s_expr_lexer_error: public arb::arbor_internal_error {
+    s_expr_lexer_error(const std::string& msg, src_location l):
+        arbor_internal_error(util::pprintf("s-expression internal error at {}: {}", l, msg))
     {}
 };
 
@@ -62,18 +87,15 @@ static std::unordered_map<std::string, tok> keyword_to_tok = {
 };
 
 class lexer {
-    const char* data_;;
-    const char* end_;;
-    const char* current_;
-    int loc_;
+    transmogrifier line_start_;
+    transmogrifier stream_;
+    unsigned line_;
     token token_;
 
 public:
 
-    lexer(const char* s):
-        data_(s),
-        end_(data_ + std::strlen(data_)),
-        current_(data_)
+    lexer(transmogrifier begin):
+        line_start_(begin), stream_(begin), line_(0)
     {
         // Prime the first token.
         parse();
@@ -91,24 +113,37 @@ public:
 
 private:
 
-    // Consume the and return the next token in the stream.
+    src_location loc() const {
+        return src_location(line_+1, stream_-line_start_+1);
+    }
+
+    bool empty() const {
+        return *stream_ == '\0';
+    }
+
+    // Consume and return the next token in the stream.
     void parse() {
         using namespace std::string_literals;
 
-        while (current_!=end_) {
-            loc_ = current_-data_;
-            switch (*current_) {
+        while (!empty()) {
+            switch (*stream_) {
                 // end of file
                 case 0      :       // end of string
-                    token_ = {loc_, tok::eof, "eof"s};
+                    token_ = {loc(), tok::eof, "eof"s};
                     return;
 
+                // new line
+                case '\n'   :
+                    line_++;
+                    ++stream_;
+                    line_start_ = stream_;
+                    continue;
+
                 // white space
                 case ' '    :
                 case '\t'   :
                 case '\v'   :
                 case '\f'   :
-                case '\n'   :
                     character();
                     continue;   // skip to next character
 
@@ -116,14 +151,14 @@ private:
                     eat_comment();
                     continue;
                 case '(':
-                    token_ = {loc_, tok::lparen, {character()}};
+                    token_ = {loc(), tok::lparen, {character()}};
                     return;
                 case ')':
-                    token_ = {loc_, tok::rparen, {character()}};
+                    token_ = {loc(), tok::rparen, {character()}};
                     return;
                 case 'a' ... 'z':
                 case 'A' ... 'Z':
-                    token_ = name();
+                    token_ = symbol();
                     return;
                 case '0' ... '9':
                     token_ = number();
@@ -133,137 +168,157 @@ private:
                     return;
                 case '-':
                 case '+':
+                case '.':
                     {
-                        char c = current_[1];
+                        if (empty()) {
+                            token_ = {loc(), tok::error, "Unexpected end of input."};
+                        }
+                        char c = stream_.peek(1);
                         if (std::isdigit(c) or c=='.') {
                             token_ = number();
                             return;
                         }
                     }
-                    token_ = {loc_, tok::error,
+                    token_ = {loc(), tok::error,
                         util::pprintf("Unexpected character '{}'.", character())};
                     return;
 
                 default:
-                    token_ = {loc_, tok::error,
+                    token_ = {loc(), tok::error,
                         util::pprintf("Unexpected character '{}'.", character())};
                     return;
             }
         }
 
-        if (current_!=end_) {
+        if (!empty()) {
             // todo: handle error: should never hit this
         }
-        token_ = {loc_, tok::eof, "eof"s};
+        token_ = {loc(), tok::eof, "eof"s};
         return;
     }
 
     // Consumes characters in the stream until end of stream or a new line.
-    // Assumes that the current_ location is the `;` that starts the comment.
+    // Assumes that the current location is the `;` that starts the comment.
     void eat_comment() {
-        while (*current_ && *current_!='\n') {
-            character();
+        while (!empty() && *stream_!='\n') {
+            ++stream_;
         }
     }
 
     // Parse alphanumeric sequence that starts with an alphabet character,
-    // and my contain alphabet, numeric or underscore '_' characters.
+    // and my contain alphabet, numeric or one of {+ -  *  /  @  $  %  ^  &  _  =  <  >  ~ .}
     //
-    // Valid names:
+    // This definition follows the symbol naming conventions of common lisp, without the
+    // use of pipes || to define symbols with arbitrary strings.
+    //
+    // Valid symbols:
     //    sub_dendrite
     //    sub-dendrite
+    //    sub-dendrite:
+    //    foo@3.2/lower
     //    temp_
     //    branch3
     //    A
-    // Invalid names:
+    // Invalid symbols:
     //    _cat          ; can't start with underscore
     //    -cat          ; can't start with hyphen
     //    2ndvar        ; can't start with numeric character
     //
-    // Returns the appropriate token kind if name is a keyword.
-    token name() {
-        std::string name;
-        char c = *current_;
+    // Returns the appropriate token kind if symbol is a keyword.
+    token symbol() {
+        auto start = loc();
+        std::string symbol;
+        char c = *stream_;
 
         // Assert that current position is at the start of an identifier
         if( !(std::isalpha(c)) ) {
-            throw parser_error(
-                "Lexer attempting to read identifier when none is available", loc_);
+            throw s_expr_lexer_error(
+                "Lexer attempting to read identifier when none is available", loc());
         }
 
-        name += c;
-        ++current_;
+        symbol += c;
+        ++stream_;
         while(1) {
-            c = *current_;
+            c = *stream_;
 
-            if(is_alphanumeric(c) || c=='_' || c=='-') {
-                name += character();
+            if(is_valid_symbol_char(c)) {
+                symbol += c;
+                ++stream_;
             }
             else {
                 break;
             }
         }
 
-        // test if the name matches a keyword
-        auto it = keyword_to_tok.find(name.c_str());
+        // test if the symbol matches a keyword
+        auto it = keyword_to_tok.find(symbol.c_str());
         if (it!=keyword_to_tok.end()) {
-            return {loc_, it->second, std::move(name)};
+            return {start, it->second, std::move(symbol)};
         }
-        return {loc_, tok::name, std::move(name)};
+        return {start, tok::symbol, std::move(symbol)};
     }
 
     token string() {
         using namespace std::string_literals;
+        if (*stream_ != '"') {
+            s_expr_lexer_error(
+                "Lexer attempting to read string without opening \"", loc());
+        }
 
-        ++current_;
-        const char* begin = current_;
-        while (current_!=end_ && character()!='"');
-
-        if (current_==end_) return {loc_, tok::error, "string missing closing \""};
+        auto start = loc();
+        ++stream_;
+        std::string str;
+        while (!empty() && *stream_!='"') {
+            str.push_back(*stream_);
+            ++stream_;
+        }
+        if (empty()) return {start, tok::error, "string missing closing \""};
+        ++stream_; // gobble the closing "
 
-        return {loc_, tok::string, std::string(begin, current_-1)};
+        return {start, tok::string, str};
     }
 
     token number() {
         using namespace std::string_literals;
 
+        auto start = loc();
         std::string str;
-        char c = *current_;
+        char c = *stream_;
 
         // Start counting the number of points in the number.
         auto num_point = (c=='.' ? 1 : 0);
         auto uses_scientific_notation = 0;
 
         str += c;
-        current_++;
+        ++stream_;
         while(1) {
-            c = *current_;
-            if(std::isdigit(c)) {
+            c = *stream_;
+            if (std::isdigit(c)) {
                 str += c;
-                current_++;
+                ++stream_;
             }
-            else if(c=='.') {
+            else if (c=='.') {
                 if (++num_point>1) {
                     // Can't have more than one '.' in a number
-                    return {int(current_-data_), tok::error, "unexpected '.'"s};
+                    return {start, tok::error, "unexpected '.'"s};
                 }
                 str += c;
-                current_++;
-                if(uses_scientific_notation) {
+                ++stream_;
+                if (uses_scientific_notation) {
                     // Can't have a '.' in the mantissa
-                    return {int(current_-data_), tok::error, "unexpected '.'"s};
+                    return {start, tok::error, "unexpected '.'"s};
                 }
             }
-            else if(!uses_scientific_notation && (c=='e' || c=='E')) {
-                if(std::isdigit(current_[1]) ||
-                   (is_plusminus(current_[1]) && std::isdigit(current_[2])))
+            else if (!uses_scientific_notation && (c=='e' || c=='E')) {
+                if ( std::isdigit(stream_.peek(1)) ||
+                    (is_plusminus(stream_.peek(1)) && std::isdigit(stream_.peek(2))))
                 {
                     uses_scientific_notation++;
                     str += c;
-                    current_++;
+                    stream_++;
                     // Consume the next char if +/-
-                    if (is_plusminus(*current_)) {
-                        str += *current_++;
+                    if (is_plusminus(*stream_)) {
+                        str += *stream_++;
                     }
                 }
                 else {
@@ -277,20 +332,14 @@ private:
         }
 
         const bool is_real = uses_scientific_notation || num_point>0;
-        return {loc_, (is_real? tok::real: tok::integer), std::move(str)};
+        return {start, (is_real? tok::real: tok::integer), std::move(str)};
     }
 
     char character() {
-        return *current_++;
+        return *stream_++;
     }
 };
 
-bool test_identifier(const char* in) {
-    lexer L(in);
-    auto x = L.current();
-    return x.kind==tok::name && x.spelling==in;
-}
-
 //
 // s expression members
 //
@@ -325,27 +374,33 @@ s_expr::operator bool() const {
 
 std::ostream& operator<<(std::ostream& o, const s_expr& x) {
     if (x.is_atom()) return o << x.atom();
-#if 0 // print full tree with terminating 'nil'
+#if 1
+    o << "(";
+    bool first = true;
+    for (auto& e: x) {
+        o << (first? "": " ") << e;
+        first = false;
+    }
+    return o << ")";
+#else
     return o << "(" << x.head() << " . " << x.tail() << ")";
-#else // print '(a . nil)' as 'a'
-    return x.tail()? o << "(" << x.head() << " . " << x.tail() << ")"
-                   : o << x.head();
 #endif
 }
 
 std::size_t length(const s_expr& l) {
+    // The length of an atom is 1.
     if (l.is_atom() && l) {
-        throw arb::arbor_internal_error(
-            util::pprintf("Internal error: can't take length of an atom in '{}'.", l));
+        return 1;
     }
-    if (!l) { // nil
+    // nil marks the end of a list.
+    if (!l) {
         return 0u;
     }
     return 1+length(l.tail());
 }
 
-int location(const s_expr& l) {
-    if (l.is_atom()) return l.atom().column;
+src_location location(const s_expr& l) {
+    if (l.is_atom()) return l.atom().loc;
     return location(l.head());
 }
 
@@ -353,6 +408,8 @@ int location(const s_expr& l) {
 // parsing s expressions
 //
 
+namespace impl {
+
 // If there is a parsing error, then an atom with kind==tok::error is returned
 // with the error string in its spelling.
 s_expr parse(lexer& L) {
@@ -366,53 +423,81 @@ s_expr parse(lexer& L) {
         s_expr* n = &node;
         while (true) {
             if (t.kind == tok::eof) {
-                return token{t.column, tok::error,
-                    "Unexpected end of input. Missing a closing parenthesis ')'."};;
+                return token{t.loc, tok::error,
+                    "Unexpected end of input. Missing a closing parenthesis ')'."};
             }
             if (t.kind == tok::error) {
                 return t;
             }
             else if (t.kind == tok::rparen) {
-                *n = token{t.column, tok::nil, "nil"};
+                *n = token{t.loc, tok::nil, "nil"};
+                t = L.next();
                 break;
             }
             else if (t.kind == tok::lparen) {
                 auto e = parse(L);
                 if (e.is_atom() && e.atom().kind==tok::error) return e;
                 *n = {std::move(e), {}};
+                t = L.current();
             }
             else {
                 *n = {s_expr(t), {}};
+                t = L.next();
             }
 
             n = &n->tail();
-            t = L.next();
         }
     }
+    else if (t.kind==tok::eof) {
+        return token{t.loc, tok::error, "Empty expression."};
+    }
+    else if (t.kind==tok::rparen) {
+        return token{t.loc, tok::error, "Missing opening parenthesis'('."};
+    }
+    // an atom or an error
     else {
-        return token{t.column, tok::error, "Missing opening parenthesis'('."};;
+        L.next(); // advance the lexer to the next token
+        return t;
     }
 
     return node;
 }
 
-s_expr parse(const char* in) {
-    lexer l(in);
-    s_expr result = parse(l);
+}
+
+s_expr parse_s_expr(transmogrifier begin) {
+    lexer l(begin);
+    s_expr result = impl::parse(l);
     const bool err = result.is_atom()? result.atom().kind==tok::error: false;
     if (!err) {
-        auto t = l.next(); // pop the last rparen token.
+        auto t = l.current();
         if (t.kind!=tok::eof) {
-            return token{t.column, tok::error,
+            return token{t.loc, tok::error,
                          util::pprintf("Unexpected '{}' at the end of input.", t)};
         }
     }
     return result;
 }
 
-s_expr parse(const std::string& in) {
-    return parse(in.c_str());
+s_expr parse_s_expr(const std::string& in) {
+    return parse_s_expr(transmogrifier{in});
+}
+
+// For parsing a file with multiple high level s expressions.
+// Returns a vector of the expressions.
+// If an error occured, terminate early and the last expression will be an error.
+std::vector<s_expr> parse_multi_s_expr(transmogrifier begin) {
+    std::vector<s_expr> result;
+    lexer l(begin);
+    bool error = false;
+    while (!error && l.current().kind!=tok::eof) {
+        result.push_back(impl::parse(l));
+        const auto& e = result.back();
+        error = e.is_atom() && e.atom().kind==tok::error;
+    }
+
+    return result;
 }
 
-} // namespace pyarb
 
+} // namespace arb
diff --git a/arbor/s_expr.hpp b/arbor/s_expr.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..edd1a7fa27a81872cdab69761fe285f625040f3c
--- /dev/null
+++ b/arbor/s_expr.hpp
@@ -0,0 +1,356 @@
+#pragma once
+
+#include <cstddef>
+#include <cstring>
+#include <iterator>
+#include <stdexcept>
+#include <string>
+#include <memory>
+#include <type_traits>
+#include <unordered_map>
+#include <variant>
+#include <vector>
+
+namespace arb {
+
+// Forward iterator that can translate a raw stream to valid s_expr input if,
+// perchance, you want to parse a half-hearted attempt at an s-expression
+// (looking at you, the guy who invented the Neurolucida .asc format).
+//
+// I am not fond of .asc files, which would be s-expressions, if they
+// didn't sometimes contain '|' and ',' characters which translate to ")(" and
+// " " respectively (not mentioning the spine syntax).
+//
+// To remedy such situations, the transmogrifier performs user-provided string
+// substitution on a characters in the input.
+//
+// For example, if you are unfortuinate enough to parse an asc file, you might want
+// to try the following:
+//
+// transmogrifier(str, {{',', " "},
+//                      {'|', ")("},
+//                      {'<', "(spine "},
+//                      {'>', ")"}});
+
+class transmogrifier {
+    using sub_map = std::unordered_map<char, std::string>;
+    using iterator_type = std::string::const_iterator;
+    using difference_type = std::string::difference_type;
+    using iterator = transmogrifier;
+
+    iterator_type pos_;
+    iterator_type end_;
+    sub_map sub_map_;
+
+    const char* sub_pos_ = nullptr;
+
+    void set_state() {
+        sub_pos_ = nullptr;
+        char next = *pos_;
+        if (auto it=sub_map_.find(next); it!=sub_map_.end()) {
+            sub_pos_ = it->second.c_str();
+        }
+    }
+
+    public:
+
+    transmogrifier(const std::string& s, sub_map map={}):
+        pos_(s.cbegin()),
+        end_(s.cend()),
+        sub_map_(std::move(map))
+    {
+        if (pos_!=end_) {
+            set_state();
+        }
+    }
+
+    char operator*() const {
+        if (pos_==end_) {
+            return '\0';
+        }
+        if (sub_pos_) {
+            return *sub_pos_;
+        }
+        return *pos_;
+    }
+
+    iterator& operator++() {
+        // If already at the end don't advance.
+        if (pos_==end_) {
+            return *this;
+        }
+
+        // If currently substituting a string, advance by one and
+        // test whether we have reached the end of the string.
+        if (sub_pos_) {
+            ++sub_pos_;
+            if (*sub_pos_=='\0') { // test for end of string
+                sub_pos_ = nullptr;
+            }
+            else {
+                return *this;
+            }
+        }
+
+        ++pos_;
+
+        set_state();
+        return *this;
+    }
+
+    iterator operator++(int) {
+        iterator it = *this;
+
+        ++(*this);
+
+        return it;
+    }
+
+    iterator operator+(unsigned n) {
+        iterator it = *this;
+
+        while (n--) ++it;
+
+        return it;
+    }
+
+    char peek(unsigned i) {
+        return *(*this+i);
+    }
+
+    bool operator==(const transmogrifier& other) const {
+        return pos_==other.pos_ && sub_pos_==other.sub_pos_;
+    }
+
+    bool operator!=(const transmogrifier& other) {
+        return !(*this==other);
+    }
+
+    operator bool() const {
+        return pos_ != end_;
+    }
+
+    difference_type operator-(const transmogrifier& rhs) const {
+        return pos_ - rhs.pos_;
+    }
+};
+
+struct src_location {
+    unsigned line = 0;
+    unsigned column = 0;
+
+    src_location() = default;
+
+    src_location(unsigned l, unsigned c):
+        line(l), column(c)
+    {}
+};
+
+std::ostream& operator<<(std::ostream& o, const src_location& l);
+
+enum class tok {
+    nil,
+    real,       // real number
+    integer,    // integer
+    symbol,     // symbol
+    lparen,     // left parenthesis '('
+    rparen,     // right parenthesis ')'
+    string,     // string, written as "spelling"
+    eof,        // end of file/input
+    error       // special error state marker
+};
+
+std::ostream& operator<<(std::ostream&, const tok&);
+
+struct token {
+    src_location loc;
+    tok kind;
+    std::string spelling;
+};
+
+std::ostream& operator<<(std::ostream&, const token&);
+
+struct s_expr {
+    template <typename U>
+    struct s_pair {
+        U head = U();
+        U tail = U();
+        s_pair(U l, U r): head(std::move(l)), tail(std::move(r)) {}
+    };
+
+    // This value_wrapper is used to wrap the shared pointer
+    template <typename T>
+    struct value_wrapper{
+        using state_t = std::unique_ptr<T>;
+        state_t state;
+
+        value_wrapper() = default;
+
+        value_wrapper(const T& v):
+            state(std::make_unique<T>(v)) {}
+
+        value_wrapper(T&& v):
+            state(std::make_unique<T>(std::move(v))) {}
+
+        value_wrapper(const value_wrapper& other):
+            state(std::make_unique<T>(other.get())) {}
+
+        value_wrapper& operator=(const value_wrapper& other) {
+            state = std::make_unique<T>(other.get());
+            return *this;
+        }
+
+        value_wrapper(value_wrapper&& other) = default;
+
+        friend std::ostream& operator<<(std::ostream& o, const value_wrapper& w) {
+            return o << *w.state;
+        }
+
+        operator T() const {
+            return *state;
+        }
+
+        const T& get() const {
+            return *state;
+        }
+
+        T& get() {
+            return *state;
+        }
+    };
+
+    template <bool Const>
+    class s_expr_iterator_impl {
+        public:
+
+        struct sentinel {};
+
+        using value_type = s_expr;
+        using iterator_category = std::forward_iterator_tag;
+        using difference_type = std::ptrdiff_t;
+        using pointer   = std::conditional_t<Const, const s_expr*, s_expr*>;
+        using reference = std::conditional_t<Const, const s_expr&, s_expr&>;
+
+        s_expr_iterator_impl(reference e):
+            inner_(&e)
+        {
+            if (inner_->is_atom()) {
+                throw std::runtime_error("Attempt to create s_expr_iterator on an atom.");
+            }
+            if (finished()) inner_ = nullptr;
+        }
+
+        s_expr_iterator_impl(const sentinel& e):
+            inner_(nullptr)
+        {}
+
+        reference operator*() const {
+            return inner_->head();
+        }
+
+        pointer operator->() const {
+            return &inner_->head();
+        }
+
+        s_expr_iterator_impl& operator++() {
+            advance();
+            return *this;
+        }
+
+        s_expr_iterator_impl operator++(int) {
+            s_expr_iterator_impl cur = *this;
+            advance();
+            return cur;
+        }
+
+        s_expr_iterator_impl operator+(difference_type i) const {
+            s_expr_iterator_impl it = *this;
+            while (i--) {
+                ++it;
+            }
+            return it;
+        }
+        bool operator==(const s_expr_iterator_impl& other) const {
+            return inner_==other.inner_;
+        }
+        bool operator!=(const s_expr_iterator_impl& other) const {
+            return !(*this==other);
+        }
+        bool operator==(const sentinel& other) const {
+            return !inner_;
+        }
+        bool operator!=(const sentinel& other) const {
+            return !(*this==other);
+        }
+
+        reference expression() const {
+            return *inner_;
+        }
+
+        private:
+
+        bool finished() const {
+            return inner_->is_atom() && inner_->atom().kind==tok::nil;
+        }
+
+        void advance() {
+            if (!inner_) return;
+            inner_ = &inner_->tail();
+            if (finished()) inner_ = nullptr;
+        }
+
+        // Pointer to the current s_expr.
+        // Set to nullptr when at the end of the range.
+        pointer inner_;
+    };
+
+    using iterator       = s_expr_iterator_impl<false>;
+    using const_iterator = s_expr_iterator_impl<true>;
+
+    // An s_expr can be one of
+    //      1. an atom
+    //      2. a pair of s_expr (head and tail)
+    // The s_expr uses a util::variant to represent these two possible states,
+    // which requires using an incomplete definition of s_expr, requiring
+    // with a std::unique_ptr via value_wrapper.
+
+    using pair_type = s_pair<value_wrapper<s_expr>>;
+    std::variant<token, pair_type> state = token{{0,0}, tok::nil, "nil"};
+
+    s_expr(const s_expr& s): state(s.state) {}
+    s_expr() = default;
+    s_expr(token t): state(std::move(t)) {}
+    s_expr(s_expr l, s_expr r):
+        state(pair_type(std::move(l), std::move(r)))
+    {}
+
+    bool is_atom() const;
+
+    const token& atom() const;
+
+    operator bool() const;
+
+    const s_expr& head() const;
+    const s_expr& tail() const;
+    s_expr& head();
+    s_expr& tail();
+
+    iterator       begin()        { return {*this}; }
+    iterator       end()          { return iterator::sentinel{}; }
+    const_iterator begin()  const { return {*this}; }
+    const_iterator end()    const { return const_iterator::sentinel{}; }
+    const_iterator cbegin() const { return {*this}; }
+    const_iterator cend()   const { return const_iterator::sentinel{}; }
+
+    friend std::ostream& operator<<(std::ostream& o, const s_expr& x);
+};
+
+std::size_t length(const s_expr& l);
+src_location location(const s_expr& l);
+
+s_expr parse_s_expr(const std::string& line);
+s_expr parse_s_expr(transmogrifier begin);
+std::vector<s_expr> parse_multi(transmogrifier begin);
+
+} // namespace arb
+
diff --git a/cmake/CompilerOptions.cmake b/cmake/CompilerOptions.cmake
index 8ff6415acc9bbc11b0d6fe6f1afc2e10ccf4a116..6734beaf84966cee8607f9f74aed469c69ae2981 100644
--- a/cmake/CompilerOptions.cmake
+++ b/cmake/CompilerOptions.cmake
@@ -8,6 +8,15 @@ if(CMAKE_CXX_COMPILER_ID MATCHES "XL")
     string(REPLACE "-qhalt=e" "" CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}")
 endif()
 
+
+if(${ARBDEV_COLOR})
+    set(colorflags
+        $<IF:$<CXX_COMPILER_ID:Clang>,-fcolor-diagnostics,>
+        $<IF:$<CXX_COMPILER_ID:AppleClang>,-fcolor-diagnostics,>
+        $<IF:$<CXX_COMPILER_ID:GNU>,-fdiagnostics-color=always,>)
+    add_compile_options("$<$<COMPILE_LANGUAGE:CXX>:${colorflags}>")
+endif()
+
 # Warning options: disable specific spurious warnings as required.
 
 set(CXXOPT_WALL
diff --git a/doc/cable_cell.rst b/doc/cable_cell.rst
index 1f0cebef76fc3f6d700c17b3b8f2d2fd098054b0..db3f2ee2365947ef00c0860429fa7dae9cc27167 100644
--- a/doc/cable_cell.rst
+++ b/doc/cable_cell.rst
@@ -120,8 +120,8 @@ for setting cell-wide defaults for properties, and the
     cell.set_properties(Vm=-70, cm=0.02, rL=30, tempK=30+273.5)
 
     # Override specific values on the soma and axon
-    cell.paint('soma', Vm=-50, cm=0.01, rL=35)
-    cell.paint('axon', Vm=-60, rL=40)
+    cell.paint('"soma"', Vm=-50, cm=0.01, rL=35)
+    cell.paint('"axon"', Vm=-60, rL=40)
 
 .. _cable-density-mechs:
 
@@ -171,9 +171,9 @@ Take for example a mechanism passive leaky dynamics:
     # Create an instance of the same mechanism, that also sets conductance (range)
     m4 = arbor.mechanism('passive/el=-45', {'g': 0.1})
 
-    cell.paint('soma', m1)
-    cell.paint('soma', m2) # error: can't place the same mechanism on overlapping regions
-    cell.paint('soma', m3) # error: technically a different mechanism?
+    cell.paint('"soma"', m1)
+    cell.paint('"soma"', m2) # error: can't place the same mechanism on overlapping regions
+    cell.paint('"soma"', m3) # error: technically a different mechanism?
 
 .. _cable-ions:
 
@@ -260,10 +260,10 @@ using the *paint* interface:
 
     # It is possible to define all of the initial condition values
     # for a ion species.
-    cell.paint('soma', arbor.ion('ca', int_con=2e-4, ext_con=2.5, rev_pot=114))
+    cell.paint('(tag 1)', arbor.ion('ca', int_con=2e-4, ext_con=2.5, rev_pot=114))
 
     # Alternatively, one can selectively overwrite the global defaults.
-    cell.paint('axon', arbor.ion('ca', rev_pot=126)
+    cell.paint('(tag 2)', arbor.ion('ca', rev_pot=126)
 
 .. _cablecell-place:
 
@@ -378,9 +378,10 @@ Creating a cable cell
         .. code-block:: Python
 
             # Specialize resistivity on soma
-            cell.paint('soma', rL=100)
-            # Specialize resistivity and capacitance on the axon
-            cell.paint('axon', cm=0.05, rL=80)
+            cell.paint('"soma"', rL=100)
+            # Specialize resistivity and capacitance on the axon, where
+            # axon is defined using a region expression.
+            cell.paint('(tag 2)', cm=0.05, rL=80)
 
 .. py:class:: ion
 
diff --git a/doc/labels.rst b/doc/labels.rst
index dce13697842b7788390c3abaa10918952882feab..a539acbb7144c1b8df9ee667703df6cf7c05d354 100644
--- a/doc/labels.rst
+++ b/doc/labels.rst
@@ -135,7 +135,7 @@ from simple expressions. For example, the expression:
 
 .. code-block:: lisp
 
-    (radius_lt (join (tag 3) (tag 4)) 0.5)
+    (radius-lt (join (tag 3) (tag 4)) 0.5)
 
 describes the region of all parts of a cell with either tag 3 or tag 4 and radius less than 0.5 μm.
 
@@ -188,9 +188,9 @@ dendritic tree where the radius first is less than or equal to 0.2 μm.
 
 .. code:: lisp
 
-    (distal_interval                   ; take subtrees that start at
+    (distal-interval                   ; take subtrees that start at
         (proximal                      ; locations closest to the soma
-            (radius_lte                ; with radius <= 0.2 um
+            (radius-le                 ; with radius <= 0.2 um
                 (join (tag 3) (tag 4)) ; on basal and apical dendrites
                 0.2)))
 
@@ -256,7 +256,7 @@ Locset Expressions
       On the left with  a seed of 0: ``(uniform (tag 3) 0 9 0)``,
       and on the right with  a seed of 1: ``(uniform (tag 3) 0 9 1)``.
 
-.. label:: (on_branches pos:double)
+.. label:: (on-branches pos:double)
 
     The set of locations ``{(location b pos) | 0 ≤ b < nbranch-1}``.
 
@@ -264,7 +264,7 @@ Locset Expressions
       :width: 300
       :align: center
 
-      The set of locations at the midpoint of every branch, expressed as ``(on_branches 0.5)``.
+      The set of locations at the midpoint of every branch, expressed as ``(on-branches 0.5)``.
 
 .. label:: (distal reg:region)
 
@@ -407,7 +407,7 @@ Region Expressions
 
     Refer to a region by its label. For example, ``(region "axon")`` would refer to a region with the label ``"axon"``.
 
-.. label:: (distal_interval start:locset extent:real)
+.. label:: (distal-interval start:locset extent:real)
 
     The distal interval of a location is the region that contains all points that are distal to the location,
     and up to ``extent`` μm from the location, measured as the distance traversed along cables between two locations.
@@ -422,9 +422,9 @@ Region Expressions
 
       .. code-block:: lisp
 
-        (distal_interval (sum (location 1 0.5) (location 2 0.7) (location 5 0.1)) 5)
+        (distal-interval (sum (location 1 0.5) (location 2 0.7) (location 5 0.1)) 5)
 
-.. label:: (distal_interval start:locset)
+.. label:: (distal-interval start:locset)
 
     When no ``extent`` distance is provided, the distal intervals are extended to all terminal
     locations that are distal to each location in ``start``.
@@ -438,10 +438,10 @@ Region Expressions
 
       .. code-block:: lisp
 
-        (distal_interval (sum (location 1 0.5) (location 2 0.7) (location 5 0.1)))
+        (distal-interval (sum (location 1 0.5) (location 2 0.7) (location 5 0.1)))
 
 
-.. label:: (proximal_interval start:locset extent:real)
+.. label:: (proximal-interval start:locset extent:real)
 
     The proximal interval of a location is the region that contains all points that are proximal to the location,
     and up to ``extent`` μm from the location, measured as the distance traversed along cables between two locations.
@@ -456,9 +456,9 @@ Region Expressions
 
       .. code-block:: lisp
 
-        (proximal_interval (sum (location 1 0.8) (location 2 0.3)) 5)
+        (proximal-interval (sum (location 1 0.8) (location 2 0.3)) 5)
 
-.. label:: (proximal_interval start:locset)
+.. label:: (proximal-interval start:locset)
 
     When no ``extent`` distance is provided, the proximal intervals are extended to the root location.
 
@@ -471,9 +471,9 @@ Region Expressions
 
       .. code-block:: lisp
 
-        (proximal_interval (sum (location 1 0.8) (location 2 0.3)))
+        (proximal-interval (sum (location 1 0.8) (location 2 0.3)))
 
-.. label:: (radius_lt reg:region radius:real)
+.. label:: (radius-lt reg:region radius:real)
 
     All parts of cable segments in the region ``reg`` with radius less than ``radius``.
 
@@ -481,16 +481,16 @@ Region Expressions
       :width: 300
       :align: center
 
-      All cable segments with radius **less than** 0.5 μm, found by applying ``radius_lt`` to all of
+      All cable segments with radius **less than** 0.5 μm, found by applying ``radius-lt`` to all of
       the cables in the morphology.
       Note that branch 2, which has a constant radius of 0.5 μm, is not in the result because its radius
       is not strictly less than 0.5 μm.
 
       .. code-block:: lisp
 
-        (radius_lt (all) 0.5)
+        (radius-lt (all) 0.5)
 
-.. label:: (radius_le reg:region radius:real)
+.. label:: (radius-le reg:region radius:real)
 
     All parts of cable segments in the region ``reg`` with radius less than or equal to ``radius``.
 
@@ -498,15 +498,15 @@ Region Expressions
       :width: 300
       :align: center
 
-      All cable segments with radius **less than or equal to** 0.5 μm, found by applying ``radius_le`` to all of
+      All cable segments with radius **less than or equal to** 0.5 μm, found by applying ``radius-le`` to all of
       the cables in the morphology.
       Note that branch 2, which has a constant radius of 0.5 μm, is in the result.
 
       .. code-block:: lisp
 
-        (radius_le (all) 0.5)
+        (radius-le (all) 0.5)
 
-.. label:: (radius_gt reg:region radius:real)
+.. label:: (radius-gt reg:region radius:real)
 
     All parts of cable segments in the region ``reg`` with radius greater than ``radius``.
 
@@ -514,16 +514,16 @@ Region Expressions
       :width: 300
       :align: center
 
-      All cable segments with radius **greater than** 0.5 μm, found by applying ``radius_ge`` to all of
+      All cable segments with radius **greater than** 0.5 μm, found by applying ``radius-ge`` to all of
       the cables in the morphology.
       Note that branch 2, which has a constant radius of 0.5 μm, is not in the result because its radius
       is not strictly greater than 0.5 μm.
 
       .. code-block:: lisp
 
-        (radius_gt (all) 0.5)
+        (radius-gt (all) 0.5)
 
-.. label:: (radius_ge reg:region radius:real)
+.. label:: (radius-ge reg:region radius:real)
 
     All parts of cable segments in the region ``reg`` with radius greater than or equal to ``radius``.
 
@@ -531,13 +531,13 @@ Region Expressions
       :width: 300
       :align: center
 
-      All cable segments with radius **greater than or equal to** 0.5 μm, found by applying ``radius_le`` to all of
+      All cable segments with radius **greater than or equal to** 0.5 μm, found by applying ``radius-le`` to all of
       the cables in the morphology.
       Note that branch 2, which has a constant radius of 0.5 μm, is in the result.
 
       .. code-block:: lisp
 
-        (radius_ge (all) 0.5)
+        (radius-ge (all) 0.5)
 
 .. label:: (join lhs:region rhs:region [...region])
 
@@ -708,18 +708,18 @@ label that has not yet been defined:
 
     d = arbor.label_dict()
     # 'reg' refers 
-    d['reg'] = '(distal_interval (locset "loc"))'
+    d['reg'] = '(distal-interval (locset "loc"))'
     d['loc'] = '(location 3 0.5)'
 
     # If d was applied to a morphology, 'reg' would refer to the region:
-    #   '(distal_interval (location 3 0.5))'
+    #   '(distal-interval (location 3 0.5))'
     # Which is the sub-tree of the matrix starting at '(location 3 0.5)'
 
     # The locset 'loc' can be redefined
     d['loc'] = '(proximal (tag 3))'
 
     # Now if d was applied to a morphology, 'reg' would refer to:
-    #   '(distal_interval (proximal (tag 3))'
+    #   '(distal-interval (proximal (tag 3))'
     # Which is the subtrees that start at the proximal locations of
     # the region '(tag 3)'
 
@@ -731,7 +731,7 @@ two labels refer to one another:
     import arbor
 
     d = arbor.label_dict()
-    d['reg'] = '(distal_interval (locset "loc"))'
+    d['reg'] = '(distal-interval (locset "loc"))'
     d['loc'] = '(proximal (region "reg"))'
 
     # Error: 'reg' needs the definition of 'loc', which in turn needs the
diff --git a/doc/mechanisms.rst b/doc/mechanisms.rst
index aa2966a3e4aa776d0525c2426977ce4559299f92..46dcf8e0964b4ad376929fcdebbca17f9ef2bce8 100644
--- a/doc/mechanisms.rst
+++ b/doc/mechanisms.rst
@@ -472,8 +472,8 @@ mechanism that is to be painted or placed on the cable cell.
 
         # Decorate the 'soma' on a cable_cell.
 
-        cell.paint('soma', m1)
-        cell.paint('soma', m2) # Error: can't place the same mechanism on overlapping regions
-        cell.paint('soma', m3) # This would be ok: m3 is a new, derived mechanism by virtue of
-                               # having a different name, i.e. 'passive/el=-45' vs. 'passive'.
+        cell.paint('"soma"', m1)
+        cell.paint('"soma"', m2) # Error: can't place the same mechanism on overlapping regions
+        cell.paint('"soma"', m3) # This would be ok: m3 is a new, derived mechanism by virtue of
+                                 # having a different name, i.e. 'passive/el=-45' vs. 'passive'.
 
diff --git a/doc/scripts/gen-labels.py b/doc/scripts/gen-labels.py
index fd1bc4c6bd3d90959a215fbf3eb86580821ee008..dac99f9863afd0a4f16a53c347f9f6655c0e6e5b 100644
--- a/doc/scripts/gen-labels.py
+++ b/doc/scripts/gen-labels.py
@@ -166,20 +166,20 @@ regions  = {
             'soma': '(region "tag1")',
             'axon': '(region "tag2")',
             'dend': '(join (region "tag3") (region "tag4"))',
-            'radlt5': '(radius_lt (all) 0.5)',
-            'radle5': '(radius_le (all) 0.5)',
-            'radgt5': '(radius_gt (all) 0.5)',
-            'radge5': '(radius_ge (all) 0.5)',
-            'rad36':  '(intersect (radius_gt (all) 0.3) (radius_lt (all) 0.6))',
+            'radlt5': '(radius-lt (all) 0.5)',
+            'radle5': '(radius-le (all) 0.5)',
+            'radgt5': '(radius-gt (all) 0.5)',
+            'radge5': '(radius-ge (all) 0.5)',
+            'rad36':  '(intersect (radius-gt (all) 0.3) (radius-lt (all) 0.6))',
             'branch0': '(branch 0)',
             'branch3': '(branch 3)',
             'cable_1_01': '(cable 1 0 1)',
             'cable_1_31': '(cable 1 0.3 1)',
             'cable_1_37': '(cable 1 0.3 0.7)',
-            'proxint':     '(proximal_interval (locset "proxint_in") 5)',
-            'proxintinf':  '(proximal_interval (locset "proxint_in"))',
-            'distint':     '(distal_interval   (locset "distint_in") 5)',
-            'distintinf':  '(distal_interval   (locset "distint_in"))',
+            'proxint':     '(proximal-interval (locset "proxint_in") 5)',
+            'proxintinf':  '(proximal-interval (locset "proxint_in"))',
+            'distint':     '(distal-interval   (locset "distint_in") 5)',
+            'distintinf':  '(distal-interval   (locset "distint_in"))',
             'lhs' : '(join (cable 0 0.5 1) (cable 1 0 0.5))',
             'rhs' : '(branch 1)',
             'and': '(intersect (region "lhs") (region "rhs"))',
@@ -192,7 +192,7 @@ locsets = {
             'loc15': '(location 1 0.5)',
             'uniform0': '(uniform (tag 3) 0 9 0)',
             'uniform1': '(uniform (tag 3) 0 9 1)',
-            'branchmid': '(on_branches 0.5)',
+            'branchmid': '(on-branches 0.5)',
             'distal':  '(distal   (region "rad36"))',
             'proximal':'(proximal (region "rad36"))',
             'distint_in': '(sum (location 1 0.5) (location 2 0.7) (location 5 0.1))',
@@ -231,12 +231,12 @@ f.write(write_morphology('ysoma_morph3',   ysoma_morph3))
 
 f.write('\n############# locsets\n\n')
 for label in locsets:
-    locs = [(l.branch, l.pos) for l in cell.locations(label)]
+    locs = [(l.branch, l.pos) for l in cell.locations('"{}"'.format(label))]
     f.write('ls_{}  = {{\'type\': \'locset\', \'value\': {}}}\n'.format(label, locs))
 
 f.write('\n############# regions\n\n')
 for label in regions:
-    comps = [(c.branch, c.prox, c.dist) for c in cell.cables(label)]
+    comps = [(c.branch, c.prox, c.dist) for c in cell.cables('"{}"'.format(label))]
     f.write('reg_{} = {{\'type\': \'region\', \'value\': {}}}\n'.format(label, comps))
 
 f.close()
diff --git a/doc/scripts/inputs.py b/doc/scripts/inputs.py
index 63b79dfe0bf6f8c10838fb8567966c0ea35f8018..d3cd78efc1b39735e0578c9b30c5ff253168658f 100644
--- a/doc/scripts/inputs.py
+++ b/doc/scripts/inputs.py
@@ -94,8 +94,8 @@ ls_loc15  = {'type': 'locset', 'value': [(1, 0.5)]}
 ls_uniform0  = {'type': 'locset', 'value': [(0, 0.5841758202819731), (1, 0.6362196581003494), (1, 0.7157318936458146), (1, 0.7464198558822775), (2, 0.6340206909931745), (2, 0.6807941995943213), (3, 0.296698184761692), (3, 0.509669134988683), (3, 0.7662305637426007), (4, 0.5701465671464049)]}
 ls_uniform1  = {'type': 'locset', 'value': [(0, 0.9778060763285382), (1, 0.19973428495790843), (1, 0.8310607916260988), (2, 0.9210229159315735), (2, 0.9244292525837472), (2, 0.9899772550845479), (3, 0.9924233395972087), (4, 0.3641426305909531), (4, 0.4787812247064867), (4, 0.5138656268861914)]}
 ls_branchmid  = {'type': 'locset', 'value': [(0, 0.5), (1, 0.5), (2, 0.5), (3, 0.5), (4, 0.5), (5, 0.5)]}
-ls_distal  = {'type': 'locset', 'value': [(1, 0.7960259763299439), (3, 0.6666666666666667), (4, 0.39052429175126996), (5, 1.0)]}
-ls_proximal  = {'type': 'locset', 'value': [(1, 0.2960259763299439), (2, 0.0), (5, 0.6124999999999999)]}
+ls_distal  = {'type': 'locset', 'value': [(1, 0.796025976329944), (3, 0.6666666666666667), (4, 0.39052429175127), (5, 1.0)]}
+ls_proximal  = {'type': 'locset', 'value': [(1, 0.29602597632994393), (2, 0.0), (5, 0.6124999999999999)]}
 ls_distint_in  = {'type': 'locset', 'value': [(1, 0.5), (2, 0.7), (5, 0.1)]}
 ls_proxint_in  = {'type': 'locset', 'value': [(1, 0.8), (2, 0.3)]}
 ls_loctest  = {'type': 'locset', 'value': [(1, 1.0), (2, 0.0), (5, 0.0)]}
@@ -112,11 +112,11 @@ reg_tag4 = {'type': 'region', 'value': []}
 reg_soma = {'type': 'region', 'value': [(0, 0.0, 0.3324708796524168)]}
 reg_axon = {'type': 'region', 'value': [(5, 0.0, 1.0)]}
 reg_dend = {'type': 'region', 'value': [(0, 0.3324708796524168, 1.0), (1, 0.0, 1.0), (2, 0.0, 1.0), (3, 0.0, 1.0), (4, 0.0, 1.0)]}
-reg_radlt5 = {'type': 'region', 'value': [(1, 0.4440389644949158, 1.0), (3, 0.0, 1.0), (4, 0.0, 1.0), (5, 0.6562500000000001, 1.0)]}
-reg_radle5 = {'type': 'region', 'value': [(1, 0.4440389644949158, 1.0), (2, 0.0, 1.0), (3, 0.0, 1.0), (4, 0.0, 1.0), (5, 0.6562500000000001, 1.0)]}
-reg_radgt5 = {'type': 'region', 'value': [(0, 0.0, 1.0), (1, 0.0, 0.4440389644949158), (5, 0.0, 0.6562500000000001)]}
-reg_radge5 = {'type': 'region', 'value': [(0, 0.0, 1.0), (1, 0.0, 0.4440389644949158), (2, 0.0, 1.0), (3, 0.0, 0.0), (4, 0.0, 0.0), (5, 0.0, 0.6562500000000001)]}
-reg_rad36 = {'type': 'region', 'value': [(1, 0.2960259763299439, 0.7960259763299439), (2, 0.0, 1.0), (3, 0.0, 0.6666666666666667), (4, 0.0, 0.39052429175126996), (5, 0.6124999999999999, 1.0)]}
+reg_radlt5 = {'type': 'region', 'value': [(1, 0.44403896449491587, 1.0), (3, 0.0, 1.0), (4, 0.0, 1.0), (5, 0.65625, 1.0)]}
+reg_radle5 = {'type': 'region', 'value': [(1, 0.44403896449491587, 1.0), (2, 0.0, 1.0), (3, 0.0, 1.0), (4, 0.0, 1.0), (5, 0.65625, 1.0)]}
+reg_radgt5 = {'type': 'region', 'value': [(0, 0.0, 1.0), (1, 0.0, 0.44403896449491587), (5, 0.0, 0.65625)]}
+reg_radge5 = {'type': 'region', 'value': [(0, 0.0, 1.0), (1, 0.0, 0.44403896449491587), (2, 0.0, 1.0), (3, 0.0, 0.0), (4, 0.0, 0.0), (5, 0.0, 0.65625)]}
+reg_rad36 = {'type': 'region', 'value': [(1, 0.29602597632994393, 0.796025976329944), (2, 0.0, 1.0), (3, 0.0, 0.6666666666666667), (4, 0.0, 0.39052429175127), (5, 0.6124999999999999, 1.0)]}
 reg_branch0 = {'type': 'region', 'value': [(0, 0.0, 1.0)]}
 reg_branch3 = {'type': 'region', 'value': [(3, 0.0, 1.0)]}
 reg_cable_1_01 = {'type': 'region', 'value': [(1, 0.0, 1.0)]}
diff --git a/example/dryrun/branch_cell.hpp b/example/dryrun/branch_cell.hpp
index 51d674f243a7dfd1160d11eec9f2b7682b6c9827..837a6569e83eaa2d5af65d18b901816d3c4e69eb 100644
--- a/example/dryrun/branch_cell.hpp
+++ b/example/dryrun/branch_cell.hpp
@@ -8,10 +8,13 @@
 #include <arbor/common_types.hpp>
 #include <arbor/cable_cell.hpp>
 #include <arbor/morph/segment_tree.hpp>
+#include <arbor/string_literals.hpp>
 
 #include <string>
 #include <sup/json_params.hpp>
 
+using namespace arb::literals;
+
 // Parameters used to generate the random cell morphologies.
 struct cell_parameters {
     cell_parameters() = default;
@@ -107,8 +110,8 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param
 
     arb::cable_cell cell(arb::morphology(tree), d);
 
-    cell.paint("soma", "hh");
-    cell.paint("dend", "pas");
+    cell.paint("soma"_lab, "hh");
+    cell.paint("dend"_lab, "pas");
     cell.default_parameters.axial_resistivity = 100; // [Ω·cm]
 
     // Add spike threshold detector at the soma.
diff --git a/example/gap_junctions/gap_junctions.cpp b/example/gap_junctions/gap_junctions.cpp
index f0c3e4b38a8d73c94a299e09730cf82e0572c52f..29dcf8731a6527cb108eb23f14e49a0a79a1732d 100644
--- a/example/gap_junctions/gap_junctions.cpp
+++ b/example/gap_junctions/gap_junctions.cpp
@@ -286,12 +286,8 @@ arb::cable_cell gj_cell(cell_gid_type gid, unsigned ncell, double stim_duration)
     double dend_rad = 3./2; // μm
     tree.append(0, {0,0,2*soma_rad, dend_rad}, {0,0,2*soma_rad+300, dend_rad}, 3);  // dendrite
 
-    // Create a label dictionary that creates a single region that covers the whole cell.
-    arb::label_dict d;
-    d.set("all",  arb::reg::all());
-
     // Create the cell and set its electrical properties.
-    arb::cable_cell cell(tree, d);
+    arb::cable_cell cell(tree);
     cell.default_parameters.axial_resistivity = 100;       // [Ω·cm]
     cell.default_parameters.membrane_capacitance = 0.018;  // [F/m²]
 
@@ -311,17 +307,17 @@ arb::cable_cell gj_cell(cell_gid_type gid, unsigned ncell, double stim_duration)
     pas["e"] =  -65;
 
     // Paint density channels on all parts of the cell
-    cell.paint("all", nax);
-    cell.paint("all", kdrmt);
-    cell.paint("all", kamt);
-    cell.paint("all", pas);
+    cell.paint("(all)", nax);
+    cell.paint("(all)", kdrmt);
+    cell.paint("(all)", kamt);
+    cell.paint("(all)", pas);
 
     // Add a spike detector to the soma.
     cell.place(arb::mlocation{0,0}, arb::threshold_detector{10});
 
     // Add two gap junction sites.
     cell.place(arb::mlocation{0, 1}, arb::gap_junction_site{});
-    cell.place(arb::mlocation{1, 1}, arb::gap_junction_site{});
+    cell.place(arb::mlocation{0, 1}, arb::gap_junction_site{});
 
     // Attach a stimulus to the second cell.
     if (!gid) {
@@ -330,7 +326,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.
-    cell.place(arb::mlocation{1, 0.5}, "expsyn");
+    cell.place(arb::mlocation{0, 0.5}, "expsyn");
 
     return cell;
 }
diff --git a/example/generators/generators.cpp b/example/generators/generators.cpp
index 48147fc064794ba181877f1956e437d99cf04711..734486369eb8fbc3f77edcda2b1b8c41287b3beb 100644
--- a/example/generators/generators.cpp
+++ b/example/generators/generators.cpp
@@ -56,7 +56,7 @@ public:
         d.set("soma", arb::reg::tagged(1));
 
         arb::cable_cell c(tree, d);
-        c.paint("soma", "pas");
+        c.paint("\"soma\"", "pas");
 
         // Add one synapse at the soma.
         // This synapse will be the target for all events, from both
diff --git a/example/ring/branch_cell.hpp b/example/ring/branch_cell.hpp
index 51d674f243a7dfd1160d11eec9f2b7682b6c9827..837a6569e83eaa2d5af65d18b901816d3c4e69eb 100644
--- a/example/ring/branch_cell.hpp
+++ b/example/ring/branch_cell.hpp
@@ -8,10 +8,13 @@
 #include <arbor/common_types.hpp>
 #include <arbor/cable_cell.hpp>
 #include <arbor/morph/segment_tree.hpp>
+#include <arbor/string_literals.hpp>
 
 #include <string>
 #include <sup/json_params.hpp>
 
+using namespace arb::literals;
+
 // Parameters used to generate the random cell morphologies.
 struct cell_parameters {
     cell_parameters() = default;
@@ -107,8 +110,8 @@ arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& param
 
     arb::cable_cell cell(arb::morphology(tree), d);
 
-    cell.paint("soma", "hh");
-    cell.paint("dend", "pas");
+    cell.paint("soma"_lab, "hh");
+    cell.paint("dend"_lab, "pas");
     cell.default_parameters.axial_resistivity = 100; // [Ω·cm]
 
     // Add spike threshold detector at the soma.
diff --git a/example/single/single.cpp b/example/single/single.cpp
index bbcb65569d364a03c0253df0e76d1a1c06a4bf61..e36740d580d50667ff5429e4afa98320af31a7ad 100644
--- a/example/single/single.cpp
+++ b/example/single/single.cpp
@@ -63,8 +63,8 @@ struct single_recipe: public arb::recipe {
         arb::cable_cell c(morpho, dict);
 
         // Add HH mechanism to soma, passive channels to dendrites.
-        c.paint("soma", "hh");
-        c.paint("dend", "pas");
+        c.paint("\"soma\"", "hh");
+        c.paint("\"dend\"", "pas");
 
         // Add synapse to last branch.
 
diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
index 64f7cb6d9979945220c84bc465e5ebc685e6bd57..5344e1f771575c7cbbd8090c38bd9d60c4a1786e 100644
--- a/python/CMakeLists.txt
+++ b/python/CMakeLists.txt
@@ -25,7 +25,6 @@ set(pyarb_source
     identifiers.cpp
     mechanism.cpp
     morphology.cpp
-    morph_parse.cpp
     mpi.cpp
     profiler.cpp
     pyarb.cpp
@@ -35,7 +34,6 @@ set(pyarb_source
     simulation.cpp
     single_cell_model.cpp
     spikes.cpp
-    s_expr.cpp
 )
 
 # compile the pyarb sources into an object library that will be
diff --git a/python/cells.cpp b/python/cells.cpp
index 6f1b67f1f227376ff2e1dbc89738bd640fe4f96e..54510c9070d5d3a359967fa45492ae3733f880a1 100644
--- a/python/cells.cpp
+++ b/python/cells.cpp
@@ -14,6 +14,7 @@
 #include <arbor/cable_cell.hpp>
 #include <arbor/lif_cell.hpp>
 #include <arbor/morph/label_dict.hpp>
+#include <arbor/morph/label_parse.hpp>
 #include <arbor/morph/locset.hpp>
 #include <arbor/morph/region.hpp>
 #include <arbor/morph/segment_tree.hpp>
@@ -26,12 +27,9 @@
 #include "cells.hpp"
 #include "conversion.hpp"
 #include "error.hpp"
-#include "morph_parse.hpp"
 #include "schedule.hpp"
 #include "strprintf.hpp"
 
-using arb::util::any_cast;
-
 namespace pyarb {
 
 // Convert a cell description inside a Python object to a cell
@@ -96,32 +94,23 @@ struct label_dict_proxy {
         // Then it parses the description, and tests whether the description
         // is a region or locset, and updates the label dictionary appropriately.
         // Errors occur when:
-        //  * the name is not a valid name.
         //  * a region is described with a name that matches an existing locset
         //    (and vice versa.)
         //  * the description is not well formed, e.g. it contains a syntax error.
         //  * the description is well-formed, but describes neither a region or locset.
         try{
-            // Test that the identifier is valid, i.e.
-            //  * only numbers, letters and underscore.
-            //  * no leading number or underscore.
-            if (!test_identifier(name)) {
-                throw std::string(util::pprintf("'{}' is not a valid label name.", name));
-            }
-            // Parse the input string into an s-expression.
-            auto parsed = parse(desc);
             // Evaluate the s-expression to build a region/locset.
-            auto result = eval(parsed);
+            auto result = arb::parse_label_expression(desc);
             if (!result) { // an error parsing / evaluating description.
-                throw std::string(result.error().message);
+                throw result.error();
             }
             else if (result->type()==typeid(arb::region)) { // describes a region.
-                dict.set(name, std::move(any_cast<arb::region&>(*result)));
+                dict.set(name, std::move(std::any_cast<arb::region&>(*result)));
                 auto it = std::lower_bound(regions.begin(), regions.end(), name);
                 if (it==regions.end() || *it!=name) regions.insert(it, name);
             }
             else if (result->type()==typeid(arb::locset)) { // describes a locset.
-                dict.set(name, std::move(any_cast<arb::locset&>(*result)));
+                dict.set(name, std::move(std::any_cast<arb::locset&>(*result)));
                 auto it = std::lower_bound(locsets.begin(), locsets.end(), name);
                 if (it==locsets.end() || *it!=name) locsets.insert(it, name);
             }
diff --git a/python/error.hpp b/python/error.hpp
index 8a650b8aed41ed2c66a0877f96ed7e8c4c00456d..91b0ddee533ee8f0e9f61e9fa7344c2bb3524967 100644
--- a/python/error.hpp
+++ b/python/error.hpp
@@ -6,11 +6,6 @@
 
 #include <pybind11/pybind11.h>
 
-#include <arbor/arbexcept.hpp>
-#include <arbor/util/expected.hpp>
-
-#include "strprintf.hpp"
-
 namespace pyarb {
 
 extern std::exception_ptr py_exception;
@@ -53,66 +48,4 @@ auto try_catch_pyexception(L func, const char* msg){
     }
 }
 
-template <typename T, typename E>
-struct hopefully {
-    using value_type = T;
-    using error_type = E;
-    arb::util::expected<value_type, error_type> state;
-
-    hopefully(const hopefully&) = default;
-
-    hopefully(value_type x): state(std::move(x)) {}
-    hopefully(error_type x): state(arb::util::unexpect, std::move(x)) {}
-
-    const value_type& operator*() const {
-        return try_get();
-    }
-    value_type& operator*() {
-        return try_get();
-    }
-    const value_type* operator->() const {
-        return &try_get();
-    }
-    value_type* operator->() {
-        return &try_get();
-    }
-
-    operator bool() const {
-        return (bool)state;
-    }
-
-    const error_type& error() const {
-        try {
-            return state.error();
-        }
-        catch(arb::util::bad_expected_access<void>& e) {
-            throw arb::arbor_internal_error(
-                "Attempt to get an error from a valid hopefully wrapper.");
-        }
-    }
-
-private:
-
-    const value_type& try_get() const {
-        try {
-            return state.value();
-        }
-        catch(arb::util::bad_expected_access<void>& e) {
-            throw arbor_internal_error(util::pprintf(
-                "Attempt to unwrap a hopefully with error state '{}'",
-                error().message));
-        }
-    }
-    value_type& try_get() {
-        try {
-            return state.value();
-        }
-        catch(arb::util::bad_expected_access<void>& e) {
-            throw arb::arbor_internal_error(util::pprintf(
-                "Attempt to unwrap a hopefully with error state '{}'",
-                error().message));
-        }
-    }
-};
-
 } // namespace pyarb
diff --git a/python/example/single_cell_builder.py b/python/example/single_cell_builder.py
index ad38d656215b51df2a1738efef7548aba265c22d..58e7054cacfa12343c85f0f390a00400cedb36bd 100644
--- a/python/example/single_cell_builder.py
+++ b/python/example/single_cell_builder.py
@@ -53,16 +53,16 @@ cell = b.build()
 # Set initial membrane potential everywhere on the cell to -40 mV.
 cell.set_properties(Vm=-40)
 # Put hh dynamics on soma, and passive properties on the dendrites.
-cell.paint('soma', 'hh')
-cell.paint('dend', 'pas')
+cell.paint('"soma"', 'hh')
+cell.paint('"dend"', 'pas')
 # Set axial resistivity in dendrite regions (Ohm.cm)
-cell.paint('dendn', rL=500)
-cell.paint('dendx', rL=10000)
+cell.paint('"dendn"', rL=500)
+cell.paint('"dendx"', rL=10000)
 # Attach stimuli with duration of 2 ms and current of 0.8 nA.
 # There are three stimuli, which activate at 10 ms, 50 ms and 80 ms.
-cell.place('stim_site', arbor.iclamp( 10, 2, 0.8))
-cell.place('stim_site', arbor.iclamp( 50, 2, 0.8))
-cell.place('stim_site', arbor.iclamp( 80, 2, 0.8))
+cell.place('"stim_site"', arbor.iclamp( 10, 2, 0.8))
+cell.place('"stim_site"', arbor.iclamp( 50, 2, 0.8))
+cell.place('"stim_site"', arbor.iclamp( 80, 2, 0.8))
 # Add a spike detector with threshold of -10 mV.
 cell.place('root', arbor.spike_detector(-10))
 
@@ -71,7 +71,7 @@ m = arbor.single_cell_model(cell)
 
 # Attach voltage probes, sampling at 10 kHz.
 m.probe('voltage', loc(0,0), 10000) # at the soma.
-m.probe('voltage', 'dtips',   10000) # at the tips of the dendrites.
+m.probe('voltage', '"dtips"',   10000) # at the tips of the dendrites.
 
 # Run simulation for 100 ms of simulated activity.
 tfinal=100
diff --git a/python/example/single_cell_swc.py b/python/example/single_cell_swc.py
index e0ac637756fb2184b4d6fc00ae9aeebdb1e9eb1b..957dc8bfa978cd745d4047cdf687102d3ecc09c3 100644
--- a/python/example/single_cell_swc.py
+++ b/python/example/single_cell_swc.py
@@ -37,17 +37,17 @@ cell.set_properties(Vm=-55)
 # Use Nernst to calculate reversal potential for calcium.
 cell.set_ion('ca', method=mech('nernst/x=ca'))
 # hh mechanism on the soma and axon.
-cell.paint('soma', 'hh')
-cell.paint('axon', 'hh')
+cell.paint('"soma"', 'hh')
+cell.paint('"axon"', 'hh')
 # pas mechanism the dendrites.
-cell.paint('dend', 'pas')
+cell.paint('"dend"', 'pas')
 # Increase resistivity on dendrites.
-cell.paint('dend', rL=500)
+cell.paint('"dend"', rL=500)
 # Attach stimuli that inject 0.8 nA currents for 1 ms, starting at 3 and 8 ms.
-cell.place('stim_site', arbor.iclamp(3, 1, current=0.5))
-cell.place('stim_site', arbor.iclamp(8, 1, current=1))
+cell.place('"stim_site"', arbor.iclamp(3, 1, current=0.5))
+cell.place('"stim_site"', arbor.iclamp(8, 1, current=1))
 # Detect spikes at the soma with a voltage threshold of -10 mV.
-cell.place('root', arbor.spike_detector(-10))
+cell.place('"root"', arbor.spike_detector(-10))
 
 # Have one compartment between each sample point.
 cell.compartments_on_segments()
diff --git a/python/flat_cell_builder.cpp b/python/flat_cell_builder.cpp
index 2a2a3ca654b06ea4a5cbab72711910a7c953f30d..512a51828789bffb00c061ba64256f6d38a60116 100644
--- a/python/flat_cell_builder.cpp
+++ b/python/flat_cell_builder.cpp
@@ -1,3 +1,4 @@
+#include <any>
 #include <mutex>
 #include <string>
 #include <unordered_map>
@@ -8,6 +9,7 @@
 #include <pybind11/stl.h>
 
 #include <arbor/cable_cell.hpp>
+#include <arbor/morph/label_parse.hpp>
 #include <arbor/morph/morphology.hpp>
 #include <arbor/morph/primitives.hpp>
 #include <arbor/morph/segment_tree.hpp>
@@ -15,12 +17,8 @@
 
 #include "conversion.hpp"
 #include "error.hpp"
-#include "morph_parse.hpp"
-#include "s_expr.hpp"
 #include "strprintf.hpp"
 
-using arb::util::any_cast;
-
 namespace pyarb {
 
 class flat_cell_builder {
@@ -65,10 +63,6 @@ public:
 
         cached_morpho_ = false;
 
-        if (!test_identifier(region)) {
-            throw pyarb_error(util::pprintf("'{}' is not a valid label name.", region));
-        }
-
         // Get tag id of region (add a new tag if region does not already exist).
         int tag = get_tag(region);
         const bool at_root = parent==mnpos;
@@ -95,17 +89,13 @@ public:
     }
 
     void add_label(const char* name, const char* description) {
-        if (!test_identifier(name)) {
-            throw pyarb_error(util::pprintf("'{}' is not a valid label name.", name));
-        }
-
-        if (auto result = eval(parse(description)) ) {
+        if (auto result = arb::parse_label_expression(description) ) {
             // The description is a region.
             if (result->type()==typeid(arb::region)) {
                 if (dict_.locset(name)) {
                     throw pyarb_error("Region name clashes with a locset.");
                 }
-                auto& reg = any_cast<arb::region&>(*result);
+                auto& reg = std::any_cast<arb::region&>(*result);
                 if (auto r = dict_.region(name)) {
                     dict_.set(name, join(std::move(reg), std::move(*r)));
                 }
@@ -113,11 +103,12 @@ public:
                     dict_.set(name, std::move(reg));
                 }
             }
+            // The description is a locset.
             else if (result->type()==typeid(arb::locset)) {
                 if (dict_.region(name)) {
                     throw pyarb_error("Locset name clashes with a region.");
                 }
-                auto& loc = any_cast<arb::locset&>(*result);
+                auto& loc = std::any_cast<arb::locset&>(*result);
                 if (auto l = dict_.locset(name)) {
                     dict_.set(name, sum(std::move(loc), std::move(*l)));
                 }
@@ -125,12 +116,13 @@ public:
                     dict_.set(name, std::move(loc));
                 }
             }
+            // Error: the description is neither.
             else {
                 throw pyarb_error("Label describes neither a region nor a locset.");
             }
         }
         else {
-            throw pyarb_error(result.error().message);
+            throw pyarb_error(result.error().what());
         }
     }
 
diff --git a/python/morph_parse.hpp b/python/morph_parse.hpp
deleted file mode 100644
index 84dd1ba30645ddb531b98a093b69cb496e094ee0..0000000000000000000000000000000000000000
--- a/python/morph_parse.hpp
+++ /dev/null
@@ -1,20 +0,0 @@
-#pragma once
-
-#include <any>
-
-#include "error.hpp"
-#include "s_expr.hpp"
-
-namespace pyarb {
-
-struct parse_error_state {
-    std::string message;
-    int location;
-};
-
-template <typename T>
-using parse_hopefully = hopefully<T, parse_error_state>;
-
-parse_hopefully<std::any> eval(const s_expr&);
-
-} // namespace pyarb
diff --git a/python/s_expr.hpp b/python/s_expr.hpp
deleted file mode 100644
index 8b041d56a3d580c1952d90218ed81364c225de1b..0000000000000000000000000000000000000000
--- a/python/s_expr.hpp
+++ /dev/null
@@ -1,123 +0,0 @@
-#pragma once
-
-#include <string>
-#include <memory>
-#include <variant>
-#include <vector>
-
-namespace pyarb {
-
-enum class tok {
-    nil,
-    real,       // real number
-    integer,    // integer
-    name,       // name
-    lparen,     // left parenthesis '('
-    rparen,     // right parenthesis ')'
-    string,     // string, written as "spelling"
-    eof,        // end of file/input
-    error       // special error state marker
-};
-
-std::ostream& operator<<(std::ostream&, const tok&);
-
-struct token {
-    int column;
-    tok kind;
-    std::string spelling;
-};
-
-std::ostream& operator<<(std::ostream&, const token&);
-
-std::vector<token> tokenize(const char* line);
-
-struct s_expr {
-    template <typename U>
-    struct s_pair {
-        U head = U();
-        U tail = U();
-        s_pair(U l, U r): head(std::move(l)), tail(std::move(r)) {}
-    };
-
-    // This value_wrapper is used to wrap the shared pointer
-    template <typename T>
-    struct value_wrapper{
-        using state_t = std::unique_ptr<T>;
-        state_t state;
-
-        value_wrapper() = default;
-
-        value_wrapper(const T& v):
-            state(std::make_unique<T>(v)) {}
-
-        value_wrapper(T&& v):
-            state(std::make_unique<T>(std::move(v))) {}
-
-        value_wrapper(const value_wrapper& other):
-            state(std::make_unique<T>(other.get())) {}
-
-        value_wrapper& operator=(const value_wrapper& other) {
-            state = std::make_unique<T>(other.get());
-            return *this;
-        }
-
-        value_wrapper(value_wrapper&& other) = default;
-
-        friend std::ostream& operator<<(std::ostream& o, const value_wrapper& w) {
-            return o << *w.state;
-        }
-
-        operator T() const {
-            return *state;
-        }
-
-        const T& get() const {
-            return *state;
-        }
-
-        T& get() {
-            return *state;
-        }
-    };
-
-    // An s_expr can be one of
-    //      1. an atom
-    //      2. a pair of s_expr (head and tail)
-    // The s_expr uses std::variant to represent these two possible states,
-    // which requires using an incomplete definition of s_expr, which is stored
-    // using a std::unique_ptr.
-
-    using pair_type = s_pair<value_wrapper<s_expr>>;
-    std::variant<token, pair_type> state = token{-1, tok::nil, "nil"};
-
-    s_expr(const s_expr& s): state(s.state) {}
-    s_expr() = default;
-    s_expr(token t): state(std::move(t)) {}
-    s_expr(s_expr l, s_expr r):
-        state(pair_type(std::move(l), std::move(r)))
-    {}
-
-    bool is_atom() const;
-
-    const token& atom() const;
-
-    operator bool() const;
-
-    const s_expr& head() const;
-    const s_expr& tail() const;
-    s_expr& head();
-    s_expr& tail();
-
-    friend std::ostream& operator<<(std::ostream& o, const s_expr& x);
-};
-
-std::size_t length(const s_expr& l);
-int location(const s_expr& l);
-
-s_expr parse(const char* line);
-s_expr parse(const std::string& line);
-
-bool test_identifier(const char* in);
-
-} // namespace pyarb
-
diff --git a/python/single_cell_model.cpp b/python/single_cell_model.cpp
index 9c12ed4ce7cc7bf74913fbe09792b95fd41c4a23..43029969127a2e45666d4911ffc45a78b80176d5 100644
--- a/python/single_cell_model.cpp
+++ b/python/single_cell_model.cpp
@@ -13,8 +13,9 @@
 #include <arbor/simulation.hpp>
 #include <arbor/util/any_cast.hpp>
 
-#include "error.hpp"
 #include "cells.hpp"
+#include "error.hpp"
+#include "strprintf.hpp"
 
 using arb::util::any_cast;
 
diff --git a/python/test/cpp/CMakeLists.txt b/python/test/cpp/CMakeLists.txt
index 291da50100446b139640fa0a2d6fe310878bba8e..98c72c7532ab342b3b603d89e6d615f9dfb7d72f 100644
--- a/python/test/cpp/CMakeLists.txt
+++ b/python/test/cpp/CMakeLists.txt
@@ -1,5 +1,6 @@
 set(py_unit_sources
-    s_expr.cpp
+    # Currently empty: the s expression tests were moved, however keep this
+    # open for the inevitable tests that will be written in the future.
 
     # unit test driver
     test.cpp
diff --git a/python/test/cpp/s_expr.cpp b/python/test/cpp/s_expr.cpp
deleted file mode 100644
index ced25e9ef2d81c29fd39ad57d4ae4141dd6ecc70..0000000000000000000000000000000000000000
--- a/python/test/cpp/s_expr.cpp
+++ /dev/null
@@ -1,117 +0,0 @@
-#include "gtest.h"
-
-#include <string>
-
-#include <arbor/morph/region.hpp>
-#include <arbor/morph/locset.hpp>
-#include <arbor/util/any_cast.hpp>
-
-#include "morph_parse.hpp"
-#include "s_expr.hpp"
-#include "strprintf.hpp"
-
-using namespace pyarb;
-using namespace std::string_literals;
-
-using arb::util::any_cast;
-
-TEST(s_expr, identifier) {
-    EXPECT_TRUE(test_identifier("foo"));
-    EXPECT_TRUE(test_identifier("f1"));
-    EXPECT_TRUE(test_identifier("f_"));
-    EXPECT_TRUE(test_identifier("f_1__"));
-    EXPECT_TRUE(test_identifier("A_1__"));
-
-    EXPECT_TRUE(test_identifier("A-1"));
-    EXPECT_TRUE(test_identifier("hello-world"));
-    EXPECT_TRUE(test_identifier("hello--world"));
-    EXPECT_TRUE(test_identifier("hello--world_"));
-
-    EXPECT_FALSE(test_identifier("_foobar"));
-    EXPECT_FALSE(test_identifier("-foobar"));
-    EXPECT_FALSE(test_identifier("2dogs"));
-    EXPECT_FALSE(test_identifier("1"));
-    EXPECT_FALSE(test_identifier("_"));
-    EXPECT_FALSE(test_identifier("-"));
-    EXPECT_FALSE(test_identifier(""));
-    EXPECT_FALSE(test_identifier(" foo"));
-    EXPECT_FALSE(test_identifier("foo "));
-    EXPECT_FALSE(test_identifier("foo bar"));
-    EXPECT_FALSE(test_identifier(""));
-}
-
-TEST(s_expr, atoms) {
-    auto get_atom = [](s_expr e) {
-        EXPECT_EQ(1u, length(e));
-        EXPECT_TRUE(e.head().is_atom());
-        return e.head().atom();
-    };
-
-    for (auto spelling: {"foo", "bar_", "car1", "car_1", "x_1__"}){
-        auto a = get_atom(parse("("s+spelling+")"));
-        EXPECT_EQ(a.kind, tok::name);
-        EXPECT_EQ(a.spelling, spelling);
-    }
-    // test parsing of integers
-    for (auto spelling: {"0", "-0", "1", "42", "-3287"}){
-        auto a = get_atom(parse("("s+spelling+")"));
-        EXPECT_EQ(a.kind, tok::integer);
-        EXPECT_EQ(a.spelling, spelling);
-    }
-    // test parsing of real numbers
-    for (auto spelling: {"0.", "-0.0", "1.21", "42.", "-3287.12"}){
-        auto a = get_atom(parse("("s+spelling+")"));
-        EXPECT_EQ(a.kind, tok::real);
-        EXPECT_EQ(a.spelling, spelling);
-    }
-    // test parsing of strings
-    for (auto spelling: {"foo", "dog cat", ""}) {
-        auto s = "(\""s+spelling+"\")";
-        auto a = get_atom(parse(s));
-        EXPECT_EQ(a.kind, tok::string);
-    }
-}
-
-TEST(s_expr, parse) {
-    auto round_trip_reg = [](const char* in) {
-        auto x = eval(parse(in));
-        return util::pprintf("{}", any_cast<arb::region>(*x));
-    };
-    auto round_trip_loc = [](const char* in) {
-        auto x = eval(parse(in));
-        return util::pprintf("{}", any_cast<arb::locset>(*x));
-    };
-
-    EXPECT_EQ("(cable 3 0 1)",      round_trip_reg("(branch 3)"));
-    EXPECT_EQ("(cable 2 0.1 0.4)",  round_trip_reg("(cable 2 0.1 0.4)"));
-    EXPECT_EQ("(all)",              round_trip_reg("(all)"));
-    EXPECT_EQ("(region \"foo\")",   round_trip_reg("(region \"foo\")"));
-
-    EXPECT_EQ("(terminal)", round_trip_loc("(terminal)"));
-    EXPECT_EQ("(root)",     round_trip_loc("(root)"));
-    EXPECT_EQ("(locset \"cat_burgler\")", round_trip_loc("(locset \"cat_burgler\")"));
-
-    auto lhs = any_cast<arb::region>(*eval(parse("(region \"dend\")")));
-    auto rhs = any_cast<arb::region>(*eval(parse("(all)")));
-
-    EXPECT_EQ(util::pprintf("{}", join(lhs,rhs)), "(join (region \"dend\") (all))");
-}
-
-TEST(s_expr, comments) {
-    auto round_trip_reg = [](const char* in) {
-        auto x = eval(parse(in));
-        return util::pprintf("{}", any_cast<arb::region>(*x));
-    };
-
-    EXPECT_EQ("(all)",  round_trip_reg("(all) ; a comment"));
-    const char *multi_line = 
-        "; comment at start\n"
-        "(radius_lt\n"
-        "    (join\n"
-        "        (tag 3) ; end of line\n"
-        " ;comment on whole line\n"
-        "        (tag 4))\n"
-        "    0.5) ; end of string";
-    EXPECT_EQ("(radius_lt (join (tag 3) (tag 4)) 0.5)",
-              round_trip_reg(multi_line));
-}
diff --git a/python/tutorial/example1.py b/python/tutorial/example1.py
index 1512803e95a029a99fa06670a6837a188eb86f21..e698363d82196e9e38ccb72bf74a012760888c3d 100644
--- a/python/tutorial/example1.py
+++ b/python/tutorial/example1.py
@@ -11,17 +11,17 @@ cell = arbor.cable_cell(tree, labels)
 # Set initial membrane potential everywhere on the cell to -40 mV.
 cell.set_properties(Vm=-40)
 # Put hh dynamics on soma, and passive properties on the dendrites.
-cell.paint('soma', 'hh')
+cell.paint('"soma"', 'hh')
 # Attach stimuli with duration of 100 ns and current of 0.8 nA.
-cell.place('center', arbor.iclamp( 10, duration=0.1, current=0.8))
+cell.place('"center"', arbor.iclamp( 10, duration=0.1, current=0.8))
 # Add a spike detector with threshold of -10 mV.
-cell.place('center', arbor.spike_detector(-10))
+cell.place('"center"', arbor.spike_detector(-10))
 
 # Make single cell model.
 m = arbor.single_cell_model(cell)
 
 # Attach voltage probes, sampling at 1 MHz.
-m.probe('voltage', 'center',  1000000)
+m.probe('voltage', '"center"',  1000000)
 
 # Run simulation for tfinal ms with time steps of 1 μs.
 tfinal=30
diff --git a/test/common_cells.cpp b/test/common_cells.cpp
index 43cbcc92da557c80afc1a31ec0705fff6c744e02..2b485360b8ba4e5a125449e3397d4ffa41fff561 100644
--- a/test/common_cells.cpp
+++ b/test/common_cells.cpp
@@ -1,3 +1,4 @@
+#include <arbor/string_literals.hpp>
 #include "common_cells.hpp"
 
 namespace arb {
@@ -173,10 +174,11 @@ cable_cell soma_cell_builder::make_cell() const {
  */
 
 cable_cell make_cell_soma_only(bool with_stim) {
+    using namespace arb::literals;
     soma_cell_builder builder(18.8/2.0);
 
     auto c = builder.make_cell();
-    c.paint("soma", "hh");
+    c.paint("soma"_lab, "hh");
     if (with_stim) {
         c.place(builder.location({0,0.5}), i_clamp{10., 100., 0.1});
     }
@@ -206,12 +208,13 @@ cable_cell make_cell_soma_only(bool with_stim) {
  */
 
 cable_cell make_cell_ball_and_stick(bool with_stim) {
+    using namespace arb::literals;
     soma_cell_builder builder(12.6157/2.0);
     builder.add_branch(0, 200, 1.0/2, 1.0/2, 4, "dend");
 
     auto c = builder.make_cell();
-    c.paint("soma", "hh");
-    c.paint("dend", "pas");
+    c.paint("soma"_lab, "hh");
+    c.paint("dend"_lab, "pas");
     if (with_stim) {
         c.place(builder.location({1,1}), i_clamp{5, 80, 0.3});
     }
@@ -244,14 +247,15 @@ cable_cell make_cell_ball_and_stick(bool with_stim) {
  */
 
 cable_cell make_cell_ball_and_3stick(bool with_stim) {
+    using namespace arb::literals;
     soma_cell_builder builder(12.6157/2.0);
     builder.add_branch(0, 100, 0.5, 0.5, 4, "dend");
     builder.add_branch(1, 100, 0.5, 0.5, 4, "dend");
     builder.add_branch(1, 100, 0.5, 0.5, 4, "dend");
 
     auto c = builder.make_cell();
-    c.paint("soma", "hh");
-    c.paint("dend", "pas");
+    c.paint("soma"_lab, "hh");
+    c.paint("dend"_lab, "pas");
     if (with_stim) {
         c.place(builder.location({2,1}), i_clamp{5.,  80., 0.45});
         c.place(builder.location({3,1}), i_clamp{40., 10.,-0.2});
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index cd9ee32d9f9d2c1080016ff0b3018573261f30e0..d51ba31461442ca41e7e84938056faaef6a713d0 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -153,6 +153,7 @@ set(unit_sources
     test_strprintf.cpp
     test_swcio.cpp
     test_synapses.cpp
+    test_s_expr.cpp
     test_thread.cpp
     test_threading_exceptions.cpp
     test_tree.cpp
diff --git a/test/unit/test_fvm_layout.cpp b/test/unit/test_fvm_layout.cpp
index 8403a30d4f829926c5d04812db44a930240bace4..3613a251cc6cbb997cccd579baf00f40e6e16051 100644
--- a/test/unit/test_fvm_layout.cpp
+++ b/test/unit/test_fvm_layout.cpp
@@ -44,8 +44,8 @@ namespace {
 
             cells.push_back(builder.make_cell());
             auto& cell = cells.back();
-            cell.paint("soma", "hh");
-            cell.paint("dend", "pas");
+            cell.paint("\"soma\"", "hh");
+            cell.paint("\"dend\"", "pas");
             cell.place(builder.location({1,1}), i_clamp{5, 80, 0.3});
 
             s.builders.push_back(std::move(builder));
@@ -87,8 +87,8 @@ namespace {
             cells.push_back(b.make_cell());
             auto& cell = cells.back();
 
-            cell.paint("soma", "hh");
-            cell.paint("dend", "pas");
+            cell.paint("\"soma\"", "hh");
+            cell.paint("\"dend\"", "pas");
 
             using ::arb::reg::branch;
             auto c1 = reg::cable(b1-1, b.location({b1, 0}).pos, 1);
@@ -542,10 +542,10 @@ TEST(fvm_layout, density_norm_area) {
     hh_3["gl"] = seg3_gl;
 
     auto cell = builder.make_cell();
-    cell.paint("soma", std::move(hh_0));
-    cell.paint("reg1", std::move(hh_1));
-    cell.paint("reg2", std::move(hh_2));
-    cell.paint("reg3", std::move(hh_3));
+    cell.paint("\"soma\"", std::move(hh_0));
+    cell.paint("\"reg1\"", std::move(hh_1));
+    cell.paint("\"reg2\"", std::move(hh_2));
+    cell.paint("\"reg3\"", std::move(hh_3));
 
     std::vector<cable_cell> cells{std::move(cell)};
 
@@ -691,7 +691,7 @@ TEST(fvm_layout, density_norm_area_partial) {
 
 TEST(fvm_layout, valence_verify) {
     auto cell = soma_cell_builder(6).make_cell();
-    cell.paint("soma", "test_cl_valence");
+    cell.paint("\"soma\"", "test_cl_valence");
     std::vector<cable_cell> cells{std::move(cell)};
 
     cable_cell_global_properties gprop;
@@ -820,9 +820,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 cell = builder.make_cell();
-    cell.paint("soma", "read_eX/c");
-    cell.paint("soma", "read_eX/a");
-    cell.paint("dend", "read_eX/a");
+    cell.paint("\"soma\"", "read_eX/c");
+    cell.paint("\"soma\"", "read_eX/a");
+    cell.paint("\"dend\"", "read_eX/a");
 
     std::vector<cable_cell> cells{cell, cell};
 
diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp
index 12e2b912f793051bb7d69626a97afe5d879be8c9..4e7027558662f9885db22c585d89b965a1e8ae2b 100644
--- a/test/unit/test_fvm_lowered.cpp
+++ b/test/unit/test_fvm_lowered.cpp
@@ -13,6 +13,7 @@
 #include <arbor/sampling.hpp>
 #include <arbor/simulation.hpp>
 #include <arbor/schedule.hpp>
+#include <arbor/string_literals.hpp>
 #include <arbor/util/any_ptr.hpp>
 
 #include <arborenv/concurrency.hpp>
@@ -542,7 +543,7 @@ TEST(fvm_lowered, read_valence) {
 
         soma_cell_builder builder(6);
         auto cell = builder.make_cell();
-        cell.paint("soma", "test_ca_read_valence");
+        cell.paint("\"soma\"", "test_ca_read_valence");
         cable1d_recipe rec({std::move(cell)});
         rec.catalogue() = make_unit_test_catalogue();
 
@@ -565,7 +566,7 @@ TEST(fvm_lowered, read_valence) {
         // Check ion renaming.
         soma_cell_builder builder(6);
         auto cell = builder.make_cell();
-        cell.paint("soma", "cr_read_valence");
+        cell.paint("\"soma\"", "cr_read_valence");
         cable1d_recipe rec({std::move(cell)});
         rec.catalogue() = make_unit_test_catalogue();
         rec.catalogue() = make_unit_test_catalogue();
@@ -652,6 +653,7 @@ TEST(fvm_lowered, ionic_concentrations) {
 }
 
 TEST(fvm_lowered, ionic_currents) {
+    using namespace arb::literals;
     arb::proc_allocation resources;
     if (auto nt = arbenv::get_env_num_threads()) {
         resources.num_threads = nt;
@@ -681,8 +683,8 @@ TEST(fvm_lowered, ionic_currents) {
     m2["coeff"] = coeff;
 
     auto c = b.make_cell();
-    c.paint("soma", m1);
-    c.paint("soma", m2);
+    c.paint("soma"_lab, m1);
+    c.paint("soma"_lab, m2);
 
     cable1d_recipe rec(std::move(c));
     rec.catalogue() = make_unit_test_catalogue();
diff --git a/test/unit/test_mc_cell_group.cpp b/test/unit/test_mc_cell_group.cpp
index 2971974752102f622f80b9cca10aa7e3f4bc0d89..327852cdc3fd0392accf37293467a68af88ef43c 100644
--- a/test/unit/test_mc_cell_group.cpp
+++ b/test/unit/test_mc_cell_group.cpp
@@ -1,6 +1,7 @@
 #include "../gtest.h"
 
 #include <arbor/common_types.hpp>
+#include <arbor/string_literals.hpp>
 
 #include "epoch.hpp"
 #include "fvm_lowered_cell.hpp"
@@ -12,6 +13,7 @@
 #include "../simple_recipes.hpp"
 
 using namespace arb;
+using namespace arb::literals;
 
 namespace {
     execution_context context;
@@ -24,8 +26,8 @@ namespace {
         soma_cell_builder builder(12.6157/2.0);
         builder.add_branch(0, 200, 0.5, 0.5, 101, "dend");
         cable_cell c = builder.make_cell();
-        c.paint("soma", "hh");
-        c.paint("dend", "pas");
+        c.paint("soma"_lab, "hh");
+        c.paint("dend"_lab, "pas");
         c.place(builder.location({1,1}), i_clamp{5, 80, 0.3});
         c.place(builder.location({0, 0}), threshold_detector{0});
         return c;
diff --git a/test/unit/test_mc_cell_group_gpu.cpp b/test/unit/test_mc_cell_group_gpu.cpp
index 5a6d11cccfa414dab265b74d8e4123db15978daf..bb7f1791e1a38ec4a4590b98467c41dfa1a15479 100644
--- a/test/unit/test_mc_cell_group_gpu.cpp
+++ b/test/unit/test_mc_cell_group_gpu.cpp
@@ -25,8 +25,8 @@ namespace {
         soma_cell_builder builder(12.6157/2.0);
         builder.add_branch(0, 200, 0.5, 0.5, 101, "dend");
         cable_cell c = builder.make_cell();
-        c.paint("soma", "hh");
-        c.paint("dend", "pas");
+        c.paint("\"soma\"", "hh");
+        c.paint("\"dend\"", "pas");
         c.place(builder.location({1, 1}), i_clamp{5, 80, 0.3});
         c.place(builder.location({0, 0}), threshold_detector{0});
         return c;
diff --git a/test/unit/test_morph_expr.cpp b/test/unit/test_morph_expr.cpp
index f3e38cd1125b8d403e4c1db362cc1a3447f42259..2ff804d23543be61cf74057ba3505fa8bc2f9cdd 100644
--- a/test/unit/test_morph_expr.cpp
+++ b/test/unit/test_morph_expr.cpp
@@ -9,6 +9,7 @@
 #include <arbor/morph/mprovider.hpp>
 #include <arbor/morph/primitives.hpp>
 #include <arbor/morph/region.hpp>
+#include <arbor/string_literals.hpp>
 
 #include "util/span.hpp"
 #include "util/strprintf.hpp"
@@ -18,6 +19,7 @@
 #include "morph_pred.hpp"
 
 using namespace arb;
+using namespace arb::literals;
 using embedding = embed_pwlin;
 
 using testing::region_eq;
@@ -104,17 +106,17 @@ TEST(locset, thingify_named) {
         dict.set("cake", cake);
 
         mprovider mp(morphology(sm), dict);
-        EXPECT_EQ(thingify(locset("cake"), mp), thingify(cake, mp));
-        EXPECT_EQ(thingify(locset("banana"), mp), thingify(banana, mp));
+        EXPECT_EQ(thingify(locset("cake"_lab), mp), thingify(cake, mp));
+        EXPECT_EQ(thingify(locset("banana"_lab), mp), thingify(banana, mp));
 
-        EXPECT_THROW(thingify(locset("durian"), mp), unbound_name);
+        EXPECT_THROW(thingify(locset("durian"_lab), mp), unbound_name);
     }
     {
         label_dict dict;
         dict.set("banana", banana);
         dict.set("cake", cake);
-        dict.set("topping", locset("fruit"));
-        dict.set("fruit", locset("strawberry"));
+        dict.set("topping", locset("fruit"_lab));
+        dict.set("fruit", locset("strawberry"_lab));
 
         EXPECT_THROW(mprovider(morphology(sm), dict), unbound_name);
     }
@@ -122,8 +124,8 @@ TEST(locset, thingify_named) {
         label_dict dict;
         dict.set("banana", banana);
         dict.set("cake", cake);
-        dict.set("topping", locset("fruit"));
-        dict.set("fruit", sum(locset("banana"), locset("topping")));
+        dict.set("topping", locset("(locset \"fruit\")"));
+        dict.set("fruit", sum(locset("banana"_lab), locset("topping"_lab)));
 
         EXPECT_THROW(mprovider(morphology(sm), dict), circular_definition);
     }
@@ -145,17 +147,17 @@ TEST(region, thingify_named) {
         dict.set("cake", cake);
 
         mprovider mp(morphology(sm), dict);
-        EXPECT_EQ(thingify(region("cake"), mp), thingify(cake, mp));
-        EXPECT_EQ(thingify(region("banana"), mp), thingify(banana, mp));
+        EXPECT_EQ(thingify(region("cake"_lab), mp), thingify(cake, mp));
+        EXPECT_EQ(thingify(region("banana"_lab), mp), thingify(banana, mp));
 
-        EXPECT_THROW(thingify(region("durian"), mp), unbound_name);
+        EXPECT_THROW(thingify(region("durian"_lab), mp), unbound_name);
     }
     {
         label_dict dict;
         dict.set("banana", banana);
         dict.set("cake", cake);
-        dict.set("topping", region("fruit"));
-        dict.set("fruit", region("strawberry"));
+        dict.set("topping", region("fruit"_lab));
+        dict.set("fruit", region("(region \"strawberry\")"));
 
         EXPECT_THROW(mprovider(morphology(sm), dict), unbound_name);
     }
@@ -163,8 +165,8 @@ TEST(region, thingify_named) {
         label_dict dict;
         dict.set("banana", banana);
         dict.set("cake", cake);
-        dict.set("topping", region("fruit"));
-        dict.set("fruit", join(region("cake"), region("topping")));
+        dict.set("topping", region("fruit"_lab));
+        dict.set("fruit", join(region("(region \"cake\")"), region("topping"_lab)));
 
         EXPECT_THROW(mprovider(morphology(sm), dict), circular_definition);
     }
diff --git a/test/unit/test_morph_stitch.cpp b/test/unit/test_morph_stitch.cpp
index 1712334b7043717d5264969342818f720b4288ea..a26b35ea58ea8e5428e37a3ed21bf312bbf91af2 100644
--- a/test/unit/test_morph_stitch.cpp
+++ b/test/unit/test_morph_stitch.cpp
@@ -8,11 +8,13 @@
 #include <arbor/morph/place_pwlin.hpp>
 #include <arbor/morph/primitives.hpp>
 #include <arbor/morph/stitch.hpp>
+#include <arbor/string_literals.hpp>
 
 #include "../test/gtest.h"
 #include "morph_pred.hpp"
 
 using namespace arb;
+using namespace arb::literals;
 using testing::region_eq;
 
 TEST(morph, stitch_none_or_one) {
@@ -33,7 +35,7 @@ TEST(morph, stitch_none_or_one) {
     EXPECT_EQ(p2, seg0.dist);
 
     mprovider p(m1, sm1.labels("stitch:"));
-    EXPECT_TRUE(region_eq(p, "stitch:first", reg::segment(0)));
+    EXPECT_TRUE(region_eq(p, "stitch:first"_lab, reg::segment(0)));
 }
 
 TEST(morph, stitch_two) {
@@ -61,8 +63,8 @@ TEST(morph, stitch_two) {
         EXPECT_EQ(p3, seg1.dist);
 
         mprovider p(m, sm.labels("stitch:"));
-        EXPECT_TRUE(region_eq(p, "stitch:0", reg::segment(0)));
-        EXPECT_TRUE(region_eq(p, "stitch:1", reg::segment(1)));
+        EXPECT_TRUE(region_eq(p, "stitch:0"_lab, reg::segment(0)));
+        EXPECT_TRUE(region_eq(p, "stitch:1"_lab, reg::segment(1)));
     }
     {
         // p1 ===== p2
@@ -95,13 +97,13 @@ TEST(morph, stitch_two) {
         mprovider p(m, sm.labels("stitch:"));
         // Branch ordering is arbitrary, so check both possibilities:
         if (seg0.dist == p2) {
-            EXPECT_TRUE(region_eq(p, "stitch:0", reg::segment(0)));
-            EXPECT_TRUE(region_eq(p, "stitch:1", reg::segment(1)));
+            EXPECT_TRUE(region_eq(p, "stitch:0"_lab, reg::segment(0)));
+            EXPECT_TRUE(region_eq(p, "stitch:1"_lab, reg::segment(1)));
         }
         else {
             ASSERT_EQ(p2, seg1.dist);
-            EXPECT_TRUE(region_eq(p, "stitch:0", reg::segment(1)));
-            EXPECT_TRUE(region_eq(p, "stitch:1", reg::segment(0)));
+            EXPECT_TRUE(region_eq(p, "stitch:0"_lab, reg::segment(1)));
+            EXPECT_TRUE(region_eq(p, "stitch:1"_lab, reg::segment(0)));
         }
     }
     {
@@ -138,13 +140,15 @@ TEST(morph, stitch_two) {
         mprovider p(m, sm.labels("stitch:"));
         // Branch ordering is arbitrary, so check both possibilities:
         if (seg2.dist == p2) {
-            EXPECT_TRUE(region_eq(p, "stitch:0", join(reg::segment(0), reg::segment(2))));
-            EXPECT_TRUE(region_eq(p, "stitch:1", reg::segment(1)));
+            EXPECT_TRUE(region_eq(p, "stitch:0"_lab, join(reg::segment(0), reg::segment(2))));
+            EXPECT_TRUE(region_eq(p, "stitch:1"_lab, reg::segment(1)));
         }
         else {
             ASSERT_EQ(p2, seg1.dist);
-            EXPECT_TRUE(region_eq(p, "stitch:0", join(reg::segment(0), reg::segment(1))));
-            EXPECT_TRUE(region_eq(p, "stitch:1", reg::segment(2)));
+            EXPECT_TRUE(region_eq(p, "stitch:0"_lab, join(reg::segment(0), reg::segment(1))));
+            EXPECT_TRUE(region_eq(p, "stitch:1"_lab, reg::segment(2)));
+            EXPECT_TRUE(region_eq(p, "(segment 2)", reg::segment(2)));
+            EXPECT_TRUE(region_eq(p, "(region \"stitch:1\")", reg::segment(2)));
         }
     }
 }
diff --git a/test/unit/test_s_expr.cpp b/test/unit/test_s_expr.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..eb050f5526bca34feb0d4d4d86a534837fbea207
--- /dev/null
+++ b/test/unit/test_s_expr.cpp
@@ -0,0 +1,230 @@
+#include <any>
+
+#include "../test/gtest.h"
+
+#include <arbor/morph/region.hpp>
+#include <arbor/morph/locset.hpp>
+#include <arbor/morph/label_parse.hpp>
+#include <typeinfo>
+
+#include "s_expr.hpp"
+#include "util/strprintf.hpp"
+
+using namespace arb;
+using namespace std::string_literals;
+
+TEST(s_expr, transmogrify) {
+    using map_t = std::unordered_map<char, std::string>;
+    auto transform = [](std::string in, map_t map) {
+        auto t = transmogrifier(in, map);
+        std::string s;
+        while (t) s.push_back(*t++);
+        return s;
+    };
+    EXPECT_EQ(transform("(42,24)", {{',', " "}}), "(42 24)");
+    EXPECT_EQ(transform("(42,24)", {{',', "hello"}}), "(42hello24)");
+    EXPECT_EQ(transform("(42,24)", {{',', " "}}), "(42 24)");
+    EXPECT_EQ(transform("(42,,24)", {{',', " "}}), "(42  24)");
+    map_t asc_map = {{',', " "},
+                     {'|', ")("},
+                     {'<', "(spine "},
+                     {'>', ")"}};
+    EXPECT_EQ(transform("(RGB 128,128,128)", asc_map), "(RGB 128 128 128)");
+    EXPECT_EQ(transform("<color blue>", asc_map), "(spine color blue)");
+    EXPECT_EQ(transform("(1 2 3 | 4 5 6)", asc_map), "(1 2 3 )( 4 5 6)");
+    EXPECT_EQ(transform("", asc_map), "");
+    EXPECT_EQ(transform("<>", asc_map), "(spine )");
+    EXPECT_EQ(transform("<32|>", asc_map), "(spine 32)()");
+}
+
+TEST(s_expr, atoms) {
+    auto get_atom = [](s_expr e) {
+        return e.atom();
+    };
+
+    for (auto spelling: {"foo", "bar_", "car1", "car_1", "x_1__", "x/bar", "x/bar@4.2"}){
+        auto a = get_atom(parse_s_expr(spelling));
+        EXPECT_EQ(a.kind, tok::symbol);
+        EXPECT_EQ(a.spelling, spelling);
+    }
+    // test parsing of integers
+    for (auto spelling: {"0", "-0", "1", "42", "-3287"}){
+        auto a = get_atom(parse_s_expr(spelling));
+        EXPECT_EQ(a.kind, tok::integer);
+        EXPECT_EQ(a.spelling, spelling);
+    }
+    // test parsing of real numbers
+    for (auto spelling: {"0.", "-0.0", "1.21", "42.", "-3287.12"}){
+        auto a = get_atom(parse_s_expr(spelling));
+        EXPECT_EQ(a.kind, tok::real);
+        EXPECT_EQ(a.spelling, spelling);
+    }
+    // test parsing of strings
+    for (auto spelling: {"foo", "dog cat", ""}) {
+        auto s = "\""s+spelling+"\"";
+        auto a = get_atom(parse_s_expr(s));
+        EXPECT_EQ(a.kind, tok::string);
+    }
+}
+
+TEST(s_expr, atoms_in_parens) {
+    auto get_atom = [](s_expr e) {
+        EXPECT_EQ(1u, length(e));
+        EXPECT_TRUE(e.head().is_atom());
+        return e.head().atom();
+    };
+
+    for (auto spelling: {"foo", "bar_", "car1", "car_1", "x_1__", "x/bar", "x/bar@4.2"}){
+        auto a = get_atom(parse_s_expr("("s+spelling+")"));
+        EXPECT_EQ(a.kind, tok::symbol);
+        EXPECT_EQ(a.spelling, spelling);
+    }
+    // test parsing of integers
+    for (auto spelling: {"0", "-0", "1", "42", "-3287"}){
+        auto a = get_atom(parse_s_expr("("s+spelling+")"));
+        EXPECT_EQ(a.kind, tok::integer);
+        EXPECT_EQ(a.spelling, spelling);
+    }
+    // test parsing of real numbers
+    for (auto spelling: {"0.", "-0.0", "1.21", "42.", "-3287.12"}){
+        auto a = get_atom(parse_s_expr("("s+spelling+")"));
+        EXPECT_EQ(a.kind, tok::real);
+        EXPECT_EQ(a.spelling, spelling);
+    }
+    // test parsing of strings
+    for (auto spelling: {"foo", "dog cat", ""}) {
+        auto s = "(\""s+spelling+"\")";
+        auto a = get_atom(parse_s_expr(s));
+        EXPECT_EQ(a.kind, tok::string);
+    }
+}
+
+template <typename L>
+std::string round_trip_label(const char* in) {
+    if (auto x = parse_label_expression(in)) {
+        return util::pprintf("{}", std::any_cast<L>(*x));
+    }
+    else {
+        return x.error().what();
+    }
+}
+
+std::string round_trip_region(const char* in) {
+    if (auto x = parse_region_expression(in)) {
+        return util::pprintf("{}", std::any_cast<arb::region>(*x));
+    }
+    else {
+        return x.error().what();
+    }
+}
+
+std::string round_trip_locset(const char* in) {
+    if (auto x = parse_locset_expression(in)) {
+        return util::pprintf("{}", std::any_cast<arb::locset>(*x));
+    }
+    else {
+        return x.error().what();
+    }
+}
+
+TEST(regloc, round_tripping) {
+    EXPECT_EQ("(cable 3 0 1)", round_trip_label<arb::region>("(branch 3)"));
+    EXPECT_EQ("(intersect (tag 1) (intersect (tag 2) (tag 3)))", round_trip_label<arb::region>("(intersect (tag 1) (tag 2) (tag 3))"));
+    auto region_literals = {
+        "(cable 2 0.1 0.4)",
+        "(region \"foo\")",
+        "(all)",
+        "(tag 42)",
+        "(distal-interval (location 3 0))",
+        "(distal-interval (location 3 0) 3.2)",
+        "(proximal-interval (location 3 0))",
+        "(proximal-interval (location 3 0) 3.2)",
+        "(join (region \"dend\") (all))",
+        "(radius-lt (tag 1) 1)",
+        "(radius-le (tag 2) 1)",
+        "(radius-gt (tag 3) 1)",
+        "(radius-ge (tag 4) 3)",
+        "(intersect (cable 2 0 0.5) (region \"axon\"))",
+    };
+    for (auto l: region_literals) {
+        EXPECT_EQ(l, round_trip_label<arb::region>(l));
+        EXPECT_EQ(l, round_trip_region(l));
+    }
+
+    auto locset_literals = {
+        "(root)",
+        "(locset \"cat man do\")",
+        "(location 3 0.2)",
+        "(terminal)",
+        "(distal (tag 2))",
+        "(proximal (join (tag 1) (tag 2)))",
+        "(uniform (tag 1) 0 100 52)",
+        "(restrict (terminal) (tag 12))",
+        "(join (terminal) (root))",
+        "(sum (terminal) (root))",
+    };
+    for (auto l: locset_literals) {
+        EXPECT_EQ(l, round_trip_label<arb::locset>(l));
+        EXPECT_EQ(l, round_trip_locset(l));
+    }
+}
+
+TEST(regloc, comments) {
+    EXPECT_EQ("(all)",  round_trip_region("(all) ; a comment"));
+    const char *multi_line = 
+        "; comment at start\n"
+        "(radius-lt\n"
+        "    (join\n"
+        "        (tag 3) ; end of line\n"
+        " ;comment on whole line\n"
+        "        (tag 4))\n"
+        "    0.5) ; end of string";
+    EXPECT_EQ("(radius-lt (join (tag 3) (tag 4)) 0.5)",
+              round_trip_region(multi_line));
+}
+
+TEST(regloc, errors) {
+    for (auto expr: {"axon",         // unquoted region name
+                     "(tag 1.2)",    // invalid argument in an otherwise valid region expression
+                     "(tag 1 2)",    // too many arguments to otherwise valid region expression
+                     "(tag 1) (tag 2)", // more than one valid description
+                     "(tag",         // syntax error in region expression
+                     "(terminal)",   // a valid locset expression
+                     "\"my region",  // unclosed quote on label
+                     })
+    {
+        // If an invalid label/expression was passed and handled correctly the parse
+        // call will return without throwing, with the error stored in the return type.
+        // So it is sufficient to assert that it evaluates to false.
+        EXPECT_FALSE(parse_region_expression(expr));
+    }
+
+    for (auto expr: {"axon",         // unquoted locset name
+                     "(location 1 \"0.5\")",  // invalid argument in an otherwise valid locset expression
+                     "(location 1 0.2 0.2)",  // too many arguments to otherwise valid locset expression
+                     "(root) (location 2 0)", // more than one valid description
+                     "(tag",         // syntax error in locset expression
+                     "(tag 3)",      // a valid region expression
+                     "\"my locset",  // unclosed quote on label
+                     })
+    {
+        // If an invalid label/expression was passed and handled correctly the parse
+        // call will return without throwing, with the error stored in the return type.
+        // So it is sufficient to assert that it evaluates to false.
+        EXPECT_FALSE(parse_locset_expression(expr));
+    }
+
+    for (auto expr: {"axon",         // unquoted locset name
+                     "(location 1 \"0.5\")",  // invalid argument in an otherwise valid locset expression
+                     "(location 1 0.2 0.2)",  // too many arguments to otherwise valid locset expression
+                     "(root) (location 2 0)", // more than one valid description
+                     "(tag",         // syntax error in locset expression
+                     "\"my locset",  // unclosed quote on label
+                     })
+    {
+        // If an invalid label/expression was passed and handled correctly the parse
+        // call will return without throwing, with the error stored in the return type.
+        // So it is sufficient to assert that it evaluates to false.
+        EXPECT_FALSE(parse_label_expression(expr));
+    }
+}