diff --git a/doc/contrib/doc.rst b/doc/contrib/doc.rst
index 180162407e5cd5b7a20cd88de95aa766f7f0561b..a77def2bd420986de2e1e845577fb5c843a2085c 100644
--- a/doc/contrib/doc.rst
+++ b/doc/contrib/doc.rst
@@ -1,16 +1,22 @@
 .. _contribdoc:
 
 Documentation
-=====================
+=============
 
 The source for `the documentation <https://docs.arbor-sim.org>`__ is
 found in the ``/doc`` subdirectory. You can add your contribution to the documentation
 in the same way you would contribute code, please see the :ref:`contribpr` section.
 
-.. todo::
-    This page is under construction.
+.. _contribdoc-tut:
 
-.. _contribdoc-update:
+Tutorials
+---------
+
+Documentation includes tutorials, which are :ref:`examples <contribexample>` with an associated
+page under ``/doc/tutorial`` explaining the goal of the example and guiding the reader through the code.
+We try to use the `literalinclude directive <https://www.sphinx-doc.org/en/master/usage/restructuredtext/directives.html#directive-literalinclude>`
+wherever possible to ensure the tutorial and example don't diverge.
+Remember to update the line numbers whenever you update the examples in the tutorial.
 
 Update policy
 -------------
diff --git a/doc/contrib/example.rst b/doc/contrib/example.rst
index 423ac49ca17d72041dede2451cd0fbb031ef2eec..5ca0601decb8544eb6102b624149aa40ad67def3 100644
--- a/doc/contrib/example.rst
+++ b/doc/contrib/example.rst
@@ -1,11 +1,19 @@
 .. _contribexample:
 
 Examples
-===============
+========
 
 The source for C++ examples are located in ``/example`` and
 Python examples in ``/python/example``. You can add yours in the same
 way you would contribute code, please see the :ref:`contribpr` section.
 
-.. todo::
-    This page is under construction.
+Especially when adding new functionality, consider adding an example to demonstrate its use.
+
+.. _contribexample-tut:
+
+Tutorials
+---------
+
+Tutorials are examples with an associated :ref:`documentation <contribdoc>` under ``/doc/tutorial``
+explaining the goal of the example and guiding the reader through the code. When you add an example,
+consider writing a tutorial page.
diff --git a/doc/tutorial/index.rst b/doc/tutorial/index.rst
index b9fbbc59296c33c805bb8e76c1c9ebb88d3480b8..c0418d9cd3bc3f3d148806fffac52c54db9466fc 100644
--- a/doc/tutorial/index.rst
+++ b/doc/tutorial/index.rst
@@ -26,4 +26,4 @@ Tutorials
    single_cell_detailed_recipe
    single_cell_cable
    network_ring
-   mpi
+   network_ring_mpi
diff --git a/doc/tutorial/network_ring.rst b/doc/tutorial/network_ring.rst
index f395d695ed7d09c94061c8b43747905a2b3b4034..ab4a341c50ad4b8061a05af6a3414fc8694bcfdc 100644
--- a/doc/tutorial/network_ring.rst
+++ b/doc/tutorial/network_ring.rst
@@ -16,7 +16,8 @@ In this example, a small *network* of cells, arranged in a ring, will be created
 The cell
 ********
 
-Step **(1)** shows how a simple cell with a dendrite is created. We construct the following :term:`morphology` and label the soma and dendrite:
+Step **(1)** shows a function that creates a simple cell with a dendrite.
+We construct the following :term:`morphology` and label the soma and dendrite:
 
 .. figure:: ../gen-images/tutorial_network_ring_morph.svg
    :width: 400
@@ -24,27 +25,9 @@ Step **(1)** shows how a simple cell with a dendrite is created. We construct th
 
    A 4-segment cell with a soma (pink) and a branched dendrite (light blue).
 
-.. code-block:: python
-
-   # (1) Build a segment tree
-   tree = arbor.segment_tree()
-
-   # Soma (tag=1) with radius 6 μm, modelled as cylinder of length 2*radius
-   s = tree.append(arbor.mnpos, arbor.mpoint(-12, 0, 0, 6), arbor.mpoint(0, 0, 0, 6), tag=1)
-
-   # Single dendrite (tag=3) of length 50 μm and radius 2 μm attached to soma.
-   b0 = tree.append(s, arbor.mpoint(0, 0, 0, 2), arbor.mpoint(50, 0, 0, 2), tag=3)
-
-   # Attach two dendrites (tag=3) of length 50 μm to the end of the first dendrite.
-   # Radius tapers from 2 to 0.5 μm over the length of the dendrite.
-   b1 = tree.append(b0, arbor.mpoint(50, 0, 0, 2), arbor.mpoint(50+50/sqrt(2), 50/sqrt(2), 0, 0.5), tag=3)
-   # Constant radius of 1 μm over the length of the dendrite.
-   b2 = tree.append(b0, arbor.mpoint(50, 0, 0, 1), arbor.mpoint(50+50/sqrt(2), -50/sqrt(2), 0, 1), tag=3)
-
-   # Associate labels to tags
-   labels = arbor.label_dict()
-   labels['soma'] = '(tag 1)'
-   labels['dend'] = '(tag 3)'
+.. literalinclude:: ../../python/example/network_ring.py
+   :language: python
+   :lines: 18-37
 
 In step **(2)** we create a :term:`label` for both the root and the site of the synapse.
 These locations will form the endpoints of the connections between the cells.
@@ -55,12 +38,9 @@ These locations will form the endpoints of the connections between the cells.
 
    We'll create labels for the root (red) and a synapse_site (black).
 
-.. code-block:: python
-
-   # (2) Mark location for synapse at the midpoint of branch 1 (the first dendrite).
-   labels['synapse_site'] = '(location 1 0.5)'
-   # Mark the root of the tree.
-   labels['root'] = '(root)'
+.. literalinclude:: ../../python/example/network_ring.py
+   :language: python
+   :lines: 39-42
 
 After we've created a basic :py:class:`arbor.decor`, step **(3)** places a synapse with an exponential decay (``'expsyn'``) on the ``'synapse_site'``.
 The synapse is given the label ``'syn'``, which is later used to form :py:class:`arbor.connection` objects terminating *at* the cell.
@@ -76,21 +56,9 @@ Step **(4)** places a spike detector at the ``'root'``. The detector is given th
 
    The same explanation applies to the number of detectors on this cell.
 
-.. code-block:: python
-
-   decor = arbor.decor()
-
-   # Put hh dynamics on soma, and passive properties on the dendrites.
-   decor.paint('"soma"', 'hh')
-   decor.paint('"dend"', 'pas')
-
-   # (3) Attach a single synapse, label it 'syn'
-   decor.place('"synapse_site"', 'expsyn', 'syn')
-
-   # (4) Attach a spike detector with threshold of -10 mV.
-   decor.place('"root"', arbor.spike_detector(-10), 'detector')
-
-   cell = arbor.cable_cell(tree, labels, decor)
+.. literalinclude:: ../../python/example/network_ring.py
+   :language: python
+   :lines: 44-59
 
 The recipe
 **********
@@ -132,57 +100,9 @@ Step **(10)** places a :term:`probe` at the ``"root"`` of each cell.
 
 Step **(11)** instantiates the recipe with 4 cells.
 
-.. code-block:: python
-
-   # (5) Create a recipe that generates a network of connected cells.
-   class ring_recipe (arbor.recipe):
-      def __init__(self, ncells):
-         # The base C++ class constructor must be called first, to ensure that
-         # all memory in the C++ class is initialized correctly.
-         arbor.recipe.__init__(self)
-         self.ncells = ncells
-         self.props = arbor.neuron_cable_properties()
-         self.cat = arbor.default_catalogue()
-         self.props.register(self.cat)
-
-      # (6) The num_cells method that returns the total number of cells in the model
-      # must be implemented.
-      def num_cells(self):
-         return self.ncells
-
-      # (7) The cell_description method returns a cell
-      def cell_description(self, gid):
-         return make_cable_cell(gid)
-
-      # The kind method returns the type of cell with gid.
-      # Note: this must agree with the type returned by cell_description.
-      def cell_kind(self, gid):
-         return arbor.cell_kind.cable
-
-      # (8) Make a ring network. For each gid, provide a list of incoming connections.
-      def connections_on(self, gid):
-         src = (gid-1)%self.ncells
-         w = 0.01
-         d = 5
-         return [arbor.connection((src,'detector'), 'syn', w, d)]
-
-      # (9) Attach a generator to the first cell in the ring.
-      def event_generators(self, gid):
-         if gid==0:
-               sched = arbor.explicit_schedule([1])
-               return [arbor.event_generator('syn', 0.1, sched)]
-         return []
-
-      # (10) Place a probe at the root of each cell.
-      def probes(self, gid):
-         return [arbor.cable_probe_membrane_voltage('"root"')]
-
-      def global_properties(self, kind):
-         return self.props
-
-   # (11) Instantiate recipe
-   ncells = 4
-   recipe = ring_recipe(ncells)
+.. literalinclude:: ../../python/example/network_ring.py
+   :language: python
+   :lines: 61-111
 
 The execution
 *************
@@ -208,34 +128,18 @@ these handles for later use.
 
 Step **(15)** executes the simulation for a duration of 100 ms.
 
