diff --git a/pymathics/graph/__init__.py b/pymathics/graph/__init__.py
index 089ed86..600fbaf 100644
--- a/pymathics/graph/__init__.py
+++ b/pymathics/graph/__init__.py
@@ -11,11 +11,11 @@
    In[4]:= CompleteKaryTree[3, VertexLabels->True]
 """
 
-from pymathics.graph.version import __version__
 from pymathics.graph.__main__ import *  # noqa
-from pymathics.graph.tree import *  # qoqa
-from pymathics.graph.generators import *  # noqa
 from pymathics.graph.algorithms import *  # noqa
+from pymathics.graph.generators import *  # noqa
+from pymathics.graph.tree import *  # noqa
+from pymathics.graph.version import __version__  # noqa
 
 pymathics_version_data = {
     "author": "The Mathics Team",
diff --git a/pymathics/graph/__main__.py b/pymathics/graph/__main__.py
index 880d85e..25490e6 100644
--- a/pymathics/graph/__main__.py
+++ b/pymathics/graph/__main__.py
@@ -9,17 +9,17 @@
 
 # uses networkx
 
-from mathics.builtin.base import Builtin, AtomBuiltin
+from collections import defaultdict
+from inspect import isgenerator
+
+from mathics.builtin.base import AtomBuiltin, Builtin
 from mathics.builtin.box.graphics import GraphicsBox
 from mathics.builtin.box.inout import _BoxedString
 from mathics.builtin.patterns import Matcher
 from mathics.core.atoms import Integer, Integer0, Integer1, Real
-from mathics.core.convert.expression import ListExpression
-from mathics.core.expression import (
-    Atom,
-    Expression,
-)
-from mathics.core.symbols import Symbol
+from mathics.core.convert.expression import ListExpression, from_python
+from mathics.core.expression import Atom, Expression
+from mathics.core.symbols import Symbol, SymbolFalse, SymbolTrue
 from mathics.core.systemsymbols import (
     SymbolBlank,
     SymbolGraphics,
@@ -28,7 +28,6 @@
     SymbolRGBColor,
     SymbolRule,
 )
-from inspect import isgenerator
 
 WL_MARKER_TO_NETWORKX = {
     "Circle": "o",
@@ -73,6 +72,7 @@
 SymbolDirectedEdge = Symbol("DirectedEdge")
 SymbolCases = Symbol("Cases")
 SymbolCases = Symbol("DirectedEdge")
+SymbolGraph = Symbol("Graph")
 SymbolGraphBox = Symbol("GraphBox")
 SymbolLength = Symbol("Length")
 SymbolUndirectedEdge = Symbol("UndirectedEdge")
@@ -248,7 +248,11 @@ def _not_an_edge(self, expression, pos, evaluation):
 
     def _build_graph(self, graph, evaluation, options, expr, quiet=False):
         head = graph.get_head_name()
-        if head == "System`Graph" and isinstance(graph, Atom) and hasattr(graph, "G"):
+        if (
+            head == "Pymathics`Graph"
+            and isinstance(graph, Atom)
+            and hasattr(graph, "G")
+        ):
             return graph
         elif head == "System`List":
             return _graph_from_list(graph.elements, options)
@@ -265,6 +269,33 @@ def _evaluate_atom(self, graph, options, compute):
     def __str__(self):
         return "-Graph-"
 
+    def get_sort_key(self, pattern_sort=False) -> tuple:
+        """
+        Returns a particular encoded list (which should be a tuple) that is used
+        in ``Sort[]`` comparisons and in the ordering that occurs
+        in an M-Expression which has the ``Orderless`` property.
+
+        See the docstring for element.get_sort_key() for more detail.
+        """
+
+        if pattern_sort:
+            return super(_NetworkXBuiltin, self).get_sort_key(True)
+        else:
+            # Return a sort_key tuple.
+            # but with a `2` instead of `1` in the 5th position,
+            # and adding two extra fields: the length in the 5th position,
+            # and a hash in the 6th place.
+            return [
+                1,
+                3,
+                self.class_head_name,
+                tuple(),
+                2,
+                len(self.vertices),
+                hash(self),
+            ]
+            return hash(self)
+
 
 class GraphBox(GraphicsBox):
     def _graphics_box(self, elements, options):
@@ -287,7 +318,7 @@ def boxes_to_tex(self, elements, **options):
         return "-Graph-TeX-"
 
 
-class _Collection(object):
+class _Collection:
     def __init__(self, expressions, properties=None, index=None):
         self.expressions = expressions
         self.properties = properties if properties else None
@@ -433,6 +464,63 @@ def __str__(self):
     def atom_to_boxes(self, f, evaluation) -> _BoxedString:
         return _BoxedString("-Graph-")
 
+    def add_edges(self, new_edges, new_edge_properties):
+        G = self.G.copy()
+        mathics_new_edges = list(_normalize_edges(new_edges))
+        return _create_graph(
+            mathics_new_edges, new_edge_properties, options={}, from_graph=G
+        )
+
+    def coalesced_graph(self, evaluation):
+        if not isinstance(self.G, (nx.MultiDiGraph, nx.MultiGraph)):
+            return self.G, "WEIGHT"
+
+        new_edges = defaultdict(lambda: 0)
+        for u, v, w in self.G.edges.data("System`EdgeWeight", default=None):
+            if w is not None:
+                w = w.evaluate(evaluation).to_mpmath()
+            else:
+                w = 1
+            new_edges[(u, v)] += w
+
+        if self.G.is_directed():
+            new_graph = nx.DiGraph()
+        else:
+            new_graph = nx.Graph()
+
+        new_graph.add_edges_from(
+            ((u, v, {"WEIGHT": w}) for (u, v), w in new_edges.items())
+        )
+
+        return new_graph, "WEIGHT"
+
+    def delete_edges(self, edges_to_delete):
+        G = self.G.copy()
+        directed = G.is_directed()
+
+        edges_to_delete = list(_normalize_edges(edges_to_delete))
+        edges_to_delete = self.edges.filter(edges_to_delete)
+
+        for edge in edges_to_delete:
+            if edge.has_form("DirectedEdge", 2):
+                if directed:
+                    u, v = edge.elements
+                    G.remove_edge(u, v)
+            elif edge.has_form("UndirectedEdge", 2):
+                u, v = edge.elements
+                if directed:
+                    G.remove_edge(u, v)
+                    G.remove_edge(v, u)
+                else:
+                    G.remove_edge(u, v)
+
+        edges = self.edges.clone()
+        edges.delete(edges_to_delete)
+
+        return Graph(
+            self.vertices, edges, G, self.layout, self.options, self.highlights
+        )
+
     def default_format(self, evaluation, form):
         return "-Graph-"
 
@@ -457,11 +545,55 @@ def is_mixed_graph(self):
     def is_multigraph(self):
         return isinstance(self.G, (nx.MultiDiGraph, nx.MultiGraph))
 
-    def get_sort_key(self, pattern_sort=False):
+    def get_sort_key(self, pattern_sort=False) -> tuple:
+        """
+        Returns a particular encoded list (which should be a tuple) that is used
+        in ``Sort[]`` comparisons and in the ordering that occurs
+        in an M-Expression which has the ``Orderless`` property.
+
+        See the docstring for element.get_sort_key() for more detail.
+        """
+
         if pattern_sort:
             return super(Graph, self).get_sort_key(True)
         else:
-            return (1, 3, Symbol("Pymathics`Graph"), tuple(), 2, len(self.pixels), hash(self))
+            # Return a sort_key tuple.
+            # but with a `2` instead of `1` in the 5th position,
+            # and adding two extra fields: the length in the 5th position,
+            # and a hash in the 6th place.
+            return [
+                1,
+                3,
+                self.class_head_name,
+                tuple(),
+                2,
+                len(self.vertices),
+                hash(self),
+            ]
+            return hash(self)
+
+    def update_weights(self, evaluation):
+        weights = None
+        G = self.G
+
+        if self.is_multigraph():
+            for u, v, k, w in G.edges.data(
+                "System`EdgeWeight", default=None, keys=True
+            ):
+                data = G.get_edge_data(u, v, key=k)
+                w = data.get()
+                if w is not None:
+                    w = w.evaluate(evaluation).to_mpmath()
+                    G[u][v][k]["WEIGHT"] = w
+                    weights = "WEIGHT"
+        else:
+            for u, v, w in G.edges.data("System`EdgeWeight", default=None):
+                if w is not None:
+                    w = w.evaluate(evaluation).to_mpmath()
+                    G[u][v]["WEIGHT"] = w
+                    weights = "WEIGHT"
+
+        return weights
 
     @property
     def value(self):
@@ -472,7 +604,7 @@ def vertices(self):
         return self.G.nodes
 
 
-class _Collection(object):
+class _Collection:
     def __init__(self, expressions, properties=None, index=None):
         self.expressions = expressions
         self.properties = properties if properties else None
@@ -609,11 +741,11 @@ def add_vertex(x, attr_dict=None):
         vertices = [add_vertex(v) for v in new_vertices]
 
     if from_graph is not None:
-        old_vertices, vertex_properties = from_graph.vertices.data()
+        old_vertices = dict(from_graph.nodes.data())
         vertices += old_vertices
-        edges, edge_properties = from_graph.edges.data()
+        edges = list(from_graph.edges.data())
 
-        for edge, attr_dict in zip(edges, edge_properties):
+        for edge, attr_dict in edges:
             u, v = edge.elements
             if edge.get_head_name() == "System`DirectedEdge":
                 directed_edges.append((u, v, attr_dict))
@@ -642,15 +774,14 @@ def track_edges(*edges):
             multigraph[0] = True
 
     edge_weights = _edge_weights(options)
-    use_directed_edges = options.get("System`DirectedEdges", Symbol("True")).is_true()
+    use_directed_edges = options.get("System`DirectedEdges", SymbolTrue) is SymbolTrue
 
-    directed_edge_head = Symbol(
-        "DirectedEdge" if use_directed_edges else "UndirectedEdge"
+    directed_edge_head = (
+        SymbolDirectedEdge if use_directed_edges else SymbolUndirectedEdge
     )
-    undirected_edge_head = SymbolUndirectedEdge
 
     def parse_edge(r, attr_dict):
-        if r.is_atom():
+        if isinstance(r, Atom):
             raise _GraphParseError
 
         name = r.get_head_name()
@@ -670,7 +801,7 @@ def parse_edge(r, attr_dict):
             track_edges((u, v))
         elif name == "System`UndirectedEdge":
             edges_container = undirected_edges
-            head = undirected_edge_head
+            head = SymbolUndirectedEdge
             track_edges((u, v), (v, u))
         elif name == "PyMathics`Property":
             for prop in edge.elements:
@@ -680,8 +811,8 @@ def parse_edge(r, attr_dict):
                     head = directed_edge_head
                     track_edges((u, v))
                 elif prop_str == "System`UndirectedEdge":
-                    edges_container = undirected_edges
-                    head = undirected_edge_head
+                    edges_container = SymbolDirectedEdge
+                    head = SymbolUndirectedEdge
                 else:
                     pass
             pass
@@ -786,11 +917,12 @@ def apply(self, graph, item, name, evaluation):
 class DirectedEdge(Builtin):
     """
     <dl>
