diff --git a/.gitignore b/.gitignore
index 3d8f66c132830ab61a84c4bbf2aec4513c1d7af3..c3bebcb1e538aa633969355770bf1e531d8cf6a5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -69,6 +69,9 @@ commit.msg
 # Mac default files
 .DS_Store
 
+# generated documentation images
+doc/gen-images
+
 # by product of generating a source distribution for pip
 *.egg-info
 dist
diff --git a/.readthedocs.yml b/.readthedocs.yml
new file mode 100644
index 0000000000000000000000000000000000000000..b8fccfee955769597aae375080063fd0d093a853
--- /dev/null
+++ b/.readthedocs.yml
@@ -0,0 +1,14 @@
+# Read the Docs configuration file
+# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
+
+version: 2
+
+sphinx:
+  configuration: doc/conf.py
+
+formats: all
+
+python:
+  version: 3.7
+  install:
+    - requirements: doc/requirements.txt
diff --git a/doc/cable_cell.rst b/doc/cable_cell.rst
new file mode 100644
index 0000000000000000000000000000000000000000..0cc1f13af5e512f0f2ab21b97eb92d24f613776c
--- /dev/null
+++ b/doc/cable_cell.rst
@@ -0,0 +1,388 @@
+.. _cablecell:
+
+Cable Cells
+===========
+
+An Arbor *cable cell* is a full description of a cell with morphology and cell
+dynamics, where cell dynamics include ion species and their properties, ion
+channels, synapses, gap junction sites, stimulii and spike detectors.
+Arbor cable cells are constructed from a morphology and a label dictionary,
+and provide a rich interface for specifying the cell's dynamics.
+
+.. note::
+    Before reading this page, it is recommended that you first read about
+    :ref:`morphology descriptions <morph-morphology>`, and also
+    :ref:`label dictionary <labels-dictionary>` that are used to describe
+    :ref:`locations <labels-locset>` and :ref:`regions <labels-region>` on a cell.
+
+.. _cablecell-decoration:
+
+Decoration
+----------------
+
+A cable cell is *decorated* by specifying the distribution and placement of dynamics
+on the cell to produce a full description
+of a cell morphology and its dynamics with all information required to build
+a standalone single-cell model, or as part of a larger network.
+
+Decoration uses region and locset descriptions, with
+their respective use for this purpose reflected in the two broad classes
+of dynamics in Arbor:
+
+* *Painted dynamics* are applied to regions of a cell, and are associated with
+  an area of the membrane or volume of the cable.
+
+  * :ref:`Cable properties <cable-properties>`.
+  * :ref:`Density mechanisms <cable-density-mechs>`.
+  * :ref:`Ion species <cable-ions>`.
+
+* *Placed dynamics* are applied to locations on the cell, and are associated
+  with entities that can be counted.
+
+  * :ref:`Synapses <cable-synapses>`.
+  * :ref:`Gap junction sites <cable-gj-sites>`.
+  * :ref:`Threshold detectors <cable-threshold-detectors>` (spike detectors).
+  * :ref:`Stimulii <cable-stimulii>`.
+  * :ref:`Probes <cable-probes>`.
+
+.. _cablecell-paint:
+
+Painted Dynamics
+''''''''''''''''
+
+Painted dynamics are applied to a subset of the surface and/or volume of cells.
+They can be specified at three different levels:
+
+* *globally*: a global default for all cells in a model.
+* *per-cell*: overide the global defaults for a specific cell.
+* *per-region*: specialize on specific cell regions.
+
+This hierarchical approach for resolving parameters and properties allows
+us to, for example, define a global default value for calcium concentration,
+then provide a different values on specific cell regions.
+
+Some dynamics, such as membrane capacitance and the initial concentration of ion species
+must be defined for all compartments. Others need only be applied where they are
+present, for example ion channels.
+The types of dynamics, and where they can be defined, are
+:ref:`tabulated <cable-painted-resolution>` below.
+
+.. _cable-painted-resolution:
+
+.. csv-table:: Painted property resolution options.
+   :widths: 20, 10, 10, 10
+
+                  ,       **region**, **cell**, **global**
+   cable properties,       ✓, ✓, ✓
+   ion initial conditions, ✓, ✓, ✓
+   density mechnism,       ✓, --, --
+   ion rev pot mechanism,  --, ✓, ✓
+   ion valence,            --, --, ✓
+
+If a property is defined at multiple levels, the most local definition will be chosen:
+a cell-local definition will override a global definition, and a definition on a region
+will override any cell-local or global definition on that region.
+
+.. warning::
+    If a property is defined on two regions that overlap, it is not possible to
+    deterministically choose the correct definition, and an error will be
+    raised during model instantiation.
+
+.. _cable-properties:
+
+Cable properties
+~~~~~~~~~~~~~~~~
+
+There are four cable properties that are defined everywhere on all cables:
+
+* *Vm*: Initial membrane voltage [mV].
+* *cm*: Membrane capacitance [F/m²].
+* *rL*: Axial resistivity of cable [Ω·cm].
+* *tempK*: Temperature [Kelvin].
+
+In Python, the :py:class:`cable_cell` interface provides the :py:func:`cable_cell.set_properties` method
+for setting cell-wide defaults for properties, and the
+:py:meth:`cable_cell.paint` interface for overriding properties on specific regions.
+
+.. code-block:: Python
+
+    import arbor
+
+    # Load a morphology from file and define basic regions.
+    tree = arbor.load_swc('granule.swc')
+    morph = arbor.morphology(tree, spherical_root=True)
+    labels = arbor.label_dict({'soma': '(tag 1)', 'axon': '(tag 2)', 'dend': '(tag 3)'})
+
+    # Create a cable cell.
+    cell = arbor.cable_cell(morph, labels)
+
+    # Set cell-wide properties that will be applied by default to # the entire cell.
+    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)
+
+.. _cable-density-mechs:
+
+Density mechanisms
+~~~~~~~~~~~~~~~~~~~~~~
+
+Regions can have density mechanisms defined over their extents.
+Density mechanisms are :ref:`NMODL mechanisms <nmodl>`
+which describe biophysical processes. These are processes
+that are distributed in space, but whose behaviour is defined purely
+by the state of the cell and the process at any given point.
+
+The most common use for density mecahnisms is to describe ion channel dynamics,
+for example the ``hh`` and ``pas`` mechanisms provided by NEURON and Arbor,
+which model classic Hodgkin-Huxley and passive leaky currents respectively.
+
+Mechanisms have two types of parameters that can be set by users
+
+* *Global* parameters are a single scalar value that is the
+  same everywhere a mechanism is defined.
+* *Range* parameters can vary spatially.
+
+Every mechanism is described by a string with its name, and
+an optional list of key-value pairs that define its range parameters.
+
+Because a global parameter is fixed over the entire spatial extent
+of a density mechanism, a new mechanism has to created for every
+combination of global parameter values.
+
+Take for example a mechanism passive leaky dynamics:
+
+* Name: ``"passive"``.
+* Global variable: reversal potential ``"el"``.
+* Range variable: conductance ``"g"``.
+
+.. code-block:: Python
+
+    # Create pas mechanism with default parameter values (set in NOMDL file).
+    m1 = arbor.mechanism('passive')
+
+    # Create default mechainsm with custom conductance (range)
+    m2 = arbor.mechanism('passive', {'g', 0.1})
+
+    # Create a new pas mechanism with that changes reversal potential (global)
+    m3 = arbor.mechanism('passive/el=-45')
+
+    # 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?
+
+.. _cable-ions:
+
+Ion species
+~~~~~~~~~~~
+
+Arbor allows arbitrary ion species to be defined, to extend the default
+calcium, potassium and sodium ion species.
+A ion species is defined globally by its name and valence, which
+can't be overriden at cell or region level.
+
+.. csv-table:: Default ion species in Arbor
+   :widths: 15, 10, 10
+
+   **Ion**,     **name**, **Valence**
+   *Calcium*,   ca,       1
+   *Potassium*,  k,       1
+   *Sodium*,    na,       2
+
+Each ion species has the following properties:
+
+1. *internal concentration*: concentration on interior of the membrane [mM].
+2. *external concentration*: concentration on exterior of the membrane [mM].
+3. *reversal potential*: reversal potential [mV].
+4. *reversal potential mechanism*:  method for calculating reversal potential.
+
+Properties 1, 2 and 3 must be defined, and are used as the initial values for
+each quantity at the start of the simulation. They are specified globally,
+then specialised at cell and region level.
+
+The reversal potential of an ion species is calculated by an
+optional *reversal potential mechanism*.
+If no reversal potential mechanism is specified for an ion species, the initial
+reversal potential values are maintained for the course of a simulation.
+Otherwise, the mechanism does the work.
+
+but it is subject to some strict restrictions.
+Specifically, a reversal potential mechanism described in NMODL:
+
+* May not maintain any STATE variables.
+* Can only write to the "eX" value associated with an ion.
+* Can not be a POINT mechanism.
+
+Essentially, reversal potential mechanisms must be pure functions of cellular
+and ionic state.
+
+.. note::
+    Arbor imposes greater restrictions on mechanisms that update ionic reversal potentials
+    than NEURON. Doing so simplifies reasoning about interactions between
+    mechanisms that share ionic species, by virtue of having one mechanism, and one
+    mechanism only, that calculates reversal potentials according to concentrations
+    that the other mechanisms use and modify.
+
+If a reversal potential mechanism that writes to multiple ions,
+it must be given for either no ions, or all of the ions it writes.
+
+Arbor's default catalogue includes a *nernst* reversal potential, which is
+parameterized over a single ion. For example, to bind it to the calcium
+ion at the cell level using the Python interface:
+
+.. code-block:: Python
+
+    cell = arbor.cable_cell(morph, labels)
+
+    # method 1: create the mechanism explicitly.
+    ca = arbor.mechanism('nernst/x=ca')
+    cell.set_ion(ion='ca', method=ca)
+
+    # method 2: set directly using a string description
+    cell.set_ion(ion='ca', method='nernst/x=ca')
+
+
+The NMODL code for the
+`Nernst mechanism  <https://github.com/arbor-sim/arbor/blob/master/mechanisms/mod/nernst.mod>`_
+can be used as a guide for how to calculate reversal potentials.
+
+While the reversal potential mechanism must be the same for a whole cell,
+the initial concentrations and reversal potential can be localised for regions
+using the *paint* interface:
+
+.. code-block:: Python
+
+    # cell is an arbor.cable_cell
+
+    # 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))
+
+    # Alternatively, one can selectively overwrite the global defaults.
+    cell.paint('axon', arbor.ion('ca', rev_pot=126)
+
+.. _cablecell-place:
+
+Placed Dynamices
+''''''''''''''''
+
+Placed dynamics are discrete countable items that affect or record the dynamics of a cell,
+and are asigned to specific locations.
+
+.. _cable-synapses:
+
+Synapses
+~~~~~~~~
+
+Synapses are instances of NMODL POINT mechanisms.
+
+.. _cable-gj-sites:
+
+Gap junction sites
+~~~~~~~~~~~~~~~~~~
+
+.. _cable-threshold-detectors:
+
+Threshold detectors (spike detectors).
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. _cable-stimulii:
+
+Stimulii
+~~~~~~~~
+
+.. _cable-probes:
+
+Probes
+~~~~~~
+
+Python API
+----------
+
+Creating a cable cell
+'''''''''''''''''''''
+
+.. py:class:: cable_cell
+
+    A cable cell is constructed from a :ref:`morphology <morph-morphology>`
+    and an optional :ref:`label dictionary <labels-dictionary>`.
+
+    .. note::
+        The regions and locsets defined in the label dictionary are
+        :ref:`concretised <labels-concretise>` when the cable cell is constructed,
+        and an exception will be thrown if an invalid label expression is found.
+
+        There are two reasons an expression might be invalid:
+
+        1. Explicitly refers to a location of cable that does not exist in the
+           morphology, for example ``(branch 12)`` on a cell with 6 branches.
+        2. Incorrect label reference: circular reference, or a label that does not exist.
+
+
+    .. code-block:: Python
+
+        import arbor
+
+        # Construct the morphology from an SWC file.
+        tree = arbor.load_swc('granule.swc')
+        morph = arbor.morphology(tree, spherical_root=True)
+
+        # Define regions using standard SWC tags
+        labels = arbor.label_dict({'soma': '(tag 1)',
+                                   'axon': '(tag 2)',
+                                   'dend': '(join (tag 3) (tag 4))'})
+
+        # Construct a cable cell.
+        cell = arbor.cable_cell(morph, labels)
+
+    .. method:: set_properties(Vm=None, cm=None, rL=None, tempK=None)
+
+        Set default values of cable properties on the whole cell.
+        Overrides the default global values, and can be overriden by painting
+        the values onto regions.
+
+        :param str region: name of the region.
+        :param Vm: Initial membrane voltage [mV].
+        :type Vm: float or None
+        :param cm: Membrane capacitance [F/m²].
+        :type cm: float or None
+        :param rL: Axial resistivity of cable [Ω·cm].
+        :type rL: float or None
+        :param tempK: Temperature [Kelvin].
+        :type tempK: float or None
+
+        .. code-block:: Python
+
+            # Set cell-wide values for properties
+            cell.set_properties(Vm=-70, cm=0.01, rL=100, tempK=280)
+
+
+    .. method:: paint(region, [Vm=None, cm=None, rL=None, tempK=None])
+
+        Set cable properties on a region.
+
+        :param str region: name of the region.
+        :param Vm: Initial membrane voltage [mV].
+        :type Vm: float or None
+        :param cm: Membrane capacitance [F/m²].
+        :type cm: float or None
+        :param rL: Axial resistivity of cable [Ω·cm].
+        :type rL: float or None
+        :param tempK: Temperature [Kelvin].
+        :type tempK: float or None
+
+        .. 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)
+
+.. py:class:: ion
+
+    properties of an ionic species.
+
diff --git a/doc/conf.py b/doc/conf.py
index 04557de4d49aec307251135981273fb9d8fb4151..ac511766a7be0c685643f4789b68d1fd631b2b7e 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -1,13 +1,18 @@
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 
+
 html_static_path = ['static']
 
 def setup(app):
     app.add_css_file('custom.css')
     app.add_object_type('generic', 'gen', 'pair: %s; generic')
+    app.add_object_type('label', 'lab', 'pair: %s; label')
 
-extensions = ['sphinx.ext.todo', 'sphinx.ext.mathjax']
+extensions = [
+    'sphinx.ext.todo',
+    'sphinx.ext.mathjax',
+]
 source_suffix = '.rst'
 master_doc = 'index'
 
@@ -22,3 +27,27 @@ html_theme = "sphinx_rtd_theme"
 html_theme_options = {
     'logo_only': True,
     'style_nav_header_background': '#dfdcdf'}
+
+# This style makes the source code pop out a bit more
+# from the background text, without being overpowering.
+pygments_style = 'perldoc'
+
+# Generate images for the documentation.
+print("--- generating images ---")
+import sys, os
+
+this_path=os.path.split(os.path.abspath(__file__))[0]
+
+# Location of scripts used to generate images
+script_path=this_path+'/scripts'
+sys.path.append(script_path)
+
+# Output path for generated images
+img_path=this_path+'/gen-images'
+if not os.path.exists(img_path):
+    os.mkdir(img_path)
+
+import make_images
+make_images.generate(img_path)
+
+print("-------------------------")
diff --git a/doc/index.rst b/doc/index.rst
index b16a15b57506a7d82c6c60b9f01f1c320f6b2073..c1319eed6ecec8ee5198226dba9994d52b26b2ac 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -43,6 +43,14 @@ Alternative citation formats for the paper can be `downloaded here <https://ieee
    python
    single_cell
 
+.. toctree::
+   :caption: Concepts:
+
+   morphology
+   labels
+   cable_cell
+   mechanisms
+
 .. toctree::
    :caption: Arbor Models:
 
@@ -80,6 +88,7 @@ Alternative citation formats for the paper can be `downloaded here <https://ieee
    :caption: Developers:
 
    library
+   nmodl
    simd_api
    profiler
    sampling_api
diff --git a/doc/labels.rst b/doc/labels.rst
new file mode 100644
index 0000000000000000000000000000000000000000..dce13697842b7788390c3abaa10918952882feab
--- /dev/null
+++ b/doc/labels.rst
@@ -0,0 +1,762 @@
+.. _labels:
+
+Labels
+=========
+
+Arbor provides a domain specific language (DSL) for describing regions and
+locations on morphologies, and a dictionary for assiciating these descriptions
+with a string label.
+
+The labels are used to refer to regions
+and locations when setting cell properties and attributes.
+For example, the membrane capacitance on a region of the cell membrane, or
+the location of synapse instances.
+
+Example Cell
+------------
+
+The following morphology is used on this page to illustrate region and location
+descriptions. It has a soma, dendritic tree and an axon with a hillock:
+
+.. _labels-morph-fig:
+
+.. figure:: gen-images/label_morph.svg
+  :width: 800
+  :align: left
+
+  Segments of the morphology are colored according to tags:
+  soma (tag 1, red), axon (tag 2, grey), dendrites (tag 3, blue) (left).
+  The 6 branches of the morphology with their branch ids (right).
+
+*Branch 0* contains the soma, which is modelled as a cylinder of length and diameter 4 μm,
+and the proximal unbranched section of the dendritic tree which has a radius of 0.75 μm,
+attached to the distal end of the soma.
+
+The other branches in the dendritic tree have the following properties:
+
+* *branch 1* tapers from 0.4 to 0.2 μm;
+* *branch 2* has a constant radius of 0.5 μm;
+* *branch 3* tapers from 0.5 to 0.2 μm;
+* *branch 4* tapers from 0.5 to 0.2 μm.
+
+*Branch 5* is the axon, composed of two cable segments: an axon hillock with a radius that
+tapers from 4 μm to 0.4 μm attached to the proximal end of the soma; and the start of the
+axon proper with constant radius 0.4 μm.
+
+Label Types
+------------
+
+.. _labels-locset:
+
+Locsets
+~~~~~~~~~~~
+
+A *locset* is a set of locations on a morphology, specifically a *multiset*,
+which may contain multiple instances of the same location, for example:
+
+* The center of the soma.
+* The locations of inhibitory synapses.
+* The tips of the dendritic tree.
+
+.. figure:: gen-images/locset_label_examples.svg
+  :width: 800
+  :align: center
+
+  Examples of locsets on the example morphology.
+  The terminal points (right).
+  Fifty random locations on the dendritic tree (left).
+  The :ref:`root <morph-segment-definitions>` of the morphology is shown with a red circle
+  for reference.
+
+
+.. _labels-region:
+
+Regions
+~~~~~~~~~~~~
+
+A *region* is a subset of a morphology's cable segments, for example:
+
+* The soma.
+* The dendritic tree.
+* An explicit reference to a specific unbranched cable, e.g. "branch 3" or "the distal half of branch 1".
+* The axon hillock.
+* The dendrites with radius less than 1 μm.
+
+It is possible for a region to be empty, for example, a region that defines
+the axon will be empty on a morphology that has no axon.
+
+Regions do not need to be complete sub-trees of a morphology, for example,
+the region of cables that have radius less than 0.5 μm
+:ref:`below <labels-region-examples>` is composed of three disjoint sub-trees.
+
+.. _labels-region-examples:
+
+.. figure:: gen-images/region_label_examples.svg
+  :width: 800
+  :align: center
+
+  Examples of regions on the example morphology. **Left**: The dendritic tree.
+  **Right**: All cables with radius less than 0.5 μm.
+
+.. _labels-expressions:
+
+Expressions
+-----------
+
+Regions and locsets are described using *expressions* written with the DSL.
+
+Examples of expressions that define regions include:
+
+* ``(all)``: the complete cell morphology.
+* ``(tag 1)``: all segments with tag 1.
+* ``(branch 2)``: branch 2.
+* ``(region "soma")``: the region with the label "soma".
+
+Examples of expressions that define locsets include:
+
+* ``(root)``: the location of the :ref:`root points <morph-segment-definitions>`.
+* ``(terminal)``: the locations of the :ref:`terminal points <morph-segment-definitions>`.
+* ``(location 3 0.5)``: the mid point of branch 3.
+* ``(locset "synapse-sites")``: the locset labelled "synapse-sites".
+
+Detailed descriptions for all of the region and locset expression types is
+given :ref:`below <labels-expr-docs>`.
+
+.. note::
+    The example expressions above will look familiar to readers who have
+    use the Lisp programming language. This is because both the DSL and Lisp use
+    *s-expressions*, which are a simple way to represent a nested list of data.
+
+    However, the DSL is not a dialect of Lisp, and has very simple semantics
+    that are only limited to describing morphology features.
+
+Expressions are *composable*, so that expressions can be constructed
+from simple expressions. For example, the expression:
+
+.. code-block:: lisp
+
+    (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.
+
+.. note:
+
+    In NEURON *prescriptive* hoc templates are typically used to calculate
+    explicit lists of sections or segments using loops and logical constructs.
+    The logic in a hoc template often makes it difficult to understand
+    what the results describe, and is error prone.
+
+    Arbor expressions are *descriptive*, in that they describe *what* a
+    region or locset is, not *how* it is to be computed.
+    As a result, label dictionaries are much more concise and easy to interpret for
+    consumers of a model than hoc templates.
+    Furthermore they are less error prone because
+    Arbor handles generation of conrete cable sections and locations when
+    expressions are applied to a morphology.
+
+.. _labels-expr-docs:
+
+Expression Syntax
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The DSL uses `s-expressions <https://en.wikipedia.org/wiki/S-expression>`_, which are composed of the following basic types:
+
+.. generic:: string
+
+    A string literal enclosed in quotes, e.g. ``"dendrites"``.
+
+.. generic:: integer
+
+    An integer. e.g: ``42``, ``-2``, ``0``.
+
+.. generic:: real
+
+    A floating point value. e.g: ``2``, ``4.3``, ``.3``, ``-2.1e3``.
+
+.. generic:: region
+
+    An expression that evaluates to a region. e.g. ``(all)``, ``(tag 3)``, ``(intersect (tag 3) (tag 4))``.
+
+.. generic:: locset
+
+    An expression that evaluates to a locset. e.g. ``(root)``, ``(location 3 0.2)``, ``(proximal (tag 2))``.
+
+Expressions can be written over multiple lines, and comments are marked with semi-colon.
+This can be used to make more complex expression easier to read, for example the
+following region that finds all the sub-trees that start at the locations on the
+dendritic tree where the radius first is less than or equal to 0.2 μm.
+
+.. code:: lisp
+
+    (distal_interval                   ; take subtrees that start at
+        (proximal                      ; locations closest to the soma
+            (radius_lte                ; with radius <= 0.2 um
+                (join (tag 3) (tag 4)) ; on basal and apical dendrites
+                0.2)))
+
+.. note::
+    If the expression above at first seems a little complex, consider how the same
+    thing could be achieved using hoc in NEURON, and whether it would be free of bugs
+    and applicable to arbitrary morphologies.
+
+
+Locset Expressions
+~~~~~~~~~~~~~~~~~~~~~
+
+.. figure:: gen-images/label_branch.svg
+  :width: 800
+  :align: center
+
+  The input morphology with branch numbers for reference in the examples below.
+
+
+.. label:: (root)
+
+    The location of the root.
+
+    Equivalent to ``(location 0 0)``.
+
+    .. figure:: gen-images/root_label.svg
+      :width: 300
+      :align: center
+
+.. _labels-location-def:
+
+.. label:: (location branch:integer pos:real)
+
+    A location on ``branch``, where ``0 ≤ pos ≤ 1`` gives the relative position
+    between the proximal and distal ends of the branch. The position is in terms
+    of branch length, so for example, on a branch of length 100 μm ``pos=0.2``
+    corresponds to 20 μm from the proximal end, or 80 μm from the distal end.
+
+    .. figure:: gen-images/location_label.svg
+      :width: 300
+      :align: center
+
+      The result of ``(location 1 0.5)``, which corresponds to the mid point of branch 1.
+
+.. label:: (terminal}
+
+    The location of terminal points, which are the most distal locations on the morphology.
+    These will typicall correspond to the tips, or end points, of dendrites and axons.
+
+    .. figure:: gen-images/term_label.svg
+      :width: 300
+      :align: center
+
+      The terminal points, generated with ``(terminal)``.
+
+.. label:: (uniform reg:region first:int last:int seed:int)
+
+    .. figure:: gen-images/uniform_label.svg
+      :width: 600
+      :align: center
+
+      Ten random locations on the dendrites drawn using different random seeds.
+      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)
+
+    The set of locations ``{(location b pos) | 0 ≤ b < nbranch-1}``.
+
+    .. figure:: gen-images/on_branches_label.svg
+      :width: 300
+      :align: center
+
+      The set of locations at the midpoint of every branch, expressed as ``(on_branches 0.5)``.
+
+.. label:: (distal reg:region)
+
+    The set of the most distal locations of a region.
+    These are defined as the locations for which there are no other locations more distal in the region.
+
+    .. figure:: gen-images/distal_label.svg
+      :width: 600
+      :align: center
+
+      On the left is the region with radius between 0.3 μm and 0.5 μm.
+      The right shows the distal set of this region.
+
+.. label:: (proximal reg:region)
+
+    The set of the most proximal locations of a region.
+    These are defined as the locations for which there are no other locations more proximal in the region.
+
+    .. figure:: gen-images/proximal_label.svg
+      :width: 600
+      :align: center
+
+      On the left is the region with radius between 0.3 μm and 0.5 μm.
+      The right shows the proximal set of this region.
+
+.. label:: (locset name:string)
+
+    Refer to a locset by its label. For example, ``(locset "synapse-sites")`` could be used in an expression to refer
+    to a locset with the name ``"synapse-sites"``.
+
+.. label:: (restrict locations:locset reg:region)
+
+    The set of locations in the locset ``loc`` that are in the region ``reg``.
+
+    .. figure:: gen-images/restrict_label.svg
+      :width: 600
+      :align: center
+
+      The result of restricting the terminal locations (left) onto the dendritic tree (middle) is the tips of the dendritic tree (right).
+
+      .. code-block:: lisp
+
+        (restrict (terminal) (tag 3))
+
+
+.. label:: (join lhs:locset rhs:locset [...locset])
+
+    Set intersection for two locsets, with duplicates removed and results sorted.
+    For example, the following:
+
+    .. code-block:: lisp
+
+        (join
+            (join (location 1 0.5) (location 2 0.1) (location 1 0.2))
+            (join (location 1 0.5) (location 4 0)))
+
+    Gives the following:
+
+    .. code-block:: lisp
+
+        (join (location 1 0.2) (location 1 0.5) (location 2 0.1) (location 4 0))
+
+    Note that ``(location 1 0.5)`` occurs in both the sets, and occurs only once in the result.
+
+.. label:: (sum lhs:locset rhs:locset [...locset])
+
+    Multiset summation of two locsets, such that ``(sum lhs rhs) = A + B``, where A and B are multisets of locations.
+    This is equivalent to contactenating the two lists, and the length of the result is the sum of
+    the lenghts of the inputs. For example:
+
+    .. code-block:: lisp
+
+        (sum
+            (join (location 1 0.5) (location 2 0.1) (location 1 0.2))
+            (join (location 1 0.5) (location 4 0)))
+
+    Gives the following:
+
+    .. code-block:: lisp
+
+        (join (location 1 0.5) (location 2 0.1) (location 1 0.2) (location 1 0.5) (location 4 0))
+
+Region Expressions
+~~~~~~~~~~~~~~~~~~~~~
+
+.. label:: (nil)
+
+    An empty region.
+
+.. label:: (all)
+
+    All branches in the morphology.
+
+    .. figure:: gen-images/nil_all_label.svg
+      :width: 600
+      :align: center
+
+      The trivial region definitions ``(nil)`` (left) and ``(all)`` (right).
+
+.. label:: (tag tag_id:integer)
+
+    All of the segments with :ref:`tag <morph-tag-definition>` ``tag_id``.
+
+    .. figure:: gen-images/tag_label.svg
+      :width: 900
+      :align: center
+
+      The soma, axon and dendritic tree, selected using ``(tag 1)``, ``(tag 2)``, and ``(tag 3)`` respectively.
+
+
+.. label:: (branch branch_id:integer)
+
+    Refer to a branch by its id.
+
+    .. figure:: gen-images/branch_label.svg
+      :width: 600
+      :align: center
+
+      Branches 0 and 3, selected using ``(branch 0)`` and ``(branch 3)`` respectively.
+
+.. _labels-cable-def:
+
+.. label:: (cable branch_id:integer prox:real dist:real)
+
+    An unbranched cable that is a subset of ``branch``.
+    The values of ``0 ≤ prox ≤ dist ≤ 1`` are the relative position
+    of the ends of the branch. The positions are in terms
+    of branch length, so for example, on a branch of length 100 μm ``prox=0.2, dist=0.8``
+    would give a cable that starts and ends 20 μm and 80 μm along the branch
+    respectively.
+
+    .. figure:: gen-images/cable_label.svg
+      :width: 600
+      :align: center
+
+      Selecting parts of branch 1, from left to right: ``(cable 1 0 1)`` to select the
+      whole branch, ``(cable 1 0.3 1)`` and ``(cable 0 0.3 0.7)`` to select part of the branch.
+
+.. label:: (region name:string)
+
+    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)
+
+    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.
+    The distal interval of the locset ``start`` is the union of the distal interval of each location in ``start``.
+
+    .. figure:: gen-images/distint_label.svg
+      :width: 600
+      :align: center
+
+      On the left is a locset of 3 locations: 1 on the axon and 2 in the dendritic tree.
+      The right shows the locset's distal interval with extent 5 μm, formed with the following expression:
+
+      .. code-block:: lisp
+
+        (distal_interval (sum (location 1 0.5) (location 2 0.7) (location 5 0.1)) 5)
+
+.. 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``.
+
+    .. figure:: gen-images/distintinf_label.svg
+      :width: 600
+      :align: center
+
+      On the left is a locset of 3 locations: 1 on the axon and 2 in the dendritic tree.
+      The right shows the locset's distal interval formed with the following expression:
+
+      .. code-block:: lisp
+
+        (distal_interval (sum (location 1 0.5) (location 2 0.7) (location 5 0.1)))
+
+
+.. 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.
+    The proximal interval of the locset ``start`` is the union of the proximal interval of each location in ``start``.
+
+    .. figure:: gen-images/proxint_label.svg
+      :width: 600
+      :align: center
+
+      On the left is a locset with two locations on separate sub-trees of the dendritic tree.
+      On the right is their proximal interval with an ``extent`` of 5 μm, formed as follows:
+
+      .. code-block:: lisp
+
+        (proximal_interval (sum (location 1 0.8) (location 2 0.3)) 5)
+
+.. label:: (proximal_interval start:locset)
+
+    When no ``extent`` distance is provided, the proximal intervals are extended to the root location.
+
+    .. figure:: gen-images/proxintinf_label.svg
+      :width: 600
+      :align: center
+
+      On the left is a locset with two locations on separate sub-trees of the dendritic tree.
+      On the right is their proximal interval formed as follows:
+
+      .. code-block:: lisp
+
+        (proximal_interval (sum (location 1 0.8) (location 2 0.3)))
+
+.. label:: (radius_lt reg:region radius:real)
+
+    All parts of cable segments in the region ``reg`` with radius less than ``radius``.
+
+    .. figure:: gen-images/radiuslt_label.svg
+      :width: 300
+      :align: center
+
+      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)
+
+.. label:: (radius_le reg:region radius:real)
+
+    All parts of cable segments in the region ``reg`` with radius less than or equal to ``radius``.
+
+    .. figure:: gen-images/radiusle_label.svg
+      :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
+      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)
+
+.. label:: (radius_gt reg:region radius:real)
+
+    All parts of cable segments in the region ``reg`` with radius greater than ``radius``.
+
+    .. figure:: gen-images/radiusgt_label.svg
+      :width: 300
+      :align: center
+
+      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)
+
+.. label:: (radius_ge reg:region radius:real)
+
+    All parts of cable segments in the region ``reg`` with radius greater than or equal to ``radius``.
+
+    .. figure:: gen-images/radiusge_label.svg
+      :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
+      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)
+
+.. label:: (join lhs:region rhs:region [...region])
+
+    The union of two or more regions.
+
+    .. figure:: gen-images/union_label.svg
+      :width: 900
+      :align: center
+
+      Two regions (left and middle) and their union (right).
+
+.. label:: (intersect lhs:region rhs:region [...region])
+
+    The intersection of two or more regions.
+
+    .. figure:: gen-images/intersect_label.svg
+      :width: 900
+      :align: center
+
+      Two regions (left and middle) and their intersection (right).
+
+.. _labels-concretise:
+
+Concretisation
+----------------
+
+When a region or locset expression is applied to a cell morphology it is
+*concretised*. Concretising a locset will return a set of *locations* on the
+morphology, and concretising a region will return a list of unbranched *cables*
+on the morphology.
+
+.. note::
+    Applying an expression to different morphologies may give different
+    concretised results.
+
+Locations
+~~~~~~~~~
+
+A *location* on a cell is described using a tuple ``(branch, pos)``, where
+``branch`` is a branch id, and ``0 ≤ pos ≤ 1`` is the relative distance along
+the branch, given that 0 and 1 are the proximal and distal ends of the branch
+respectively.
+
+Examples of locations, :ref:`expressed using the DSL <labels-location-def>`, include:
+
+* The root ``(location 0 0)``.
+* The start of branch 5 ``(location 5 0)``.
+* The end of branch 5 ``(location 5 1)``.
+* One quarter of the way along branch 5 ``(location 5 0.25)``.
+
+Cables
+~~~~~~~~~
+
+An unbranched *cable* is a tuple of the form ``(branch, prox, dist)``,
+where ``branch`` is the branch id, and ``0 ≤ prox ≤ dist ≤ 1`` define the relative position
+of the end points of the section on the branch.
+
+Examples of cables, :ref:`expressed using the DSL <labels-cable-def>`, include:
+
+* All of branch 2 ``(cable 2 0 1)``.
+* The middle third of branch 2 ``(cable 2 0.333 0.667)``.
+* A zero length cable in the middle of branch 2 ``(cable 2 0.5 0.5)``.
+
+.. note::
+    Zero length cables are permitted.
+    They are not useful for defining membrane properties, which are applied to
+    the surface of a region.
+    However, they can occur as the result of sub-expressions in larger
+    expressions that define non-trivial regions and locsets.
+
+.. _labels-dictionary:
+
+Label Dictionaries
+------------------
+
+*Labels* can be assigned to expressions, and used to refer to the expression or the
+concrete region or locset generated when the expression is applied to a morphology.
+A label is a string with the following rules:
+
+* may contain alpha-numeric values, ``{a-z}[A-z][0-9]``, underscore ``_`` and hyphen ``-``.
+* no leading underscore, hyphen or numeric values: for example ``_myregion``,
+  ``-samples``, and ``2ndpoint`` are invalid labels.
+
+labels are stored with their associated expressions as key-value pairs in *label dictionaries*.
+
+Python API
+----------
+
+The ``arbor.label_dict`` type is used for creating and manipulating label dictionaries,
+which can be initialised with a dictionary that defines (label, expression)
+pairs. For example, a dictionary that uses tags that correspond to SWC
+`structure identifiers <http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html>`_
+to label soma, axon, dendrite and apical dendrites is:
+
+
+.. code-block:: python
+
+    import arbor
+
+    labels = {'soma': '(tag 1)',
+              'axon': '(tag 2)',
+              'dend': '(tag 3)',
+              'apic': '(tag 4)'}
+
+    d = arbor.label_dict(labels)
+
+Alternatively, start with an empty label dictionary and add the labels and
+their definitions one by one:
+
+.. code-block:: python
+
+    import arbor
+
+    d = arbor.label_dict()
+
+    d['soma'] = '(tag 1)'
+    d['axon'] = '(tag 2)'
+    d['dend'] = '(tag 3)'
+    d['apic'] = '(tag 4)'
+
+The square bracket operator is used above to add label definitions. It can
+be used to modify existing definitions, so long as the new new definition has the
+same type (region or locset):
+
+.. code-block:: python
+
+    import arbor
+
+    # A label dictionary that defines the label "dend" that defines a region.
+    d = arbor.label_dict({'dend': '(tag 3)')
+
+    # The definition of a label can be overwritten with a definition of the
+    # same type, in this case a region.
+    d['dend'] = '(join (tag 3) (tag 4))'
+
+    # However, a region can't be overwritten by a locset, or vice-versa.
+    d['dend'] = '(terminal)' # error: '(terminal)' defines a locset.
+
+    # New labels can be added to the dictionary.
+    d['soma'] = '(tag 1)'
+    d['axon'] = '(tag 2)'
+
+    # Square brackets can also be used to get a label's definition.
+    assert(d['soma'] == '(tag 1)')
+
+Expressions can refer to other regions and locsets in a label dictionary.
+In the example below, we define a region labeled *'tree'* that is the union
+of both the *'dend'* and *'apic'* regions.
+
+.. code-block:: python
+
+    import arbor
+
+    d = arbor.label_dict({
+            'soma': '(tag 1)',
+            'axon': '(tag 2)',
+            'dend': '(tag 3)',
+            'apic': '(tag 4)',
+            # equivalent to (join (tag 3) (tag 4))
+            'tree': '(join (region "dend") (region "apic"))'})
+
+The order that labels are defined does not matter, so an expression can refer to a
+label that has not yet been defined:
+
+.. code-block:: python
+
+    import arbor
+
+    d = arbor.label_dict()
+    # 'reg' refers 
+    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))'
+    # 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))'
+    # Which is the subtrees that start at the proximal locations of
+    # the region '(tag 3)'
+
+Cyclic dependencies are not permitted, as in the following example where
+two labels refer to one another:
+
+.. code-block:: python
+
+    import arbor
+
+    d = arbor.label_dict()
+    d['reg'] = '(distal_interval (locset "loc"))'
+    d['loc'] = '(proximal (region "reg"))'
+
+    # Error: 'reg' needs the definition of 'loc', which in turn needs the
+    # definition of 'reg'.
+
+.. note::
+    In the example above there will be no error when the label dictionary is defined.
+    Instead, there will be an error later when the label dictionary is applied to
+    a morphology, and the cyclic dependency is detected when concretising the locations
+    in the locsets and the cable segments in the regions.
+
+
+The type of an expression, locset or region, is inferred automatically when it is
+input into a label dictionary.
+Lists of the labels for regions and locsets are available as attributes:
+
+.. code-block:: python
+
+    import arbor
+
+    d = arbor.label_dict({
+            'soma': '(tag 1)',
+            'axon': '(tag 2)',
+            'dend': '(tag 3)',
+            'apic': '(tag 4)',
+            'site': '(location 2 0.5)',
+            'term': '(terminal)'})
+
diff --git a/doc/mechanisms.rst b/doc/mechanisms.rst
new file mode 100644
index 0000000000000000000000000000000000000000..aa2966a3e4aa776d0525c2426977ce4559299f92
--- /dev/null
+++ b/doc/mechanisms.rst
@@ -0,0 +1,479 @@
+.. _mechanisms:
+
+Mechanisms
+===========
+
+Mechanisms describe biophysical processes such as ion channels and synapses.
+Mechanisms are assigned to regions and locations on a cell morphology
+through a process that called :ref:`decoration <cablecell-decoration>`.
+Mechanisms are described using a dialect of the :ref:`NMODL <nmodl>` domain
+specific language that is similarly used in `NEURON <https://neuron.yale.edu/neuron/>`_.
+
+Mechanism Catalogues
+----------------------
+
+A *mechanism catalogue* is a collection of mechanisms that maintains:
+
+1. Collection of mechanism metadata indexed by name.
+2. A further hierarchy of *derived* mechanisms, that allow specialization of
+   global parameters, ion bindings, and implementations.
+3. A map for looking up a concrete mechanism implementation on a target hardware back end.
+
+A derived mechanism can be given a new name. Alternatively derived mechanisms can
+be created implitlitly.
+When a mechanism name of the form ``"mech/param=value,..."`` is :ref:`requested <mechanisms-name-note>`,
+if the mechanism of that name does not already exist in the catalogue, it will be
+implicitly derived from an existing mechanism ``"mech"``, with global parameters
+and ion bindings overridden by the supplied assignments that follow the slash.
+If the mechanism in question has a single ion dependence, then that ion name
+can be omitted in the assignments; ``"mech/oldion=newion"`` will make the same
+derived mechanism as simply ``"mech/newion"``.
+
+In additional being able to derive new mechanisms, catalogtues provide and interface
+for looking up a mechanism by name, and querying the following properties:
+
+* Global parameter: name, units and default value.
+* Range parameters: name, units and default value.
+* State variables: name, units and default value.
+* Ion dependency: for each ion whether it writes concentrations or reversal potential, and
+  whether the mechanism reads the reversal potential.
+
+Default Mechanisms
+''''''''''''''''''
+
+Arbor provides a default catalogue with the following mechanisms:
+
+* *pas*: Leaky current (:ref:`density mechanism <mechanisms-density>`).
+* *hh*:  Classic Hodgkin-Huxley dynamics (:ref:`density mechanism <mechanisms-density>`).
+* *nernst*: Calculate reversal potential for an ionic species using the Nernst equation (:ref:`reversal potential mechanism <mechanisms-revpot>`)
+* *expsyn*: Synapse with discontinuous change in conductance at an event followed by an exponential decay (:ref:`point mechanism <mechanisms-point>`).
+* *exp2syn*: Two state kinetic scheme synapse described by two time constants: rise and decay (:ref:`point mechanism <mechanisms-point>`).
+
+With the exception of *nernst*, these mechanisms are the same as those available in NEURON.
+
+Parameters
+''''''''''
+
+Mechanism behavior can be tuned using parameters and ion channel dependencies,
+prescribed in the NMODL description.
+Parameters and ion species are set initially before a simulation starts, and remain
+unchanged thereafter, for the duration of the simulation.
+There are two types of parameters that can be set by users:
+
+* *Global* parameters are a single scalar value that is the same everywhere a mechanism is defined.
+* *Range* parameters can vary spatially.
+
+Every mechanism is described with a *mechanism description*, a
+``(name, range_parameters)`` tuple, where ``name`` is a string,
+and ``range_parameters`` is an optional dictionary of key-value pairs
+that specifies values for range parameters.
+For example, consider a mechanism that models passive leaky dynamics with
+the following parameters:
+
+* *Name*: ``"passive"``.
+* *Global parameter*: reversal potential ``el``, default -65 mV.
+* *Range parameter*: conductance ``g``, default 0.001 S⋅cm⁻².
+
+The following example mechanism descriptions for our passive mechanism show that parameters and
+ion species dependencies only need to be specified when they differ from their defaults:
+
+* ``("passive")``: the passive mechanism with default parameters.
+* ``("passive/el=-80")``: derive a new passive mechanism with a non-default value for global parameter.
+* ``("passive", {"gl": 0.005})``: passive mechanism with a new a non-default range parameter value.
+* ``("passive/el=-80", {"gl": 0.005})``: derive a new passive mechanism that overrides both 
+
+Similarly to global parameters, ion species can be renamed in the mechanism name.
+This allows the use of generic mechanisms that can be adapted to a specific species
+during model instantiation.
+For example, the ``nernst`` mechanism in Arbor's default mechanism catalogue calculates
+the reversal potential of a generic ionic species ``x`` according to its internal
+and external concentrations and valence. To specialize ``nersnt`` for calcium name it
+``("nernst/x=ca")``, or if there is only one ions species in the mechanism the following
+shorthand ``("nernst/ca")`` can be used unambiguously.
+
+.. _mechanisms-name-note:
+
+.. note::
+    Global parameter values and ionic dependencies are the same for each instance of
+    a mechanism, so when these are redifeind a new mechanism is created, derived from
+    the parent mechanism.
+    For this reason, new global parameter values and ion renaming are part of the name of
+    the new mechanism, or a mechanism with a new unique name must be defined.
+
+
+Mechanism Types
+---------------
+
+There are two broad categories of mechanism, density mechanisms and
+point mechanisms, and a third special density mechanism for
+computing ionic reversal potentials.
+
+.. _mechanisms-density:
+
+Density mechanisms
+''''''''''''''''''''''
+
+Density mechanisms are :ref:`NMODL mechanisms <nmodl>`
+which describe biophysical processes that are distributed in space, but whose behaviour
+is defined purely by the state of the cell and the process at any given point.
+
+Density mechanisms are commonly used to describe ion channel dynamics,
+for example the ``hh`` and ``pas`` mechanisms provided by NEURON and Arbor,
+which model classic Hodgkin-Huxley and passive leaky currents respectively.
+
+.. _mechanisms-revpot:
+
+Ion reversal potential mechanisms
+'''''''''''''''''''''''''''''''''
+
+These mechanisms, which describe ionic reversal potential
+behaviour, can be specified for cells or the whole model.
+
+The reversal potential of an ion species is calculated by an
+optional *reversal potential mechanism*.
+If no such mechanism is specified for an ion species, the initial
+reversal potential values are maintained for the course of a simulation.
+Otherwise, the mechanism does the work.
+
+Reversal potential mechanisms are density mechanisms subject to some strict restrictions.
+Specifically, a reversal potential mechanism described in NMODL:
+
+* May not maintain any state variables.
+* Can only write to the reversal potential (``eX``) value of the ion species.
+* Can not be a :ref:`point mechanism <mechanisms-point>`.
+
+Essentially, reversal potential mechanisms must be pure functions of cellular
+and ionic state.
+
+.. note::
+    Arbor imposes greater restrictions on mechanisms that update ionic reversal potentials
+    than NEURON. Doing so simplifies reasoning about interactions between
+    mechanisms that share ionic species, by virtue of having one mechanism, and one
+    mechanism only, that calculates reversal potentials according to concentrations
+    that the other mechanisms use and modify.
+
+.. _mechanisms-point:
+
+Point mechanisms
+'''''''''''''''''''''''''''''''''
+
+*Point mechanisms*, which are associated with connection end points on a
+cable cell, are placed at discrete locations on the cell.
+Unlike density mechanisms, whose behaviour is defined purely by the state of the cell and the process,
+their behavior is additionally governed by the timing and weight of events delivered
+via incoming connections.
+
+Python API
+----------
+
+Mechanism Catalogues
+''''''''''''''''''''
+
+.. py:class:: catalogue
+
+    A *mechanism catalogue* is a collection of mechanisms that maintains:
+
+    1. Collection of mechanism metadata indexed by name.
+    2. A further hierarchy of *derived* mechanisms, that allow specialization of
+       global parameters, ion bindings, and implementations.
+
+    .. py:method:: has(name)
+
+        Test if mechanism with *name* is in the catalogue.
+
+        :param name: name of mechanism.
+        :type name: str
+        :return: bool
+
+    .. py:method:: is_derived(name)
+
+        Is *name* a derived mechanism or can it be implicitly derived?
+
+        :param name: name of mechanism.
+        :type name: str
+        :return: bool
+
+    .. py:method:: __getitem__(name)
+
+        Look up mechanism meta data with *name*.
+
+        .. code-block:: Python
+
+            import arbor
+
+            cat = arbor.default_catalogue()
+
+            # Print default value and units for gnabar parameter of hh.
+            print(cat['hh'].parameters['gnabar'])
+
+        :param name: name of mechanism.
+        :type name: str
+        :return: mechanism metadata
+        :rtype: :class:`mechanism_info`
+
+    .. py:method:: derive(name, parent, globals={}, ions={})
+
+        Derive a new mechanism with *name* from the mechanism *parent*.
+
+        If no parameters or ion renaming are specified with *globals* or *ions*,
+        the method will attempt to implicitly derive a new mechanism from parent by parsing global and
+        ions from the parent string.
+
+        .. code-block:: Python
+
+            import arbor
+
+            cat = arbor.default_catalogue()
+
+            # Use the value of the Faraday constant as published by CODATA in 1986,
+            # and bind to pottasium ion species.
+            cat.derive('krev',  'nernst', globals={'F': 96485.309}, ions={'x': 'k'})
+
+            # Derive a reversal potential mechanism for sodium from the one we defined
+            # for potasium, which will inherit the redefined Faraday constant.
+            cat.derive('narev', 'krev', ions={'k': 'na'})
+
+            # Alternatively, we can derive a mechanism with global parameters and ion renaming
+            # specified in the parent name string.
+            cat.derive('krev_imp', 'nernst/F=96485.309,k')
+            cat.derive('carev', 'krev_imp/ca')
+
+        :param name: name of new derived mechanism.
+        :type name: str
+        :param parent: name of parent mechanism.
+        :type parent: str
+        :param globals: a dictionary mapping global parameter names to their values, if any.
+        :type globals: dict[str, float]
+        :param ions: a dictionary renaming ion species, if any.
+        :type ions: dict[str, str]
+
+.. py:class:: mechanism_info
+
+    Meta data about the fields and ion dependencies of a mechanism.
+    The data is presented as read-only attributes.
+
+    .. code-block:: Python
+
+        import arbor
+        cat = arbor.default_catalogue()
+
+        # Get mechanism_info for the 'expsyn' mechanism.
+        mech = cat['expsyn']
+
+        # Query the mechanism_info for information about parameters.
+
+        print(mech.parameters.keys())
+        # dict_keys(['e', 'tau'])
+
+        print(mech.parameters['tau'].units)
+        # 'ms'
+
+        print(mech.parameters['tau'].default)
+        # 2.0
+
+    .. py:attribute:: globals
+        :type: dict[str, mechanism_field]
+
+        Global fields have one value common to an instance of a mechanism, are constant in time and set at instantiation.
+
+    .. py:attribute:: parameters
+        :type: dict[str, mechanism_field]
+
+        Parameter fields may vary across the extent of a mechanism, but are constant in time and set at instantiation.
+
+    .. py:attribute:: state
+        :type: dict[str, mechanism_field]
+
+        State fields vary in time and across the extent of a mechanism, and potentially can be sampled at run-time.
+
+    .. py:attribute:: ions
+        :type: dict[str, ion_dependency]
+
+        Ion dependencies.
+
+    .. py:attribute:: linear
+        :type: bool
+
+        True if a synapse mechanism has linear current contributions so that multiple instances on the same compartment can be coalesced.
+
+
+.. py:class:: ion_dependency
+
+    Meta data about a mechanism's dependence on an ion species,
+    presented as read-only attributes.
+
+    .. code-block:: Python
+
+        import arbor
+        cat = arbor.default_catalogue()
+
+        # Get ion_dependency for the 'hh' mechanism.
+        ions = cat['hh'].ions
+
+        # Query the ion_dependency.
+
+        print(ions.keys())
+        # dict_keys(['k', 'na'])
+
+        print(ions['k'].write_rev_pot)
+        # False
+
+        print(ions['k'].read_rev_pot)
+        # True
+
+    .. py:attribute:: write_int_con
+        :type: bool
+
+        If the mechanism contributes to the internal concentration of the ion species.
+
+    .. py:attribute:: write_ext_con
+        :type: bool
+
+        If the mechanism contributes to the external concentration of the ion species.
+
+    .. py:attribute:: write_rev_pot
+        :type: bool
+
+        If the mechanism calculates the reversal potential of the ion species.
+
+    .. py:attribute:: read_rev_pot
+        :type: bool
+
+        If the mechanism depends on the reversal potential of the ion species.
+
+
+.. py:class:: mechanism_field
+
+    Meta data about a specific field of a mechanism, presented as read-only attributes.
+
+    .. py:attribute:: units
+        :type: string
+
+        The units of the field.
+
+    .. py:attribute:: default
+        :type: float
+
+        The default value of the field.
+
+    .. py:attribute:: min
+        :type: float
+
+        The minimum permissible value of the field.
+
+    .. py:attribute:: max
+        :type: float
+
+        The maximum permissible value of the field.
+
+The :py:class:`mechanism_info` type above presents read-only information about a mechanism that is available in a catalogue.
+
+When :ref:`decorating <cablecell-decoration>` a cable cell, we use a :py:class:`mechanism` type to describe a
+mechanism that is to be painted or placed on the cable cell.
+
+.. py:class:: mechanism
+
+    Mechanisms describe physical processes, distributed over the membrane of the cell.
+    *Density mechanisms* are associated with regions of the cell, whose dynamics are
+    a function of the cell state and their own state where they are present.
+    *Point mechanisms* are defined at discrete locations on the cell, which receive
+    events from the network.
+    A third, specific type of density mechanism, which describes ionic reversal potential
+    behaviour, can be specified for cells or the whole model.
+
+    The :class:`mechanism` type is a simple wrapper around a mechanism
+    :attr:`mechanism.name` and a dictionary of named parameters.
+
+    Mechanisms have two types of parameters:
+
+    * global parameters: a scalar value that is the same for all instances
+      of a mechanism.
+    * range parameters: the value of range parameters is defined for each instance
+      of the mechanism on a cell. For density mechanisms, this means one value for
+      each compartment on which it is present.
+
+    The method for setting a parameter depends on its type.
+    If global parameters change, we are effectively defining a new type
+    of mechanism, so global parameter information is encoded in the
+    name.
+    Range parameters are set using a dictionary of name-value pairs.
+
+    .. code-block:: Python
+
+        import arbor
+
+        # hh dynamics with default parameters.
+        hh = arbor.mechanism('hh')
+
+        # A passive leaky channel with custom parameters
+        pas = arbor.mechanism('pas', {'e': -55, 'gl': 0.02})
+
+        # Reversal potential using Nernst equation with GLOBAL parameter values
+        # for Faraday's constant and the target ion species, set with a '/' followed
+        # by comma-separated list of parameter after the base mechanism name.
+        rev = arbor.mechanism('nernst/F=96485,x=ca')
+
+    .. method:: mechanism(name, params)
+
+        constructor for mechanism with *name* and range parameter overrides *params*,
+        for example: ``arbor.mechanism(name='pas', params={'g': 0.01})``.
+
+        :param name: name of mechanism.
+        :type name: str
+        :param params: A dictionary of parameter values, with parameter name as key.
+        :type params: dict[str, double]
+
+    .. method:: mechanism(name)
+        :noindex:
+
+        constructor for mechanism.
+        The *name* can be either the name of a mechanism in the catalogue,
+        e.g.  ``arbor.mechanism('pas')``, or an implicitly derived mechanism,
+        e.g. ``arbor.mechanism('nernst/k')``.
+
+    .. method:: set(name, value)
+
+        Set new value for a parameter.
+
+        :param name: name of the parameter.
+        :type name: str
+        :param value: value of the parameter.
+        :type value: float
+
+    .. py:attribute:: name
+        :type: str
+
+        The name of the mechanism.
+
+    .. py:attribute:: values
+        :type: dict
+
+        A dictionary of key-value pairs for the parameters.
+
+    .. code-block:: Python
+
+        import arbor
+
+        # Create pas mechanism with default parameter values (set in NOMDL file).
+        m1 = arbor.mechanism('passive')
+
+        # Create default mechainsm with custom conductance (range).
+        m2 = arbor.mechanism('passive', {'g', 0.1})
+
+        # Create a new pas mechanism with that changes reversal potential (global).
+        m3 = arbor.mechanism('passive/el=-45')
+
+        # Create an instance of the same mechanism, that also sets conductance (range).
+        m4 = arbor.mechanism('passive/el=-45', {'g', 0.1})
+
+        # This is an equivalent to m4, using set method to specify range parameters.
+        m5 = arbor.mechanism('passive/el=-45')
+        m5.set('g', 0.1)
+
+        # 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'.
+
diff --git a/doc/morphology.rst b/doc/morphology.rst
new file mode 100644
index 0000000000000000000000000000000000000000..a0338c23b6d31df8ee3fe879c7133d6265e7de9a
--- /dev/null
+++ b/doc/morphology.rst
@@ -0,0 +1,855 @@
+.. _morphology:
+
+Morphology
+==========
+
+A cell's *morphology* describes both its geometry and branching structure.
+Morphologies in Arbor are modeled as a set of one dimensional cables of variable radius,
+joined together to form a tree.
+
+The building blocks of morpholgies tree are points and segments.
+A *point* is a three-dimensional location and a radius, used to mark the centre and radius
+of the cable.
+
+.. csv-table::
+   :widths: 10, 10, 30
+
+   **Field**,   **Type**, **Description**
+   ``x``,       real, x coordinate of centre of cable (μm).
+   ``y``,       real, y coordinate of centre of cable (μm).
+   ``z``,       real, z coordinate of centre of cable (μm).
+   ``radius``,  real, cross sectional radius of cable (μm).
+
+
+A *segment* is a frustum (cylinder or truncated cone), with the centre and radius at each
+end defined by a pair of points.
+
+.. csv-table::
+   :widths: 10, 10, 30
+
+   **Field**,   **Type**, **Description**
+   ``prox``,       :py:class:`point <arbor.mpoint>`,   the center and radius of the proximal end.
+   ``dist``,       :py:class:`point <arbor.mpoint>`,   the center and radius of the distal end.
+   ``tag``,        integer, tag meta-data.
+
+.. _morph-tag-definition:
+
+A *tag* is an integer label on every segment, which can be used to define disjoint
+regions on cells.
+The meaning of tag values are not fixed in Arbor, however we typically use tag values that correspond
+to SWC `structure identifiers <http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html>`_.
+
+.. _morph-segment_tree:
+
+Segment Trees
+--------------
+
+A *segment tree* describes a morphology as a set of segments and their connections,
+designed to support both the diverse descriptions
+of cell morphologies (e.g. SWC, NeuroLicida, NeuroML), and tools that
+iteratively construct cell morphologies (e.g. L-system generators, interactive cell-builders).
+
+Segment trees comprise a sequence of segments starting from at lease one *root* segment,
+together with a parent-child adjacency relationship where a child segment is
+distal to its parent.
+Branches in the tree occur where a segment has more than one child.
+Furthermore, a segment can not have more than one parent.
+In this manner, neuron morphologies are modeled as a *tree*, where cables that
+represent dendrites and axons can branch, but branches can not rejoin.
+
+.. _morph-segment-definitions:
+
+The following definitions are used to refer to segments in a segment tree:
+
+* *root*: segments at the root or start of the tree. A non-empty tree must have at least one root segment,
+  and the first segment will always be a root.
+
+* *parent*: Each segment has one parent, except for root segments which have :data:`mnpos <arbor.mnpos>` as their parent.
+
+  * The id of a segment is always greater than the id of its parent.
+  * The ids of segments on the same unbranched sequence of segments do not need to be contiguous.
+
+* *child*: A segment's children are the segments that have the segment as their parent.
+* *terminal*: A segment with no children. Terminals lie at the end of dendritic trees or axons.
+* *fork*: A segment with more than one child. The distal end of a fork segment are *fork points*,
+  where a cable splits into two or more branches.
+
+  * Arbor allows more than two branches at a fork point.
+
+The following segment tree models a soma as a cylinder, a branching dendritic tree and
+an axon with an axonal hillock. The segments are colored according to their tag, which
+in this case are SWC structure identifiers: tag 1 colored pink for soma;
+tag 2 colored grey for axon; tag 3 colored blue for basal dendrites.
+
+.. _morph-label-seg-fig:
+
+.. figure:: gen-images/label_seg.svg
+  :width: 600
+  :align: center
+
+  Example Python code to generate this morphology is in the :class:`segment_tree<arbor.segment_tree>`
+  documentation :ref:`below <morph-label-seg-code>`.
+
+We can apply the following labels to the segments:
+
+* The tree is composed of 11 segments (1 soma, 2 axon, 8 dendrite).
+* The proximal ends of segments 0 and 9 (the soma and axon hillock respectively) are attached to the root of the tree.
+* Segment 2 is a fork, with segments 3 and 5 as children.
+* Segment 5 is a fork, with segments 6 and 7 as children.
+* There is also a fork at the root, whith segments 0 and 9 as children.
+* Segments 4, 6, 8 and 10 are terminal segments.
+
+In the example above there are no gaps between segments, however
+it is possible for segments to be detached, where the proximal end of a segment is not coincident
+with the distal end of its parent. The following morphology has gaps between the start of the
+axon and dendritic tree and the soma segment to which they attach.
+
+.. _morph-detached-seg-fig:
+
+.. figure:: gen-images/detached_seg.svg
+  :width: 600
+  :align: center
+
+.. note::
+    In Arbor, segments are always treated as though they are connected directly
+    to their parents, regardless of whether ends where they attached are collocated.
+
+    Gaps are frequently the result of simplifying the soma,
+    whereby the complex geometry of a soma is represented using a cylinder or sphere
+    (spheres are represented by a cylinder with length and diameter equal to that of
+    the sphere in simulation tools like Arbor and NEURON).
+
+    A gap between a cylindrical soma and segments attached to it does not mean
+    that the segmentation is invalid.
+    To illustrate why this can occur, consider a potato-shaped soma modeled with a
+    cylinder of the same surface area.
+    If the cell description places the first segment of a dendritic tree where it attaches to
+    the "potato soma", it is unlikely to be collocated with an end of the simplified soma.
+    The cell model will correctly represent the location and dimension of the dendritic tree,
+    while preserving the soma surface area with a simplified cylindrical model.
+
+Because Arbor supports tapered segments (where radius varies linearly along a segment) it is possible to
+represent more complex soma shapes using multiple segments, for example the segmentation below
+uses 4 segments to model the soma.
+
+.. _morph-stacked-seg-fig:
+
+.. figure:: gen-images/stacked_seg.svg
+  :width: 600
+  :align: center
+
+.. _morph-morphology:
+
+Morphology
+----------
+
+A *morphology* describes the geometry of a cell as unbranched cables with variable radius
+, and their associated tree structure. 
+Every segment tree can be used to generate a unique morphology, which derives and enumerates
+*branches* from the segments.
+The branches of a morphology are unbranched cables, composed of one or more segments, where:
+
+  * the first (proximal) segment of the branch is either a root or the child of fork segment;
+  * the last (distal) segment of the branch is either a fork or terminal segment;
+  * branches are enumerated according to the ids of their proximal segments in the segment trree.
+
+When constructed in this manner, the following statements are true for the branches and
+their enumeration:
+
+  * Because a branch must have root, fork or terminal ends, a branch can not be sub-divided
+    into two or more branches, and hence there is only one possible set of branches that
+    can be derived from a segment tree.
+  * Because branches are enumerated according to the id of their proximal segments,
+    there is only one branch enumeration representation for a segment tree.
+  * However, it is possible for two topologically equivalent morphologies to be
+    derived from different segment trees (e.g. two trees with the same segments, however
+    different valid segment enumerations), and potentially have different branch numbers.
+  * Every valid segment tree can be used to construct a valid morphology.
+
+.. Note::
+
+    Because two topologically-equivalent morphologies may have different segment and
+    branch numbering, it is important that model descriptions should avoid refering to
+    branches or segments by id.
+    This should be relaxed only in well-understood situations, for example when working with
+    models that always represent to soma with a single segment at the root of the tree,
+    which will always have segment id 0.
+
+To illustrate branch generation, consider the first segment tree example on this page,
+which is illustrated along with its branches below.
+
+.. _morph-label-morph-fig:
+
+.. figure:: gen-images/label_morph.svg
+  :width: 800
+  :align: center
+
+  The code used to generate this morphology is in the :class:`segment_tree<arbor.segment_tree>`
+  documentation :ref:`below <morph-label-seg-code>`.
+
+The first branch contains the soma and the first two segments of the dendritic tree.
+There are four more branches in the dendritic tree, and one representing the two
+segments of the axon.
+
+Note, that though it is possible to create an unbranched sequence of segments composed
+of the axon, soma and first two segements in the dendritic tree, this sequence is decomposed
+as two branches because segments 0 (soma) and 9 (first segment in axon) are at the
+root of the tree.
+
+Similarly to segments, the branches in a morphology have a parent child relationship.
+Every branch has one parent, with branches at the root of the tree having the placeholder
+parent index :data:`mnpos <arbor.mnpos>`. Segments can have any non-negative number of children,
+however by nature of their construction, no branch can have only one child: a branch has
+either no children, or two or more children.
+The parent-child information and segments for the morphology are summarised:
+
+.. csv-table::
+   :widths: 10, 10, 10, 10
+
+   **Branch**, **Parent**, **Children**, **Segments**
+   0,          ``mnpos``,  "[1, 2]",       "[0, 1, 2]"
+   1,          0,          "[]",           "[3, 4]"
+   2,          0,          "[3, 4]",       "[5]"
+   3,          2,          "[]",           "[6]"
+   4,          2,          "[]",           "[7, 8]"
+   5,          ``mnpos``,  "[]",           "[9, 10]"
+
+Gaps between segments do not influence branch creation, hence branches
+can contain gaps between segments. Take the example of a morphology with
+a gap between the soma and the axona and dendritic trees:
+
+.. figure:: gen-images/detached_morph.svg
+  :width: 800
+  :align: center
+
+The soma is part of branch 0, despite the gap:
+
+.. csv-table::
+   :widths: 10, 10, 10, 10
+
+   **Branch**, **Parent**, **Children**, **Segments**
+   0,          ``mnpos``,  "[1, 2]",       "[0, 1, 2]"
+   1,          0,          "[]",           "[3, 4]"
+   2,          0,          "[3, 4]",       "[5]"
+   3,          2,          "[]",           "[6]"
+   4,          2,          "[]",           "[7, 8]"
+   5,          ``mnpos``,  "[]",           "[9]"
+
+Tag information is not used when creating branches, so that a branch can
+contain segments with different tags, which in our examples gives branches
+that contain both soma and dendrite segments. For example, when building the
+soma from multiple segments:
+
+.. figure:: gen-images/stacked_morph.svg
+  :width: 800
+  :align: center
+
+The morphology has the same number of branches as the other examples, with
+multiple soma and dendrite segments in branch 0.
+
+.. csv-table::
+   :widths: 10, 10, 10, 10
+
+   **Branch**, **Parent**, **Children**, **Segments**
+   0,          ``mnpos``,  "[1, 2]",       "[0, 1, 2, 3, 4, 5]"
+   1,          0,          "[]",           "[6, 7]"
+   2,          0,          "[3, 4]",       "[8]"
+   3,          2,          "[]",           "[9]"
+   4,          2,          "[]",           "[10, 11]"
+   5,          ``mnpos``,  "[]",           "[12, 13]"
+
+.. Note::
+    Arbor provides a consistent representation of morphologies with no
+    special cases for concepts like magical soma branches, in order to
+    build reproducable and consistent model descriptions.
+
+    Users of NEURON who are used to creating a separate soma section
+    that is always the first section in a morphology should not
+    worry that the soma is not treated as a special branch
+    in the examples above.
+
+    The soma in the examples above can be referred to in later model
+    building phases, for example when describing the distribution of
+    ion channels, by using refering to all parts of the cell with
+    :ref:`tag 1 <labels-expressions>`.
+
+
+Examples
+~~~~~~~~~~~~~~~
+
+Here we present a series of morphology examples of increasing complexity.
+The examples use the Python API are two-dimensional, with the z-dimension set to zero.
+
+.. _morph-tree1:
+
+Example 1: Spherical cell
+""""""""""""""""""""""""""""""
+
+A simple model of a cell as a sphere can be modeled using a cylinder with length
+and diameter equal to the diameter of the sphere, which will have the same
+surface area (disregarding the area of the cylinder's circular ends).
+
+Here a cylinder of length and diameter 5 μm is used to represent a *spherical cell*
+with a radius of 2 μm, centered at the origin.
+
+.. code:: Python
+
+    tree = arbor.segment_tree()
+    tree.append(mnpos, mpoint(-2, 0, 0, 2), mpoint(2, 0, 0, 2), tag=1)
+    morph = arbor.morphology(tree)
+
+.. figure:: gen-images/sphere_morph.svg
+  :width: 400
+  :align: center
+
+  The morphology is a single cylinder segment (left) that forms branch 0 (right).
+
+.. _morph-tree2:
+
+Example 2: Unbranched cable
+""""""""""""""""""""""""""""""
+
+Consider a cable of length 10 μm, with a radius that tapers from 1 μm to 0.5 μm
+at the proximal and distal ends respectively.
+This can be described using a single segment.
+
+.. code:: Python
+
+    tree = arbor.segment_tree()
+    tree.append(mnpos, mpoint(0, 0, 0, 1), mpoint(10, 0, 0, 0.5), tag=3)
+    morph = arbor.morphology(tree)
+
+.. figure:: gen-images/branch_morph1.svg
+  :width: 600
+  :align: center
+
+  A tapered cable with one cable segment (left), generates a morphology with one branch (right).
+
+The radius of a cable segment varies lineary between its end points. To define an unbranched cable
+with irregular radius and "squiggly" shape, use multiple segments to build a piecewise linear reconstruction
+of the cable geometry.
+This example starts and ends at the same locations as the previous, however it is constructed from 4
+distinct cable segments:
+
+.. code:: Python
+
+    tree = arbor.segment_tree()
+    tree.append(mnpos, mpoint( 0.0,  0.0,  0.0, 1.0), mpoint( 3.0,  0.2,  0.0, 0.8), tag=1)
+    tree.append(0,     mpoint( 3.0,  0.2,  0.0, 0.8), mpoint( 5.0, -0.1,  0.0, 0.7), tag=2)
+    tree.append(1,     mpoint( 5.0, -0.1,  0.0, 0.7), mpoint( 8.0,  0.0,  0.0, 0.6), tag=2)
+    tree.append(2,     mpoint( 8.0,  0.0,  0.0, 0.6), mpoint(10.0,  0.0,  0.0, 0.5), tag=3)
+    morph = arbor.morphology(tree)
+
+    morph = arbor.morphology(tree)
+
+.. figure:: gen-images/branch_morph2.svg
+  :width: 600
+  :align: center
+
+  The morphology is an ubranched cable comprised of 4 cable segments,
+  colored according to their tags: tag 1 red; tag 2 gree; tag 3 blue (left).
+  The four segments form one branch (right).
+
+Gaps are possible between two segments. The example below inserts a 1 μm gap between the second
+and third segments of the previous morphology. Note that Arbor will ignore the gap, effectively
+joining the segments together, such that the morphology with the gap is the same as that without.
+
+.. code:: Python
+
+    tree = arbor.segment_tree()
+    tree.append(mnpos, mpoint( 0.0,  0.0,  0.0, 1.0), mpoint(3.0,  0.2,  0.0, 0.8), tag=1)
+    tree.append(0,     mpoint( 3.0,  0.2,  0.0, 0.8), mpoint(5.0, -0.1,  0.0, 0.7), tag=2)
+    tree.append(1,     mpoint( 7.0, -0.1,  0.0, 0.7), mpoint(10.0, 0.0,  0.0, 0.6), tag=2)
+    tree.append(2,     mpoint(10.0,  0.0,  0.0, 0.6), mpoint(12.0, 0.0,  0.0, 0.5), tag=3)
+    morph = arbor.morphology(tree)
+
+.. figure:: gen-images/branch_morph3.svg
+  :width: 600
+  :align: center
+
+  There is a gap between segment 1 and segment 2 (left), and there is a single branch (right).
+
+The radius of a cable is piecewise linear, with discontinuities permited at the
+interface between segments.
+The next example adds a discontinuity to the previous example between segments
+3 and 4, where the radius changes from 0.5 μm to 0.3 μm:
+
+.. code:: Python
+
+    tree = arbor.segment_tree()
+    tree.append(mnpos, mpoint( 0.0,  0.0,  0.0, 1.0), mpoint( 3.0,  0.2,  0.0, 0.8), tag=1)
+    tree.append(0,     mpoint( 3.0,  0.2,  0.0, 0.8), mpoint( 5.0, -0.1,  0.0, 0.7), tag=2)
+    tree.append(1,     mpoint( 5.0, -0.1,  0.0, 0.7), mpoint( 8.0,  0.0,  0.0, 0.5), tag=2)
+    tree.append(2,     mpoint( 8.0,  0.0,  0.0, 0.3), mpoint(10.0,  0.0,  0.0, 0.5), tag=3)
+    morph = arbor.morphology(tree)
+
+.. figure:: gen-images/branch_morph4.svg
+  :width: 600
+  :align: center
+
+  The resulting morphology has a step discontinuity in radius.
+
+.. _morph-example4:
+
+Example 3: Y-shaped cell
+""""""""""""""""""""""""""""""
+
+The simplest branching morphology is a cable that bifurcates into two branches,
+which we will call a *y-shaped cell*.
+In the example below, the first branch of the tree is a cable of length 10 μm with a
+a radius that tapers from 1 μm to 0.5 μm.
+The two child branches are attached to the end of the first branch, and taper from from 0.5 μ m
+to 0.2 μm.
+
+Note that only the distal point is required to describe the child segments,
+because the proximal end of each child segment has the same location and
+radius as the distal end of the parent.
+
+.. code:: Python
+
+    tree = arbor.segment_tree()
+    tree.append(mnpos, mpoint( 0.0, 0.0, 0.0, 1.0), mpoint(10.0, 0.0, 0.0, 0.5), tag= 3)
+    tree.append(0,     mpoint(15.0, 3.0, 0.0, 0.2), tag= 3)
+    tree.append(0,     mpoint(15.0,-3.0, 0.0, 0.2), tag= 3)
+    morph = arbor.morphology(tree)
+
+.. figure:: gen-images/yshaped_morph.svg
+  :width: 800
+  :align: center
+
+Example 4: Soma with branches
+""""""""""""""""""""""""""""""
+
+Now let's look at cell with a simple dendritic tree attached to a spherical soma.
+The spherical soma of radius 3 μm is modelled with a cylinder with length and
+diameter equal to 6 μm, which has the same surface area as the sphere.
+
+.. code:: Python
+
+    tree = arbor.segment_tree()
+    tree.append(mnpos, mpoint(-3.0, 0.0, 0.0, 3.0), mpoint( 3.0, 0.0, 0.0, 3.0), tag=1)
+    tree.append(0, mpoint( 4.0, -1.0,  0.0, 0.6), mpoint(10.0,  -2.0,  0.0, 0.5), tag=3)
+    tree.append(1, mpoint(15.0, -1.0,  0.0, 0.5), tag=3)
+    tree.append(2, mpoint(18.0, -5.0,  0.0, 0.3), tag=3)
+    tree.append(2, mpoint(20.0,  2.0,  0.0, 0.3), tag=3)
+    morph = arbor.morphology(tree)
+
+
+.. figure:: gen-images/ysoma_morph1.svg
+  :width: 900
+  :align: center
+
+  Note that branch 0 (right) is composed of segments 0, 1, and 2 (left).
+
+The soma is the first segment, labeled with tag 1. The dendritic tree is a simple
+y-shaped tree composed of 4 segments, each labeled with tag 3.
+The first branch is composed of 3 segments: the soma segment and the first two segments
+in the dendritic tree because the segments have parent child ordering and no fork points.
+
+.. note::
+    The first branch is derived directly from the topological relationship between the segments,
+    and no special treatment is given to the soma.
+    There is no need to treat segments with different tags (e.g. tags that we might associate
+    with soma, axon, basal dendrite and apical dendrite) when defining geometric primitives like
+    segments and branches, because they can later be referenced later using
+    :ref:`region expressions <labels-expressions>`.
+
+Now we can attach another dendrite and an axon to the soma, to make a total of three cables
+attached to the soma (two dendrites and an axon).
+The dendrites are attached to the distal end of the soma (segment 0), so they have the
+0 as their parent.
+The axon is attached to the proximal end of the soma, which is at the root of the tree,
+so it has :data:`mnpos` as its parent.
+There are 7 branches generated from 10 segments, and soma segment is its own branch,
+because it has two children: the dendrites attached to its distal end.
+
+.. figure:: gen-images/ysoma_morph2.svg
+  :width: 900
+  :align: center
+
+
+.. note::
+    The discretization process, which converts segments and branches into compartments,
+    will ignore gaps between segments in the input. The cell below, in which the dendrites
+    and axon have been translated to remove any gaps, is equivalent to the previous example
+    for the back end simulator.
+
+    Note that the dendrites are children of the soma segment, so they are coincident with
+    the distal end of the soma, and the axon is translated to the proximal end of the
+    soma segment because both it and the soma have :py:data:`mnpos <arbor.mnpos>` as a parent.
+    More generally, segments at the root of the tree are connected electrically at their
+    proximal ends.
+
+    .. figure:: gen-images/ysoma_morph3.svg
+      :width: 900
+      :align: center
+
+Python API
+----------
+
+.. currentmodule:: arbor
+
+.. data:: mnpos
+    :type: int
+
+    Value used to indicate "no parent" in :class:`segment_tree` and :class:`morphology`
+    trees of segments and branches respectively.
+
+    .. code-block:: python
+
+        import arbor
+
+        tree = arbor.segment_tree()
+
+        # mnpos can be used to explicitly specify that a segment
+        # is at the root of the tree. More than one segment can
+        # be at the root, and they will all be joined electrically
+        # at their proximal ends.
+        tree.append(parent=arbor.mnpos, # attach segment to root.
+                    prox=arbor.mpoint(0, 0,-5, 5),
+                    dist=arbor.mpoint(0, 0, 5, 5),
+                    tag=1)
+        tree.append(parent=0,
+                    prox=arbor.mpoint(0, 0, 5, 0.5),
+                    dist=arbor.mpoint(0, 0,50, 0.2),
+                    tag=3)
+
+        # mnpos can also be used when querying a sample_tree or morphology,
+        # for example the following snippet that finds all branches in the
+        # morphology that are attached to the root of the morphology.
+        m = arbor.morphology(tree)
+        base_branches = [i for i in range(m.num_branches)
+                            if m.branch_parent(i) == arbor.mnpos]
+
+        print(base_branches)
+
+
+
+.. class:: location
+
+    A location on :attr:`branch`, where :attr:`pos`, in the range ``0 ≤ pos ≤ 1``,
+    gives the relative position
+    between the proximal and distal ends of the branch. The position is in terms
+    of branch path length, so for example, on a branch of path length 100 μm ``pos=0.2``
+    corresponds to 20 μm and 80 μm from the proximal and distal ends of the
+    branch respectively.
+
+    .. function:: location(branch, pos)
+
+        Constructor.
+
+    .. attribute:: branch
+        :type: int
+
+        The branch id of the location.
+
+    .. attribute:: pos
+        :type: float
+
+        The relative position of the location on the branch.
+
+.. class:: cable
+
+    An unbranched cable that is a subset of a branch.
+    The values of ``0 ≤ prox ≤ dist ≤ 1`` are the relative position
+    of the cable's end points on the branch, in terms
+    of branch path length. For example, on a branch of path length 100 μm, the values
+    :attr:`prox` =0.2, :attr:`dist` =0.8 describe a cable that starts and
+    ends 20 μm and 80 μm along the branch respectively.
+
+    .. function:: cable(branch, prox, dist)
+
+        Constructor.
+
+    .. attribute:: branch
+        :type: int
+
+        The branch id of the cable.
+
+    .. attribute:: prox
+        :type: float
+
+        The relative position of the proximal end of the cable on the branch.
+
+    .. attribute:: dist
+        :type: float
+
+        The relative position of the distal end of the cable on the branch.
+
+.. class:: mpoint
+
+    A location of a cell morphology at a fixed location in space. Describes the location
+    of the as three-dimensional coordinates (:attr:`x`, :attr:`y`, :attr:`z`) and
+    the :attr:`radius` of the cable.
+
+    .. attribute:: x
+        :type: real
+
+        X coordinate (μm)
+
+    .. attribute:: y
+        :type: real
+
+        Y coordinate (μm)
+
+    .. attribute:: z
+        :type: real
+
+        x coordinate (μm)
+
+    .. attribute:: radius
+        :type: real
+
+        Radius of the cable (μm)
+
+.. class:: segment
+
+    .. attribute:: prox
+        :type: mpoint
+
+        The location and radius at the proximal end of the segment.
+
+    .. attribute:: dist
+        :type: mpoint
+
+        The location and radius at the distal end of the segment.
+
+    .. attribute:: tag
+        :type: int
+
+        Integer tag meta-data associated with the segment.
+        Typically the tag would correspond to the SWC structure identifier:
+        soma=1, axon=2, dendrite=3, apical dendrite=4, however arbitrary
+        tags, including zero and negative values, can be used.
+
+.. class:: segment_tree
+
+    A segment tree is a description of a the segments and their connections
+    Segment trees comprise a sequence of segments starting from at lease one *root* segment,
+    together with a parent-child adjacency relationship where a child segment is
+    distal to its parent.
+    Branches in the tree occur where a segment has more than one child.
+    Furthermore, a segment can not have more than one parent.
+    In this manner, neuron morphologies are modeled as a *tree*, where cables that
+    represent dendrites and axons can branch, but branches can not rejoin.
+    A segment tree is a segment-based description of a cell's morphology.
+
+    .. function:: segment_tree()
+
+        Construct an empty segment tree.
+
+    The tree is constructed by *appending* segments to the tree.
+    Segments are numbered starting at 0 in the order that they are added,
+    with the first segment getting id 0, the second segment id 1, and so forth.
+
+    A segment can not be added before its parent, hence the first segment
+    is always at the root. In this manner, a segment tree is
+    always guarenteed to be in a correct state, with consistent parent-child
+    indexing, and with *n* segments numbered from *0* to *n-1*.
+
+    To illustrate how a segment tree is constructed by appending segments,
+    take the segment tree used in the :ref:`documentation above <morph-label-seg-fig>`.
+
+
+    .. figure:: gen-images/label_seg.svg
+
+
+    Which is constructed as follows.
+
+    .. _morph-label-seg-code:
+
+    .. code-block:: Python
+
+        import arbor
+        from arbor import mpoint
+        from arbor import mpos
+
+        tree = arbor.segment_tree()
+        # Start with a cylinder segment for the soma (with tag 1)
+        tree.append(mnpos, mpoint(0,   0.0, 0, 2.0), mpoint( 4,  0.0, 0, 2.0), tag=1)
+        # Construct the first section of the dendritic tree,
+        # comprised of segments 1 and 2, attached to soma segment 0.
+        tree.append(0,     mpoint(4,   0.0, 0, 0.8), mpoint( 8,  0.0, 0, 0.8), tag=3)
+        tree.append(1,     mpoint(8,   0.0, 0, 0.8), mpoint(12, -0.5, 0, 0.8), tag=3)
+        # Construct the rest of the dendritic tree.
+        tree.append(2,     mpoint(12, -0.5, 0, 0.8), mpoint(20,  4.0, 0, 0.4), tag=3)
+        tree.append(3,     mpoint(20,  4.0, 0, 0.4), mpoint(26,  6.0, 0, 0.2), tag=3)
+        tree.append(2,     mpoint(12, -0.5, 0, 0.5), mpoint(19, -3.0, 0, 0.5), tag=3)
+        tree.append(5,     mpoint(19, -3.0, 0, 0.5), mpoint(24, -7.0, 0, 0.2), tag=3)
+        tree.append(5,     mpoint(19, -3.0, 0, 0.5), mpoint(23, -1.0, 0, 0.2), tag=3)
+        tree.append(7,     mpoint(23, -1.0, 0, 0.2), mpoint(26, -2.0, 0, 0.2), tag=3)
+        # Two segments that define the axon, with the first at the root, where its proximal
+        # end will be connected with the proximal end of the soma segment.
+        tree.append(mnpos, mpoint(0,   0.0, 0, 2.0), mpoint(-7,  0.0, 0, 0.4), tag=2)
+        tree.append(9,     mpoint(-7,  0.0, 0, 0.4), mpoint(-10, 0.0, 0, 0.4), tag=2)
+
+    .. method:: append(parent, prox, dist, tag)
+
+        Append a segment to the tree.
+
+        :return: index of the new segment
+        :param int parent: index of segment
+        :param mpoint prox: proximal end of the segment
+        :param mpoint dist: distal end of the segment
+        :param int tag: tag meta data of segment
+
+    .. method:: append(parent, dist, tag)
+        :noindex:
+
+        Append a segment to the tree whose proximal end has the location and
+        radius of the distal end of the parent segment.
+
+        This version of append can't be used for a segment at the root of the
+        tree, that is, when ``parent`` is :data:`mnpos`, in which case both proximal
+        and distal ends of the segment must be specified.
+
+        :return: index of the new segment
+        :param int parent: index of segment
+        :param mpoint dist: distal end of the segment
+        :param int tag: tag meta data of segment
+
+    .. method:: append(parent, x, y, z, radius, tag)
+        :noindex:
+
+        Append a segment to the tree whose proximal end has the location and
+        radius of the distal end of the parent segment.
+
+        This version of append can't be used for a segment at the root of the
+        tree, that is, when ``parent`` is :data:`mnpos`, in which case both proximal
+        and distal ends of the segment must be specified.
+
+        :return: index of the new segment
+        :param int parent: index of segment
+        :param float x: distal x coordinate (μm)
+        :param float y: distal y coordinate (μm)
+        :param float z: distal z coordinate (μm)
+        :param float radius: distal radius (μm)
+        :param int tag: tag meta data of segment
+
+    .. attribute:: empty
+        :type: bool
+
+        If the tree is empty (i.e. whether it has size 0)
+
+    .. attribute:: size
+        :type: int
+
+        The number of segments.
+
+    .. attribute:: parents
+        :type: list
+
+        A list of parent indexes of the segments.
+
+    .. attribute:: segments
+        :type: list
+
+        A list of the segments.
+
+.. py:function:: load_swc(filename)
+
+    Loads the morphology in an SWC file as a :class:`segment_tree`.
+
+    The samples in the SWC files are treated as the end points of segments, where a
+    sample and its parent form a segment.
+    The :py:attr:`tag <segment.tag>` of each segment is the
+    `structure identifier <http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html>`_
+    of the distal sample.
+    The structure identifier of the first (root) sample is ignored, as it can only be the
+    proximal end of any segment.
+
+    .. note::
+        This method does not interpret the first sample, typically associated with the soma,
+        as a sphere. SWCs with single point somas are, unfortunately, reasonably common, for example
+        `SONATA <https://github.com/AllenInstitute/sonata/blob/master/docs/SONATA_DEVELOPER_GUIDE.md#representing-biophysical-neuron-morphologies>`_
+        model descriptions.
+
+        Such representations are unfortunate because simulation tools like Arbor and NEURON require
+        the use of cylinders or fustrums to describe morphologies, and it is not possible to
+        infer how branches attached to the soma should be connected.
+
+        The :func:`load_swc_allen` function provides support for interpreting
+        such SWC files.
+
+
+    :param str filename: the name of the SWC file.
+    :rtype: segment_tree
+
+.. py:function:: load_swc_allen(filename, no_gaps=False)
+
+    Generate a segment tree from an SWC file following the rules prescribed by
+    AllenDB and Sonata. Specifically:
+
+        * The first sample (the root) is treated as the center of the soma.
+        * The morphology is translated such that the soma is centered at (0,0,0).
+        * The first sample has tag 1 (soma).
+        * All other samples have tags 2, 3 or 4 (axon, apic and dend respectively)
+
+    SONATA prescribes that there should be no gaps, however some models in AllenDB
+    have gaps between the start of sections and the soma. The ``no_gaps`` flag can be
+    used to enforce this requirement.
+
+    Arbor does not support modelling the soma as a sphere, so a cylinder with length
+    equal to the soma diameter is used. The cylinder is centered on the origin, and
+    aligned along the z axis.
+    Axons and apical dendrites are attached to the proximal end of the cylinder, and
+    dendrites to the distal end, with a gap between the start of each branch and the
+    end of the soma cylinder to which it is attached.
+
+    :param str filename: the name of the SWC file.
+    :param bool no_gaps: enforce that distance between soma center and branches attached to soma is the soma radius.
+    :rtype: segment_tree
+
+.. py:class:: morphology
+
+    A *morphology* describes the geometry of a cell as unbranched cables
+    with variable radius and their associated tree structure.
+
+    .. note::
+        A morphology takes a segment tree and construct the cable branches.
+        Meta data about branches and their properties that may be expensive to calculate
+        is stored for fast look up during later stages of model building, and
+        querying by users.
+
+        For this reason, morpholgies are read only. To change a morphology, a new
+        morphology should be created using a new segment tree.
+
+    There is one *constructor* for a morphology:
+
+    .. function:: morphology(segment_tree)
+
+        Construct from a segment tree.
+
+    The morphology provides an interface for querying morphology properties:
+
+    .. attribute:: empty
+            :type: bool
+
+            Indicates if the morphology is empty.
+
+    .. attribute:: num_branches
+            :type: int
+
+            The number of branches in the morphology.
+
+    .. method:: branch_parent(i)
+
+            The parent branch of a branch.
+
+            :param int i: branch index
+            :rtype: int
+
+    .. method:: branch_children(i)
+
+            The child branches of a branch.
+
+            :param int i: branch index
+            :rtype: list
+
+    .. method:: branch_segments(i)
+
+            A list of the segments in a branch, ordered from proximal to distal.
+
+            :param int i: branch index
+            :rtype: list
+
diff --git a/doc/nmodl.rst b/doc/nmodl.rst
new file mode 100644
index 0000000000000000000000000000000000000000..84474d31efcf4a64bc3b61a00b2c8f8603f3f791
--- /dev/null
+++ b/doc/nmodl.rst
@@ -0,0 +1,79 @@
+.. _nmodl:
+
+NMODL
+======
+
+*NMODL* is a `DSL <https://www.neuron.yale.edu/neuron/static/py_doc/modelspec/programmatic/mechanisms/nmodl.html>`_
+for describing ion channel and synapse dynamics that is used by NEURON,
+which provides the mod2c compiler parses dynamics described in NMODL to
+generate C code that is called from NEURON.
+
+Arbor has an NMODL compiler, *modcc*, that generates
+optimized code in C++ and CUDA, which is optimzed for
+the target architecture. NMODL does not have a formal specification,
+and its semantis are often
+ambiguous. To manage this, Arbor uses its own dialect of NMODL that
+does not allow some constructions used in NEURON's NMODL.
+
+.. note::
+    We hope to replace NMODL with a DSL that is well defined, and easier
+    for both users and the Arbor developers to work with in the long term.
+    Until then, please write issues on our GitHub with any questions
+    that you have about getting your NMODL files to work in Arbor.
+
+This page is a collection of NMODL rules for Arbor. It assumes that the reader
+alreay has a working knowledge of NMODL.
+
+Ions
+-----
+
+* Arbor recognizes ``na``, ``ca`` and ``k`` ions by default. Any new ions
+  must be added explicitly in Arbor along with their default properties and
+  valence (this can be done in the recipe or on a single cell model).
+  Simply specifying them in NMODL will not work.
+* The parameters and variabnles of each ion referenced in a ``USEION`` statement
+  are available automatically to the mechanism. The exposed variables are:
+  internal concentration ``Xi``, external concentration ``Xo``, reversal potential
+  ``eX`` and current ``iX``. It is an error to also mark these as
+  ``PARAMETER``, ``ASSIGNED`` or ``CONSTANT``.
+* ``READ`` and ``WRITE`` permissions of ``Xi``, ``Xo``, ``eX`` and ``iX`` can be set
+  in NMODL in the ``NEURON`` block. If a parameter is writable it is automatically
+  readable and doesn't need to be specified as both.
+* If ``Xi``, ``Xo``, ``eX``, ``iX`` are used in a ``PROCEDURE`` or ``FUNCTION``,
+  they need to be passed as arguments.
+* If ``Xi`` or ``Xo`` (internal and external concentrations) are written in the
+  NMODL mechanism they need to be specified as ``STATE`` variables.
+
+Special variables
+-----------------
+
+* Arbor exposes some parameters from the simulation to the NMODL mechanisms.
+  These include ``v``, ``diam``, ``celsius`` in addition to the previously
+  mentioned ion parameters.
+* Special variables should not be ``ASSIGNED`` or ``CONSTANT``,
+  they are ``PARAMETER``.
+* ``diam`` and ``celsius`` can be set from the simulation side.
+* ``v`` is a reserved varible name and can be written in NMODL.
+* If Special variables are used in a ``PROCEDURE`` or ``FUNCTION``, they need
+  to be passed as arguments.
+* ``dt`` is not exposed to NMODL mechanisms.
+
+Functions and blocks
+---------------------
+
+* ``SOLVE`` statements should be the first statement in the ``BREAKPOINT`` block.
+* The return variable of ``FUNCTION`` has to always be set. ``if`` without associated
+  ``else`` can break that if users are not careful.
+
+Unsupported features
+--------------------
+
+* Unit conversion is not supported in Arbor (there is limited support for parsing
+  units, which are just ignored).
+* Unit declaration is not supported (ex: ``FARADAY = (faraday)  (10000 coulomb)``).
+  They can be replaced by declaring them and setting their values in ``CONSTANT``.
+* ``FROM`` - ``TO`` clamping of variables is not supported. The tokens are parsed and ignored.
+  However, ``CONSERVE`` statements are supported.
+* ``TABLE`` is not supported, calculations are exact.
+* ``derivimplicit`` solving method is not supported, use ``cnexp`` instead.
+
diff --git a/doc/py_cable_cell.rst b/doc/py_cable_cell.rst
index ab437312c32dc3eb82ed5690de615341f1c438a3..7e783c2f45f79cc3451e673259de644cb4c28fab 100644
--- a/doc/py_cable_cell.rst
+++ b/doc/py_cable_cell.rst
@@ -3,73 +3,4 @@
 Python Cable Cells
 ====================
 
-The interface for specifying cell morphologies with the distribution of ion channels
-and synapses is a key part of the user interface. Arbor will have an advanced and user-friendly
-interface for this, which is currently under construction.
 
-To allow users to experiment will multi-compartment cells, we provide some helpers
-for generating cells with random morphologies, which are documented here.
-
-.. Warning::
-
-    These features will be deprecated once the morphology interface has been implemented.
-
-.. function:: make_cable_cell(seed, params)
-
-    Construct a branching :class:`cable_cell` with a random morphology (via parameter ``seed``) and
-    synapse end points locations described by parameter ``params``.
-
-    The soma has an area of 500 μm², a bulk resistivity of 100 Ω·cm,
-    and the ion channel and synapse dynamics are described by a Hodgkin-Huxley (HH) mechanism.
-    The default parameters of HH mechanisms are:
-
-    - Na-conductance        0.12 S⋅cm⁻²,
-    - K-conductance         0.036 S⋅cm⁻²,
-    - passive conductance   0.0003 S⋅cm⁻², and
-    - passive potential     -54.3 mV
-
-    Each cable has a diameter of 1 μm, a bulk resistivity of 100 Ω·cm,
-    and the ion channel and synapse dynamics are described by a passive/ leaky integrate-and-fire model with parameters:
-
-    - passive conductance   0.001 S⋅cm⁻², and
-    - resting potential     -65 mV
-
-    Further, a spike detector is added at the soma with threshold 10 mV,
-    and a synapse is added to the mid point of the first dendrite with an exponential synapse model:
-
-    - time decaying constant    2 ms
-    - resting potential         0 mV
-
-    Additional synapses are added based on the number of randomly generated :attr:`cell_parameters.synapses` on the cell.
-
-    :param seed: The seed is an integral value used to seed the random number generator, for which the :attr:`arbor.cell_member.gid` of the cell is a good default.
-
-    :param params: By default set to :class:`cell_parameters()`.
-
-.. class:: cell_parameters
-
-        Parameters used to generate random cell morphologies.
-        Where parameters must be given as ranges, the first value is at the soma,
-        and the last value is used on the last level.
-        Values at levels in between are found by linear interpolation.
-
-    .. attribute:: depth
-
-        The maximum depth of the branch structure
-        (i.e., maximum number of levels in the cell (not including the soma)).
-
-    .. attribute:: lengths
-
-        The length of the branch [μm], given as a range ``[l1, l2]``.
-
-    .. attribute:: synapses
-
-        The number of randomly generated synapses on the cell.
-
-    .. attribute:: branch_probs
-
-        The probability of a branch occuring, given as a range ``[p1, p2]``.
-
-    .. attribute:: compartments
-
-        The compartment count on a branch, given as a range ``[n1, n2]``.
diff --git a/doc/py_common.rst b/doc/py_common.rst
index 11cf1cb5882705581dd1aa255e160a4b7c178fa3..8fc9edf9f953d7b9890de3f86e4d87bd9de06b5a 100644
--- a/doc/py_common.rst
+++ b/doc/py_common.rst
@@ -79,18 +79,3 @@ The types defined below are used as identifiers for cells and members of cell-lo
 
             kind = arbor.cell_kind.cable
 
-Morphology
-----------
-
-.. class:: location
-
-    Construct a location specification with the :attr:`branch` id and the relative :attr:`position` on the branch ∈ [0.,1.], where 0. means proximal and 1. distal.
-
-    .. attribute:: branch
-
-        The id of the branch.
-
-    .. attribute:: position
-
-        The relative position (from 0., proximal, to 1., distal) on the branch.
-
diff --git a/doc/py_hardware.rst b/doc/py_hardware.rst
index 0a5e5cd168a15f469fa53c5098c7f785673f33bc..94ba7e8d6bea1ee0806086dae6a2794882d43a4e 100644
--- a/doc/py_hardware.rst
+++ b/doc/py_hardware.rst
@@ -74,11 +74,7 @@ The Python wrapper provides an API for:
     Enumerates the computational resources on a node to be used for a simulation,
     specifically the number of threads and identifier of a GPU if available.
 
-    .. function:: proc_allocation()
-
-        By default selects one thread and no GPU.
-
-    .. function:: proc_allocation(threads, gpu_id)
+    .. function:: proc_allocation([threads=1, gpu_id=None])
 
         Constructor that sets the number of :attr:`threads` and the id :attr:`gpu_id` of the available GPU.
 
diff --git a/doc/py_recipe.rst b/doc/py_recipe.rst
index 2d7f63cfe2c287180b79db99b3ca7d77a6100747..9f6577ceb8a616640bcc5c8b0593699cfa132b2e 100644
--- a/doc/py_recipe.rst
+++ b/doc/py_recipe.rst
@@ -298,10 +298,6 @@ An example of an event generator reads as follows:
 Cells
 ------
 
-.. class:: cable_cell
-
-   See :ref:`pycable_cell`.
-
 .. class:: lif_cell
 
     A benchmarking cell (leaky integrate-and-fire), used by Arbor developers to test communication performance,
diff --git a/doc/py_simulation.rst b/doc/py_simulation.rst
index d2a6c5d2e12d8cef9d86dc4feca822ef1ec3bf92..1b60539ea4c87907180aed8f94cece853e1ac33c 100644
--- a/doc/py_simulation.rst
+++ b/doc/py_simulation.rst
@@ -198,7 +198,7 @@ In order to analyze the data collected from an :class:`arbor.probe` the samples
 
 **Types**:
 
-.. class:: sample
+.. class:: trace_sample
 
     .. attribute:: time
 
diff --git a/doc/requirements.txt b/doc/requirements.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c4e363475fc7c5aa3fadbe0113a1f914ccc12ce0
--- /dev/null
+++ b/doc/requirements.txt
@@ -0,0 +1,2 @@
+sphinx>=2.3.0
+svgwrite
diff --git a/doc/scripts/gen-labels.py b/doc/scripts/gen-labels.py
new file mode 100644
index 0000000000000000000000000000000000000000..fd1bc4c6bd3d90959a215fbf3eb86580821ee008
--- /dev/null
+++ b/doc/scripts/gen-labels.py
@@ -0,0 +1,243 @@
+import arbor
+from arbor import mpoint
+
+def is_collocated(l, r):
+    return l[0]==r[0] and l[1]==r[1]
+
+def write_morphology(name, morph):
+    string = 'tmp = ['.format(name)
+    for i in range(morph.num_branches):
+        first = True
+        sections = '['
+        for seg in morph.branch_segments(i):
+            if not first:
+                if is_collocated((seg.prox.x, seg.prox.y), (last_dist.x, last_dist.y)):
+                    sections += ', '
+                else:
+                    sections += '], ['
+
+            first = False
+            p = seg.prox
+            d = seg.dist
+            sections += 'Segment(({}, {}, {}), ({}, {}, {}), {})'.format(p.x, p.y, p.radius, d.x, d.y, d.radius, seg.tag)
+            last_dist = seg.dist
+        sections += ']'
+
+        string += '\n    [{}],'.format(sections)
+    string += ']\n'
+    string += '{} = representation.make_morph(tmp)\n\n'.format(name)
+    return string
+
+# Describe the morphologies
+
+mnpos = arbor.mnpos
+
+# The morphology used for all of the region/locset illustrations
+label_tree = arbor.segment_tree()
+label_tree.append(mnpos, mpoint(0,   0.0, 0, 2.0), mpoint( 4,  0.0, 0, 2.0), tag=1)
+label_tree.append(0,     mpoint(4,   0.0, 0, 0.8), mpoint( 8,  0.0, 0, 0.8), tag=3)
+label_tree.append(1,     mpoint(8,   0.0, 0, 0.8), mpoint(12, -0.5, 0, 0.8), tag=3)
+label_tree.append(2,     mpoint(12, -0.5, 0, 0.8), mpoint(20,  4.0, 0, 0.4), tag=3)
+label_tree.append(3,     mpoint(20,  4.0, 0, 0.4), mpoint(26,  6.0, 0, 0.2), tag=3)
+label_tree.append(2,     mpoint(12, -0.5, 0, 0.5), mpoint(19, -3.0, 0, 0.5), tag=3)
+label_tree.append(5,     mpoint(19, -3.0, 0, 0.5), mpoint(24, -7.0, 0, 0.2), tag=3)
+label_tree.append(5,     mpoint(19, -3.0, 0, 0.5), mpoint(23, -1.0, 0, 0.2), tag=3)
+label_tree.append(7,     mpoint(23, -1.0, 0, 0.2), mpoint(26, -2.0, 0, 0.2), tag=3)
+label_tree.append(mnpos, mpoint(0,   0.0, 0, 2.0), mpoint(-7,  0.0, 0, 0.4), tag=2)
+label_tree.append(9,     mpoint(-7,  0.0, 0, 0.4), mpoint(-10, 0.0, 0, 0.4), tag=2)
+
+label_morph = arbor.morphology(label_tree)
+
+# The label morphology with some gaps (at start of dendritic tree and remove the axon hillock)
+label_tree = arbor.segment_tree()
+label_tree.append(mnpos, mpoint(0,   0.0, 0, 2.0), mpoint( 4,  0.0, 0, 2.0), tag=1)
+label_tree.append(0,    mpoint(5,   0.0, 0, 0.8), mpoint( 8,  0.0, 0, 0.8), tag=3)
+label_tree.append(1,    mpoint(8,   0.0, 0, 0.8), mpoint(12, -0.5, 0, 0.8), tag=3)
+label_tree.append(2,    mpoint(12, -0.5, 0, 0.8), mpoint(20,  4.0, 0, 0.4), tag=3)
+label_tree.append(3,    mpoint(20,  4.0, 0, 0.4), mpoint(26,  6.0, 0, 0.2), tag=3)
+label_tree.append(2,    mpoint(12, -0.5, 0, 0.5), mpoint(19, -3.0, 0, 0.5), tag=3)
+label_tree.append(5,    mpoint(19, -3.0, 0, 0.5), mpoint(24, -7.0, 0, 0.2), tag=3)
+label_tree.append(5,    mpoint(19, -3.0, 0, 0.5), mpoint(23, -1.0, 0, 0.2), tag=3)
+label_tree.append(7,    mpoint(23, -1.0, 0, 0.2), mpoint(26, -2.0, 0, 0.2), tag=3)
+label_tree.append(mnpos, mpoint(-2,  0.0, 0, 0.4), mpoint(-10, 0.0, 0, 0.4), tag=2)
+
+detached_morph = arbor.morphology(label_tree)
+
+# soma with "stacked cylinders"
+stacked_tree = arbor.segment_tree()
+stacked_tree.append(mnpos, mpoint(0,   0.0, 0, 0.5), mpoint( 1,  0.0, 0, 1.5), tag=1)
+stacked_tree.append(0,    mpoint(1,   0.0, 0, 1.5), mpoint( 2,  0.0, 0, 2.5), tag=1)
+stacked_tree.append(1,    mpoint(2,   0.0, 0, 2.5), mpoint( 3,  0.0, 0, 2.5), tag=1)
+stacked_tree.append(2,    mpoint(3,   0.0, 0, 2.5), mpoint( 4,  0.0, 0, 1.2), tag=1)
+stacked_tree.append(3,    mpoint(4,   0.0, 0, 0.8), mpoint( 8,  0.0, 0, 0.8), tag=3)
+stacked_tree.append(4,    mpoint(8,   0.0, 0, 0.8), mpoint(12, -0.5, 0, 0.8), tag=3)
+stacked_tree.append(5,    mpoint(12, -0.5, 0, 0.8), mpoint(20,  4.0, 0, 0.4), tag=3)
+stacked_tree.append(6,    mpoint(20,  4.0, 0, 0.4), mpoint(26,  6.0, 0, 0.2), tag=3)
+stacked_tree.append(5,    mpoint(12, -0.5, 0, 0.5), mpoint(19, -3.0, 0, 0.5), tag=3)
+stacked_tree.append(8,    mpoint(19, -3.0, 0, 0.5), mpoint(24, -7.0, 0, 0.2), tag=3)
+stacked_tree.append(8,    mpoint(19, -3.0, 0, 0.5), mpoint(23, -1.0, 0, 0.2), tag=3)
+stacked_tree.append(10,   mpoint(23, -1.0, 0, 0.2), mpoint(26, -2.0, 0, 0.2), tag=3)
+stacked_tree.append(mnpos, mpoint(0,   0.0, 0, 0.4), mpoint(-7,  0.0, 0, 0.4), tag=2)
+stacked_tree.append(12,   mpoint(-7,  0.0, 0, 0.4), mpoint(-10, 0.0, 0, 0.4), tag=2)
+
+stacked_morph = arbor.morphology(stacked_tree)
+
+# spherical cell with radius 2 μm
+tree = arbor.segment_tree()
+tree.append(mnpos, mpoint(-2, 0, 0, 2), mpoint(2, 0, 0, 2), tag=1)
+sphere_morph = arbor.morphology(tree)
+
+# single branch: one tapered segment
+tree = arbor.segment_tree()
+tree.append(mnpos, mpoint(0, 0, 0, 1), mpoint(10, 0, 0, 0.5), tag=3)
+branch_morph1 = arbor.morphology(tree)
+
+# single branch: multiple segments, continuous radius
+tree = arbor.segment_tree()
+tree.append(mnpos, mpoint( 0.0,  0.0,  0.0, 1.0), mpoint( 3.0,  0.2,  0.0, 0.8), tag=1)
+tree.append(0,     mpoint( 3.0,  0.2,  0.0, 0.8), mpoint( 5.0, -0.1,  0.0, 0.7), tag=2)
+tree.append(1,     mpoint( 5.0, -0.1,  0.0, 0.7), mpoint( 8.0,  0.0,  0.0, 0.6), tag=2)
+tree.append(2,     mpoint( 8.0,  0.0,  0.0, 0.6), mpoint(10.0,  0.0,  0.0, 0.5), tag=3)
+branch_morph2 = arbor.morphology(tree)
+
+# single branch: multiple segments, gaps
+tree = arbor.segment_tree()
+tree.append(mnpos, mpoint( 0.0,  0.0,  0.0, 1.0), mpoint(3.0,  0.2,  0.0, 0.8), tag=1)
+tree.append(0,     mpoint( 3.0,  0.2,  0.0, 0.8), mpoint(5.0, -0.1,  0.0, 0.7), tag=2)
+tree.append(1,     mpoint( 6.0, -0.1,  0.0, 0.7), mpoint(9.0, 0.0,  0.0, 0.6), tag=2)
+tree.append(2,     mpoint( 9.0,  0.0,  0.0, 0.6), mpoint(11.0, 0.0,  0.0, 0.5), tag=3)
+branch_morph3 = arbor.morphology(tree)
+
+# single branch: multiple segments, discontinuous radius
+tree = arbor.segment_tree()
+tree.append(mnpos, mpoint( 0.0,  0.0,  0.0, 1.0), mpoint( 3.0,  0.2,  0.0, 0.8), tag=1)
+tree.append(0,     mpoint( 3.0,  0.2,  0.0, 0.8), mpoint( 5.0, -0.1,  0.0, 0.7), tag=2)
+tree.append(1,     mpoint( 5.0, -0.1,  0.0, 0.7), mpoint( 8.0,  0.0,  0.0, 0.5), tag=2)
+tree.append(2,     mpoint( 8.0,  0.0,  0.0, 0.3), mpoint(10.0,  0.0,  0.0, 0.5), tag=3)
+branch_morph4 = arbor.morphology(tree)
+
+tree = arbor.segment_tree()
+tree.append(mnpos, mpoint( 0.0, 0.0, 0.0, 1.0), mpoint(10.0, 0.0, 0.0, 0.5), tag= 3)
+tree.append(0,     mpoint(15.0, 3.0, 0.0, 0.2), tag= 3)
+tree.append(0,     mpoint(15.0,-3.0, 0.0, 0.2), tag= 3)
+yshaped_morph = arbor.morphology(tree)
+
+tree = arbor.segment_tree()
+tree.append(mnpos, mpoint(-3.0, 0.0, 0.0, 3.0), mpoint( 3.0, 0.0, 0.0, 3.0), tag=1)
+tree.append(0, mpoint( 4.0, -1.0,  0.0, 0.6), mpoint(10.0,  -2.0,  0.0, 0.5), tag=3)
+tree.append(1, mpoint(15.0, -1.0,  0.0, 0.5), tag=3)
+tree.append(2, mpoint(18.0, -5.0,  0.0, 0.3), tag=3)
+tree.append(2, mpoint(20.0,  2.0,  0.0, 0.3), tag=3)
+ysoma_morph1 = arbor.morphology(tree)
+
+tree = arbor.segment_tree()
+tree.append(mnpos, mpoint(-3.0, 0.0, 0.0, 3.0),   mpoint( 3.0, 0.0, 0.0, 3.0), tag=1)
+tree.append(0,     mpoint( 4.0, -1.0,  0.0, 0.6), mpoint(10.0,-2.0,  0.0, 0.5), tag=3)
+tree.append(1,     mpoint(15.0, -1.0,  0.0, 0.5), tag= 3)
+tree.append(2,     mpoint(18.0, -5.0,  0.0, 0.3), tag= 3)
+tree.append(2,     mpoint(20.0,  2.0,  0.0, 0.3), tag= 3)
+tree.append(0,     mpoint( 2.0,  1.0,  0.0, 0.6), mpoint(12.0, 4.0,  0.0, 0.5), tag=3)
+tree.append(5,     mpoint(18.0,  4.0,  0.0, 0.3), tag= 3)
+tree.append(5,     mpoint(16.0,  9.0,  0.0, 0.1), tag= 3)
+tree.append(mnpos, mpoint(-3.5,  0.0,  0.0, 1.5), mpoint(-6.0,-0.2,  0.0, 0.5), tag=2)
+tree.append(8,     mpoint(-15.0,-0.1,  0.0, 0.5), tag=2)
+ysoma_morph2 = arbor.morphology(tree)
+
+tree = arbor.segment_tree()
+tree.append(mnpos, mpoint(-3.0, 0.0, 0.0, 3.0),   mpoint( 3.0, 0.0, 0.0, 3.0), tag=1)
+tree.append(0,     mpoint( 3.0,  0.0,  0.0, 0.6), mpoint(9.0,-1.0,  0.0, 0.5), tag=3)
+tree.append(1,     mpoint(14.0,  0.0,  0.0, 0.5), tag= 3)
+tree.append(2,     mpoint(17.0, -4.0,  0.0, 0.3), tag= 3)
+tree.append(2,     mpoint(19.0,  3.0,  0.0, 0.3), tag= 3)
+tree.append(0,     mpoint( 3.0,  0.0,  0.0, 0.6), mpoint(13.0, 3.0,  0.0, 0.5), tag=3)
+tree.append(5,     mpoint(19.0,  3.0,  0.0, 0.3), tag= 3)
+tree.append(5,     mpoint(17.0,  8.0,  0.0, 0.1), tag= 3)
+tree.append(mnpos, mpoint(-3.0,  0.0,  0.0, 1.5), mpoint(-5.5,-0.2,  0.0, 0.5), tag=2)
+tree.append(8,     mpoint(-14.5,-0.1,  0.0, 0.5), tag=2)
+ysoma_morph3 = arbor.morphology(tree)
+
+regions  = {
+            'empty': '(nil)',
+            'all': '(all)',
+            'tag1': '(tag 1)',
+            'tag2': '(tag 2)',
+            'tag3': '(tag 3)',
+            'tag4': '(tag 4)',
+            '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))',
+            '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"))',
+            'lhs' : '(join (cable 0 0.5 1) (cable 1 0 0.5))',
+            'rhs' : '(branch 1)',
+            'and': '(intersect (region "lhs") (region "rhs"))',
+            'or':  '(join      (region "lhs") (region "rhs"))',
+          }
+locsets = {
+            'root': '(root)',
+            'term': '(terminal)',
+            'rand_dend': '(uniform (region "dend") 0 50 0)',
+            'loc15': '(location 1 0.5)',
+            'uniform0': '(uniform (tag 3) 0 9 0)',
+            'uniform1': '(uniform (tag 3) 0 9 1)',
+            '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))',
+            'proxint_in': '(sum (location 1 0.8) (location 2 0.3))',
+            'loctest' : '(distal (complete (join (branch 1) (branch 0))))',
+            'restrict': '(restrict  (terminal) (tag 3))',
+          }
+
+labels = {**regions, **locsets}
+d = arbor.label_dict(labels)
+
+# Create a cell to concretise the region and locset definitions
+cell = arbor.cable_cell(label_morph, d)
+
+################################################################################
+# Output all of the morphologies and reion/locset definitions to a Python script
+# that can be run during the documentation build to generate images.
+################################################################################
+f = open('inputs.py', 'w')
+f.write('import representation\n')
+f.write('from representation import Segment\n')
+
+f.write('\n############# morphologies\n\n')
+f.write(write_morphology('label_morph',    label_morph))
+f.write(write_morphology('detached_morph', detached_morph))
+f.write(write_morphology('stacked_morph',  stacked_morph))
+f.write(write_morphology('sphere_morph',   sphere_morph))
+f.write(write_morphology('branch_morph1',  branch_morph1))
+f.write(write_morphology('branch_morph2',  branch_morph2))
+f.write(write_morphology('branch_morph3',  branch_morph3))
+f.write(write_morphology('branch_morph4',  branch_morph4))
+f.write(write_morphology('yshaped_morph',  yshaped_morph))
+f.write(write_morphology('ysoma_morph1',   ysoma_morph1))
+f.write(write_morphology('ysoma_morph2',   ysoma_morph2))
+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)]
+    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)]
+    f.write('reg_{} = {{\'type\': \'region\', \'value\': {}}}\n'.format(label, comps))
+
+f.close()
+
diff --git a/doc/scripts/inputs.py b/doc/scripts/inputs.py
new file mode 100644
index 0000000000000000000000000000000000000000..63b79dfe0bf6f8c10838fb8567966c0ea35f8018
--- /dev/null
+++ b/doc/scripts/inputs.py
@@ -0,0 +1,132 @@
+import representation
+from representation import Segment
+
+############# morphologies
+
+tmp = [
+    [[Segment((0.0, 0.0, 2.0), (4.0, 0.0, 2.0), 1), Segment((4.0, 0.0, 0.8), (8.0, 0.0, 0.8), 3), Segment((8.0, 0.0, 0.8), (12.0, -0.5, 0.8), 3)]],
+    [[Segment((12.0, -0.5, 0.8), (20.0, 4.0, 0.4), 3), Segment((20.0, 4.0, 0.4), (26.0, 6.0, 0.2), 3)]],
+    [[Segment((12.0, -0.5, 0.5), (19.0, -3.0, 0.5), 3)]],
+    [[Segment((19.0, -3.0, 0.5), (24.0, -7.0, 0.2), 3)]],
+    [[Segment((19.0, -3.0, 0.5), (23.0, -1.0, 0.2), 3), Segment((23.0, -1.0, 0.2), (26.0, -2.0, 0.2), 3)]],
+    [[Segment((0.0, 0.0, 2.0), (-7.0, 0.0, 0.4), 2), Segment((-7.0, 0.0, 0.4), (-10.0, 0.0, 0.4), 2)]],]
+label_morph = representation.make_morph(tmp)
+
+tmp = [
+    [[Segment((0.0, 0.0, 2.0), (4.0, 0.0, 2.0), 1)], [Segment((5.0, 0.0, 0.8), (8.0, 0.0, 0.8), 3), Segment((8.0, 0.0, 0.8), (12.0, -0.5, 0.8), 3)]],
+    [[Segment((12.0, -0.5, 0.8), (20.0, 4.0, 0.4), 3), Segment((20.0, 4.0, 0.4), (26.0, 6.0, 0.2), 3)]],
+    [[Segment((12.0, -0.5, 0.5), (19.0, -3.0, 0.5), 3)]],
+    [[Segment((19.0, -3.0, 0.5), (24.0, -7.0, 0.2), 3)]],
+    [[Segment((19.0, -3.0, 0.5), (23.0, -1.0, 0.2), 3), Segment((23.0, -1.0, 0.2), (26.0, -2.0, 0.2), 3)]],
+    [[Segment((-2.0, 0.0, 0.4), (-10.0, 0.0, 0.4), 2)]],]
+detached_morph = representation.make_morph(tmp)
+
+tmp = [
+    [[Segment((0.0, 0.0, 0.5), (1.0, 0.0, 1.5), 1), Segment((1.0, 0.0, 1.5), (2.0, 0.0, 2.5), 1), Segment((2.0, 0.0, 2.5), (3.0, 0.0, 2.5), 1), Segment((3.0, 0.0, 2.5), (4.0, 0.0, 1.2), 1), Segment((4.0, 0.0, 0.8), (8.0, 0.0, 0.8), 3), Segment((8.0, 0.0, 0.8), (12.0, -0.5, 0.8), 3)]],
+    [[Segment((12.0, -0.5, 0.8), (20.0, 4.0, 0.4), 3), Segment((20.0, 4.0, 0.4), (26.0, 6.0, 0.2), 3)]],
+    [[Segment((12.0, -0.5, 0.5), (19.0, -3.0, 0.5), 3)]],
+    [[Segment((19.0, -3.0, 0.5), (24.0, -7.0, 0.2), 3)]],
+    [[Segment((19.0, -3.0, 0.5), (23.0, -1.0, 0.2), 3), Segment((23.0, -1.0, 0.2), (26.0, -2.0, 0.2), 3)]],
+    [[Segment((0.0, 0.0, 0.4), (-7.0, 0.0, 0.4), 2), Segment((-7.0, 0.0, 0.4), (-10.0, 0.0, 0.4), 2)]],]
+stacked_morph = representation.make_morph(tmp)
+
+tmp = [
+    [[Segment((-2.0, 0.0, 2.0), (2.0, 0.0, 2.0), 1)]],]
+sphere_morph = representation.make_morph(tmp)
+
+tmp = [
+    [[Segment((0.0, 0.0, 1.0), (10.0, 0.0, 0.5), 3)]],]
+branch_morph1 = representation.make_morph(tmp)
+
+tmp = [
+    [[Segment((0.0, 0.0, 1.0), (3.0, 0.2, 0.8), 1), Segment((3.0, 0.2, 0.8), (5.0, -0.1, 0.7), 2), Segment((5.0, -0.1, 0.7), (8.0, 0.0, 0.6), 2), Segment((8.0, 0.0, 0.6), (10.0, 0.0, 0.5), 3)]],]
+branch_morph2 = representation.make_morph(tmp)
+
+tmp = [
+    [[Segment((0.0, 0.0, 1.0), (3.0, 0.2, 0.8), 1), Segment((3.0, 0.2, 0.8), (5.0, -0.1, 0.7), 2)], [Segment((6.0, -0.1, 0.7), (9.0, 0.0, 0.6), 2), Segment((9.0, 0.0, 0.6), (11.0, 0.0, 0.5), 3)]],]
+branch_morph3 = representation.make_morph(tmp)
+
+tmp = [
+    [[Segment((0.0, 0.0, 1.0), (3.0, 0.2, 0.8), 1), Segment((3.0, 0.2, 0.8), (5.0, -0.1, 0.7), 2), Segment((5.0, -0.1, 0.7), (8.0, 0.0, 0.5), 2), Segment((8.0, 0.0, 0.3), (10.0, 0.0, 0.5), 3)]],]
+branch_morph4 = representation.make_morph(tmp)
+
+tmp = [
+    [[Segment((0.0, 0.0, 1.0), (10.0, 0.0, 0.5), 3)]],
+    [[Segment((10.0, 0.0, 0.5), (15.0, 3.0, 0.2), 3)]],
+    [[Segment((10.0, 0.0, 0.5), (15.0, -3.0, 0.2), 3)]],]
+yshaped_morph = representation.make_morph(tmp)
+
+tmp = [
+    [[Segment((-3.0, 0.0, 3.0), (3.0, 0.0, 3.0), 1)], [Segment((4.0, -1.0, 0.6), (10.0, -2.0, 0.5), 3), Segment((10.0, -2.0, 0.5), (15.0, -1.0, 0.5), 3)]],
+    [[Segment((15.0, -1.0, 0.5), (18.0, -5.0, 0.3), 3)]],
+    [[Segment((15.0, -1.0, 0.5), (20.0, 2.0, 0.3), 3)]],]
+ysoma_morph1 = representation.make_morph(tmp)
+
+tmp = [
+    [[Segment((-3.0, 0.0, 3.0), (3.0, 0.0, 3.0), 1)]],
+    [[Segment((4.0, -1.0, 0.6), (10.0, -2.0, 0.5), 3), Segment((10.0, -2.0, 0.5), (15.0, -1.0, 0.5), 3)]],
+    [[Segment((15.0, -1.0, 0.5), (18.0, -5.0, 0.3), 3)]],
+    [[Segment((15.0, -1.0, 0.5), (20.0, 2.0, 0.3), 3)]],
+    [[Segment((2.0, 1.0, 0.6), (12.0, 4.0, 0.5), 3)]],
+    [[Segment((12.0, 4.0, 0.5), (18.0, 4.0, 0.3), 3)]],
+    [[Segment((12.0, 4.0, 0.5), (16.0, 9.0, 0.1), 3)]],
+    [[Segment((-3.5, 0.0, 1.5), (-6.0, -0.2, 0.5), 2), Segment((-6.0, -0.2, 0.5), (-15.0, -0.1, 0.5), 2)]],]
+ysoma_morph2 = representation.make_morph(tmp)
+
+tmp = [
+    [[Segment((-3.0, 0.0, 3.0), (3.0, 0.0, 3.0), 1)]],
+    [[Segment((3.0, 0.0, 0.6), (9.0, -1.0, 0.5), 3), Segment((9.0, -1.0, 0.5), (14.0, 0.0, 0.5), 3)]],
+    [[Segment((14.0, 0.0, 0.5), (17.0, -4.0, 0.3), 3)]],
+    [[Segment((14.0, 0.0, 0.5), (19.0, 3.0, 0.3), 3)]],
+    [[Segment((3.0, 0.0, 0.6), (13.0, 3.0, 0.5), 3)]],
+    [[Segment((13.0, 3.0, 0.5), (19.0, 3.0, 0.3), 3)]],
+    [[Segment((13.0, 3.0, 0.5), (17.0, 8.0, 0.1), 3)]],
+    [[Segment((-3.0, 0.0, 1.5), (-5.5, -0.2, 0.5), 2), Segment((-5.5, -0.2, 0.5), (-14.5, -0.1, 0.5), 2)]],]
+ysoma_morph3 = representation.make_morph(tmp)
+
+
+############# locsets
+
+ls_root  = {'type': 'locset', 'value': [(0, 0.0)]}
+ls_term  = {'type': 'locset', 'value': [(1, 1.0), (3, 1.0), (4, 1.0), (5, 1.0)]}
+ls_rand_dend  = {'type': 'locset', 'value': [(0, 0.5547193370156588), (0, 0.5841758202819731), (0, 0.607192003545501), (0, 0.6181091003428546), (0, 0.6190845627201184), (0, 0.7027325639263277), (0, 0.7616129092226993), (0, 0.9645150497869694), (1, 0.15382287505908834), (1, 0.2594719824047551), (1, 0.28087652335178354), (1, 0.3729681478609085), (1, 0.3959560134241004), (1, 0.4629424550242548), (1, 0.47346867377446744), (1, 0.5493486883630476), (1, 0.6227685370674116), (1, 0.6362196581003494), (1, 0.6646511214508091), (1, 0.7157318936458146), (1, 0.7464198558822775), (1, 0.77074507802833), (1, 0.7860238136304932), (1, 0.8988928261704698), (1, 0.9581259332943499), (2, 0.12773985425987294), (2, 0.3365926476076694), (2, 0.44454300804769703), (2, 0.5409466695719178), (2, 0.5767511435223905), (2, 0.6340206909931745), (2, 0.6354772583375223), (2, 0.6807941995943213), (2, 0.774655947503608), (3, 0.05020708596877571), (3, 0.25581431877212274), (3, 0.2958305460715556), (3, 0.296698184761692), (3, 0.509669134988683), (3, 0.7662305637426007), (3, 0.8565839889923518), (3, 0.8889077221517746), (4, 0.24311286693286885), (4, 0.4354361205546333), (4, 0.4467752481260171), (4, 0.5308169153994543), (4, 0.5701465671464049), (4, 0.670081739879954), (4, 0.6995486862583797), (4, 0.8186709628604206), (4, 0.9141224600171143)]}
+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_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)]}
+ls_restrict  = {'type': 'locset', 'value': [(1, 1.0), (3, 1.0), (4, 1.0)]}
+
+############# regions
+
+reg_empty = {'type': 'region', 'value': []}
+reg_all = {'type': 'region', 'value': [(0, 0.0, 1.0), (1, 0.0, 1.0), (2, 0.0, 1.0), (3, 0.0, 1.0), (4, 0.0, 1.0), (5, 0.0, 1.0)]}
+reg_tag1 = {'type': 'region', 'value': [(0, 0.0, 0.3324708796524168)]}
+reg_tag2 = {'type': 'region', 'value': [(5, 0.0, 1.0)]}
+reg_tag3 = {'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_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_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)]}
+reg_cable_1_31 = {'type': 'region', 'value': [(1, 0.3, 1.0)]}
+reg_cable_1_37 = {'type': 'region', 'value': [(1, 0.3, 0.7)]}
+reg_proxint = {'type': 'region', 'value': [(0, 0.7697564611867647, 1.0), (1, 0.4774887508467626, 0.8), (2, 0.0, 0.3)]}
+reg_proxintinf = {'type': 'region', 'value': [(0, 0.0, 1.0), (1, 0.0, 0.8), (2, 0.0, 0.3)]}
+reg_distint = {'type': 'region', 'value': [(1, 0.5, 0.8225112491532374), (2, 0.7, 1.0), (3, 0.0, 0.432615327328525), (4, 0.0, 0.3628424955125098), (5, 0.1, 0.6)]}
+reg_distintinf = {'type': 'region', 'value': [(1, 0.5, 1.0), (2, 0.7, 1.0), (3, 0.0, 1.0), (4, 0.0, 1.0), (5, 0.1, 1.0)]}
+reg_lhs = {'type': 'region', 'value': [(0, 0.5, 1.0), (1, 0.0, 0.5)]}
+reg_rhs = {'type': 'region', 'value': [(1, 0.0, 1.0)]}
+reg_and = {'type': 'region', 'value': [(1, 0.0, 0.5)]}
+reg_or = {'type': 'region', 'value': [(0, 0.5, 1.0), (1, 0.0, 1.0)]}
diff --git a/doc/scripts/make_images.py b/doc/scripts/make_images.py
new file mode 100644
index 0000000000000000000000000000000000000000..0df83945c9063a7b80ce3fd821eae19b240c6b1b
--- /dev/null
+++ b/doc/scripts/make_images.py
@@ -0,0 +1,285 @@
+import copy
+import svgwrite
+import math
+import inputs
+
+tag_colors = ['white', '#ffc2c2', 'gray', '#c2caff']
+
+#
+# ############################################
+#
+
+def translate(x, f, xshift):
+    return (f*x[0]+xshift, -f*x[1])
+
+def translate_all(points, f, xshift):
+    return [translate(x, f, xshift) for x in points]
+
+# Draw one or more morphologies, side by side.
+# Each morphology can be drawn as segments or branches.
+def morph_image(morphs, methods, filename, sc=20):
+    assert(len(morphs)==len(methods))
+
+    print('generating:', filename)
+    dwg = svgwrite.Drawing(filename=filename, debug=True)
+
+    # Width of lines and circle strokes.
+    line_width=0.1*sc
+
+    # Padding around image.
+    fudge=1.5*sc
+
+    linecolor='black'
+    pointcolor='red'
+    lines = dwg.add(dwg.g(id='lines',
+                          stroke=linecolor,
+                          fill='white',
+                          stroke_width=line_width))
+    points = dwg.add(dwg.g(id='points',
+                           stroke=pointcolor,
+                           fill=pointcolor,
+                           stroke_width=line_width))
+    numbers = dwg.add(dwg.g(id='numbers',
+                             text_anchor='middle',
+                             alignment_baseline='middle'))
+
+    minx = math.inf
+    miny = math.inf
+    maxx = -math.inf
+    maxy = -math.inf
+
+    offset = 0
+
+    bcolor = 'mediumslateblue'
+    branchfillcolor = 'lightgray'
+
+    nmorph = len(morphs)
+
+    for l in range(nmorph):
+        morph = morphs[l]
+        method = methods[l]
+
+        nbranches = len(morph)
+
+        segid = 0
+        for i in range(nbranches):
+            branch = morph[i]
+
+            lx, ux, ly, uy = branch.minmax()
+            minx = min(minx,  sc*lx+offset)
+            miny = min(miny,  sc*ly)
+            maxx = max(maxx,  sc*ux+offset)
+            maxy = max(maxy,  sc*uy)
+
+            if method=='segments':
+                for sec in branch.sections:
+                    for seg in sec:
+                        if seg.length>0.00001: # only draw nonzero length segments
+                            line = translate_all(seg.corners(), sc, offset)
+                            lines.add(dwg.polygon(points=line, fill=tag_colors[seg.tag]))
+
+                            pos = translate(seg.location(0.5), sc, offset)
+                            points.add(dwg.circle(center=pos,
+                                                  stroke='black',
+                                                  r=sc*0.55,
+                                                  fill='white'))
+                            # The svg alignment_baseline attribute:
+                            #   - works on Chrome/Chromium
+                            #   - doesn't work on Firefox
+                            # so for now we just shift the relative position by sc/3
+                            label_pos = (pos[0], pos[1]+sc/3)
+                            numbers.add(dwg.text(str(segid),
+                                                  insert=label_pos,
+                                                  stroke='black',
+                                                  fill='black'))
+                        segid += 1
+
+            elif method=='branches':
+                for line in branch.outline():
+                    lines.add(dwg.polygon(points=translate_all(line, sc, offset),
+                                          fill=branchfillcolor))
+
+                pos = translate(branch.location(0.5), sc, offset)
+                points.add(dwg.circle(center=pos,
+                                      stroke=bcolor,
+                                      r=sc*0.55,
+                                      fill=bcolor))
+                # The svg alignment_baseline attribute:
+                #   - works on Chrome/Chromium
+                #   - doesn't work on Firefox
+                # so for now we just shift the relative position by sc/3
+                label_pos = (pos[0], pos[1]+sc/3)
+                numbers.add(dwg.text(str(i),
+                                      insert=label_pos,
+                                      stroke='white',
+                                      fill='white'))
+        offset = maxx - minx + sc
+
+
+    # Find extent of image.
+    minx -= fudge
+    miny -= fudge
+    maxx += fudge
+    maxy += fudge
+    width = maxx-minx
+    height = maxy-miny
+    dwg.viewbox(minx, -maxy, width, height)
+
+    # Write the image to file.
+    dwg.save()
+
+# Generate an image that illustrates regions and locsets on a morphology.
+#
+# Can't handle morpholgies with gaps, where segemnts with a parent-child
+# ordering don't have collocated distal-proximal locations respectively.
+# Handling this case would make rendering regions more complex, but would
+# not bee too hard to support.
+def label_image(morphology, labels, filename, sc=20):
+    morph = morphology
+    print('generating:', filename)
+    dwg = svgwrite.Drawing(filename=filename, debug=True)
+
+    # Width of lines and circle strokes.
+    line_width=0.2*sc
+
+    # Padding around image.
+    fudge=1.5*sc
+
+    linecolor='black'
+    pointcolor='red'
+    lines = dwg.add(dwg.g(id='lines',
+                          stroke=linecolor,
+                          fill='white',
+                          stroke_width=line_width))
+    points = dwg.add(dwg.g(id='points',
+                           stroke=pointcolor,
+                           fill=pointcolor,
+                           stroke_width=line_width))
+
+    minx = math.inf
+    miny = math.inf
+    maxx = -math.inf
+    maxy = -math.inf
+
+    offset = 0
+
+    branchfillcolor = 'lightgray'
+
+    nimage = len(labels)
+    for l in range(nimage):
+        lab = labels[l]
+
+        nbranches = len(morph)
+
+        # Draw the outline of the cell
+        for i in range(nbranches):
+            branch = morph[i]
+
+            lx, ux, ly, uy = branch.minmax()
+            minx = min(minx,  sc*lx+offset)
+            miny = min(miny,  sc*ly)
+            maxx = max(maxx,  sc*ux+offset)
+            maxy = max(maxy,  sc*uy)
+
+            for line in branch.outline():
+                lines.add(dwg.polygon(points=translate_all(line, sc, offset),
+                                      fill=branchfillcolor,
+                                      stroke=branchfillcolor))
+
+        # Draw the root
+        root = translate(morph[0].location(0), sc, offset)
+        points.add(dwg.circle(center=root, stroke='red', r=sc/2.5, fill='white'))
+
+        if lab['type'] == 'locset':
+            for loc in lab['value']:
+                bid = loc[0]
+                pos = loc[1]
+
+                loc = translate(morph[bid].location(pos), sc, offset)
+                points.add(dwg.circle(center=loc, stroke='black', r=sc/3, fill='white'))
+
+        if lab['type'] == 'region':
+            for cab in lab['value']:
+                # skip zero length cables
+                bid  = cab[0]
+                ppos = cab[1]
+                dpos = cab[2]
+
+                # Don't draw zero-length cables
+                # How should these be drawn: with a line or a circle?
+                if ppos==dpos: continue
+
+                for line in morph[bid].outline(ppos, dpos):
+                    lines.add(dwg.polygon(points=translate_all(line, sc, offset),
+                                          fill='black',
+                                          stroke=branchfillcolor))
+
+        offset = maxx - minx + sc
+
+    # Find extent of image.
+    minx -= fudge
+    miny -= fudge
+    maxx += fudge
+    maxy += fudge
+    width = maxx-minx
+    height = maxy-miny
+    dwg.viewbox(minx, -maxy, width, height)
+
+    # Write the image to file.
+    dwg.save()
+
+def generate(path=''):
+
+    morph_image([inputs.label_morph],    ['branches'], path+'/label_branch.svg')
+
+    morph_image([inputs.label_morph],    ['segments'], path+'/label_seg.svg')
+    morph_image([inputs.detached_morph], ['segments'], path+'/detached_seg.svg')
+    morph_image([inputs.stacked_morph],  ['segments'], path+'/stacked_seg.svg')
+
+    morph_image([inputs.label_morph, inputs.label_morph], ['segments', 'branches'], path+'/label_morph.svg')
+    morph_image([inputs.detached_morph, inputs.detached_morph], ['segments', 'branches'], path+'/detached_morph.svg')
+    morph_image([inputs.stacked_morph, inputs.stacked_morph], ['segments', 'branches'], path+'/stacked_morph.svg')
+    morph_image([inputs.sphere_morph, inputs.sphere_morph], ['segments', 'branches'], path+'/sphere_morph.svg')
+    morph_image([inputs.branch_morph1, inputs.branch_morph1], ['segments', 'branches'], path+'/branch_morph1.svg')
+    morph_image([inputs.branch_morph2, inputs.branch_morph2], ['segments', 'branches'], path+'/branch_morph2.svg')
+    morph_image([inputs.branch_morph3, inputs.branch_morph3], ['segments', 'branches'], path+'/branch_morph3.svg')
+    morph_image([inputs.branch_morph4, inputs.branch_morph4], ['segments', 'branches'], path+'/branch_morph4.svg')
+    morph_image([inputs.yshaped_morph, inputs.yshaped_morph], ['segments', 'branches'], path+'/yshaped_morph.svg')
+    morph_image([inputs.ysoma_morph1,  inputs.ysoma_morph1],  ['segments', 'branches'], path+'/ysoma_morph1.svg')
+    morph_image([inputs.ysoma_morph2,  inputs.ysoma_morph2],  ['segments', 'branches'], path+'/ysoma_morph2.svg')
+    morph_image([inputs.ysoma_morph3,  inputs.ysoma_morph3],  ['segments', 'branches'], path+'/ysoma_morph3.svg')
+
+    ####################### locsets
+
+    label_image(inputs.label_morph, [inputs.ls_term, inputs.ls_rand_dend], path+'/locset_label_examples.svg')
+    label_image(inputs.label_morph, [inputs.reg_dend, inputs.reg_radlt5], path+'/region_label_examples.svg')
+    label_image(inputs.label_morph, [inputs.ls_root], path+'/root_label.svg')
+    label_image(inputs.label_morph, [inputs.ls_term], path+'/term_label.svg')
+    label_image(inputs.label_morph, [inputs.ls_loc15], path+'/location_label.svg')
+    label_image(inputs.label_morph, [inputs.reg_rad36, inputs.ls_distal], path+'/distal_label.svg')
+    label_image(inputs.label_morph, [inputs.reg_rad36, inputs.ls_proximal], path+'/proximal_label.svg')
+    label_image(inputs.label_morph, [inputs.ls_uniform0, inputs.ls_uniform1], path+'/uniform_label.svg')
+    label_image(inputs.label_morph, [inputs.ls_branchmid], path+'/on_branches_label.svg')
+    label_image(inputs.label_morph, [inputs.ls_term, inputs.reg_tag3, inputs.ls_restrict], path+'/restrict_label.svg')
+
+    ####################### regions
+
+    label_image(inputs.label_morph, [inputs.reg_empty, inputs.reg_all], path+'/nil_all_label.svg')
+    label_image(inputs.label_morph, [inputs.reg_tag1, inputs.reg_tag2, inputs.reg_tag3], path+'/tag_label.svg')
+    label_image(inputs.label_morph, [inputs.reg_tag1, inputs.reg_tag3], path+'/tag_label.svg')
+    label_image(inputs.label_morph, [inputs.reg_branch0, inputs.reg_branch3], path+'/branch_label.svg')
+    label_image(inputs.label_morph, [inputs.reg_cable_1_01, inputs.reg_cable_1_31, inputs.reg_cable_1_37], path+'/cable_label.svg')
+    label_image(inputs.label_morph, [inputs.ls_proxint_in, inputs.reg_proxint],    path+'/proxint_label.svg')
+    label_image(inputs.label_morph, [inputs.ls_proxint_in, inputs.reg_proxintinf], path+'/proxintinf_label.svg')
+    label_image(inputs.label_morph, [inputs.ls_distint_in, inputs.reg_distint],    path+'/distint_label.svg')
+    label_image(inputs.label_morph, [inputs.ls_distint_in, inputs.reg_distintinf], path+'/distintinf_label.svg')
+    label_image(inputs.label_morph, [inputs.reg_lhs, inputs.reg_rhs, inputs.reg_or],  path+'/union_label.svg')
+    label_image(inputs.label_morph, [inputs.reg_lhs, inputs.reg_rhs, inputs.reg_and], path+'/intersect_label.svg')
+    label_image(inputs.label_morph, [inputs.reg_radlt5],  path+'/radiuslt_label.svg')
+    label_image(inputs.label_morph, [inputs.reg_radle5],  path+'/radiusle_label.svg')
+    label_image(inputs.label_morph, [inputs.reg_radgt5],  path+'/radiusgt_label.svg')
+    label_image(inputs.label_morph, [inputs.reg_radge5],  path+'/radiusge_label.svg')
+
+
+if __name__ == '__main__':
+    generate('.')
diff --git a/doc/scripts/representation.py b/doc/scripts/representation.py
new file mode 100644
index 0000000000000000000000000000000000000000..e29e9de5afc9cb0ce4b0297dc480775b3ca46b52
--- /dev/null
+++ b/doc/scripts/representation.py
@@ -0,0 +1,184 @@
+import math
+
+def add_vec(u, v):
+    return (u[0]+v[0], u[1]+v[1])
+
+def norm_vec(v):
+    return math.sqrt(v[0]*v[0] + v[1]*v[1])
+
+def sub_vec(u, v):
+    return (u[0]-v[0], u[1]-v[1])
+
+def rot90_vec(v):
+    return (v[1], -v[0])
+
+def scal_vec(alpha, v):
+    return (alpha*v[0], alpha*v[1])
+
+def unit_vec(v):
+    L = norm_vec(v)
+    return (v[0]/L, v[1]/L)
+
+def is_collocated(x, y):
+    return x[0]==y[0] and x[1]==y[1]
+
+class Segment:
+    def __init__(self, prox, dist, tag):
+        self.prox = prox
+        self.dist = dist
+        self.tag = tag
+        self.length = norm_vec(sub_vec(prox, dist))
+
+    def location(self, pos):
+        b = self.prox
+        e = self.dist
+        return (b[0]+pos*(e[0]-b[0]), b[1]+pos*(e[1]-b[1]), b[2]+pos*(e[2]-b[2]))
+
+    def corners(self, prox=0, dist=1):
+        b = self.location(prox)
+        e = self.location(dist)
+        rb = b[2]
+        re = e[2]
+        d = sub_vec(self.dist, self.prox)
+        o = rot90_vec(unit_vec(d))
+        p1 = add_vec(b, scal_vec(rb, o))
+        p2 = add_vec(e, scal_vec(re, o))
+        p3 = sub_vec(e, scal_vec(re, o))
+        p4 = sub_vec(b, scal_vec(rb, o))
+        return [p1,p2,p3,p4]
+
+    def __str__(self):
+        return 'seg({}, {})'.format(self.prox, self.dist)
+
+    def __repr__(self):
+        return 'seg({}, {})'.format(self.prox, self.dist)
+
+
+# Represent and query a cable cell branch for rendering.
+# The branch is composed of segments, which are grouped into "sections".
+# A section is a fully connected sequence of segments, and a branch
+# will have more than one section if there are jumps/gaps between
+# segments.
+# The section representation makes things a bit messy, but it is required
+# to be able to draw morphologies with gaps.
+
+class Branch:
+
+    def __init__(self, sections):
+        self.sections = sections
+        length = 0
+        for sec in self.sections:
+            for seg in sec:
+                length += seg.length
+        self.length = length
+
+    def minmax(self):
+        minx = math.inf
+        miny = math.inf
+        maxx = -math.inf
+        maxy = -math.inf
+        for sec in self.sections:
+            for seg in sec:
+                px, py, pr = seg.prox
+                dx, dy, dr = seg.dist
+                minx = min(minx, px-pr, dx-dr)
+                miny = min(miny, py-pr, dy-dr)
+                maxx = max(maxx, px+pr, dx+dr)
+                maxy = max(maxy, py+pr, dy+dr)
+
+        return (minx, maxx, miny, maxy)
+
+    # Return the segment location that contains the location
+    # that is 0 ≤ pos ≤ 1 along the branch
+    #
+    # return tuple: (sec, seg, pos)
+    #       sec: id of the section containing the segment
+    #       seg: index of the segment in the section
+    #       pos: relative position of the location inside the segment
+    def segment_id(self, pos):
+        assert(pos>=0 and pos<=1)
+        if pos==0:
+            return (0, 0, 0.0)
+        if pos==1:
+            return (len(self.sections)-1, len(self.sections[-1])-1, 1.0)
+        l = pos * self.length
+
+        part = 0
+        for secid in range(len(self.sections)):
+            sec = self.sections[secid]
+            for segid in range(len(sec)):
+                seg = sec[segid]
+                if part+seg.length >= l:
+                    segpos = (l-part)/seg.length
+                    return (secid, segid, segpos)
+
+                part += seg.length
+
+    def location(self, pos):
+        assert(pos>=0 and pos<=1)
+
+        secid, segid, segpos = self.segment_id(pos)
+        return self.sections[secid][segid].location(segpos)
+
+    def sec_outline(self, secid, pseg, ppos, dseg, dpos):
+        sec = self.sections[secid]
+
+        # Handle the case where the cable is in one segment
+        if pseg==dseg:
+            assert(ppos<=dpos)
+            return sec[pseg].corners(ppos, dpos)
+
+        left = []
+        right = []
+
+        # Handle the partial partial proximal segment
+        p1, p2, p3, p4 = sec[pseg].corners(ppos, 1)
+        left += [p1, p2]
+        right += [p4, p3]
+
+        # Handle the full segments in the middle
+        for segid in range(pseg+1, dseg):
+            p1, p2, p3, p4 = sec[segid].corners()
+            left += [p1, p2]
+            right += [p4, p3]
+
+        # Handle the partial distal segment
+        p1, p2, p3, p4 = sec[dseg].corners(0, dpos)
+        left += [p1, p2]
+        right += [p4, p3]
+
+        right.reverse()
+        return left + right
+
+
+
+    # Return outline of all (sub)sections in the branch between the relative
+    # locations: 0 ≤ prox ≤ dist ≤ 1
+    def outline(self, prox=0, dist=1):
+        psec, pseg, ppos = self.segment_id(prox)
+        dsec, dseg, dpos = self.segment_id(dist)
+
+        if psec==dsec and pseg==dseg:
+            return [self.sections[psec][pseg].corners(ppos, dpos)]
+        if psec==dsec:
+            return [self.sec_outline(psec, pseg, ppos, dseg, dpos)]
+
+        outlines = [self.sec_outline(psec, pseg, ppos, len(self.sections[psec])-1, 1)]
+        for secid in range(psec+1,dsec):
+            outlines.append(self.sec_outline(secid, 0, 0, len(self.sections[secid])-1, 1))
+        outlines.append(self.sec_outline(dsec, 0, 0, dseg, dpos))
+
+        return outlines
+
+# A morphology for rendering is a flat list of branches, with no
+# parent-child information for the branches.
+# Each branch is itself a list of sections, where each section
+# represents a sequence of segments with no gaps.
+
+def make_morph(branches):
+    m = []
+    for branch_sections in branches:
+        m.append(Branch(branch_sections))
+
+    return m
+
diff --git a/doc/single_cell.rst b/doc/single_cell.rst
index 94de6c7306a6eb8d8d33335ba1a3ad263861fd64..a92cab5b22f1eb5b98599b94e459f0f65c9f908f 100644
--- a/doc/single_cell.rst
+++ b/doc/single_cell.rst
@@ -1,31 +1,98 @@
-.. _single_cell:
+.. _single:
 
 Single Cell Models
 ==================
 
-Arbor breaks down building *single cell models* into the following steps:
+Building and testing detailed models of individual cells, then optimizing their parameters
+is usually the first step in building models with multi-compartment cells.
+Arbor supports a *single cell model* workflow for this purpose, which is a good way to
+introduce Arbor's cell modeling concepts and approach.
 
-1. Defining the `morphology <single_morpho_>`_ of the cell.
-2. Labeling regions and locations on the morphology.
-3. Defining the mechanisms that will be applied to the cell.
-4. Adding ion channels and synapses (mechanisms) to labeled regions and locations.
-5. Attaching stimuli, spike detectors, event generators and probes to locations (inputs & outputs).
+This guide will walk through a series of single cell models of increasing
+complexity that the reader is encouraged to follow along with in Python. Links
+are provide to separate documentation that covers relevant topics in more detail.
 
-This guide will provide a simple example workflow of this kind that the reader
-can follow along in Python, with links to separate documentation that covers
-relevant topics in more detail.
+.. _single_soma:
 
-.. note::
-    Most readers will be familiar with NEURON. Boxes like this
-    will be used to highlight differences between NEURON and Arbor
-    throughout the guide.
+Example 1: Single compartment cell with HH dynamics
+----------------------------------------------------
 
-    NEURON users will recognise that Arbor uses many similar concepts, and
-    an effort has been made to use the same nomenclature wherever possible.
+The most trivial representation of a cell in Arbor is to model the entire cell as a sphere.
+The following model shows the steps required to construct a model of a spherical cell with
+radius 3 μm, Hodgkin–Huxley dynamics and a current clamp stimulus, then run the model for
+30 ms.
 
-    Arbor takes a more structured approach to model building,
-    from morphology descriptions up to network connectivity, to allow model
-    descriptions that are more scalable and portable.
+The first step is to construct the cell. In Arbor, the abstract representation used to define
+a cell with branching "cable" morphology is a ``cable_cell``, which holds a description
+of the cell's morphology, named regions and locations on the morphology, and descriptions of
+ion channels, synapses, spike detectors and electrical properties.
+
+Our "single-compartment HH cell" has a trivial morphology and dynamics, so the steps to
+create the ``cable_cell`` that represents it are quite straightforward:
+
+.. code-block:: python
+
+    import arbor
+
+    # (1) Create a morphology with a single segment of length=diameter=6 μm
+    tree = arbor.segment_tree()
+    tree.append(arbor.mnpos, arbor.mpoint(-3, 0, 0, 3), arbor.mpoint(3, 0, 0, 3), tag=1)
+
+    # (2) Define the soma and its center
+    labels = arbor.label_dict({'soma':   '(tag 1)',
+                               'center': '(location 0 0.5)'})
+
+    # (3) Create cell and set properties
+    cell = arbor.cable_cell(tree, labels)
+    cell.set_properties(Vm=-40)
+    cell.paint('soma', 'hh')
+    cell.place('center', arbor.iclamp( 10, 2, 0.8))
+    cell.place('center', arbor.spike_detector(-10))
+
+Arbor's cell morphologies are constructed from a :ref:`segment tree<morph-segment_tree>`,
+which is a list of segments, which are tapered cones with a *tag*.
+Step **(1)** above shows how the spherical cell is represented using a single segment.
+
+Cell builders need to refer to *regions* and *locations* on a cell morphology.
+Arbor uses a domains specific language (DSL) to describe regions and locations,
+which are given labels. In step **(2)** a dictionary of labels is created
+with two labels:
+
+* ``soma`` defines a *region* with ``(tag  2)``. Note that this corresponds to the ``tag`` parameter that was used to define the single segment in step (1).
+* ``center`` defines a *location* at ``(location 0 0.5)``, which is the mid point ``0.5`` of branch ``0``, which corresponds to the center of the soma on the morphology defined in Step (1).
+
+In step **(3)** the cable cell is constructed by combining the segment tree with
+the named regions and locations.
+
+* Set initial membrane potential everywhere on the cell to -40 mV.
+* Use HH dynamics on soma.
+* Attach stimuli with duration of 2 ms and current of 0.8 nA.
+* Add a spike detector with threshold of -10 mV.
+
+Arbor can simulate networks with multiple individual cells, connected together in a network.
+Single cell models do not require the full *recipe* interface used to describing such
+network models, with many unique cells, network and gap junctions.
+Arbor provides a ``single_cell_model`` helper that wraps a cell description, and provides
+an interface for recording 
+
+
+.. code-block:: python
+
+    # (4) Make single cell model.
+    m = arbor.single_cell_model(cell)
+
+    # (5) Attach voltage probe sampling at 10 kHz (every 0.1 ms).
+    m.probe('voltage', 'center', frequency=10000)
+
+    # (6) Run simulation for 100 ms of simulated activity.
+    m.run(tfinal=100)
+
+Step (4) instantiates the single cell model using our single-compartment cell.
+To record variables from the model three pieces of information are prov
+* 
+
+
+**Everything below here is to be discarded/moved**
 
 .. _single_morpho:
 
@@ -35,22 +102,8 @@ Morphology
 The first step in building a cell model is to define the cell's *morphology*.
 Conceptually, Arbor treats morphologies as a tree of truncated frustums, with
 an optional spherical segment at the root of the tree.
-These are represented as a tree of sample points, where each sample has a 3D location,
-a radius, and a tag, and a parent sample.
-
-.. note::
-    NEURON represents morphologies as a tree of cylindrical *segments*, whereas
-    in Arbor the radius can vary linearly between two sample locations.
-
-    A cylinder with equal diameter and length is used to model spherical somata
-    in NEURON, which coincidently has the same surface area as a sphere of the same diameter.
-    Arbor allows the user to optionally use a spherical section at the root
-    of the tree to represent spherical somata.
-
-.. note::
-    In NEURON cell morphologies are constructed by creating individual sections,
-    then connecting them together. In Arbor we start with an "empty"
-    sample tree, to which samples are appended to build a connected morphology.
+These are represented as a tree of segments, where each segment is defined
+by two end points with radius, and a tag, and a parent segment to which it is attached.
 
 Let's start with a simple "ball and stick" model cell.
 
@@ -74,7 +127,7 @@ Let's start with a simple "ball and stick" model cell.
 
 
 Building the morphology there are two approaches: construct it manually using
-``sample_tree`` or ``flat_cell_builder``, or load from swc file.
+``segment_tree`` or ``flat_cell_builder``, or load from swc file.
 
 TODO: cover all methods here?
     - we could just ``flat_cell_builder`` because it is most comfortable for
@@ -84,7 +137,43 @@ TODO: cover all methods here?
       complicated, well-defined, morphology, and illustrate how to build
       it using all of the different methods.
 
-Labeling Regions and Locations
+NEURON erratum
 ------------------------------
 
+These should probably be combined into a single section that describes the differences
+between Arbor and NEURON, because the alternative is to keep stopping the
+narative to point out the difference with NEURON, instead of explaining what
+Arbor is from a fresh start.
+
+.. note::
+    Most readers will be familiar with NEURON. Boxes like this
+    will be used to highlight differences between NEURON and Arbor
+    throughout the guide.
+
+    NEURON users will recognise that Arbor uses many similar concepts, and
+    an effort has been made to use the same nomenclature wherever possible.
+
+    Arbor takes a more structured approach to model building,
+    from morphology descriptions up to network connectivity, to allow model
+    descriptions that are more scalable and portable.
+
+.. note::
+    NEURON represents morphologies as a tree of cylindrical *segments*, whereas
+    in Arbor the radius can vary linearly along segments.
+
+    A cylinder with equal diameter and length is used to model spherical somata
+    in NEURON, which coincidently has the same surface area as a sphere of the same diameter.
+    Arbor allows the user to optionally use a spherical section at the root
+    of the tree to represent spherical somata.
+
+.. note::
+    In NEURON cell morphologies are constructed by creating individual sections,
+    then connecting them together. In Arbor we start with an "empty"
+    segment tree, to which segments are appended to build a connected morphology.
+
+1. Defining the `morphology <single_morpho_>`_ of the cell.
+2. Labeling regions and locations on the morphology.
+3. Defining the mechanisms that will be applied to the cell.
+4. Adding ion channels and synapses (mechanisms) to labeled regions and locations.
+5. Attaching stimuli, spike detectors, event generators and probes to locations (inputs & outputs).
 
diff --git a/doc/static/custom.css b/doc/static/custom.css
index 5f9ebf31bc22939d67e9488cd1ad1d6bfb1c3027..2430843e70c6d2a3043d20a48cf268cfa46b8c43 100644
--- a/doc/static/custom.css
+++ b/doc/static/custom.css
@@ -29,6 +29,11 @@ table.docutils td {
     vertical-align: top !important;
 }
 
+em {
+  font-style: bold;
+  color: #000099;
+}
+
 @media screen and (min-width: 40em) {
     .wy-table-responsive table td,
     .wy-table-responsive table th {
diff --git a/python/morphology.cpp b/python/morphology.cpp
index c5b3d20920f0d5f313145bed8b534d1e74a53637..30be0372bea136e8eec2f7783e67bd529503aea8 100644
--- a/python/morphology.cpp
+++ b/python/morphology.cpp
@@ -218,13 +218,13 @@ void register_morphology(pybind11::module& m) {
                 "Append a segment to the tree, using the distal location of the parent segment as the proximal end.")
         // properties
         .def_property_readonly("empty", [](const arb::segment_tree& st){return st.empty();},
-                "Indicates whether the sample tree is empty (i.e. whether it has size 0)")
+                "Indicates whether the tree is empty (i.e. whether it has size 0)")
         .def_property_readonly("size", [](const arb::segment_tree& st){return st.size();},
-                "The number of samples in the sample tree.")
+                "The number of segments in the tree.")
         .def_property_readonly("parents", [](const arb::segment_tree& st){return st.parents();},
-                "A list with the parent index of each sample.")
+                "A list with the parent index of each segment.")
         .def_property_readonly("segments", [](const arb::segment_tree& st){return st.segments();},
-                "A list of the samples.")
+                "A list of the segments.")
         .def("__str__", [](const arb::segment_tree& s) {
                 return util::pprintf("<arbor.segment_tree:\n{}>", s);});
 
@@ -247,7 +247,7 @@ void register_morphology(pybind11::module& m) {
                                   e.line_number, fname, e.what()));
             }
         },
-        "Load an swc file and convert to a segment_tree.");
+        "Load an swc file and as a segment_tree.");
 
     m.def("load_swc_allen", &load_swc_allen,
             "filename"_a, "no_gaps"_a=false,
diff --git a/python/tutorial/example1.py b/python/tutorial/example1.py
index 5e9b49cdb2767e42993326b3bf49ea4d17719baf..1512803e95a029a99fa06670a6837a188eb86f21 100644
--- a/python/tutorial/example1.py
+++ b/python/tutorial/example1.py
@@ -1,10 +1,10 @@
 import arbor
 
-# Create a sample tree with a single sample of radius 3 μm
-tree = arbor.sample_tree()
-tree.append(arbor.sample(x=0, y=0, z=0, radius=3, tag=2))
+# Create a morphology with a single segment of length=diameter=5 μm
+tree = arbor.segment_tree()
+tree.append(arbor.mnpos, arbor.mpoint(-3, 0, 0, 3), arbor.mpoint(3, 0, 0, 3), tag=1)
 
-labels = arbor.label_dict({'soma': '(tag 2)', 'center': '(location 0 0.5)'})
+labels = arbor.label_dict({'soma': '(tag 1)', 'center': '(location 0 0.5)'})
 
 cell = arbor.cable_cell(tree, labels)
 
@@ -12,8 +12,8 @@ cell = arbor.cable_cell(tree, labels)
 cell.set_properties(Vm=-40)
 # Put hh dynamics on soma, and passive properties on the dendrites.
 cell.paint('soma', 'hh')
-# Attach stimuli with duration of 1 ms and current of 0.8 nA.
-cell.place('center', arbor.iclamp( 10, duration=1, current=0.8))
+# 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))
 # Add a spike detector with threshold of -10 mV.
 cell.place('center', arbor.spike_detector(-10))
 
@@ -45,7 +45,7 @@ legend_labels = ['{}: {}'.format(s.variable, s.location) for s in m.traces]
 ax.legend(legend_labels)
 ax.set(xlabel='time (ms)', ylabel='voltage (mV)', title='cell builder demo')
 plt.xlim(0,tfinal)
-plt.ylim(-80,80)
+plt.ylim(-80,50)
 ax.grid()
 
 # Set to True to save the image to file instead of opening a plot window.