-.. code-block:: python
-
-   # (12) Create a default execution context, domain decomposition and simulation
-   context = arbor.context()
-   decomp = arbor.partition_load_balance(recipe, context)
-   sim = arbor.simulation(recipe, decomp, context)
-
-   # (13) Set spike generators to record
-   sim.record(arbor.spike_recording.all)
-
-   # (14) Attach a sampler to the voltage probe on cell 0. Sample rate of 10 sample every ms.
-   handles = [sim.sample((gid, 0), arbor.regular_schedule(0.1)) for gid in range(ncells)]
-
-   # (15) Run simulation
-   sim.run(100)
-   print('Simulation finished')
+.. literalinclude:: ../../python/example/network_ring.py
+   :language: python
+   :lines: 113-126
 
 The results
 ***********
 
 Step **(16)** prints the timestamps of the spikes:
 
-.. code-block:: python
-
-   # Print spike times
-   print('spikes:')
-   for sp in sim.spikes():
-      print(' ', sp)
+.. literalinclude:: ../../python/example/network_ring.py
+   :language: python
+   :lines: 128-131
 
 Step **(17)** generates a plot of the sampling data.
 :py:func:`arbor.simulation.samples` takes a ``handle`` of the probe we wish to examine. It returns a list
@@ -245,18 +149,9 @@ discrete locations pointed to by the handle, which in this case is 1, so we can
 (Recall that in step **(10)** we attached a probe to the ``"root"``, which describes one location.
 It could have described a :term:`locset`.)
 
-.. code-block:: python
-
-   # Plot the recorded voltages over time.
-   print("Plotting results ...")
-   df_list = []
-   for gid in range(ncells):
-      samples, meta = sim.samples(handles[gid])[0]
-      df_list.append(pandas.DataFrame({'t/ms': samples[:, 0], 'U/mV': samples[:, 1], 'Cell': f"cell {gid}"}))
-
-   df = pandas.concat(df_list)
-   seaborn.relplot(data=df, kind="line", x="t/ms", y="U/mV",hue="Cell",ci=None).savefig('network_ring_result.svg')
-
+.. literalinclude:: ../../python/example/network_ring.py
+   :language: python
+   :lines: 133-141
 
 Since we have created ``ncells`` cells, we have ``ncells`` traces. We should be seeing phase shifted traces, as the action potential propagated through the network.
 
@@ -266,7 +161,6 @@ We plot the results using pandas and seaborn:
     :width: 400
     :align: center
 
-
 The full code
 *************
 
diff --git a/doc/tutorial/mpi.rst b/doc/tutorial/network_ring_mpi.rst
similarity index 69%
rename from doc/tutorial/mpi.rst
rename to doc/tutorial/network_ring_mpi.rst
index 236a5af07d9c523b28b0c5ec49fa4ad1fcfd076a..91093e1004f2baacff08ff0abf339fdbb741eb83 100644
--- a/doc/tutorial/mpi.rst
+++ b/doc/tutorial/network_ring_mpi.rst
@@ -19,11 +19,9 @@ The recipe
 
 Step **(11)** is changed to generate a network with five hundred cells.
 
-.. code-block:: python
-
-   # (11) Instantiate recipe
-   ncells = 500
-   recipe = ring_recipe(ncells)
+.. literalinclude:: ../../python/example/network_ring_mpi.py
+   :language: python
+   :lines: 111-113
 
 The hardware context
 ********************
@@ -33,25 +31,18 @@ Step **(12)** uses the Arbor-built-in :py:class:`MPI communicator <arbor.mpi_com
 communicator for its ``mpi`` parameter. Note that you can also pass in communicators created with ``mpi4py``.
 We print both the communicator and context to observe how Arbor configures their defaults.
 
-.. code-block:: python
-
-   # (12) Create an MPI communicator, and use it to create a hardware context
-   arbor.mpi_init()
-   comm = arbor.mpi_comm()
-   print(comm)
-   context = arbor.context(mpi=comm)
-   print(context)
+.. literalinclude:: ../../python/example/network_ring_mpi.py
+   :language: python
+   :lines: 115-120
 
 The execution
 *************
 
-Step **(16)** runs the simulation. Since we have more cells this time, which are connected in series, it will take some time for the action potential to propagate. In the :ref:`ring network <tutorialnetworkring>` we could see it takes about 5 ms for the signal to propgate through one cell, so let's set the runtime to ``5*ncells``.
-
-.. code-block:: python
+Step **(16)** runs the simulation. Since we have more cells this time, which are connected in series, it will take some time for the action potential to propagate. In the :ref:`ring network <tutorialnetworkring>` we could see it takes about 5 ms for the signal to propagate through one cell, so let's set the runtime to ``5*ncells``.
 
-   # (16) Run simulation
-   sim.run(ncells*5)
-   print('Simulation finished')
+.. literalinclude:: ../../python/example/network_ring_mpi.py
+   :language: python
+   :lines: 133-135
 
 An important change in the execution is how the script is run. Whereas normally you run the Python script by passing
 it as an argument to the ``python`` command, you need to use ``srun`` or ``mpirun`` (depending on your MPI
@@ -83,44 +74,22 @@ length. The effect is that we collect the results generated on this particular n
 instances of our script, and we can't access the results between nodes, we have to write the results to disk and
 analyse them later. We query :py:attr:`arbor.context.rank` to generate a unique filename for the result.
 
-.. code-block:: python
-
-   # (18) Store the recorded voltages
-   print("Storing results ...")
-   df_list = []
-   for gid in range(ncells):
-      if len(sim.samples(handles[gid])):
-         samples, meta = sim.samples(handles[gid])[0]
-         df_list.append(pandas.DataFrame({'t/ms': samples[:, 0], 'U/mV': samples[:, 1], 'Cell': f"cell {gid}"}))
-
-   if len(df_list):
-      df = pandas.concat(df_list)
-      df.to_csv(f"result_mpi_{context.rank}.csv", float_format='%g')
+.. literalinclude:: ../../python/example/network_ring_mpi.py
+   :language: python
+   :lines: 137-147
 
-In a second script, ``mpi_plot.py``, we load the results stored to disk into a pandas table, and plot the concatenated table as before:
+In a second script, ``network_ring_mpi_plot.py``, we load the results stored to disk into a pandas table, and plot the concatenated table as before:
 
-.. code-block:: python
-
-   import glob
-   import pandas, seaborn
-
-   results = glob.glob("result_mpi_*.csv")
-
-   df_list = []
-   for result in results:
-      df_list.append(pandas.read_csv(result))
-
-   df = pandas.concat(df_list)
-   seaborn.relplot(data=df, kind="line", x="t/ms", y="U/mV",hue="Cell",ci=None).savefig('mpi_result.svg')
+.. literalinclude:: ../../python/example/network_ring_mpi_plot.py
+   :language: python
 
 To avoid an overcrowded plot, this plot was generated with just 50 cells:
 
-.. figure:: mpi_result.svg
+.. figure:: network_ring_mpi_result.svg
     :width: 400
     :align: center
 
-
 The full code
 *************
 
-You can find the full code of the example at ``python/examples/mpi.py`` and ``python/examples/mpi_plot.py``.
+You can find the full code of the example at ``python/examples/network_ring_mpi.py`` and ``python/examples/network_ring_mpi_plot.py``.
diff --git a/doc/tutorial/mpi_result.svg b/doc/tutorial/network_ring_mpi_result.svg
similarity index 100%
rename from doc/tutorial/mpi_result.svg
rename to doc/tutorial/network_ring_mpi_result.svg
diff --git a/doc/tutorial/single_cell_cable.rst b/doc/tutorial/single_cell_cable.rst
index ccd6e9ef2aa2159a285c3976799ba84a7d833b3e..2926652f759bf75d70a2bb1a4a2a8d58dcd798fb 100644
--- a/doc/tutorial/single_cell_cable.rst
+++ b/doc/tutorial/single_cell_cable.rst
@@ -25,7 +25,7 @@ Executing the script will run the simulation with default parameters.
 
 
 Walk-through
-**********
+************
 
 We set up a recipe for the simulation of a single dendrite with some parameters.
 
diff --git a/doc/tutorial/single_cell_detailed.rst b/doc/tutorial/single_cell_detailed.rst
index 52207ef06bca47b9788408eb2c883f64dfad844e..7d24460b96f288665d5f7da2cf5fa6c10f64da2a 100644
--- a/doc/tutorial/single_cell_detailed.rst
+++ b/doc/tutorial/single_cell_detailed.rst
@@ -18,8 +18,13 @@ complex single cell model, and go through the process in more detail.
    6. Setting and overriding model and cell parameters.
    7. Running a simulation and visualising the results using a :class:`arbor.single_cell_model`.
 
+.. _tutorialsinglecellswc-cell:
+
+The cell
+********
+
 We start by building the cell. This will be a :ref:`cable cell <cablecell>` with complex
-geometry and dynamics which can be constructed from 3 components:
+geometry and dynamics which is constructed from 3 components:
 
 1. A **morphology** defining the geometry of the cell.
 2. A **label dictionary** storing labelled expressions which define regions and locations of
@@ -28,31 +33,12 @@ geometry and dynamics which can be constructed from 3 components:
    The decor also includes hints about how the cell is to be modelled under the hood, by
    splitting it into discrete control volumes (CV).
 
-Next, we construct a :class:`arbor.single_cell_model`. This model takes care of a lot of details
-behind the scenes: it sets up a recipe (more on recipes :ref:`here <modelrecipe>`), creates
-a simulation object, manages the hardware etc. These details become more important when modelling
-a network of cells, but can be abstracted away when working with single cell networks.
-
-The single cell model has 4 main functions:
-
-1. It holds the **global properties** of the model
-2. It registers **probes** on specific locations on the cell to measure the voltage.
-3. It **runs** the simulation.
-4. It collects **spikes** from spike detectors and voltage **traces** from registered probes.
-
-.. _tutorialsinglecellswc-cell:
-
-The cell
-********
-
-Before creating the actual cell object, we have to create its components.
-
 The morphology
 ^^^^^^^^^^^^^^^
 We begin by constructing the following morphology:
 
 .. figure:: ../gen-images/tutorial_morph.svg
-   :width: 400
+   :width: 600
    :align: center
 
 This can be done by manually building a segment tree:
@@ -94,23 +80,8 @@ The same morphology can be represented using an SWC file (interpreted according
 to :ref:`Arbor's specifications <morph-formats>`). We can save the following in
 ``single_cell_detailed.swc``.
 