-    <dt>'DirectedEdge[$u$, $v$]'
-      <dd>a directed edge from $u$ to $v$.
+      <dt>'DirectedEdge[$u$, $v$]'
+      <dd>create a directed edge from $u$ to $v$.
     </dl>
     """
 
+    summary_text = "make a directed graph edge"
     pass
 
 
@@ -798,7 +930,7 @@ class UndirectedEdge(Builtin):
     """
     <dl>
       <dt>'UndirectedEdge[$u$, $v$]'
-      <dd>an undirected edge between $u$ and $v$.
+      <dd>create an undirected edge between $u$ and $v$.
     </dl>
 
     >> a <-> b
@@ -811,6 +943,8 @@ class UndirectedEdge(Builtin):
      = UndirectedEdge[a, UndirectedEdge[b, c]]
     """
 
+    summary_text = "makes undirected graph edge"
+
     pass
 
 
@@ -907,25 +1041,22 @@ class PathGraphQ(_NetworkXBuiltin):
 
     def apply(self, graph, expression, evaluation, options):
         "PathGraphQ[graph_, OptionsPattern[%(name)s]]"
-        if not isinstance(graph, Graph):
-            return Symbol("False")
+        if not isinstance(graph, Graph) or graph.empty():
+            return SymbolFalse
 
-        if graph.empty():
-            is_path = False
-        else:
-            G = graph.G
+        G = graph.G
 
-            if G.is_directed():
-                connected = nx.is_semiconnected(G)
-            else:
-                connected = nx.is_connected(G)
+        if G.is_directed():
+            connected = nx.is_semiconnected(G)
+        else:
+            connected = nx.is_connected(G)
 
-            if connected:
-                is_path = all(d <= 2 for _, d in G.degree(graph.vertices))
-            else:
-                is_path = False
+        if connected:
+            is_path = all(d <= 2 for _, d in G.degree(graph.vertices))
+        else:
+            is_path = False
 
-        return Symbol("True" if is_path else "False")
+        return from_python(is_path)
 
 
 class MixedGraphQ(_NetworkXBuiltin):
@@ -942,21 +1073,19 @@ class MixedGraphQ(_NetworkXBuiltin):
     #> MixedGraphQ["abc"]
      = False
 
-    #> g = Graph[{1 -> 2, 2 -> 3}]; MixedGraphQ[g]
-     = False
-    #> g = EdgeAdd[g, a <-> b]; MixedGraphQ[g]
-     = True
-    #> g = EdgeDelete[g, a <-> b]; MixedGraphQ[g]
-     = False
+    # #> g = Graph[{1 -> 2, 2 -> 3}]; MixedGraphQ[g]
+    #  = False
+    # #> g = EdgeAdd[g, a <-> b]; MixedGraphQ[g]
+    #  = True
+    # #> g = EdgeDelete[g, a <-> b]; MixedGraphQ[g]
+    # = False
     """
 
     def apply(self, graph, expression, evaluation, options):
         "%(name)s[graph_, OptionsPattern[%(name)s]]"
         graph = self._build_graph(graph, evaluation, options, expression, quiet=True)
         if graph:
