diff --git a/examples/density_plot.py b/examples/density_plot.py
index fd03e77d..c831ae00 100644
--- a/examples/density_plot.py
+++ b/examples/density_plot.py
@@ -42,15 +42,15 @@
 PACKAGE_DIR = Path(__file__).resolve().parent.parent
 
 
-def extract_density(population, plane='xy', bins=100, neurite_type=NeuriteType.basal_dendrite):
+def extract_density(population, plane="xy", bins=100, neurite_type=NeuriteType.basal_dendrite):
     """Extracts the 2d histogram of the center
     coordinates of segments in the selected plane.
     """
     segment_midpoints = np.array(
-        get_feat('segment_midpoints', population, neurite_type=neurite_type)
+        get_feat("segment_midpoints", population, neurite_type=neurite_type)
     )
-    horiz = segment_midpoints[:, 'xyz'.index(plane[0])]
-    vert = segment_midpoints[:, 'xyz'.index(plane[1])]
+    horiz = segment_midpoints[:, "xyz".index(plane[0])]
+    vert = segment_midpoints[:, "xyz".index(plane[1])]
     return np.histogram2d(np.array(horiz), np.array(vert), bins=(bins, bins))
 
 
@@ -60,10 +60,10 @@ def plot_density(
     new_fig=True,
     subplot=111,
     levels=None,
-    plane='xy',
-    colorlabel='Nodes per unit area',
+    plane="xy",
+    colorlabel="Nodes per unit area",
     labelfontsize=16,
-    color_map='Reds',
+    color_map="Reds",
     no_colorbar=False,
     threshold=0.01,
     neurite_type=NeuriteType.basal_dendrite,
@@ -82,7 +82,7 @@ def plot_density(
     H2 = np.ma.masked_array(H1, mask)
 
     colormap = mpl.cm.get_cmap(color_map).copy()
-    colormap.set_bad(color='white', alpha=None)
+    colormap.set_bad(color="white", alpha=None)
 
     plots = ax.contourf(
         (xedges1[:-1] + xedges1[1:]) / 2,
@@ -96,9 +96,9 @@ def plot_density(
         cbar = plt.colorbar(plots)
         cbar.ax.set_ylabel(colorlabel, fontsize=labelfontsize)
 
-    kwargs['title'] = kwargs.get('title', '')
-    kwargs['xlabel'] = kwargs.get('xlabel', plane[0])
-    kwargs['ylabel'] = kwargs.get('ylabel', plane[1])
+    kwargs["title"] = kwargs.get("title", "")
+    kwargs["xlabel"] = kwargs.get("xlabel", plane[0])
+    kwargs["ylabel"] = kwargs.get("ylabel", plane[1])
 
     return matplotlib_utils.plot_style(fig=fig, ax=ax, **kwargs)
 
@@ -109,10 +109,10 @@ def plot_neuron_on_density(
     new_fig=True,
     subplot=111,
     levels=None,
-    plane='xy',
-    colorlabel='Nodes per unit area',
+    plane="xy",
+    colorlabel="Nodes per unit area",
     labelfontsize=16,
-    color_map='Reds',
+    color_map="Reds",
     no_colorbar=False,
     threshold=0.01,
     neurite_type=NeuriteType.basal_dendrite,
diff --git a/examples/end_to_end_distance.py b/examples/end_to_end_distance.py
index bf2460ae..2b22a6ec 100755
--- a/examples/end_to_end_distance.py
+++ b/examples/end_to_end_distance.py
@@ -55,8 +55,8 @@ def make_end_to_end_distance_plot(nb_segments, end_to_end_distance, neurite_type
     plt.figure()
     plt.plot(nb_segments, end_to_end_distance)
     plt.title(neurite_type)
-    plt.xlabel('Number of segments')
-    plt.ylabel('End-to-end distance')
+    plt.xlabel("Number of segments")
+    plt.ylabel("End-to-end distance")
     # uncomment to show
     # plt.show()
 
@@ -79,24 +79,24 @@ def _dist(seg):
 
 def main():
     #  load a neuron from an SWC file
-    filename = Path(PACKAGE_DIR, 'tests/data/swc/Neuron_3_random_walker_branches.swc')
+    filename = Path(PACKAGE_DIR, "tests/data/swc/Neuron_3_random_walker_branches.swc")
     m = nm.load_morphology(filename)
 
     # print mean end-to-end distance per neurite type
     print(
-        'Mean end-to-end distance for axons: ',
+        "Mean end-to-end distance for axons: ",
         mean_end_to_end_dist(n for n in m.neurites if n.type == nm.AXON),
     )
     print(
-        'Mean end-to-end distance for basal dendrites: ',
+        "Mean end-to-end distance for basal dendrites: ",
         mean_end_to_end_dist(n for n in m.neurites if n.type == nm.BASAL_DENDRITE),
     )
     print(
-        'Mean end-to-end distance for apical dendrites: ',
+        "Mean end-to-end distance for apical dendrites: ",
         mean_end_to_end_dist(n for n in m.neurites if n.type == nm.APICAL_DENDRITE),
     )
 
-    print('End-to-end distance per neurite (nb segments, end-to-end distance, neurite type):')
+    print("End-to-end distance per neurite (nb segments, end-to-end distance, neurite type):")
     for nrte in m.neurites:
         # plot end-to-end distance for increasingly larger parts of neurite
         calculate_and_plot_end_to_end_distance(nrte)
@@ -108,5 +108,5 @@ def main():
         )
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     main()
diff --git a/examples/extract_distribution.py b/examples/extract_distribution.py
index 3821efa7..9595b8bc 100755
--- a/examples/extract_distribution.py
+++ b/examples/extract_distribution.py
@@ -66,8 +66,8 @@ def main():
         find_optimal_distribution(population_directory, "section_lengths")
     )
 
-    print(json.dumps(result, indent=2, separators=(',', ': '), cls=NeuromJSON))
+    print(json.dumps(result, indent=2, separators=(",", ": "), cls=NeuromJSON))
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     main()
diff --git a/examples/features_graph_table.py b/examples/features_graph_table.py
index e488ce42..9de66ec2 100755
--- a/examples/features_graph_table.py
+++ b/examples/features_graph_table.py
@@ -41,7 +41,7 @@
 def stylize(ax, name, feature):
     """Stylization modifications to the plots"""
     ax.set_ylabel(feature)
-    ax.set_title(name, fontsize='small')
+    ax.set_title(name, fontsize="small")
 
 
 def histogram(neuron, feature, ax, bins=15, normed=True, cumulative=False):
@@ -104,5 +104,5 @@ def main():
     )
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     main()
diff --git a/examples/get_features.py b/examples/get_features.py
index e9529fce..a3a9739a 100755
--- a/examples/get_features.py
+++ b/examples/get_features.py
@@ -51,12 +51,12 @@ def stats(data):
             min and max of data
     """
     return {
-        'len': len(data),
-        'mean': np.mean(data),
-        'sum': np.sum(data),
-        'std': np.std(data),
-        'min': np.min(data),
-        'max': np.max(data),
+        "len": len(data),
+        "mean": np.mean(data),
+        "sum": np.sum(data),
+        "std": np.std(data),
+        "min": np.min(data),
+        "max": np.max(data),
     }
 
 
@@ -66,15 +66,15 @@ def pprint_stats(data):
 
 
 def main():
-    filename = Path(PACKAGE_DIR, 'tests/data/swc/Neuron.swc')
+    filename = Path(PACKAGE_DIR, "tests/data/swc/Neuron.swc")
 
     #  load a neuron from an SWC file
     m = nm.load_morphology(filename)
 
     # Get some soma information
     # Soma radius and surface area
-    print("Soma radius", nm.get('soma_radius', m))
-    print("Soma surface area", nm.get('soma_surface_area', m))
+    print("Soma radius", nm.get("soma_radius", m))
+    print("Soma surface area", nm.get("soma_surface_area", m))
 
     # Get information about neurites
     # Most neurite data can be queried for a particular type of neurite.
@@ -85,48 +85,48 @@ def main():
     # to warm up...
 
     # number of neurites
-    print('Number of neurites (all):', nm.get('number_of_neurites', m))
+    print("Number of neurites (all):", nm.get("number_of_neurites", m))
     print(
-        'Number of neurites (axons):',
-        nm.get('number_of_neurites', m, neurite_type=nm.NeuriteType.axon),
+        "Number of neurites (axons):",
+        nm.get("number_of_neurites", m, neurite_type=nm.NeuriteType.axon),
     )
     print(
-        'Number of neurites (apical dendrites):',
-        nm.get('number_of_neurites', m, neurite_type=nm.NeuriteType.apical_dendrite),
+        "Number of neurites (apical dendrites):",
+        nm.get("number_of_neurites", m, neurite_type=nm.NeuriteType.apical_dendrite),
     )
     print(
-        'Number of neurites (basal dendrites):',
-        nm.get('number_of_neurites', m, neurite_type=nm.NeuriteType.basal_dendrite),
+        "Number of neurites (basal dendrites):",
+        nm.get("number_of_neurites", m, neurite_type=nm.NeuriteType.basal_dendrite),
     )
 
     # number of sections
-    print('Number of sections:', nm.get('number_of_sections', m))
+    print("Number of sections:", nm.get("number_of_sections", m))
     print(
-        'Number of sections (axons):',
-        nm.get('number_of_sections', m, neurite_type=nm.NeuriteType.axon),
+        "Number of sections (axons):",
+        nm.get("number_of_sections", m, neurite_type=nm.NeuriteType.axon),
     )
     print(
-        'Number of sections (apical dendrites):',
-        nm.get('number_of_sections', m, neurite_type=nm.NeuriteType.apical_dendrite),
+        "Number of sections (apical dendrites):",
+        nm.get("number_of_sections", m, neurite_type=nm.NeuriteType.apical_dendrite),
     )
     print(
-        'Number of sections (basal dendrites):',
-        nm.get('number_of_sections', m, neurite_type=nm.NeuriteType.basal_dendrite),
+        "Number of sections (basal dendrites):",
+        nm.get("number_of_sections", m, neurite_type=nm.NeuriteType.basal_dendrite),
     )
 
     # number of sections per neurite
-    print('Number of sections per neurite:', nm.get('number_of_sections_per_neurite', m))
+    print("Number of sections per neurite:", nm.get("number_of_sections_per_neurite", m))
     print(
-        'Number of sections per neurite (axons):',
-        nm.get('number_of_sections_per_neurite', m, neurite_type=nm.NeuriteType.axon),
+        "Number of sections per neurite (axons):",
+        nm.get("number_of_sections_per_neurite", m, neurite_type=nm.NeuriteType.axon),
     )
     print(
-        'Number of sections per neurite (apical dendrites):',
-        nm.get('number_of_sections_per_neurite', m, neurite_type=nm.NeuriteType.apical_dendrite),
+        "Number of sections per neurite (apical dendrites):",
+        nm.get("number_of_sections_per_neurite", m, neurite_type=nm.NeuriteType.apical_dendrite),
     )
     print(
-        'Number of sections per neurite (basal dendrites):',
-        nm.get('number_of_sections_per_neurite', m, neurite_type=nm.NeuriteType.apical_dendrite),
+        "Number of sections per neurite (basal dendrites):",
+        nm.get("number_of_sections_per_neurite", m, neurite_type=nm.NeuriteType.apical_dendrite),
     )
 
     # OK, this is getting repetitive, so lets loop over valid neurite types.
@@ -135,42 +135,42 @@ def main():
 
     # Section lengths for all and different types of neurite
     for ttype in nm.NEURITE_TYPES:
-        sec_len = nm.get('section_lengths', m, neurite_type=ttype)
-        print('Section lengths (', ttype, '):', sep='')
+        sec_len = nm.get("section_lengths", m, neurite_type=ttype)
+        print("Section lengths (", ttype, "):", sep="")
         pprint_stats(sec_len)
 
     # Segment lengths for all and different types of neurite
     for ttype in nm.NEURITE_TYPES:
-        seg_len = nm.get('segment_lengths', m, neurite_type=ttype)
-        print('Segment lengths (', ttype, '):', sep='')
+        seg_len = nm.get("segment_lengths", m, neurite_type=ttype)
+        print("Segment lengths (", ttype, "):", sep="")
         pprint_stats(seg_len)
 
     # Section radial distances for all and different types of neurite
     # Careful! Here we need to pass tree type as a named argument
     for ttype in nm.NEURITE_TYPES:
-        sec_rad_dist = nm.get('section_radial_distances', m, neurite_type=ttype)
-        print('Section radial distance (', ttype, '):', sep='')
+        sec_rad_dist = nm.get("section_radial_distances", m, neurite_type=ttype)
+        print("Section radial distance (", ttype, "):", sep="")
         pprint_stats(sec_rad_dist)
 
     # Section path distances for all and different types of neurite
     # Careful! Here we need to pass tree type as a named argument
     for ttype in nm.NEURITE_TYPES:
-        sec_path_dist = nm.get('section_path_distances', m, neurite_type=ttype)
-        print('Section path distance (', ttype, '):', sep='')
+        sec_path_dist = nm.get("section_path_distances", m, neurite_type=ttype)
+        print("Section path distance (", ttype, "):", sep="")
         pprint_stats(sec_path_dist)
 
     # Local bifurcation angles for all and different types of neurite
     for ttype in nm.NEURITE_TYPES:
-        local_bifangles = nm.get('local_bifurcation_angles', m, neurite_type=ttype)
-        print('Local bifurcation angles (', ttype, '):', sep='')
+        local_bifangles = nm.get("local_bifurcation_angles", m, neurite_type=ttype)
+        print("Local bifurcation angles (", ttype, "):", sep="")
         pprint_stats(local_bifangles)
 
     # Remote bifurcation angles for all and different types of neurite
     for ttype in nm.NEURITE_TYPES:
-        rem_bifangles = nm.get('remote_bifurcation_angles', m, neurite_type=ttype)
-        print('Local bifurcation angles (', ttype, '):', sep='')
+        rem_bifangles = nm.get("remote_bifurcation_angles", m, neurite_type=ttype)
+        print("Local bifurcation angles (", ttype, "):", sep="")
         pprint_stats(rem_bifangles)
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     main()
diff --git a/examples/histogram.py b/examples/histogram.py
index 1ad439b1..4c9223ed 100644
--- a/examples/histogram.py
+++ b/examples/histogram.py
@@ -76,16 +76,16 @@ def histogram(neurons, feature, new_fig=True, subplot=111, normed=False, **kwarg
             figure file.
     """
 
-    bins = kwargs.get('bins', 25)
-    cumulative = kwargs.get('cumulative', False)
+    bins = kwargs.get("bins", 25)
+    cumulative = kwargs.get("cumulative", False)
 
     fig, ax = matplotlib_utils.get_figure(new_fig=new_fig, subplot=subplot)
 
-    kwargs['xlabel'] = kwargs.get('xlabel', feature)
+    kwargs["xlabel"] = kwargs.get("xlabel", feature)
 
-    kwargs['ylabel'] = kwargs.get('ylabel', feature + ' fraction')
+    kwargs["ylabel"] = kwargs.get("ylabel", feature + " fraction")
 
-    kwargs['title'] = kwargs.get('title', feature + ' histogram')
+    kwargs["title"] = kwargs.get("title", feature + " histogram")
 
     feature_values = [neurom.features.get(feature, neu) for neu in neurons]
 
@@ -93,7 +93,7 @@ def histogram(neurons, feature, new_fig=True, subplot=111, normed=False, **kwarg
 
     ax.hist(feature_values, bins=bins, cumulative=cumulative, label=neu_labels, density=normed)
 
-    kwargs['no_legend'] = len(neu_labels) == 1
+    kwargs["no_legend"] = len(neu_labels) == 1
 
     return matplotlib_utils.plot_style(fig=fig, ax=ax, **kwargs)
 
@@ -149,16 +149,16 @@ def population_histogram(pops, feature, new_fig=True, normed=False, subplot=111,
             figure file.
     """
 
-    bins = kwargs.get('bins', 25)
-    cumulative = kwargs.get('cumulative', False)
+    bins = kwargs.get("bins", 25)
+    cumulative = kwargs.get("cumulative", False)
 
     fig, ax = matplotlib_utils.get_figure(new_fig=new_fig, subplot=subplot)
 
-    kwargs['xlabel'] = kwargs.get('xlabel', feature)
+    kwargs["xlabel"] = kwargs.get("xlabel", feature)
 
-    kwargs['ylabel'] = kwargs.get('ylabel', feature + ' fraction')
+    kwargs["ylabel"] = kwargs.get("ylabel", feature + " fraction")
 
-    kwargs['title'] = kwargs.get('title', feature + ' histogram')
+    kwargs["title"] = kwargs.get("title", feature + " histogram")
 
     pops_feature_values = population_feature_values(pops, feature)
 
@@ -168,7 +168,7 @@ def population_histogram(pops, feature, new_fig=True, normed=False, subplot=111,
         pops_feature_values, bins=bins, cumulative=cumulative, label=pops_labels, density=normed
     )
 
-    kwargs['no_legend'] = len(pops_labels) == 1
+    kwargs["no_legend"] = len(pops_labels) == 1
 
     return matplotlib_utils.plot_style(fig=fig, ax=ax, **kwargs)
 
diff --git a/examples/iteration_analysis.py b/examples/iteration_analysis.py
index 240e8cb0..020edbb9 100755
--- a/examples/iteration_analysis.py
+++ b/examples/iteration_analysis.py
@@ -48,7 +48,7 @@
 
 
 def main():
-    filename = Path(PACKAGE_DIR, 'tests/data/swc/Neuron.swc')
+    filename = Path(PACKAGE_DIR, "tests/data/swc/Neuron.swc")
 
     #  load a neuron from an SWC file
     m = nm.load_morphology(filename)
@@ -64,22 +64,22 @@ def sec_len(sec):
         """Return the length of a section."""
         return mm.section_length(sec.points)
 
-    print('Total neurite length (sections):', sum(sec_len(s) for s in nm.iter_sections(m)))
+    print("Total neurite length (sections):", sum(sec_len(s) for s in nm.iter_sections(m)))
 
     # Get length of all neurites in cell by iterating over segments,
     # and summing the segment lengths.
     # This should yield the same result as iterating over sections.
     print(
-        'Total neurite length (segments):', sum(mm.segment_length(s) for s in nm.iter_segments(m))
+        "Total neurite length (segments):", sum(mm.segment_length(s) for s in nm.iter_segments(m))
     )
 
     # get volume of all neurites in cell by summing over segment
     # volumes
-    print('Total neurite volume:', sum(mm.segment_volume(s) for s in nm.iter_segments(m)))
+    print("Total neurite volume:", sum(mm.segment_volume(s) for s in nm.iter_segments(m)))
 
     # get area of all neurites in cell by summing over segment
     # areas
-    print('Total neurite surface area:', sum(mm.segment_area(s) for s in nm.iter_segments(m)))
+    print("Total neurite surface area:", sum(mm.segment_area(s) for s in nm.iter_segments(m)))
 
     # get total number of neurite points in cell.
     def n_points(sec):
@@ -88,24 +88,24 @@ def n_points(sec):
         # Non-root sections have duplicate first point
         return n if sec.parent is None else n - 1
 
-    print('Total number of points:', sum(n_points(s) for s in nm.iter_sections(m)))
+    print("Total number of points:", sum(n_points(s) for s in nm.iter_sections(m)))
 
     # get mean radius of neurite points in cell.
     # p[COLS.R] yields the radius for point p.
     # Note: this includes duplicated points at beginning of
     # non-trunk sections
-    print('Mean radius of points:', np.mean([s.points[:, COLS.R] for s in nm.iter_sections(m)]))
+    print("Mean radius of points:", np.mean([s.points[:, COLS.R] for s in nm.iter_sections(m)]))
 
     # get mean radius of neurite points in cell.
     # p[COLS.R] yields the radius for point p.
     # Note: this includes duplicated points at beginning of
     # non-trunk sections
     pts = [p[COLS.R] for s in m.sections[1:] for p in s.points]
-    print('Mean radius of points:', np.mean(pts))
+    print("Mean radius of points:", np.mean(pts))
 
     # get mean radius of segments
     print(
-        'Mean radius of segments:', np.mean(list(mm.segment_radius(s) for s in nm.iter_segments(m)))
+        "Mean radius of segments:", np.mean(list(mm.segment_radius(s) for s in nm.iter_segments(m)))
     )
 
     # get stats for the segment taper rate, for different types of neurite
@@ -117,28 +117,28 @@ def n_points(sec):
         ]
 
         print(
-            'Segment taper rate (',
+            "Segment taper rate (",
             ttype,
-            '):\n  mean=',
+            "):\n  mean=",
             np.mean(seg_taper_rate),
-            ', std=',
+            ", std=",
             np.std(seg_taper_rate),
-            ', min=',
+            ", min=",
             np.min(seg_taper_rate),
-            ', max=',
+            ", max=",
             np.max(seg_taper_rate),
-            sep='',
+            sep="",
         )
 
     # Number of bifurcation points.
     print(
-        'Number of bifurcation points:',
+        "Number of bifurcation points:",
         sum(1 for _ in nm.iter_sections(m, iterator_type=Section.ibifurcation_point)),
     )
 
     # Number of bifurcation points for apical dendrites
     print(
-        'Number of bifurcation points (apical dendrites):',
+        "Number of bifurcation points (apical dendrites):",
         sum(
             1
             for _ in nm.iter_sections(
@@ -150,12 +150,12 @@ def n_points(sec):
     )
 
     # Maximum branch order
-    print('Maximum branch order:', max(section.branch_order(s) for s in nm.iter_sections(m)))
+    print("Maximum branch order:", max(section.branch_order(s) for s in nm.iter_sections(m)))
 
     # Morphology's bounding box
     # Note: does not account for soma radius
-    print('Bounding box ((min x, y, z), (max x, y, z))', geom.bounding_box(m))
+    print("Bounding box ((min x, y, z), (max x, y, z))", geom.bounding_box(m))
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     main()
diff --git a/examples/nl_fst_compat.py b/examples/nl_fst_compat.py
index 671a3b9c..36623e64 100755
--- a/examples/nl_fst_compat.py
+++ b/examples/nl_fst_compat.py
@@ -40,19 +40,19 @@
 
 
 def main():
-    m_h5 = nm.load_morphology(Path(PACKAGE_DIR, 'tests/data/h5/v1/bio_neuron-001.h5'))
-    m_asc = nm.load_morphology(Path(PACKAGE_DIR, 'tests/data/neurolucida/bio_neuron-001.asc'))
+    m_h5 = nm.load_morphology(Path(PACKAGE_DIR, "tests/data/h5/v1/bio_neuron-001.h5"))
+    m_asc = nm.load_morphology(Path(PACKAGE_DIR, "tests/data/neurolucida/bio_neuron-001.asc"))
 
-    print('h5 number of sections:', nm.get('number_of_sections', m_h5))
-    print('nl number of sections:', nm.get('number_of_sections', m_asc))
-    print('h5 number of segments:', nm.get('number_of_segments', m_h5))
-    print('nl number of segments:', nm.get('number_of_segments', m_asc))
-    print('h5 total neurite length:', np.sum(nm.get('section_lengths', m_h5)))
-    print('nl total neurite length:', np.sum(nm.get('section_lengths', m_asc)))
-    print('h5 principal direction extents:', nm.get('principal_direction_extents', m_h5))
-    print('nl principal direction extents:', nm.get('principal_direction_extents', m_asc))
+    print("h5 number of sections:", nm.get("number_of_sections", m_h5))
+    print("nl number of sections:", nm.get("number_of_sections", m_asc))
+    print("h5 number of segments:", nm.get("number_of_segments", m_h5))
+    print("nl number of segments:", nm.get("number_of_segments", m_asc))
+    print("h5 total neurite length:", np.sum(nm.get("section_lengths", m_h5)))
+    print("nl total neurite length:", np.sum(nm.get("section_lengths", m_asc)))
+    print("h5 principal direction extents:", nm.get("principal_direction_extents", m_h5))
+    print("nl principal direction extents:", nm.get("principal_direction_extents", m_asc))
 
-    print('\nNumber of neurites:')
+    print("\nNumber of neurites:")
     for nt in iter(nm.NeuriteType):
         print(
             nt,
@@ -60,7 +60,7 @@ def main():
             mf.number_of_neurites(m_asc, neurite_type=nt),
         )
 
-    print('\nNumber of segments:')
+    print("\nNumber of segments:")
     for nt in iter(nm.NeuriteType):
         print(nt, nf.number_of_segments(m_h5.neurites[0]), nf.number_of_segments(m_asc.neurites[0]))
 
diff --git a/examples/plot_somas.py b/examples/plot_somas.py
index ea4580d0..d3e9d596 100755
--- a/examples/plot_somas.py
+++ b/examples/plot_somas.py
@@ -36,7 +36,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-DATA_PATH = Path(__file__).resolve().parent.parent / 'tests/data/swc'
+DATA_PATH = Path(__file__).resolve().parent.parent / "tests/data/swc"
 
 
 def random_color():
@@ -47,7 +47,7 @@ def random_color():
 def plot_somas(somas):
     """Plot set of somas on same figure as spheres, each with different color."""
     _, ax = matplotlib_utils.get_figure(
-        new_fig=True, subplot=111, params={'projection': '3d', 'aspect': 'auto'}
+        new_fig=True, subplot=111, params={"projection": "3d", "aspect": "auto"}
     )
     for s in somas:
         matplotlib_utils.plot_sphere(ax, s.center, s.radius, color=random_color(), alpha=1)
@@ -60,7 +60,7 @@ def main():
     #  define set of files containing relevant morphs
     file_nms = [
         Path(DATA_PATH, file_nm)
-        for file_nm in ['Soma_origin.swc', 'Soma_translated_1.swc', 'Soma_translated_2.swc']
+        for file_nm in ["Soma_origin.swc", "Soma_translated_1.swc", "Soma_translated_2.swc"]
     ]
 
     # load from file and plot
@@ -68,5 +68,5 @@ def main():
     plot_somas(sms)
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     main()
diff --git a/examples/radius_of_gyration.py b/examples/radius_of_gyration.py
index 7f442904..8054203d 100755
--- a/examples/radius_of_gyration.py
+++ b/examples/radius_of_gyration.py
@@ -93,7 +93,7 @@ def mean_rad_of_gyration(neurites):
 
 def main():
     #  load a neuron from an SWC file
-    filename = Path(PACKAGE_DIR, 'tests/data/swc/Neuron.swc')
+    filename = Path(PACKAGE_DIR, "tests/data/swc/Neuron.swc")
     m = nm.load_morphology(filename)
 
     # for every neurite, print (number of segments, radius of gyration, neurite type)
@@ -110,18 +110,18 @@ def main():
 
     # print mean radius of gyration per neurite type
     print(
-        'Mean radius of gyration for axons: ',
+        "Mean radius of gyration for axons: ",
         mean_rad_of_gyration(n for n in m.neurites if n.type == nm.AXON),
     )
     print(
-        'Mean radius of gyration for basal dendrites: ',
+        "Mean radius of gyration for basal dendrites: ",
         mean_rad_of_gyration(n for n in m.neurites if n.type == nm.BASAL_DENDRITE),
     )
     print(
-        'Mean radius of gyration for apical dendrites: ',
+        "Mean radius of gyration for apical dendrites: ",
         mean_rad_of_gyration(n for n in m.neurites if n.type == nm.APICAL_DENDRITE),
     )
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     main()
diff --git a/examples/section_ids.py b/examples/section_ids.py
index 8a2178d6..e615bb7f 100755
--- a/examples/section_ids.py
+++ b/examples/section_ids.py
@@ -49,13 +49,13 @@ def get_segment(neuron, section_id, segment_id):
 
 
 def main():
-    m = nm.load_morphology(Path(PACKAGE_DIR, 'tests/data/h5/v1/Neuron.h5'))
+    m = nm.load_morphology(Path(PACKAGE_DIR, "tests/data/h5/v1/Neuron.h5"))
 
     seg = get_segment(m, 3, 2)
-    print('Segment:\n', seg)
-    print('Mid-point (x, y, z):\n', mm.linear_interpolate(seg[0], seg[1], 0.5))
-    print('Mid-point R:\n', mm.interpolate_radius(seg[0][COLS.R], seg[1][COLS.R], 0.5))
+    print("Segment:\n", seg)
+    print("Mid-point (x, y, z):\n", mm.linear_interpolate(seg[0], seg[1], 0.5))
+    print("Mid-point R:\n", mm.interpolate_radius(seg[0][COLS.R], seg[1][COLS.R], 0.5))
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     main()
diff --git a/examples/soma_radius_fit.py b/examples/soma_radius_fit.py
index 09a6142a..7169d6c6 100755
--- a/examples/soma_radius_fit.py
+++ b/examples/soma_radius_fit.py
@@ -48,10 +48,10 @@ def test_multiple_distr(filepath):
     population = nm.load_morphologies(filepath)
 
     # Create a list of basic distributions
-    distr_to_check = ('norm', 'expon', 'uniform')
+    distr_to_check = ("norm", "expon", "uniform")
 
     # Get the soma radii of a population of morphs
-    soma_size = nm.get('soma_radius', population)
+    soma_size = nm.get("soma_radius", population)
 
     # Find the best fit distribution
     return st.optimal_distribution(soma_size, distr_to_check)
@@ -66,5 +66,5 @@ def main():
     )
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     main()
diff --git a/pyproject.toml b/pyproject.toml
index a3324bf8..3f3c5b87 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -83,6 +83,7 @@ target-version = [
     'py311',
 ]
 skip-string-normalization = true
+include = 'neurom\/.*\.py$|tests\/.*\.py$|doc\/source\/conf\.py$|setup\.py$|examples\/.*\.py$'
 
 [tool.isort]
 profile = "black"
diff --git a/tox.ini b/tox.ini
index 102aed15..ca60be8e 100644
--- a/tox.ini
+++ b/tox.ini
@@ -55,8 +55,16 @@ extras = docs
 commands =
     # remove autosummary output and cleanup
     rm -rf {toxinidir}/doc/source/_neurom_build {toxinidir}/doc/build/*
-    sphinx-build -W -b doctest {toxinidir}/doc/source {toxinidir}/doc/build/doctest
-    sphinx-build -W -b html {toxinidir}/doc/source {toxinidir}/doc/build/html
+    sphinx-build -b doctest \
+        {toxinidir}/doc/source \
+        {toxinidir}/doc/build/doctest \
+        -d {toxinidir}/doc/build/doctrees \
+        -W
+    sphinx-build -b html \
+        {toxinidir}/doc/source \
+        {toxinidir}/doc/build/html \
+        -d {toxinidir}/doc/build/doctrees \
+        -W
 allowlist_externals =
     rm