-.. code-block:: python
-
-   # id,  tag,      x,      y,      z,      r,    parent
-       1     1     0.0     0.0     0.0     2.0        -1  # seg0 prox / seg9 prox
-       2     1    40.0     0.0     0.0     2.0         1  # seg0 dist
-       3     3    40.0     0.0     0.0     0.8         2  # seg1 prox
-       4     3    80.0     0.0     0.0     0.8         3  # seg1 dist / seg2 prox
-       5     3   120.0    -5.0     0.0     0.8         4  # seg2 dist / seg3 prox
-       6     3   200.0    40.0     0.0     0.4         5  # seg3 dist / seg4 prox
-       7     3   260.0    60.0     0.0     0.2         6  # seg4 dist
-       8     3   120.0    -5.0     0.0     0.5         5  # seg5 prox
-       9     3   190.0   -30.0     0.0     0.5         8  # seg5 dist / seg6 prox / seg7 prox
-      10     4   240.0   -70.0     0.0     0.2         9  # seg6 dist
-      11     4   230.0   -10.0     0.0     0.2         9  # seg7 dist / seg8 prox
-      12     4   360.0   -20.0     0.0     0.2        11  # seg8 dist
-      13     2   -70.0     0.0     0.0     0.4         1  # seg9 dist / seg10 prox
-      14     2  -100.0     0.0     0.0     0.4        13  # seg10 dist
+.. literalinclude:: ../../python/example/single_cell_detailed.swc
+   :language: python
 
 .. Note::
 
@@ -128,7 +99,7 @@ The morphology can then be loaded from ``single_cell_detailed.swc`` in the follo
 
     import arbor
 
-    # Read the morphology from an SWC file
+    # (1) Read the morphology from an SWC file
 
     morph = arbor.load_swc_arbor("single_cell_detailed.swc")
 
@@ -155,17 +126,9 @@ They can correspond to full segments or parts of segments. Our morphology alread
 pre-established regions determined by the ``tag`` parameter of the segments. They are
 defined as follows:
 
-.. code-block:: python
-
-    # Create a label dictionary
-
-    labels = arbor.label_dict()
-
-    # Add labels for tag 1, 2, 3, 4
-    labels['soma'] = '(tag 1)'
-    labels['axon'] = '(tag 2)'
-    labels['dend'] = '(tag 3)'
-    labels['last'] = '(tag 4)'
+.. literalinclude:: ../../python/example/single_cell_detailed.py
+   :language: python
+   :lines: 22-32
 
 This will generate the following regions when applied to the previously defined morphology:
 
@@ -179,13 +142,9 @@ We can also define a region that represents the whole cell; and to make things a
 a region that includes the parts of the morphology that have a radius greater than 1.5 μm. This is done
 in the following way:
 
-.. code-block:: python
-
-    # Add a label for a region that includes the whole morphology
-    labels['all'] = '(all)'
-
-    # Add a label for the parts of the morphology with radius greater than 1.5 μm.
-    labels['gt_1.5'] = '(radius-gt (region "all") 1.5)'
+.. literalinclude:: ../../python/example/single_cell_detailed.py
+   :language: python
+   :lines: 33-36
 
 This will generate the following regions when applied to the previously defined morphology:
 
@@ -201,10 +160,9 @@ segment 9.
 Finally, let's define a region that includes two already defined regions: "last" and "gt_1.5". This can
 be done as follows:
 
-.. code-block:: python
-
-    # Join regions "last" and "gt_1.5"
-    labels['custom'] = '(join (region "last") (region "gt_1.5"))'
+.. literalinclude:: ../../python/example/single_cell_detailed.py
+   :language: python
+   :lines: 37-38
 
 This will generate the following region when applied to the previously defined morphology:
 
@@ -216,11 +174,9 @@ Our label dictionary so far only contains regions. We can also add some **locati
 with a location that is the root of the morphology, and the set of locations that represent all the
 terminal points of the morphology.
 
-.. code-block:: python
-
-    # Add a labels for the root of the morphology and all the terminal points
-    labels['root'] = '(root)'
-    labels['terminal'] = '(terminal)'
+.. literalinclude:: ../../python/example/single_cell_detailed.py
+   :language: python
+   :lines: 42-44
 
 This will generate the following **locsets** (sets of one or more locations) when applied to the
 previously defined morphology:
@@ -235,13 +191,9 @@ To make things more interesting, let's select only the terminal points which bel
 previously defined "custom" region; and, separately, the terminal points which belong to the
 "axon" region:
 
-.. code-block:: python
-
-    # Add a label for the terminal locations in the "custom" region:
-    labels['custom_terminal'] = '(restrict (locset "terminal") (region "custom"))'
-
-    # Add a label for the terminal locations in the "axon" region:
-    labels['axon_terminal'] = '(restrict (locset "terminal") (region "axon"))'
+.. literalinclude:: ../../python/example/single_cell_detailed.py
+   :language: python
+   :lines: 45-48
 
 This will generate the following 2 locsets when applied to the previously defined morphology:
 
@@ -272,15 +224,9 @@ regions. It is in general better to explicitly set all the default properties of
 to avoid the confusion to having simulator-specific default values. This will therefore be our first
 step:
 
-.. code-block:: python
-
-    # Create a decor object
-    decor = arbor.decor()
-
-    # Set the default properties
-    decor.set_property(Vm =-55, tempK=300, rL=35.4, cm=0.01)
-    decor.set_ion('na', int_con=10,   ext_con=140, rev_pot=50, method='nernst/na')
-    decor.set_ion('k',  int_con=54.4, ext_con=2.5, rev_pot=-77)
+.. literalinclude:: ../../python/example/single_cell_detailed.py
+   :language: python
+   :lines: 50-57
 
 We have set the default initial membrane voltage to -55 mV; the default initial
 temperature to 300 K; the default axial resistivity to 35.4 Ω·cm; and the default membrane
@@ -298,28 +244,18 @@ dictionary earlier to be colder, and the initial voltage of the "soma" region to
 We can override the default properties by *painting* new values on the relevant regions using
 :meth:`arbor.decor.paint`.
 
-.. code-block:: python
-
-    # Override default parameters on certain regions
-
-   decor.paint('"custom"', tempK=270)
-   decor.paint('"soma"', Vm=-50)
+.. literalinclude:: ../../python/example/single_cell_detailed.py
+   :language: python
+   :lines: 59-61
 
 With the default and initial values taken care of, we now add some density mechanisms. Let's *paint*
 a *pas* mechanism everywhere on the cell using the previously defined "all" region; an *hh* mechanism
 on the "custom" region; and an *Ih* mechanism on the "dend" region. The *Ih* mechanism is explicitly
 constructed in order to change the default values of its 'gbar' parameter.
 
-
-.. code-block:: python
-
-   # Paint density mechanisms on certain regions
-
-   from arbor import mechanism as mech
-
-   decor.paint('"all"', 'pas')
-   decor.paint('"custom"', 'hh')
-   decor.paint('"dend"',  mech('Ih', {'gbar': 0.001}))
+.. literalinclude:: ../../python/example/single_cell_detailed.py
+   :language: python
+   :lines: 8,62-66
 
 The decor object is also used to *place* stimuli and spike detectors on the cell using :meth:`arbor.decor.place`.
 We place 3 current clamps of 2 nA on the "root" locset defined earlier, starting at time = 10, 30, 50 ms and
@@ -327,14 +263,9 @@ lasting 1ms each. As well as spike detectors on the "axon_terminal" locset for v
 Every placement gets a label. The labels of detectors and synapses are used to form connection from and to them
 in the recipe.
 
-.. code-block:: python
-
-   # Place stimuli and spike detectors on certain locsets
-
-   decor.place('"root"', arbor.iclamp(10, 1, current=2), 'iclamp0')
-   decor.place('"root"', arbor.iclamp(30, 1, current=2), 'iclamp1')
-   decor.place('"root"', arbor.iclamp(50, 1, current=2), 'iclamp2')
-   decor.place('"axon_terminal"', arbor.spike_detector(-10), 'detector')
+.. literalinclude:: ../../python/example/single_cell_detailed.py
+   :language: python
+   :lines: 68-72
 
 .. Note::
 
@@ -352,30 +283,24 @@ The user controls the discretisation using a :class:`arbor.cv_policy`. There are
 choose from, and they can be composed with one another. In this example, we would like the "soma" region
 to be a single CV, and the rest of the morphology to be comprised of CVs with a maximum length of 1 μm:
 