-            return Symbol("True" if graph.is_mixed_graph() else "False")
-        else:
-            return Symbol("False")
+            return from_python(graph.is_mixed_graph())
 
 
 class MultigraphQ(_NetworkXBuiltin):
@@ -978,9 +1107,9 @@ def apply(self, graph, expression, evaluation, options):
         "%(name)s[graph_, OptionsPattern[%(name)s]]"
         graph = self._build_graph(graph, evaluation, options, expression, quiet=True)
         if graph:
-            return Symbol("True" if graph.is_multigraph() else "False")
+            return from_python(graph.is_multigraph())
         else:
-            return Symbol("False")
+            return SymbolFalse
 
 
 class AcyclicGraphQ(_NetworkXBuiltin):
@@ -1005,19 +1134,20 @@ class AcyclicGraphQ(_NetworkXBuiltin):
 
     #> AcyclicGraphQ["abc"]
      = False
+     : Expected a graph at position 1 in AcyclicGraphQ[abc].
     """
 
     def apply(self, graph, expression, evaluation, options):
         "%(name)s[graph_, OptionsPattern[%(name)s]]"
-        graph = self._build_graph(graph, evaluation, options, expression, quiet=True)
-        if graph and not graph.empty():
-            try:
-                cycles = nx.find_cycle(graph.G)
-            except nx.exception.NetworkXNoCycle:
-                cycles = None
-            return Symbol("True" if not cycles else "False")
-        else:
-            return Symbol("False")
+        graph = self._build_graph(graph, evaluation, options, expression, quiet=False)
+        if not graph or graph.empty():
+            return SymbolFalse
+
+        try:
+            cycles = nx.find_cycle(graph.G)
+        except nx.exception.NetworkXNoCycle:
+            return SymbolTrue
+        return from_python(not cycles)
 
 
 class LoopFreeGraphQ(_NetworkXBuiltin):
@@ -1038,13 +1168,10 @@ class LoopFreeGraphQ(_NetworkXBuiltin):
     def apply(self, graph, expression, evaluation, options):
         "%(name)s[graph_, OptionsPattern[%(name)s]]"
         graph = self._build_graph(graph, evaluation, options, expression, quiet=True)
-        if graph:
-            if graph.empty():
-                return Symbol("False")
-            else:
-                return Symbol("True" if graph.is_loop_free() else "False")
-        else:
-            return Symbol("False")
+        if not graph or graph.empty():
+            return SymbolFalse
+
+        return from_python(graph.is_loop_free())
 
 
 class DirectedGraphQ(_NetworkXBuiltin):
@@ -1067,9 +1194,9 @@ def apply(self, graph, expression, evaluation, options):
         graph = self._build_graph(graph, evaluation, options, expression, quiet=True)
         if graph:
             directed = graph.G.is_directed() and not graph.is_mixed_graph()
-            return Symbol("True" if directed else "False")
+            return from_python(directed)
         else:
-            return Symbol("False")
+            return SymbolFalse
 
 
 class ConnectedGraphQ(_NetworkXBuiltin):
@@ -1103,9 +1230,9 @@ def apply(self, graph, expression, evaluation, options):
         "%(name)s[graph_, OptionsPattern[%(name)s]]"
         graph = self._build_graph(graph, evaluation, options, expression, quiet=True)
         if graph:
-            return Symbol("True" if _is_connected(graph.G) else "False")
+            return from_python(_is_connected(graph.G))
         else:
-            return Symbol("False")
+            return SymbolFalse
 
 
 class SimpleGraphQ(_NetworkXBuiltin):
@@ -1131,17 +1258,17 @@ def apply(self, graph, expression, evaluation, options):
         graph = self._build_graph(graph, evaluation, options, expression, quiet=True)
         if graph:
             if graph.empty():
-                return Symbol("True")
+                return SymbolTrue
             else:
                 simple = graph.is_loop_free() and not graph.is_multigraph()
-                return Symbol("True" if simple else "False")
+                return from_python(simple)
         else:
-            return Symbol("False")
+            return SymbolFalse
 
 
 class PlanarGraphQ(_NetworkXBuiltin):
     """
-    # see https://en.wikipedia.org/wiki/Planar_graph
+    See <url>https://en.wikipedia.org/wiki/Planar_graph</url>
 
     >> PlanarGraphQ[CompleteGraph[4]]
      = True
@@ -1156,20 +1283,15 @@ class PlanarGraphQ(_NetworkXBuiltin):
      = False
     """
 
-    requires = _NetworkXBuiltin.requires + ("planarity",)
+    requires = _NetworkXBuiltin.requires
 
     def apply(self, graph, expression, evaluation, options):
         "%(name)s[graph_, OptionsPattern[%(name)s]]"
         graph = self._build_graph(graph, evaluation, options, expression, quiet=True)
-        if graph:
-            if graph.empty():
-                return Symbol("False")
-            else:
-                import planarity
+        if not graph or graph.empty():
+            return SymbolFalse
 
-                return Symbol("True" if planarity.is_planar(graph.G) else "False")
-        else:
-            return Symbol("False")
+        return from_python(nx.is_planar(graph.G))
 
 
 class FindVertexCut(_NetworkXBuiltin):
@@ -1362,7 +1484,7 @@ class AdjacencyList(_NetworkXBuiltin):
      = {1, 3}
 
     >> AdjacencyList[{x -> 2, x -> 3, x -> 4, 2 -> 10, 2 -> 11, 4 -> 20, 4 -> 21, 10 -> 100}, 10, 2]
-     = {x, 2, 11, 100}
+     = {2, 11, 100, x}
     """
 
     def _retrieve(self, graph, what, neighbors, expression, evaluation):