-.. code-block:: python
-
-   # Single CV for the "soma" region
-   soma_policy = arbor.cv_policy_single('"soma"')
-
-   # CVs with max length = 1 μm as default
-   dflt_policy = arbor.cv_policy_max_extent(1.0)
+.. literalinclude:: ../../python/example/single_cell_detailed.py
+   :language: python
+   :lines: 74-81
 
-   # default policy everywhere except the soma
-   policy = dflt_policy | soma_policy
-
-   decor.discretization(policy)
+Finally, we create the cell.
 
+.. literalinclude:: ../../python/example/single_cell_detailed.py
+   :language: python
+   :lines: 83-85
 
 The model
 *********
 
-We begin by constructing an :class:`arbor.single_cell_model` of the cell we just created.
+Having created the cell, we construct an :class:`arbor.single_cell_model`.
 
-.. code-block:: python
-
-   # Construct the model
-
-   model = arbor.single_cell_model(cell)
+.. literalinclude:: ../../python/example/single_cell_detailed.py
+   :language: python
+   :lines: 87-89
 
 The global properties
 ^^^^^^^^^^^^^^^^^^^^^
@@ -415,13 +340,9 @@ model:
 
 .. _tutorialsinglecellswc-gprop:
 
-.. code-block:: python
-
-   # Set the model default properties
-
-   model.properties.set_property(Vm =-65, tempK=300, rL=35.4, cm=0.01)
-   model.properties.set_ion('na', int_con=10,   ext_con=140, rev_pot=50, method='nernst/na')
-   model.properties.set_ion('k',  int_con=54.4, ext_con=2.5, rev_pot=-77)
+.. literalinclude:: ../../python/example/single_cell_detailed.py
+   :language: python
+   :lines: 91-95
 
 We set the same properties as we did earlier when we were creating the *decor* of the cell, except
 for the initial membrane voltage, which is -65 mV as opposed to -55 mV.
@@ -430,14 +351,9 @@ During the decoration step, we also made use of 3 mechanisms: *pas*, *hh* and *I
 the *pas* and *hh* mechanisms are in the default Arbor catalogue, whereas the *Ih* mechanism is in
 the "allen" catalogue. We can extend the default catalogue as follow:
 
-.. code-block:: python
-
-   # Extend the default catalogue with the allen catalogue.
-   # The function takes a second string parameter that can prefix
-   # the name of the mechanisms to avoid collisions between catalogues
-   # in this case we have no collisions so we use an empty prefix string.
-
-   model.catalogue.extend(arbor.allen_catalogue(), "")
+.. literalinclude:: ../../python/example/single_cell_detailed.py
+   :language: python
+   :lines: 97-101
 
 Now all three mechanisms in the *decor* object have been made available to the model.
 
@@ -450,22 +366,18 @@ measure at this point is the spikes from the spike detectors placed in the decor
 The :class:`arbor.single_cell_model` can also measure the voltage on specific locations of the cell.
 We can indicate the location we would like to probe using labels from the :class:`label_dict`:
 
-.. code-block:: python
-
-   # Add voltage probes on the "custom_terminal" locset
-   # which sample the voltage at 50 kHz
-
-   model.probe('voltage', where='"custom_terminal"', frequency=50)
+.. literalinclude:: ../../python/example/single_cell_detailed.py
+   :language: python
+   :lines: 103-107
 
 The simulation
 ^^^^^^^^^^^^^^
 
 The cell and model descriptions are now complete and we can run the simulation:
 
-.. code-block:: python
-
-   # Run the simulation for 100 ms, with a dt of 0.025 ms
-   model.run(tfinal=100, dt=0.025)
+.. literalinclude:: ../../python/example/single_cell_detailed.py
+   :language: python
+   :lines: 109-111
 
 The results
 ^^^^^^^^^^^
@@ -474,32 +386,18 @@ Finally we move on to the data collection segment of the example. We have added
 on the "axon_terminal" locset. The :class:`arbor.single_cell_model` automatically registers all
 spikes on the cell from all spike detectors on the cell and saves the times at which they occurred.
 
-.. code-block:: python
-
-   # Print the number of spikes.
-   print(len(model.spikes), 'spikes recorded:')
-
-   # Print the spike times.
-   for s in model.spikes:
-       print(s)
-
+.. literalinclude:: ../../python/example/single_cell_detailed.py
+   :language: python
+   :lines: 113-117
 
 A more interesting result of the simulation is perhaps the output of the voltage probe previously
 placed on the "custom_terminal" locset. The model saves the output of the probes as [time, value]
 pairs which can then be plotted. We use `pandas` and `seaborn` for the plotting, but the user can
 choose the any other library:
 
-.. code-block:: python
-
-   import pandas
-   import seaborn
-
-   # Plot the output of the probes
-   df = pandas.DataFrame()
-   for t in model.traces:
-      df=df.append(pandas.DataFrame({'t/ms': t.time, 'U/mV': t.value, 'Location': str(t.location), 'Variable': t.variable}))
-
-   seaborn.relplot(data=df, kind="line", x="t/ms", y="U/mV",hue="Location",col="Variable",ci=None).savefig('single_cell_detailed_result.svg')
+.. literalinclude:: ../../python/example/single_cell_detailed.py
+   :language: python
+   :lines: 5,6,118-125
 
 The following plot is generated. The orange line is slightly delayed from the blue line, which is
 what we'd expect because branch 4 is longer than branch 3 of the morphology. We also see 3 spikes,
diff --git a/doc/tutorial/single_cell_detailed_recipe.rst b/doc/tutorial/single_cell_detailed_recipe.rst
index a05849cf5b7f49791ae7043619ad7203faa139c8..73af662c1599c2df6ccb4dda972334ff8372aef3 100644
--- a/doc/tutorial/single_cell_detailed_recipe.rst
+++ b/doc/tutorial/single_cell_detailed_recipe.rst
@@ -37,86 +37,18 @@ We outline the following steps of this example:
 The cell
 ********
 
-We can immediately paste the cell description code from the
+We can copy the cell description code or reuse ``single_cell_detailed.swc`` from the
 :ref:`previous example <tutorialsinglecellswc-cell>` where it is explained in detail.
 
-.. code-block:: python
-
-   import arbor
-   from arbor import mechanism as mech
-
-   #(1) Read the morphology from an SWC file.
-
-   morph = arbor.load_swc_arbor("single_cell_detailed.swc")
-
-   #(2) Create and populate the label dictionary.
-
-   labels = arbor.label_dict()
-
-   # Regions:
-
-   labels['soma'] = '(tag 1)'
-   labels['axon'] = '(tag 2)'
-   labels['dend'] = '(tag 3)'
-   labels['last'] = '(tag 4)'
-
-   labels['all'] = '(all)'
-
-   labels['gt_1.5'] = '(radius-ge (region "all") 1.5)'
-   labels['custom'] = '(join (region "last") (region "gt_1.5"))'
-
-   # Locsets:
-
-   labels['root']     = '(root)'
-   labels['terminal'] = '(terminal)'
-   labels['custom_terminal'] = '(restrict (locset "terminal") (region "custom"))'
-   labels['axon_terminal'] = '(restrict (locset "terminal") (region "axon"))'
-
-   # (3) Create and populate the decor.
-
-   decor = arbor.decor()
-
-   # Set the default properties of the cell (this overrides the model defaults)
-
-   decor.set_property(Vm =-55)
-
-   # Override the cell defaults.
-
-   decor.paint('"custom"', tempK=270)
-   decor.paint('"soma"',   Vm=-50)
-
-   # Paint density mechanisms.
-
-   decor.paint('"all"', 'pas')
-   decor.paint('"custom"', 'hh')
-   decor.paint('"dend"',  mech('Ih', {'gbar': 0.001}))
-
-   # Place stimuli and spike detectors.
-
-   decor.place('"root"', arbor.iclamp(10, 1, current=2), 'iclamp0')
-   decor.place('"root"', arbor.iclamp(30, 1, current=2), 'iclamp1')
-   decor.place('"root"', arbor.iclamp(50, 1, current=2), 'iclamp2')
-   decor.place('"axon_terminal"', arbor.spike_detector(-10), 'detector')
-
-   # Set cv_policy
-
-   soma_policy = arbor.cv_policy_single('"soma"')
-   dflt_policy = arbor.cv_policy_max_extent(1.0)
-   policy = dflt_policy | soma_policy
-   decor.discretization(policy)
-
-   # (4) Create the cell.
-
-   cell = arbor.cable_cell(morph, labels, decor)
-
-We will add one more thing to this section. We will create the voltage probe at the "custom_terminal" locset.
+We will need to add one more thing to the cell. We will create the voltage probe at the "custom_terminal" locset.
 In the previous example, this probe was registered directly using the :class:`arbor.single_cell_model` object.
 Now it has to be explicitly created and registered in the recipe.
 
 .. _tutorialsinglecellswcrecipe-probe:
-.. code-block:: python
 
-   probe = arbor.cable_probe_membrane_voltage('"custom_terminal"')
+.. literalinclude:: ../../python/example/single_cell_detailed_recipe.py
+   :language: python
+   :lines: 87-89
 
 The recipe
 **********