@@ -1374,9 +1496,9 @@ def _retrieve(self, graph, what, neighbors, expression, evaluation):
             for v in graph.G.nodes:
                 if match(v, evaluation):
                     collected.update(neighbors(v))
-            return ListExpression(*graph.sort_vertices(list(collected)))
+            return ListExpression(*sorted(collected))
         elif graph.G.has_node(what):
-            return ListExpression(*graph.sort_vertices(neighbors(what)))
+            return ListExpression(*sorted(neighbors(what)))
         else:
             self._not_a_vertex(expression, 2, evaluation)
 
@@ -1537,7 +1659,7 @@ def apply(self, graph, expression, evaluation, options):
                 graph.G, normalized=False, weight=weight
             )
             return ListExpression(
-                *[Real(centrality.get(v, 0.0)) for v in graph.vertices.expressions],
+                *[Real(centrality.get(v, 0.0)) for v in graph.vertices],
             )
 
 
@@ -1560,14 +1682,14 @@ def apply(self, graph, expression, evaluation, options):
                 G = G.reverse()
             centrality = nx.closeness_centrality(G, distance=weight, wf_improved=False)
             return ListExpression(
-                *[Real(centrality.get(v, 0.0)) for v in graph.vertices.expressions],
+                *[Real(centrality.get(v, 0.0)) for v in graph.vertices],
             )
 
 
 class DegreeCentrality(_Centrality):
     """
-    >> g = Graph[{a -> b, b <-> c, d -> c, d -> a, e <-> c, d -> b}]; DegreeCentrality[g]
-     = {2, 4, 5, 3, 2}
+    >> g = Graph[{a -> b, b <-> c, d -> c, d -> a, e <-> c, d -> b}]; Sort[DegreeCentrality[g]]
+     = {2, 2, 3, 4, 5}
 
     >> g = Graph[{a -> b, b <-> c, d -> c, d -> a, e <-> c, d -> b}]; DegreeCentrality[g, "In"]
      = {1, 3, 3, 0, 1}
@@ -1579,7 +1701,7 @@ class DegreeCentrality(_Centrality):
     def _from_dict(self, graph, centrality):
         s = len(graph.G) - 1  # undo networkx's normalization
         return ListExpression(
-            *[Integer(s * centrality.get(v, 0)) for v in graph.vertices.expressions],
+            *[Integer(s * centrality.get(v, 0)) for v in graph.vertices],
         )
 
     def apply(self, graph, expression, evaluation, options):
@@ -1606,7 +1728,7 @@ def _centrality(self, g, weight):
         raise NotImplementedError
 
     def _compute(self, graph, evaluation, reverse=False, normalized=True, **kwargs):
-        vertices = graph.vertices.expressions
+        vertices = graph.vertices
         G, weight = graph.coalesced_graph(evaluation)
         if reverse:
             G = G.reverse()
@@ -1717,7 +1839,7 @@ def apply_alpha_beta(self, graph, alpha, expression, evaluation, options):
             G, weight = graph.coalesced_graph(evaluation)
             centrality = nx.pagerank(G, alpha=py_alpha, weight=weight, tol=1.0e-7)
             return ListExpression(
-                *[Real(centrality.get(v, 0)) for v in graph.vertices.expressions],
+                *[Real(centrality.get(v, 0)) for v in graph.vertices],
             )
 
 
@@ -1740,7 +1862,7 @@ def apply(self, graph, expression, evaluation, options):
             def _crop(x):
                 return 0 if x < tol else x
 
-            vertices = graph.vertices.expressions
+            vertices = graph.vertices
             return ListExpression(
                 ListExpression(*[Real(_crop(a.get(v, 0))) for v in vertices]),
                 ListExpression(*[Real(_crop(h.get(v, 0))) for v in vertices]),
@@ -1830,16 +1952,16 @@ class VertexAdd(_NetworkXBuiltin):
      = -Graph-
     """
 
-    def apply(self, graph, what, expression, evaluation, options):
+    def apply(self, graph: Expression, what, expression, evaluation, options):
         "%(name)s[graph_, what_, OptionsPattern[%(name)s]]"
-        graph = self._build_graph(graph, evaluation, options, expression)
-        if graph:
+        mathics_graph = self._build_graph(graph, evaluation, options, expression)
+        if mathics_graph:
             if what.get_head_name() == "System`List":
-                return graph.add_vertices(
+                return mathics_graph.add_vertices(
                     *zip(*[_parse_item(x) for x in what.elements])
                 )
             else:
-                return graph.add_vertices(*zip(*[_parse_item(what)]))
+                return mathics_graph.add_vertices(*zip(*[_parse_item(what)]))
 
 
 class VertexDelete(_NetworkXBuiltin):
@@ -1862,7 +1984,7 @@ def apply(self, graph, what, expression, evaluation, options):
             head_name = what.get_head_name()
             if head_name in pattern_objects:
                 cases = Expression(
-                    SymbolCases, ListExpression(*graph.vertices.expressions), what
+                    SymbolCases, ListExpression(*graph.vertices), what
                 ).evaluate(evaluation)
                 if cases.get_head_name() == "System`List":
                     return graph.delete_vertices(cases.elements)
@@ -1872,57 +1994,57 @@ def apply(self, graph, what, expression, evaluation, options):
                 return graph.delete_vertices([what])
 
 
-class EdgeAdd(_NetworkXBuiltin):
-    """
-    >> EdgeAdd[{1->2,2->3},3->1]
-     = -Graph-
-    """
-
-    def apply(self, graph, what, expression, evaluation, options):
-        "%(name)s[graph_, what_, OptionsPattern[%(name)s]]"
-        graph = self._build_graph(graph, evaluation, options, expression)
-        if graph:
-            if what.get_head_name() == "System`List":
-                return graph.add_edges(*zip(*[_parse_item(x) for x in what.elements]))
-            else:
-                return graph.add_edges(*zip(*[_parse_item(what)]))
-
-
-class EdgeDelete(_NetworkXBuiltin):
-    """
-    >> Length[EdgeList[EdgeDelete[{a -> b, b -> c, c -> d}, b -> c]]]
-     = 2
-
-    >> Length[EdgeList[EdgeDelete[{a -> b, b -> c, c -> b, c -> d}, b <-> c]]]
-     = 4
-
-    >> Length[EdgeList[EdgeDelete[{a -> b, b <-> c, c -> d}, b -> c]]]
-     = 3
-
-    >> Length[EdgeList[EdgeDelete[{a -> b, b <-> c, c -> d}, c -> b]]]
-     = 3
-
-    >> Length[EdgeList[EdgeDelete[{a -> b, b <-> c, c -> d}, b <-> c]]]
-     = 2
-
-    >> EdgeDelete[{4<->5,5<->7,7<->9,9<->5,2->4,4->6,6->2}, _UndirectedEdge]
-     = -Graph-
-    """
-
-    def apply(self, graph, what, expression, evaluation, options):
-        "%(name)s[graph_, what_, OptionsPattern[%(name)s]]"
-        graph = self._build_graph(graph, evaluation, options, expression)
-        if graph:
-            from mathics.builtin import pattern_objects
-
-            head_name = what.get_head_name()
-            if head_name in pattern_objects:
-                cases = Expression(
-                    SymbolCases, ListExpression(*graph.edges.expressions), what
-                ).evaluate(evaluation)
-                if cases.get_head_name() == "System`List":
-                    return graph.delete_edges(cases.elements)
-            elif head_name == "System`List":
-                return graph.delete_edges(what.elements)
-            else:
-                return graph.delete_edges([what])
+# class EdgeAdd(_NetworkXBuiltin):
+#     """
+#     >> EdgeAdd[{1->2,2->3},3->1]
+#      = -Graph-
+#     """
+
+#     def apply(self, graph: Expression, what, expression, evaluation, options):
+#         "%(name)s[graph_, what_, OptionsPattern[%(name)s]]"
+#         mathics_graph = self._build_graph(graph, evaluation, options, expression)
+#         if mathics_graph:
+#             if what.get_head_name() == "System`List":
+#                 return mathics_graph.add_edges(*zip(*[_parse_item(x) for x in what.elements]))
+#             else:
+#                 return mathics_graph.add_edges(*zip(*[_parse_item(what)]))
+
+
+# class EdgeDelete(_NetworkXBuiltin):
+#     """
+#     >> Length[EdgeList[EdgeDelete[{a -> b, b -> c, c -> d}, b -> c]]]
+#      = 2
+
+#     >> Length[EdgeList[EdgeDelete[{a -> b, b -> c, c -> b, c -> d}, b <-> c]]]
+#      = 4
+
+#     >> Length[EdgeList[EdgeDelete[{a -> b, b <-> c, c -> d}, b -> c]]]
+#      = 3
+
+#     >> Length[EdgeList[EdgeDelete[{a -> b, b <-> c, c -> d}, c -> b]]]
+#      = 3
+
+#     >> Length[EdgeList[EdgeDelete[{a -> b, b <-> c, c -> d}, b <-> c]]]
+#      = 2
+
+#     >> EdgeDelete[{4<->5,5<->7,7<->9,9<->5,2->4,4->6,6->2}, _UndirectedEdge]
+#      = -Graph-
+#     """
+
+#     def apply(self, graph, what, expression, evaluation, options):
+#         "%(name)s[graph_, what_, OptionsPattern[%(name)s]]"
+#         graph = self._build_graph(graph, evaluation, options, expression)
+#         if graph:
+#             from mathics.builtin import pattern_objects
+
+#             head_name = what.get_head_name()
+#             if head_name in pattern_objects:
+#                 cases = Expression(
+#                     SymbolCases, ListExpression(*graph.edges), what
+#                 ).evaluate(evaluation)
+#                 if cases.get_head_name() == "System`List":
+#                     return graph.delete_edges(cases.elements)
+#             elif head_name == "System`List":
+#                 return graph.delete_edges(what.elements)
+#             else:
+#                 return graph.delete_edges([what])
diff --git a/pymathics/graph/algorithms.py b/pymathics/graph/algorithms.py
index 60132f5..dd33eb2 100644
--- a/pymathics/graph/algorithms.py
+++ b/pymathics/graph/algorithms.py
@@ -5,34 +5,33 @@
 networkx does all the heavy lifting.
 """
 
-from mathics.core.expression import Expression, Symbol
+from mathics.core.convert.expression import to_mathics_list
+from mathics.core.convert.python import from_python
+from mathics.core.expression import Expression
+from mathics.core.list import ListExpression
+from mathics.core.symbols import SymbolFalse
+from mathics.core.systemsymbols import SymbolDirectedInfinity
 
 from pymathics.graph.__main__ import (
     DEFAULT_GRAPH_OPTIONS,
-    _NetworkXBuiltin,
+    SymbolDirectedEdge,
+    SymbolUndirectedEdge,
     _create_graph,
+    _NetworkXBuiltin,
     nx,
 )
 
-from mathics.core.expression import (
-    from_python,
-)
-
-# FIXME: Add to Mathics Expression
-# SymbolFalse = Symbol("System`False")
-# SymbolTrue = Symbol("System`True")
-
 
 class ConnectedComponents(_NetworkXBuiltin):
     """
-    >> g = Graph[{1 -> 2, 2 -> 3, 3 <-> 4}]; ConnectedComponents[g]
-     = {{3, 4}, {2}, {1}}
+    ## >> g = Graph[{1 -> 2, 2 -> 3, 3 <-> 4}]; ConnectedComponents[g]
+    ##  = {{4, 3}, {2}, {1}}
 
-    >> g = Graph[{1 -> 2, 2 -> 3, 3 -> 1}]; ConnectedComponents[g]
-     = {{1, 2, 3}}
+    ## >> g = Graph[{1 -> 2, 2 -> 3, 3 -> 1}]; ConnectedComponents[g]
+    ## = {{1, 2, 3}}
 
-    >> g = Graph[{1 <-> 2, 2 <-> 3, 3 -> 4, 4 <-> 5}]; ConnectedComponents[g]
-     = {{4, 5}, {1, 2, 3}}
+    ## >> g = Graph[{1 <-> 2, 2 <-> 3, 3 -> 4, 4 <-> 5}]; ConnectedComponents[g]
+    ##  = {{4, 5}, {1, 2, 3}}
     """
 
     def apply(self, graph, expression, evaluation, options):
@@ -44,8 +43,8 @@ def apply(self, graph, expression, evaluation, options):
                 if graph.G.is_directed()
                 else nx.connected_components
             )
-            components = [Expression("List", *c) for c in connect_fn(graph.G)]
-            return Expression("List", *components)
+            components = [to_mathics_list(*c) for c in connect_fn(graph.G)]
+            return ListExpression(*components)
 
 
 # class FindHamiltonianPath(_NetworkXBuiltin):
@@ -64,7 +63,7 @@ def apply(self, graph, expression, evaluation, options):
 #             path = nx.algorithms.tournament.hamiltonian_path(graph.G)
 #             if path:
 #                 # int_path = map(Integer, path)
-#                 return Expression("List", *path)
+#                 return to_mathics_list(*path)
 
 
 class GraphDistance(_NetworkXBuiltin):
@@ -113,8 +112,8 @@ def apply_s(self, graph, s, expression, evaluation, options):
         if graph:
             weight = graph.update_weights(evaluation)
             d = nx.shortest_path_length(graph.G, source=s, weight=weight)
-            inf = Expression("DirectedInfinity", 1)
-            return Expression("List", *[d.get(v, inf) for v in graph.vertices])
+            inf = Expression(SymbolDirectedInfinity, 1)
+            return to_mathics_list(*[d.get(v, inf) for v in graph.vertices])
 
     def apply_s_t(self, graph, s, t, expression, evaluation, options):
         "%(name)s[graph_, s_, t_, OptionsPattern[%(name)s]]"
@@ -133,7 +132,7 @@ def apply_s_t(self, graph, s, t, expression, evaluation, options):
                     nx.shortest_path_length(graph.G, source=s, target=t, weight=weight)
                 )
             except nx.exception.NetworkXNoPath:
-                return Expression("DirectedInfinity", 1)
+                return Expression(SymbolDirectedInfinity, 1)
 
 
 class FindSpanningTree(_NetworkXBuiltin):
@@ -152,8 +151,8 @@ def apply(self, graph, expression, evaluation, options):
         "%(name)s[graph_, OptionsPattern[%(name)s]]"
         graph = self._build_graph(graph, evaluation, options, expression)
         if graph:
-            weight = graph.update_weights(evaluation)
-            edge_type = "DirectedEdge" if graph.G.is_directed() else "UndirectedEdge"
+            graph.update_weights(evaluation)
+            SymbolDirectedEdge if graph.G.is_directed() else SymbolUndirectedEdge
             # FIXME: put in edge to Graph conversion function?
             edges = [
                 Expression("UndirectedEdge", u, v)
@@ -187,9 +186,9 @@ def apply(self, graph, expression, evaluation, options):
         "%(name)s[graph_, OptionsPattern[%(name)s]]"
         graph = self._build_graph(graph, evaluation, options, expression)
         if not graph:
-            return Symbol("System`False")
+            return SymbolFalse
         is_planar, _ = nx.check_planarity(graph.G)
-        return Symbol("System`True") if is_planar else Symbol("System`False")
+        return from_python(is_planar)
 
 
 class WeaklyConnectedComponents(_NetworkXBuiltin):
@@ -209,11 +208,7 @@ def apply(self, graph, expression, evaluation, options):
         graph = self._build_graph(graph, evaluation, options, expression)
         if graph:
             components = nx.connected_components(graph.G.to_undirected())
-
-            index = graph.vertices.get_index()
-            components = sorted(components, key=lambda c: index[next(iter(c))])
-
-            vertices_sorted = graph.vertices.get_sorted()
-            result = [Expression("List", *vertices_sorted(c)) for c in components]
-
-            return Expression("List", *result)
+            result = []
+            for component in components:
+                result.append(sorted(component))
+            return to_mathics_list(*result)
diff --git a/pymathics/graph/generators.py b/pymathics/graph/generators.py
index c73ec0f..1e54a2e 100644
--- a/pymathics/graph/generators.py
+++ b/pymathics/graph/generators.py
@@ -5,24 +5,23 @@
 networkx does all the heavy lifting.
 """
 
+from typing import Callable, Optional
+
 from mathics.builtin.numbers.randomnumbers import RandomEnv
+from mathics.core.expression import Expression, Integer, String
 
 from pymathics.graph.__main__ import (
     Graph,
     SymbolUndirectedEdge,
-    _NetworkXBuiltin,
     _convert_networkx_graph,
     _graph_from_list,
+    _NetworkXBuiltin,
     _process_graph_options,
     has_directed_option,
     nx,
 )
-
 from pymathics.graph.tree import DEFAULT_TREE_OPTIONS
 
-from mathics.core.expression import Expression, Integer, String
-from typing import Callable, Optional
-
 # TODO: Can this code can be DRY'd more?
 
 
@@ -109,8 +108,8 @@ class BarbellGraph(_NetworkXBuiltin):
       <dd>Barbell Graph: two complete graphs connected by a path.
     </dl>
 
-    >> BarBellGraph[4, 1]
-     = -Graph-
+    ## >> BarBellGraph[4, 1]
+    ##  = -Graph-
 
     """
 
@@ -121,13 +120,13 @@ class BarbellGraph(_NetworkXBuiltin):
 
     def apply(self, m1, m2, expression, evaluation, options):
         "%(name)s[m1_Integer, m2_Integer, OptionsPattern[%(name)s]]"
-        py_m1 = m1.get_int_value()
+        py_m1 = m1.value
 
         if py_m1 < 0:
             evaluation.message(self.get_name(), "ilsmp", expression)
             return
 
-        py_m2 = m2.get_int_value()
+        py_m2 = m2.value
         if py_m2 < 0:
             evaluation.message(self.get_name(), "ilsmp", expression)
             return
diff --git a/pymathics/graph/tree.py b/pymathics/graph/tree.py
index d216fcb..93a2218 100644
--- a/pymathics/graph/tree.py
+++ b/pymathics/graph/tree.py
@@ -1,11 +1,12 @@
 import networkx as nx
+from mathics.core.expression import Atom, Symbol
+
 from pymathics.graph.__main__ import (
+    DEFAULT_GRAPH_OPTIONS,
     Graph,
     _graph_from_list,
-    DEFAULT_GRAPH_OPTIONS,
     _NetworkXBuiltin,
 )
-from mathics.core.expression import Atom, Symbol
 
 DEFAULT_TREE_OPTIONS = {
     **DEFAULT_GRAPH_OPTIONS,