@@ -125,72 +57,20 @@ The :class:`arbor.single_cell_model` of the previous example created a :class:`a
 the hood, and abstracted away the details so we were unaware of its existence. In this example, we will
 examine the recipe in detail: how to create one, and why it is needed.
 
-.. code-block:: python
-
-   # (1) Create a class that inherits from arbor.recipe
-   class single_recipe (arbor.recipe):
-
-       # (2) Define the class constructor
-       def __init__(self, cell, probes):
-           # The base C++ class constructor must be called first, to ensure that
-           # all memory in the C++ class is initialized correctly.
-           arbor.recipe.__init__(self)
-           self.the_cell = cell
-           self.the_probes = probes
-
-           self.the_cat = arbor.default_catalogue()
-           self.the_cat.extend(arbor.allen_catalogue(), "")
-
-           self.the_props = arbor.cable_global_properties()
-           self.the_props.set_property(Vm=-65, tempK=300, rL=35.4, cm=0.01)
-           self.the_props.set_ion(ion='na', int_con=10,   ext_con=140, rev_pot=50, method='nernst/na')
-           self.the_props.set_ion(ion='k',  int_con=54.4, ext_con=2.5, rev_pot=-77)
-           self.the_props.set_ion(ion='ca', int_con=5e-5, ext_con=2, rev_pot=132.5)
-
-           self.the_props.register(self.the_cat)
-
-       # (3) Override the num_cells method
-       def num_cells(self):
-           return 1
-
-       # (4) Override the cell_kind method
-       def cell_kind(self, gid):
-           return arbor.cell_kind.cable
-
-       # (5) Override the cell_description method
-       def cell_description(self, gid):
-           return self.the_cell
-
-       # (6) Override the probes method
-       def probes(self, gid):
-           return self.the_probes
-
-       # (7) Override the connections_on method
-       def connections_on(self, gid):
-           return []
-
-       # (8) Override the gap_junction_on method
-       def gap_junction_on(self, gid):
-           return []
-
-       # (9) Override the event_generators method
-       def event_generators(self, gid):
-           return []
-
-       # (10) Overrode the global_properties method
-       def global_properties(self, gid):
-          return self.the_props
+.. literalinclude:: ../../python/example/single_cell_detailed_recipe.py
+   :language: python
+   :lines: 91-143
 
 Let's go through the recipe point by point.
 
-Step **(1)** creates a ``single_recipe`` class that inherits from :class:`arbor.recipe`. The base recipe
+Step **(6)** creates a ``single_recipe`` class that inherits from :class:`arbor.recipe`. The base recipe
 implements all the methods defined above with default values except :meth:`arbor.recipe.num_cells`,
 :meth:`arbor.recipe.cell_kind` and :meth:`arbor.recipe.cell_description` which always have to be implemented
 by the user. The :meth:`arbor.recipe.global_properties` also needs to be implemented for
 :class:`arbor.cell_kind.cable` cells. The inherited recipe can implement any number of additional methods and
 have any number of instance or class variables.
 
-Step **(2)** defines the class constructor. In this case, we pass a ``cell`` and a set of ``probes`` as
+Step **(6.1)** defines the class constructor. In this case, we pass a ``cell`` and a set of ``probes`` as
 arguments. These will be used to initialize the instance variables ``self.the_cell`` and ``self.the_probes``,
 which will be used in the overloaded ``cell_description`` and ``get_probes`` methods. Before variable
 initialization, we call the base C++ class constructor ``arbor.recipe.__init__(self)``. This ensures correct
@@ -210,49 +90,57 @@ what we did in the :ref:`previous example <tutorialsinglecellswc-gprop>`. One la
    The mechanism catalogue needs to live in the recipe as an instance variable. Its lifetime needs to extend
    to the entire duration of the simulation.
 
-Step **(3)** overrides the :meth:`arbor.recipe.num_cells` method. It takes no arguments. We simply return 1,
+Step **(6.2)** overrides the :meth:`arbor.recipe.num_cells` method. It takes no arguments. We simply return 1,
 as we are only simulating one cell in this example.
 
-Step **(4)** overrides the :meth:`arbor.recipe.cell_kind` method. It takes one argument: ``gid``.
+Step **(6.3)** overrides the :meth:`arbor.recipe.cell_kind` method. It takes one argument: ``gid``.
 Given the gid, this method returns the kind of the cell. Our defined cell is a
 :class:`arbor.cell_kind.cable`, so we simply return that.
 
-Step **(5)** overrides the :meth:`arbor.recipe.cell_description` method. It takes one argument: ``gid``.
+Step **(6.4)** overrides the :meth:`arbor.recipe.cell_description` method. It takes one argument: ``gid``.
 Given the gid, this method returns the cell description which is the cell object passed to the constructor
 of the recipe. We return ``self.the_cell``.
 
-Step **(6)** overrides the :meth:`arbor.recipe.get_probes` method. It takes one argument: ``gid``.
+Step **(6.5)** overrides the :meth:`arbor.recipe.get_probes` method. It takes one argument: ``gid``.
 Given the gid, this method returns all the probes on the cell. The probes can be of many different kinds
 measuring different quantities on different locations of the cell. We pass these probes explicitly to the recipe
 and they are stored in ``self.the_probes``, so we return that variable.
 
-Step **(7)** overrides the :meth:`arbor.recipe.connections_on` method. It takes one argument: ``gid``.
+Step **(6.6)** overrides the :meth:`arbor.recipe.connections_on` method. It takes one argument: ``gid``.
 Given the gid, this method returns all the connections ending on that cell. These are typically synapse
 connections from other cell *sources* to specific *targets* on the cell with id ``gid``. Since we are
 simulating a single cell, and self-connections are not possible, we return an empty list.
 
-Step **(8)** overrides the :meth:`arbor.recipe.gap_junctions_on` method. It takes one argument: ``gid``.
+Step **(6.7)** overrides the :meth:`arbor.recipe.gap_junctions_on` method. It takes one argument: ``gid``.
 Given the gid, this method returns all the gap junctions on that cell. Gap junctions require 2 separate cells.
 Since we are simulating a single cell, we return an empty list.
 
-Step **(9)** overrides the :meth:`arbor.recipe.event_generators` method. It takes one argument: ``gid``.
+Step **(6.8)** overrides the :meth:`arbor.recipe.event_generators` method. It takes one argument: ``gid``.
 Given the gid, this method returns *event generators* on that cell. These generators trigger events (or
 spikes) on specific *targets* on the cell. They can be used to simulate spikes from other cells, to kick-start
 a simulation for example. Our cell uses a current clamp as a stimulus, and has no targets, so we return
 an empty list.
 
-Step **(10)** overrides the :meth:`arbor.recipe.global_properties` method. It takes one argument: ``kind``.
+Step **(6.9)** overrides the :meth:`arbor.recipe.global_properties` method. It takes one argument: ``kind``.
 This method returns the default global properties of the model which apply to all cells in the network of
 that kind. We return ``self.the_props`` which we defined in step **(1)**.
 
+.. Note::
+
+   You may wonder why the method :meth:`arbor.recipe.cell_kind` is required, since it can be inferred by examining the cell description.
+   The recipe was designed to allow building simulations efficiently in a distributed system with minimum
+   communication. Some parts of the model initialization require only the cell kind,
+   not the full cell description which can be quite expensive to build. Providing these
+   descriptions separately saves time and resources for the user.
+
+   More information on the recipe can be found :ref:`here <modelrecipe>`.
+
 Now we can instantiate a ``single_recipe`` object using the ``cell`` and ``probe`` we created in the
 previous section:
 
-.. code-block:: python
-
-   # Instantiate recipe
-   # Pass the probe in a list because that it what single_recipe expects.
-   recipe = single_recipe(cell, [probe])
+.. literalinclude:: ../../python/example/single_cell_detailed_recipe.py
+   :language: python
+   :lines: 145-147
 
 The execution context
 *********************
@@ -267,9 +155,9 @@ The details of the execution context can be customized by the user. We may speci
 in the thread pool; determine the id of the GPU to be used; or create our own MPI communicator. However,
 the ideal settings can usually be inferred from the system, and Arbor can do that with a simple command.
 
-.. code-block:: python
-
-   context = arbor.context()
+.. literalinclude:: ../../python/example/single_cell_detailed_recipe.py
+   :language: python
+   :lines: 149-151
 
 The domain decomposition
 ************************
@@ -286,18 +174,18 @@ context, creates the :class:`arbor.domain_decomposition` object for us.
 Our example is a simple one, with just one cell. We don't need any sophisticated partitioning algorithms, so
 we can use the load balancer, which does a good job distributing simple networks.
 
-.. code-block:: python
-
-   domains = arbor.partition_load_balance(recipe, context)
+.. literalinclude:: ../../python/example/single_cell_detailed_recipe.py
+   :language: python
+   :lines: 153-155
 
 The simulation
 **************
 
 Finally we have the 3 components needed to create a :class:`arbor.simulation` object.
 
-.. code-block:: python
-
-   sim = arbor.simulation(recipe, domains, context)
+.. literalinclude:: ../../python/example/single_cell_detailed_recipe.py
+   :language: python
+   :lines: 157-159
 
 Before we run the simulation, however, we need to register what results we expect once execution is over.
 This was handled by the :class:`arbor.single_cell_model` object in the previous example.
@@ -305,15 +193,9 @@ This was handled by the :class:`arbor.single_cell_model` object in the previous
 We would like to get a list of the spikes on the cell during the runtime of the simulation, and we would like
 to plot the voltage registered by the probe on the "custom_terminal" locset.
 
-.. code-block:: python
-
-   # Instruct the simulation to record the spikes
-   sim.record(arbor.spike_recording.all)
-
-   # Instruct the simulation to sample the probe (0, 0)
-   # at a regular schedule with period = 0.02 ms (50 kHz)
-   probe_id = arbor.cell_member(0,0)
-   handle = sim.sample(probe_id, arbor.regular_schedule(0.02))
+.. literalinclude:: ../../python/example/single_cell_detailed_recipe.py
+   :language: python
+   :lines: 161-166
 
 The lines handling probe sampling warrant a second look. First, we declared ``probe_id`` to be a
 :class:`arbor.cell_member`, with :class:`arbor.cell_member.gid` = 0 and :class:`arbor.cell_member.index` = 0.
@@ -329,10 +211,9 @@ The execution
 
 We can now run the simulation we just instantiated for a duration of 100 ms with a time step of 0.025 ms.
 
-.. code-block:: python
-
-   sim.run(tfinal=100, dt=0.025)
-
+.. literalinclude:: ../../python/example/single_cell_detailed_recipe.py
+   :language: python
+   :lines: 168-170
 
 The results
 ***********
@@ -342,26 +223,15 @@ to sample the probe.
 
 We can print the times of the spikes:
 
-.. code-block:: python
-
-   spikes = sim.spikes()
-
-   # Print the number of spikes.
-   print(len(spikes), 'spikes recorded:')
-
-   # Print the spike times.
-   for s in spikes:
-       print(s)
+.. literalinclude:: ../../python/example/single_cell_detailed_recipe.py
+   :language: python
+   :lines: 172-177
 
 The probe results, again, warrant some more explanation:
 
-.. code-block:: python
-
-   data = []
-   meta = []
-   for d, m in sim.samples(handle):
-      data.append(d)
-      meta.append(m)
+.. literalinclude:: ../../python/example/single_cell_detailed_recipe.py
+   :language: python
+   :lines: 179-183
 
 ``sim.samples()`` takes a ``handle`` of the probe we wish to examine. It returns a list
 of ``(data, meta)`` terms: ``data`` being the time and value series of the probed quantity; and
@@ -372,12 +242,9 @@ to be 2.
 
 We plot the results using pandas and seaborn as we did in the previous example, and expect the same results:
 
-.. code-block:: python
-
-   df = pandas.DataFrame()
-   for i in range(len(data)):
-       df = df.append(pandas.DataFrame({'t/ms': data[i][:, 0], 'U/mV': data[i][:, 1], 'Location': str(meta[i])}))
-   seaborn.relplot(data=df, kind="line", x="t/ms", y="U/mV",hue="Location",col="Variable",ci=None).savefig('single_cell_detailed_recipe_result.svg')
+.. literalinclude:: ../../python/example/single_cell_detailed_recipe.py
+   :language: python
+   :lines: 185-188
 
 The following plot is generated. Identical to the plot of the previous example.
 
diff --git a/doc/tutorial/single_cell_model.rst b/doc/tutorial/single_cell_model.rst
index d659d0dc2eec6f7b2db3cea346b2e2eb038e70de..4ecd04f08c1930cd2e4bf38eb47ce5f1d835a656 100644
--- a/doc/tutorial/single_cell_model.rst
+++ b/doc/tutorial/single_cell_model.rst
@@ -36,27 +36,9 @@ descriptions of ion channels, synapses, spike detectors and electrical propertie
 
 Our *single-segment HH cell* has a simple morphology and dynamics, constructed as follows:
 
-.. code-block:: python
-
-    import arbor
-
-    # (1) Create a morphology with a single (cylindrical) 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 midpoint
-    labels = arbor.label_dict({'soma':   '(tag 1)',
-                               'midpoint': '(location 0 0.5)'})
-
-    # (3) Create and set up a decor object
-    decor = arbor.decor()
-    decor.set_property(Vm=-40)
-    decor.paint('"soma"', 'hh')
-    decor.place('"midpoint"', arbor.iclamp( 10, 2, 0.8), 'iclamp')
-    decor.place('"midpoint"', arbor.spike_detector(-10), 'detector')
-
-    # (4) Create cell
-    cell = arbor.cable_cell(tree, labels, decor)
+.. literalinclude:: ../../python/example/single_cell_model.py
+   :language: python
+   :lines: 4,6-23
 
 Step **(1)** constructs a :class:`arbor.segment_tree` (see also :ref:`segment tree<morph-segment_tree>`).
 The segment tree is the representation used to construct the morphology of a cell. A segment is
@@ -106,16 +88,18 @@ Arbor provides an interface for constructing single cell models with the
 :class:`arbor.single_cell_model` helper that creates a model from a cell description, with
 an interface for recording outputs and running the simulation.
 
-.. code-block:: python
+The single cell model has 4 main functions:
 
-    # (5) Make single cell model.
-    m = arbor.single_cell_model(cell)
+1. It holds the **global properties** of the model
+2. It registers **probes** on specific locations on the cell to measure the voltage.
+3. It **runs** the simulation.
+4. It collects **spikes** from spike detectors and voltage **traces** from registered probes.
 
-    # (6) Attach voltage probe sampling at 10 kHz (every 0.1 ms).
-    m.probe('voltage', '"midpoint"', frequency=10)
+Right now, we'll only set a probe and run the simulation.
 
-    # (7) Run simulation for 30 ms of simulated activity.
-    m.run(tfinal=30)
+.. literalinclude:: ../../python/example/single_cell_model.py
+   :language: python
+   :lines: 25-32
 
 Step **(5)** instantiates the :class:`arbor.single_cell_model`
 with our single-compartment cell.
@@ -133,21 +117,9 @@ The results
 Our cell and model have been defined and we have run our simulation. Now we can look at what
 the spike detector and a voltage probes from our model have produced.
 
-.. code-block:: python
-
-    # (8) Print spike times.
-    if len(m.spikes)>0:
-        print('{} spikes:'.format(len(m.spikes)))
-        for s in m.spikes:
-            print('{:3.3f}'.format(s))
-    else:
-        print('no spikes')
-
-    # (9) Plot the recorded voltages over time.
-    import pandas, seaborn # You may have to pip install these.
-    seaborn.set_theme() # Apply some styling to the plot
-    df = pandas.DataFrame({'t/ms': m.traces[0].time, 'U/mV': m.traces[0].value})
-    seaborn.relplot(data=df, kind="line", x="t/ms", y="U/mV",ci=None).savefig('single_cell_model_result.svg')
+.. literalinclude:: ../../python/example/single_cell_model.py
+   :language: python
+   :lines: 34-46
 
 Step **(8)** accesses :meth:`arbor.single_cell_model.spikes`
 to print the spike times. A single spike should be generated at around the same time the stimulus
diff --git a/doc/tutorial/single_cell_recipe.rst b/doc/tutorial/single_cell_recipe.rst
index 0012e2f099eb21b1e28583b1e2ff6150bd363ab6..bdc4b0dc739c1e8b389fb9ad600fbdc9975dec6f 100644
--- a/doc/tutorial/single_cell_recipe.rst
+++ b/doc/tutorial/single_cell_recipe.rst
@@ -12,86 +12,35 @@ and :class:`arbor.simulation` instead of a :class:`arbor.single_cell_model`.
    **Concepts covered in this example:**
 
    1. Building a :class:`arbor.recipe`.
-   2. Using the recipe, context and domain decomposition to create a :class:`arbor.simulation`
+   2. Using the recipe, default context and domain decomposition to create an :class:`arbor.simulation`
    3. Running the simulation and visualizing the results.
 
 The cell
 --------
 
-We can immediately paste the cell description code from the
-:ref:`previous example <tutorialsinglecell-cell>` where it is explained in detail.
+Let's copy the cell description from the :ref:`previous example <tutorialsinglecell-cell>`,
+where construction of the cell is explained in detail.
 
-.. code-block:: python
-
-    import arbor
-
-    # (1) Create a morphology with a single (cylindrical) 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 midpoint
-    labels = arbor.label_dict({'soma':   '(tag 1)',
-                              'midpoint': '(location 0 0.5)'})
-
-    # (3) Create cell and set properties
-    decor = arbor.decor()
-    decor.set_property(Vm=-40)
-    decor.paint('"soma"', 'hh')
-    decor.place('"midpoint"', arbor.iclamp( 10, 2, 0.8), 'iclamp')
-    decor.place('"midpoint"', arbor.spike_detector(-10), 'detector')
-    cell = arbor.cable_cell(tree, labels, decor)
+.. literalinclude:: ../../python/example/single_cell_recipe.py
+   :language: python
+   :lines: 4,8-23
 
 The recipe
 ----------
 
-The :class:`arbor.single_cell_model` of the previous example created a :class:`arbor.recipe` under
-the hood, and abstracted away the details so we were unaware of its existence.
+In the :ref:`previous example <tutorialsinglecell-cell>`, the :class:`arbor.single_cell_model` creates
+a :class:`arbor.recipe` under the hood, and abstracts away a few details that you may want control over
+in more complex simulations. Let's go into those abstractions and create an analogous :class:`arbor.recipe`
+manually.
 
-Creating an analogous recipe starts with creating a class that inherits from :class:`arbor.recipe`
+Creating a recipe starts with creating a class that inherits from :class:`arbor.recipe`
 and overrides and implements some of :class:`arbor.recipe` methods. Not all methods
 have to be overridden, but some will always have to be, such as :meth:`arbor.recipe.num_cells`.
 It returns `0` by default and models without cells are quite boring!
 
-.. code-block:: python
-
-    # (4) Define a recipe for a single cell and set of probes upon it.
-    # This constitutes the corresponding generic recipe version of
-    # `single_cell_model.py`.
-
-    class single_recipe (arbor.recipe):
-        def __init__(self, cell, probes):
-            # (4.1) The base C++ class constructor must be called first, to ensure that
-            # all memory in the C++ class is initialized correctly.
-            arbor.recipe.__init__(self)
-            self.the_cell = cell
-            self.the_probes = probes
-            self.the_props = arbor.neuron_cable_properties()
-            self.the_cat = arbor.default_catalogue()
-            self.the_props.register(self.the_cat)
-
-        def num_cells(self):
-            # (4.2) Override the num_cells method
-            return 1
-
-        def cell_kind(self, gid):
-            # (4.3) Override the cell_kind method
-            return arbor.cell_kind.cable
-
-        def cell_description(self, gid):
-            # (4.4) Override the cell_description method
-            return self.the_cell
-
-        def probes(self, gid):
-            # (4.5) Override the probes method
-            return self.the_probes
-
-        def global_properties(self, kind):
-            # (4.6) Override the global_properties method
-            return self.the_props
-
-    # (5) Instantiate recipe with a voltage probe located on "midpoint".
-
-    recipe = single_recipe(cell, [arbor.cable_probe_membrane_voltage('"midpoint"')])
+.. literalinclude:: ../../python/example/single_cell_recipe.py
+   :language: python
+   :lines: 25-62
 
 Step **(4)** describes the recipe that will reflect our single cell model.
 
@@ -134,27 +83,18 @@ and a domain decomposition. Fortunately, the default constructors of
 this model, and is what :class:`arbor.single_cell_model` does under the hood! We'll
 leave the details of this subject for another tutorial.
 
-.. code-block:: python
-
-    # (6) Create a default execution context and a default domain decomposition.
-
-    context = arbor.context()
-    domains = arbor.partition_load_balance(recipe, context)
+.. literalinclude:: ../../python/example/single_cell_recipe.py
+   :language: python
+   :lines: 64-67
 
 Step **(6)** sets up a default context and domains.
 
 The simulation
 --------------
 
-.. code-block:: python
-
-    # (7) Create and run simulation and set up 10 kHz (every 0.1 ms) sampling on the probe.
-    # The probe is located on cell 0, and is the 0th probe on that cell, thus has probe_id (0, 0).
-
-    sim = arbor.simulation(recipe, domains, context)
-    sim.record(arbor.spike_recording.all)
-    handle = sim.sample((0, 0), arbor.regular_schedule(0.1))
-    sim.run(tfinal=30)
+.. literalinclude:: ../../python/example/single_cell_recipe.py
+   :language: python
+   :lines: 69-75
 
 Step **(7)** instantiates the simulation and sets up the probe added in step 5. In the
 :class:`arbor.single_cell_model` version of this example, the probe frequency and
@@ -169,26 +109,9 @@ Apart from creating :class:`arbor.recipe` ourselves, we have changed nothing
 about this simulation compared to :ref:`the previous tutorial <tutorialsinglecell>`.
 If we create the same analysis of the results we therefore expect the same results.
 
-.. code-block:: python
-
-    # (8) Collect results.
-
-    spikes = sim.spikes()
-    data, meta = sim.samples(handle)[0]
-
-    if len(spikes)>0:
-        print('{} spikes:'.format(len(spikes)))
-        for t in spikes['time']:
-            print('{:3.3f}'.format(t))
-    else:
-        print('no spikes')
-
-    print("Plotting results ...")
-    seaborn.set_theme() # Apply some styling to the plot
-    df = pandas.DataFrame({'t/ms': data[:, 0], 'U/mV': data[:, 1]})
-    seaborn.relplot(data=df, kind="line", x="t/ms", y="U/mV", ci=None).savefig('single_cell_recipe_result.svg')
-
-    df.to_csv('single_cell_recipe_result.csv', float_format='%g')
+.. literalinclude:: ../../python/example/single_cell_recipe.py
+   :language: python
+   :lines: 77-94
 
 Step **(8)** plots the measured potentials during the runtime of the simulation.
 Retrieving the sampled quantities is a little different, these have to be accessed
diff --git a/python/example/network_ring.py b/python/example/network_ring.py
index 09b3a75cd322690a9408f5f3b42528aa76e94a37..2617a2de136ab92a4058ed5d779319202c298369 100755
--- a/python/example/network_ring.py
+++ b/python/example/network_ring.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# This script is included in documentation. Adapt line numbers if touched.
 
 import arbor
 import pandas, seaborn
diff --git a/python/example/network_ring_mpi.py b/python/example/network_ring_mpi.py
index 98a53f69f17e08be3fc7df14ef9d1ff14b81e209..51bde9d46e16f9acb27c72eae17e545b15db259d 100644
--- a/python/example/network_ring_mpi.py
+++ b/python/example/network_ring_mpi.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# This script is included in documentation. Adapt line numbers if touched.
 
 import arbor
 import pandas, seaborn
@@ -108,7 +109,7 @@ class ring_recipe (arbor.recipe):
         return self.props
 
 # (11) Instantiate recipe
-ncells = 50
+ncells = 500
 recipe = ring_recipe(ncells)
 
 # (12) Create an MPI communicator, and use it to create a hardware context
diff --git a/python/example/network_ring_mpi_plot.py b/python/example/network_ring_mpi_plot.py
index ba6a4f1fcaf7c6c6a58beac45dffac45a2721462..c5ec346c38529f9a4a5daab9a9d2cc3c4e30618d 100644
--- a/python/example/network_ring_mpi_plot.py
+++ b/python/example/network_ring_mpi_plot.py
@@ -1,3 +1,6 @@
+#!/usr/bin/env python3
+# This script is included in documentation. Adapt line numbers if touched.
+
 import glob
 import pandas, seaborn
 
diff --git a/python/example/single_cell_detailed.py b/python/example/single_cell_detailed.py
index fb3cd6fac4ae3bf95729c64d7c3d815263869600..b19627f946c090d1215d245bad9f94a267a79a96 100755
--- a/python/example/single_cell_detailed.py
+++ b/python/example/single_cell_detailed.py
@@ -1,10 +1,13 @@
+#!/usr/bin/env python3
+# This script is included in documentation. Adapt line numbers if touched.
+
 import arbor
 import pandas
 import seaborn
 import sys
 from arbor import mechanism as mech
 
-#(1) Read the morphology from an SWC file.
+# (1) Read the morphology from an SWC file.
 
 # Read the SWC filename from input
 # Example from docs: single_cell_detailed.swc
@@ -16,27 +19,32 @@ if len(sys.argv) < 2:
 filename = sys.argv[1]
 morph = arbor.load_swc_arbor(filename)
 
-#(2) Create and populate the label dictionary.
+# (2) Create and populate the label dictionary.
 
 labels = arbor.label_dict()
 
 # Regions:
 
+# Add labels for tag 1, 2, 3, 4
 labels['soma'] = '(tag 1)'
 labels['axon'] = '(tag 2)'
 labels['dend'] = '(tag 3)'
 labels['last'] = '(tag 4)'
-
+# Add a label for a region that includes the whole morphology
 labels['all'] = '(all)'
-
+# Add a label for the parts of the morphology with radius greater than 1.5 μm.
 labels['gt_1.5'] = '(radius-ge (region "all") 1.5)'
+# Join regions "last" and "gt_1.5"
 labels['custom'] = '(join (region "last") (region "gt_1.5"))'
 
 # Locsets:
 
+# Add a labels for the root of the morphology and all the terminal points
 labels['root']     = '(root)'
 labels['terminal'] = '(terminal)'
+# Add a label for the terminal locations in the "custom" region:
 labels['custom_terminal'] = '(restrict (locset "terminal") (region "custom"))'
+# Add a label for the terminal locations in the "axon" region:
 labels['axon_terminal'] = '(restrict (locset "terminal") (region "axon"))'
 
 # (3) Create and populate the decor.
@@ -44,32 +52,32 @@ labels['axon_terminal'] = '(restrict (locset "terminal") (region "axon"))'
 decor = arbor.decor()
 
 # Set the default properties of the cell (this overrides the model defaults).
-
 decor.set_property(Vm =-55)
+decor.set_ion('na', int_con=10,   ext_con=140, rev_pot=50, method='nernst/na')
+decor.set_ion('k',  int_con=54.4, ext_con=2.5, rev_pot=-77)
 
 # Override the cell defaults.
-
 decor.paint('"custom"', tempK=270)
 decor.paint('"soma"',   Vm=-50)
 
 # Paint density mechanisms.
-
 decor.paint('"all"', 'pas')
 decor.paint('"custom"', 'hh')
 decor.paint('"dend"',  mech('Ih', {'gbar': 0.001}))
 
 # Place stimuli and spike detectors.
+decor.place('"root"', arbor.iclamp(10, 1, current=2), 'iclamp0')
+decor.place('"root"', arbor.iclamp(30, 1, current=2), 'iclamp1')
+decor.place('"root"', arbor.iclamp(50, 1, current=2), 'iclamp2')
+decor.place('"axon_terminal"', arbor.spike_detector(-10), 'detector')
 
-decor.place('"root"', arbor.iclamp(10, 1, current=2), "iclamp0")
-decor.place('"root"', arbor.iclamp(30, 1, current=2), "iclamp1")
-decor.place('"root"', arbor.iclamp(50, 1, current=2), "iclamp2")
-decor.place('"axon_terminal"', arbor.spike_detector(-10), "detector")
-
-# Set cv_policy
-
+# Single CV for the "soma" region
 soma_policy = arbor.cv_policy_single('"soma"')
+# Single CV for the "soma" region
 dflt_policy = arbor.cv_policy_max_extent(1.0)
+# default policy everywhere except the soma
 policy = dflt_policy | soma_policy
+# Set cv_policy
 decor.discretization(policy)
 
 # (4) Create the cell.
@@ -86,24 +94,25 @@ model.properties.set_property(Vm =-65, tempK=300, rL=35.4, cm=0.01)
 model.properties.set_ion('na', int_con=10,   ext_con=140, rev_pot=50, method='nernst/na')
 model.properties.set_ion('k',  int_con=54.4, ext_con=2.5, rev_pot=-77)
 
-# Extend the default catalogue with the allen catalogue.
-
+# Extend the default catalogue with the Allen catalogue.
+# The function takes a second string parameter that can prefix
+# the name of the mechanisms to avoid collisions between catalogues
+# in this case we have no collisions so we use an empty prefix string.
 model.catalogue.extend(arbor.allen_catalogue(), "")
 
 # (7) Add probes.
 
+# Add voltage probes on the "custom_terminal" locset
+# which sample the voltage at 50 kHz
 model.probe('voltage', where='"custom_terminal"',  frequency=50)
 
-# (8) Run the simulation.
+# (8) Run the simulation for 100 ms, with a dt of 0.025 ms
 
 model.run(tfinal=100, dt=0.025)
 
 # (9) Print the spikes.
 
 print(len(model.spikes), 'spikes recorded:')
-
-# Print the spike times.
-
 for s in model.spikes:
     print(s)
 
diff --git a/python/example/single_cell_detailed_recipe.py b/python/example/single_cell_detailed_recipe.py
index 83ff1f056ccb6decc08926f52e822241aa6aa0d6..281da877c3b492af7312440eb68593cce80ecc39 100644
--- a/python/example/single_cell_detailed_recipe.py
+++ b/python/example/single_cell_detailed_recipe.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# This script is included in documentation. Adapt line numbers if touched.
 
 import arbor
 import pandas
@@ -6,9 +7,7 @@ import seaborn
 import sys
 from arbor import mechanism as mech
 
-#(1) Creat a cell.
-
-# Create the morphology
+# (1) Read the morphology from an SWC file.
 
 # Read the SWC filename from input
 # Example from docs: single_cell_detailed.swc
@@ -20,74 +19,79 @@ if len(sys.argv) < 2:
 filename = sys.argv[1]
 morph = arbor.load_swc_arbor(filename)
 
-# Create and populate the label dictionary.
+# (2) Create and populate the label dictionary.
 
 labels = arbor.label_dict()
 
 # Regions:
 
+# Add labels for tag 1, 2, 3, 4
 labels['soma'] = '(tag 1)'
 labels['axon'] = '(tag 2)'
 labels['dend'] = '(tag 3)'
 labels['last'] = '(tag 4)'
-
+# Add a label for a region that includes the whole morphology
 labels['all'] = '(all)'
-
+# Add a label for the parts of the morphology with radius greater than 1.5 μm.
 labels['gt_1.5'] = '(radius-ge (region "all") 1.5)'
+# Join regions "last" and "gt_1.5"
 labels['custom'] = '(join (region "last") (region "gt_1.5"))'
 
 # Locsets:
 
+# Add a labels for the root of the morphology and all the terminal points
 labels['root']     = '(root)'
 labels['terminal'] = '(terminal)'
+# Add a label for the terminal locations in the "custom" region:
 labels['custom_terminal'] = '(restrict (locset "terminal") (region "custom"))'
+# Add a label for the terminal locations in the "axon" region:
 labels['axon_terminal'] = '(restrict (locset "terminal") (region "axon"))'
 
-# Create and populate the decor.
+# (3) Create and populate the decor.
 
 decor = arbor.decor()
 
-# Set the default properties.
-
+# Set the default properties of the cell (this overrides the model defaults).
 decor.set_property(Vm =-55)
+decor.set_ion('na', int_con=10,   ext_con=140, rev_pot=50, method='nernst/na')
+decor.set_ion('k',  int_con=54.4, ext_con=2.5, rev_pot=-77)
 
-# Override the defaults.
-
+# Override the cell defaults.
 decor.paint('"custom"', tempK=270)
 decor.paint('"soma"',   Vm=-50)
 
 # Paint density mechanisms.
-
 decor.paint('"all"', 'pas')
 decor.paint('"custom"', 'hh')
 decor.paint('"dend"',  mech('Ih', {'gbar': 0.001}))
 
 # Place stimuli and spike detectors.
+decor.place('"root"', arbor.iclamp(10, 1, current=2), 'iclamp0')
+decor.place('"root"', arbor.iclamp(30, 1, current=2), 'iclamp1')
+decor.place('"root"', arbor.iclamp(50, 1, current=2), 'iclamp2')
+decor.place('"axon_terminal"', arbor.spike_detector(-10), 'detector')
 
-decor.place('"root"', arbor.iclamp(10, 1, current=2), "iclamp0")
-decor.place('"root"', arbor.iclamp(30, 1, current=2), "iclamp1")
-decor.place('"root"', arbor.iclamp(50, 1, current=2), "iclamp2")
-decor.place('"axon_terminal"', arbor.spike_detector(-10), "detector")
-
-# Set cv_policy
-
+# Single CV for the "soma" region
 soma_policy = arbor.cv_policy_single('"soma"')
+# Single CV for the "soma" region
 dflt_policy = arbor.cv_policy_max_extent(1.0)
+# default policy everywhere except the soma
 policy = dflt_policy | soma_policy
+# Set cv_policy
 decor.discretization(policy)
 
-# Create a cell
+# (4) Create the cell.
 
 cell = arbor.cable_cell(morph, labels, decor)
 
-# (2) Declare a probe.
+# (5) Declare a probe.
 
 probe = arbor.cable_probe_membrane_voltage('"custom_terminal"')
 
-# (3) Create a recipe class and instantiate a recipe
-
+# (6) Create a class that inherits from arbor.recipe
 class single_recipe (arbor.recipe):
 
+    # (6.1) Define the class constructor
     def __init__(self, cell, probes):
         # The base C++ class constructor must be called first, to ensure that
         # all memory in the C++ class is initialized correctly.
@@ -106,30 +110,40 @@ class single_recipe (arbor.recipe):
 
         self.the_props.register(self.the_cat)
 
+    # (6.2) Override the num_cells method
     def num_cells(self):
         return 1
 
+    # (6.3) Override the num_targets method
     def cell_kind(self, gid):
         return arbor.cell_kind.cable
 
+    # (6.4) Override the cell_description method
     def cell_description(self, gid):
         return self.the_cell
 
+    # (6.5) Override the probes method
     def probes(self, gid):
         return self.the_probes
 
+    # (6.6) Override the connections_on method
     def connections_on(self, gid):
         return []
 
+    # (6.7) Override the gap_junction_on method
     def gap_junction_on(self, gid):
         return []
 
+    # (6.8) Override the event_generators method
     def event_generators(self, gid):
         return []
 
+    # (6.9) Overrode the global_properties method
     def global_properties(self, gid):
         return self.the_props
 
+# Instantiate recipe
+# Pass the probe in a list because that it what single_recipe expects.
 recipe = single_recipe(cell, [probe])
 
 # (4) Create an execution context
diff --git a/python/example/single_cell_model.py b/python/example/single_cell_model.py
index d0e83d1851ae94019dc3a46d4d92b9f5a3a09314..51431cc63754ecd366c228da2464ff592d73dca5 100755
--- a/python/example/single_cell_model.py
+++ b/python/example/single_cell_model.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# This script is included in documentation. Adapt line numbers if touched.
 
 import arbor
 import pandas, seaborn # You may have to pip install these.
diff --git a/python/example/single_cell_recipe.py b/python/example/single_cell_recipe.py
index f3d51e1dc304b90e889399102a85c7430830aff8..01af5d469dbeaf7f094428804a6f4250492896f7 100644
--- a/python/example/single_cell_recipe.py
+++ b/python/example/single_cell_recipe.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# This script is included in documentation. Adapt line numbers if touched.
 
 import arbor
 import pandas, seaborn # You may have to pip install these.
@@ -37,18 +38,23 @@ class single_recipe (arbor.recipe):
         self.the_props.register(self.the_cat)
 
     def num_cells(self):
+        # (4.2) Override the num_cells method
         return 1
 
     def cell_kind(self, gid):
+        # (4.3) Override the cell_kind method
         return arbor.cell_kind.cable
 
     def cell_description(self, gid):
+        # (4.4) Override the cell_description method
         return self.the_cell
 
     def probes(self, gid):
+        # (4.5) Override the probes method
         return self.the_probes
 
     def global_properties(self, kind):
+        # (4.6) Override the global_properties method
         return self.the_props
 
 # (5) Instantiate recipe with a voltage probe located on "midpoint".