diff --git a/src/pathpyG/__init__.py b/src/pathpyG/__init__.py
index 4c365455..c4c23c9a 100644
--- a/src/pathpyG/__init__.py
+++ b/src/pathpyG/__init__.py
@@ -8,6 +8,7 @@
 __version__ = get_version("pathpyG")
 
 from pathpyG.core.graph import Graph
+from pathpyG.core.higher_order_graph import HigherOrderGraph
 from pathpyG.core.index_map import IndexMap
 from pathpyG.core.multi_order_model import MultiOrderModel
 from pathpyG.core.path_data import PathData
@@ -21,6 +22,7 @@
 
 __all__ = [
     "Graph",
+    "HigherOrderGraph",
     "TemporalGraph",
     "EventGraph",
     "PathData",
diff --git a/src/pathpyG/algorithms/lift_order.py b/src/pathpyG/algorithms/lift_order.py
index 1b2269b9..1d95e4b2 100644
--- a/src/pathpyG/algorithms/lift_order.py
+++ b/src/pathpyG/algorithms/lift_order.py
@@ -106,6 +106,57 @@ def lift_order_edge_index_weighted(
     return ho_index, ho_edge_weight
 
 
+def lift_node_sequence(edge_index: torch.Tensor, node_sequence: torch.Tensor) -> torch.Tensor:
+    """Extend node sequences by one order along an edge index.
+
+    Each edge `(u, v)` of the (k-1)-th order graph becomes a node of the k-th order graph,
+    representing the path of `u` followed by the last first-order node of `v`.
+
+    Args:
+        edge_index: A **sorted** edge index tensor of shape (2, num_edges).
+        node_sequence: The node sequences of the (k-1)-th order graph, of shape (num_nodes, k-1).
+
+    Returns:
+        The node sequences of the k-th order graph, of shape (num_edges, k).
+    """
+    return torch.cat([node_sequence[edge_index[0]], node_sequence[edge_index[1]][:, -1:]], dim=1)
+
+
+def lift_order_step(
+    edge_index: torch.Tensor,
+    node_sequence: torch.Tensor,
+    edge_weight: torch.Tensor | None = None,
+    aggr: str = "src",
+) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor | None]:
+    """Lift an edge index together with its node sequences by one order.
+
+    Combines the line-graph transformation of the edge index with the corresponding
+    extension of the node sequences, so that the result again describes a graph whose
+    nodes are paths of first-order nodes. The result is **not** aggregated: duplicate
+    node sequences are left for [`aggregate_edge_index`][pathpyG.algorithms.lift_order.aggregate_edge_index]
+    (or [`HigherOrderGraph.from_aggregated`][pathpyG.HigherOrderGraph.from_aggregated]) to merge.
+
+    Args:
+        edge_index: A **sorted** edge index tensor of shape (2, num_edges).
+        node_sequence: The node sequences of the (k-1)-th order graph.
+        edge_weight: The edge weights of the (k-1)-th order graph. If None, the lifted
+            graph is returned without weights.
+        aggr: The aggregation method for the edge weights. One of "src", "dst", "max",
+            "mul" or "add". Ignored if `edge_weight` is None.
+
+    Returns:
+        A tuple of the lifted edge index, the lifted node sequences and the aggregated
+        edge weights (None if `edge_weight` was None).
+    """
+    if edge_weight is None:
+        ho_index = lift_order_edge_index(edge_index, num_nodes=node_sequence.size(0))
+    else:
+        ho_index, edge_weight = lift_order_edge_index_weighted(
+            edge_index, edge_weight=edge_weight, num_nodes=node_sequence.size(0), aggr=aggr
+        )
+    return ho_index, lift_node_sequence(edge_index, node_sequence), edge_weight
+
+
 def aggregate_edge_index(
     edge_index: torch.Tensor, node_sequence: torch.Tensor, edge_weight: torch.Tensor | None = None, aggr: str = "sum"
 ) -> Graph:
diff --git a/src/pathpyG/core/higher_order_graph.py b/src/pathpyG/core/higher_order_graph.py
new file mode 100644
index 00000000..378f708d
--- /dev/null
+++ b/src/pathpyG/core/higher_order_graph.py
@@ -0,0 +1,490 @@
+"""Higher-order De Bruijn graph representation and related operations."""
+
+from __future__ import annotations
+
+import logging
+from typing import Optional, Union
+
+import torch
+from torch_geometric.data import Data
+from torch_geometric.utils import coalesce
+
+from pathpyG.algorithms.lift_order import aggregate_edge_index, lift_order_step
+from pathpyG.core.event_graph import EventGraph
+from pathpyG.core.graph import Graph
+from pathpyG.core.index_map import IndexMap
+from pathpyG.core.path_data import PathData
+from pathpyG.core.temporal_graph import TemporalGraph
+
+logger = logging.getLogger("root")
+
+
+class HigherOrderGraph(Graph):
+    """A De Bruijn graph of order `k`, whose nodes are paths of `k` first-order nodes.
+
+    A `HigherOrderGraph` has one node per distinct path of length `k`
+    in the underlying first-order graph. Repeated observations of the same path are
+    aggregated into an `edge_weight`. Timestamps are not represented: this
+    is a model of how paths flow rather than a record of what happened.
+
+    Order 1 is the degenerate case and is simply the weighted first-order graph, with
+    plain node IDs rather than tuples.
+
+    Info:
+        In addition to the attributes of [`Graph`][pathpyG.Graph], the `data` object holds:
+
+        - `node_sequence`: [Tensor][torch.Tensor] of shape `(num_nodes, order)`, the
+            first-order node indices making up the path each higher-order node represents.
+        - `edge_weight`: [Tensor][torch.Tensor] with the aggregated weight of each transition.
+        - `inverse_idx`: [Tensor][torch.Tensor] mapping each row of the *pre-aggregation*
+            node sequence to the index of the higher-order node it was merged into.
+
+    Attributes:
+        data (Data): PyG Data object containing edges and attributes.
+        mapping (IndexMap): Mapping from higher-order node IDs (tuples, for order > 1) to indices.
+        first_order_mapping (IndexMap): Mapping of the underlying first-order node IDs to indices.
+        n_first_order (int): Number of first-order nodes the higher-order nodes are built from.
+
+    Examples:
+        >>> import pathpyG as pp
+        >>> from pathpyG.core.higher_order_graph import HigherOrderGraph
+        >>> g = pp.Graph.from_edge_list([("a", "c"), ("c", "d")])
+        >>> h = HigherOrderGraph.from_graph(g)
+        >>> print(h.order, h.nodes)
+        1 ['a', 'c', 'd']
+    """
+
+    def __init__(
+        self,
+        data: Data,
+        order: Optional[int] = None,
+        first_order_mapping: Optional[IndexMap] = None,
+        n_first_order: Optional[int] = None,
+        mapping: Optional[IndexMap] = None,
+    ) -> None:
+        """Create a HigherOrderGraph from a `Data` object carrying a `node_sequence`.
+
+        Args:
+            data: PyG `Data` object with an `edge_index` and a `node_sequence` of shape
+                `(num_nodes, order)`. For order 1, the `node_sequence` may be omitted and
+                is then taken to be the identity.
+            order: Expected order `k`. If given, it is validated against the width of the
+                node sequence; if omitted, the order is inferred from it.
+            first_order_mapping: Mapping of the underlying first-order node IDs. Defaults
+                to an empty mapping.
+            n_first_order: Number of first-order nodes. Defaults to the number of IDs in
+                `first_order_mapping`, or the largest index in the node sequence plus one.
+            mapping: Mapping of higher-order node IDs to indices. For order > 1 this must
+                use tuple IDs; for order 1 it must not.
+
+        Raises:
+            ValueError: If the order, the node sequence, and the mapping disagree, or if
+                the node sequence refers to first-order nodes that do not exist.
+        """
+        if "node_sequence" not in data and order not in (None, 1):
+            raise ValueError(f"A HigherOrderGraph of order {order} requires a `node_sequence` node attribute.")
+
+        super().__init__(data, mapping=mapping)
+
+        # `Graph` creates an identity node sequence if none is given, so `self.order`
+        # (inherited: the width of the node sequence) is now well-defined.
+        if order is not None and order != self.order:
+            raise ValueError(f"order={order} does not match node sequence of width {self.order}")
+
+        if first_order_mapping is not None:
+            self.first_order_mapping = first_order_mapping
+        elif self.order == 1:
+            # For order 1 the higher-order nodes *are* the first-order nodes.
+            self.first_order_mapping = self.mapping
+        else:
+            self.first_order_mapping = IndexMap()
+
+        if n_first_order is not None:
+            self._n_first_order = int(n_first_order)
+        elif self.first_order_mapping.has_ids:
+            self._n_first_order = self.first_order_mapping.num_ids()
+        elif self.data.node_sequence.numel() > 0:
+            self._n_first_order = int(self.data.node_sequence.max().item()) + 1
+        else:
+            self._n_first_order = 0
+
+        self._validate()
+
+    def _validate(self) -> None:
+        """Check that order, node sequence, mapping and first-order node set agree."""
+        if self.data.node_sequence.numel() > 0:
+            max_idx = int(self.data.node_sequence.max().item())
+            if max_idx >= self._n_first_order:
+                raise ValueError(
+                    f"node sequence refers to first-order node {max_idx}, "
+                    f"but there are only {self._n_first_order} first-order nodes"
+                )
+
+        if self.mapping.has_ids:
+            # Higher-order nodes are paths and are identified by tuples; first-order
+            # nodes are entities and are identified by plain IDs.
+            if self.mapping.has_tuple_ids != (self.order > 1):
+                raise ValueError(
+                    f"a mapping for a graph of order {self.order} must "
+                    f"{'use' if self.order > 1 else 'not use'} tuple IDs"
+                )
+            if self.mapping.num_ids() != self.n:
+                logger.warning(
+                    "mapping has %s IDs but graph has %s nodes", self.mapping.num_ids(), self.n
+                )
+
+    @staticmethod
+    def _validate_order(order: int) -> None:
+        """Reject orders for which no De Bruijn graph is defined."""
+        if order < 1:
+            logger.error("order must be at least 1, got %s", order)
+            raise ValueError(f"order must be at least 1, got {order}")
+
+    @staticmethod
+    def _build_mapping(node_sequence: torch.Tensor, first_order_mapping: IndexMap) -> IndexMap:
+        """Build the higher-order `IndexMap` naming each node by the path it represents."""
+        # TODO: Is it better to have a single HigherOrderMapping class?
+        order = node_sequence.size(1)
+        if node_sequence.size(0) == 0:
+            # An order beyond the longest observed path yields a graph without nodes,
+            # and `IndexMap` cannot be built from an empty list of IDs.
+            return IndexMap()
+        if order == 1:
+            # Order-1 node indices are first-order node indices, so the mapping carries over.
+            return first_order_mapping
+        if first_order_mapping.has_ids:
+            return IndexMap([tuple(first_order_mapping.to_ids(v.cpu())) for v in node_sequence])
+        return IndexMap([tuple(v.tolist()) for v in node_sequence])
+
+    @classmethod
+    def aggregate(
+        cls,
+        edge_index: torch.Tensor,
+        node_sequence: torch.Tensor,
+        first_order_mapping: Optional[IndexMap] = None,
+        edge_weight: Optional[torch.Tensor] = None,
+        n_first_order: Optional[int] = None,
+        aggr: str = "sum",
+    ) -> HigherOrderGraph:
+        """Aggregate a (possibly duplicated) higher-order edge index into a De Bruijn graph.
+
+        Args:
+            edge_index: Edge index whose nodes are indices into `node_sequence`.
+            node_sequence: Tensor of shape `(num_nodes, order)` with the first-order path
+                each (not yet aggregated) node represents.
+            first_order_mapping: Mapping of the underlying first-order node IDs.
+            edge_weight: Weight of each edge prior to aggregation. Defaults to ones.
+            n_first_order: Number of first-order nodes, including isolated ones.
+            aggr: Reduction used for the edge weights. One of "sum", "mean", "min", "max".
+
+        Returns:
+            HigherOrderGraph: The aggregated higher-order graph.
+        """
+        order = node_sequence.size(1)
+        if first_order_mapping is None:
+            first_order_mapping = IndexMap()
+        if n_first_order is None:
+            if first_order_mapping.has_ids:
+                n_first_order = first_order_mapping.num_ids()
+            else:
+                n_first_order = int(node_sequence.max().item()) + 1 if node_sequence.numel() > 0 else 0
+
+        data = aggregate_edge_index(edge_index, node_sequence, edge_weight, aggr=aggr).data
+
+        if order == 1 and n_first_order > data.num_nodes:
+            # Order-1 indices are first-order indices, so first-order nodes that are not
+            # traversed by any path are simply isolated nodes of the order-1 graph.
+            data.num_nodes = n_first_order
+            data.node_sequence = torch.arange(n_first_order, device=edge_index.device).unsqueeze(1)
+
+        return cls(
+            data,
+            order=order,
+            first_order_mapping=first_order_mapping,
+            n_first_order=n_first_order,
+            mapping=cls._build_mapping(data.node_sequence, first_order_mapping),
+        )
+
+    @classmethod
+    def from_aggregated_graph(
+        cls,
+        g: Graph,
+        first_order_mapping: Optional[IndexMap] = None,
+        n_first_order: Optional[int] = None,
+    ) -> HigherOrderGraph:
+        """Adopt an already-aggregated [`Graph`][pathpyG.Graph] as a higher-order graph.
+
+        Args:
+            g: Aggregated graph carrying a `node_sequence` of shape `(num_nodes, order)`.
+            first_order_mapping: Mapping of the underlying first-order node IDs.
+            n_first_order: Number of first-order nodes.
+
+        Returns:
+            HigherOrderGraph: The same graph, typed as a higher-order graph.
+        """
+        if isinstance(g, HigherOrderGraph):
+            return g
+        return cls(
+            g.data,
+            first_order_mapping=first_order_mapping,
+            n_first_order=n_first_order,
+            mapping=g.mapping,
+        )
+
+    @classmethod
+    def from_graph(cls, g: Graph, weight: str = "edge_weight") -> HigherOrderGraph:
+        """Create the order-1 graph corresponding to a first-order graph.
+
+        Multi-edges are coalesced into a single weighted edge.
+
+        Args:
+            g: First-order graph.
+            weight: Name of the edge attribute to use as edge weight. If absent, each
+                edge counts once.
+
+        Returns:
+            HigherOrderGraph: A higher-order graph of order 1.
+        """
+        edge_index = g.data.edge_index.as_tensor()
+        if weight in g.data:
+            edge_weight = g.data[weight]
+        else:
+            edge_weight = torch.ones(edge_index.size(1), device=edge_index.device)
+        node_sequence = torch.arange(g.n, device=edge_index.device).unsqueeze(1)
+
+        return cls.aggregate(
+            edge_index,
+            node_sequence,
+            first_order_mapping=g.mapping,
+            edge_weight=edge_weight,
+            n_first_order=g.n,
+        )
+
+    @classmethod
+    def from_temporal_graph(
+        cls,
+        g: TemporalGraph,
+        order: int = 1,
+        delta: float | int = 1,
+        weight: str = "edge_weight",
+    ) -> HigherOrderGraph:
+        """Create the De Bruijn graph of order `k` for time-respecting paths in a temporal graph.
+
+        Order 1 is simply the weighted static graph and ignores `delta`; for higher orders
+        the nodes are the time-respecting paths of `k` nodes, i.e. those whose consecutive
+        interactions are at most `delta` apart. Orders above 2 are reached by repeatedly
+        lifting the unaggregated data.
+
+        Args:
+            g: The temporal graph.
+            order: The order `k` of the graph to compute.
+            delta: The maximum time difference between two consecutive interactions of a path.
+            weight: The edge attribute of `g` to use as edge weight.
+
+        Returns:
+            HigherOrderGraph: A higher-order graph of order `order`. It has no nodes if
+            there is no time-respecting path of that length.
+
+        Note:
+            Each call rebuilds the whole chain of lifts from order 1. To obtain several
+            orders, build a [`MultiOrderModel`][pathpyG.MultiOrderModel] with
+            `cached=True` once and read its `layers` instead.
+
+        Examples:
+            >>> import pathpyG as pp
+            >>> t = pp.TemporalGraph.from_edge_list([("a", "c", 1), ("c", "d", 2)])
+            >>> print(pp.HigherOrderGraph.from_temporal_graph(t, order=2, delta=1).nodes)
+            [('a', 'c'), ('c', 'd')]
+        """
+        cls._validate_order(order)
+        from pathpyG.core.multi_order_model import MultiOrderModel
+
+        return MultiOrderModel.from_temporal_graph(
+            g, delta=delta, max_order=order, weight=weight, cached=False
+        ).layers[order]
+
+    @classmethod
+    def from_path_data(cls, path_data: PathData, order: int = 1, mode: str = "propagation") -> HigherOrderGraph:
+        """Create the De Bruijn graph of order `k` modelling paths in [`PathData`][pathpyG.PathData].
+
+        Args:
+            path_data: The observed paths.
+            order: The order `k` of the graph to compute.
+            mode: The process that we assume. Either "diffusion" or "propagation".
+
+        Returns:
+            HigherOrderGraph: A higher-order graph of order `order`. It has no nodes if
+            no observed path is that long.
+
+        Note:
+            Each call rebuilds the whole chain of lifts from order 1. To obtain several
+            orders, build a [`MultiOrderModel`][pathpyG.MultiOrderModel] with
+            `cached=True` once and read its `layers` instead.
+
+        Examples:
+            >>> import pathpyG as pp
+            >>> paths = pp.PathData(pp.IndexMap(list("acd")))
+            >>> paths.append_walk(("a", "c", "d"), weight=2)
+            >>> print(pp.HigherOrderGraph.from_path_data(paths, order=2).nodes)
+            [('a', 'c'), ('c', 'd')]
+        """
+        cls._validate_order(order)
+        from pathpyG.core.multi_order_model import MultiOrderModel
+
+        return MultiOrderModel.from_path_data(path_data, max_order=order, mode=mode, cached=False).layers[order]
+
+    @classmethod
+    def from_event_graph(cls, eg: EventGraph, order: int = 2) -> HigherOrderGraph:
+        """Aggregate an [`EventGraph`][pathpyG.core.event_graph.EventGraph] into an order-`k` graph.
+
+        Equivalent to [`from_temporal_graph`][pathpyG.HigherOrderGraph.from_temporal_graph]
+        on the underlying temporal graph with the event graph's `delta`.
+
+        Args:
+            eg: The second-order temporal event graph to aggregate.
+            order: The order `k` of the graph to compute.
+
+        Returns:
+            HigherOrderGraph: A higher-order graph of order `order`. It has no nodes if
+            there is no time-respecting path of that length.
+
+        Note:
+            Each call rebuilds the whole chain of lifts from order 1. To obtain several
+            orders, build a [`MultiOrderModel`][pathpyG.MultiOrderModel] with
+            `cached=True` once and read its `layers` instead.
+        """
+        cls._validate_order(order)
+        from pathpyG.core.multi_order_model import MultiOrderModel
+
+        return MultiOrderModel.from_event_graph(eg, max_order=order, cached=False).layers[order]
+
+    def lift(self, aggr: str = "src") -> HigherOrderGraph:
+        """Return the De Bruijn graph of order `k + 1` obtained by lifting this graph.
+
+        Nodes of the result are the edges of this graph, i.e. the paths of length `k + 1`
+        that exist in this graph's topology.
+
+        Args:
+            aggr: Aggregation used for the lifted edge weights. One of "src", "dst",
+                "max", "mul" or "add".
+
+        Returns:
+            HigherOrderGraph: A higher-order graph of order `k + 1`.
+        """
+        edge_index = self.data.edge_index.as_tensor()
+        if "edge_weight" in self.data:
+            edge_weight = self.data.edge_weight
+        else:
+            edge_weight = torch.ones(edge_index.size(1), device=edge_index.device)
+
+        ho_index, node_sequence, ho_weight = lift_order_step(
+            edge_index, self.data.node_sequence, edge_weight=edge_weight, aggr=aggr
+        )
+
+        return HigherOrderGraph.aggregate(
+            ho_index,
+            node_sequence,
+            first_order_mapping=self.first_order_mapping,
+            edge_weight=ho_weight,
+            n_first_order=self.n_first_order,
+        )
+
+    def to_first_order(self, mode: str = "last") -> Graph:
+        """Project the higher-order graph back onto the first-order nodes.
+
+        Each higher-order node is replaced by one of the first-order nodes of its path,
+        and the weights of higher-order edges mapping to the same first-order edge are
+        summed. First-order nodes not traversed by any path remain as isolated nodes.
+        
+        Warning: This is a projection, not an inverse transformation
+            This method does not reconstruct the original first-order graph from
+            which this higher-order graph was built. Instead, it maps each higher-
+            order node to either the first or last first-order node in its represented path.
+            
+            Consequently, the result preserves flow encoded by the higher-order model
+            under the selected projection, but may differ from the original graph in its
+            edge set. In particular, isolated first-order edges cannot be recovered.
+
+        Args:
+            mode: Which first-order node of the path represents it. Either "last" or "first".
+
+        Returns:
+            Graph: A weighted first-order graph.
+        """
+        if mode == "last":
+            projection = self.data.node_sequence[:, -1]
+        elif mode == "first":
+            projection = self.data.node_sequence[:, 0]
+        else:
+            raise ValueError(f"Unknown mode {mode}. Only 'last' and 'first' are accepted.")
+
+        edge_index = projection[self.data.edge_index.as_tensor()]
+        if "edge_weight" in self.data:
+            edge_weight = self.data.edge_weight
+        else:
+            edge_weight = torch.ones(edge_index.size(1), device=edge_index.device)
+        edge_index, edge_weight = coalesce(
+            edge_index, edge_attr=edge_weight, num_nodes=self.n_first_order, reduce="sum"
+        )
+
+        return Graph(
+            Data(edge_index=edge_index, edge_weight=edge_weight, num_nodes=self.n_first_order),
+            mapping=self.first_order_mapping,
+        )
+
+    def bipartite_edge_index(
+        self,
+        first_order_graph: Optional[Graph] = None,
+        mapping: str = "last",
+        device: Optional[torch.device] = None,
+    ) -> torch.Tensor:
+        """Return the edge index connecting higher-order nodes to first-order nodes.
+
+        Args:
+            first_order_graph: The first-order graph. Optional; accepted so that call
+                sites read symmetrically, and used only for its device.
+            mapping: Which first-order nodes to connect to. One of "last", "first" or "both".
+            device: Device on which to create the tensor.
+
+        Returns:
+            torch.Tensor: Edge index of shape `(2, ·)`, higher-order nodes in the first row.
+        """
+        if device is None:
+            device = first_order_graph.device if first_order_graph is not None else self.device
+
+        node_sequence = self.data.node_sequence
+        ho_idx = torch.arange(self.n, device=device)
+
+        if mapping == "last":
+            fo_idx = node_sequence[:, -1].to(device)
+        elif mapping == "first":
+            fo_idx = node_sequence[:, 0].to(device)
+        elif mapping == "both":
+            fo_idx = torch.cat([node_sequence[:, 0], node_sequence[:, -1]]).to(device)
+            ho_idx = torch.cat([ho_idx, ho_idx])
+        else:
+            raise ValueError(f"Unknown mapping {mapping}. Only 'last', 'first' and 'both' are accepted.")
+
+        return torch.stack([ho_idx, fo_idx])
+
+    @property
+    def n_first_order(self) -> int:
+        """Number of first-order nodes underlying the higher-order nodes."""
+        return self._n_first_order
+
+    def node_id(self, idx: int) -> Union[str, int, tuple]:
+        """Return the first-order path represented by the higher-order node `idx`."""
+        seq = self.data.node_sequence[idx]
+        if self.order == 1:
+            return self.first_order_mapping.to_id(int(seq[0].item()))
+        if self.first_order_mapping.has_ids:
+            return tuple(self.first_order_mapping.to_ids(seq.cpu()).tolist())
+        return tuple(seq.tolist())
+
+    def __str__(self) -> str:
+        """Return a human-readable summary of the higher-order graph."""
+        s = (
+            f"Higher-order graph of order {self.order} with {self.n} nodes and {self.m} edges\n"
+            f"(over {self.n_first_order} first-order nodes)\n"
+        )
+        return s + "\n".join(super().__str__().split("\n")[1:])
diff --git a/src/pathpyG/core/multi_order_model.py b/src/pathpyG/core/multi_order_model.py
index 542156ac..30d02a1d 100644
--- a/src/pathpyG/core/multi_order_model.py
+++ b/src/pathpyG/core/multi_order_model.py
@@ -13,15 +13,15 @@
 from pathpyG.algorithms.lift_order import (
     aggregate_edge_index,
     aggregate_node_attributes,
+    lift_node_sequence,
     lift_order_edge_index,
-    lift_order_edge_index_weighted,
+    lift_order_step,
 )
 from pathpyG.core.event_graph import EventGraph
-from pathpyG.core.graph import Graph
+from pathpyG.core.higher_order_graph import HigherOrderGraph
 from pathpyG.core.index_map import IndexMap
 from pathpyG.core.path_data import PathData
 from pathpyG.core.temporal_graph import TemporalGraph
-from pathpyG.utils.dbgnn import generate_bipartite_edge_index
 
 logger = logging.getLogger("root")
 
@@ -31,13 +31,16 @@ class MultiOrderModel:
 
     This class stores multiple higher-order De Bruijn graphs as layers in a dictionary.
     Each layer corresponds to a De Bruijn graph of order k, where k is the key in the dictionary.
-    Each graph layer is represented as a [pathpyG.Graph][] object.
+    Each graph layer is represented as a
+    [HigherOrderGraph][pathpyG.core.higher_order_graph.HigherOrderGraph] object, layer 1
+    included. Each layer therefore knows its own order and the first-order nodes it was
+    built from.
     This class provides methods to search for the optimal order of the model based on likelihood ratio tests,
     as well as methods to compute the log-likelihood of observed paths given the model.
 
     Attributes:
-        layers (dict[int, Graph]): A dictionary mapping the order k to the corresponding
-            higher-order De Bruijn graph of order k.
+        layers (dict[int, HigherOrderGraph]): A dictionary mapping the order k to the
+            corresponding higher-order De Bruijn graph of order k.
 
     Examples:
         Example where the optimal order is 1:
@@ -56,11 +59,15 @@ class MultiOrderModel:
         >>> m = MultiOrderModel.from_path_data(paths, max_order=2)
         >>> print(m.estimate_order(paths, max_order=2))
         2
+
+        Each layer knows the first-order path that each of its nodes represents:
+        >>> print(m.layers[2].order, m.layers[2].nodes)
+        2 [('a', 'c'), ('b', 'c'), ('c', 'd'), ('c', 'e')]
     """
 
     def __init__(self) -> None:
         """Initialize an empty MultiOrderModel."""
-        self.layers: dict[int, Graph] = {}
+        self.layers: dict[int, HigherOrderGraph] = {}
 
     def __str__(self) -> str:
         """Return a string representation of the higher-order graph."""
@@ -88,7 +95,8 @@ def iterate_lift_order(
         edge_weight: torch.Tensor | None = None,
         aggr: str = "src",
         save: bool = True,
-    ) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor | None, Graph | None]:
+        n_first_order: Optional[int] = None,
+    ) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor | None, HigherOrderGraph | None]:
         """Lift order by one and save the result in the layers dictionary of the object.
 
         This is a helper function that should not be called directly.
@@ -103,20 +111,23 @@ def iterate_lift_order(
             k: The order of the graph that should be computed.
             aggr: The aggregation method to use. One of "src", "dst", "max", "mul".
             save: Whether to compute the aggregated graph and later save it in the layers dictionary.
+            n_first_order: The number of first-order nodes the node sequences refer to.
+                Defaults to the number of IDs in `mapping`.
         """
         # Lift order
-        if edge_weight is None:
-            ho_index = lift_order_edge_index(edge_index, num_nodes=node_sequence.size(0))
-        else:
-            ho_index, edge_weight = lift_order_edge_index_weighted(
-                edge_index, edge_weight=edge_weight, num_nodes=node_sequence.size(0), aggr=aggr
-            )
-        node_sequence = torch.cat([node_sequence[edge_index[0]], node_sequence[edge_index[1]][:, -1:]], dim=1)
+        ho_index, node_sequence, edge_weight = lift_order_step(
+            edge_index, node_sequence, edge_weight=edge_weight, aggr=aggr
+        )
 
         # Aggregate
         if save:
-            gk = aggregate_edge_index(ho_index, node_sequence, edge_weight)
-            gk.mapping = IndexMap([tuple(mapping.to_ids(v.cpu())) for v in gk.data.node_sequence])
+            gk = HigherOrderGraph.aggregate(
+                ho_index,
+                node_sequence,
+                first_order_mapping=mapping,
+                edge_weight=edge_weight,
+                n_first_order=n_first_order,
+            )
         else:
             gk = None
         return ho_index, node_sequence, edge_weight, gk
@@ -156,13 +167,16 @@ def from_temporal_graph(
         else:
             edge_weight = torch.ones(edge_index.size(1), device=edge_index.device)
         if cached or max_order == 1:
-            m.layers[1] = aggregate_edge_index(
-                edge_index=edge_index, node_sequence=node_sequence, edge_weight=edge_weight
+            m.layers[1] = HigherOrderGraph.aggregate(
+                edge_index=edge_index,
+                node_sequence=node_sequence,
+                edge_weight=edge_weight,
+                first_order_mapping=g.mapping,
+                n_first_order=g.n,
             )
-            m.layers[1].mapping = g.mapping
 
         if max_order > 1:
-            node_sequence = torch.cat([node_sequence[edge_index[0]], node_sequence[edge_index[1]][:, -1:]], dim=1)
+            node_sequence = lift_node_sequence(edge_index, node_sequence)
             if event_graph is None:
                 edge_index = EventGraph.build_edge_index(g, delta)
             else:
@@ -171,11 +185,12 @@ def from_temporal_graph(
 
             # Aggregate
             if cached or max_order == 2:
-                m.layers[2] = aggregate_edge_index(
-                    edge_index=edge_index, node_sequence=node_sequence, edge_weight=edge_weight
-                )
-                m.layers[2].mapping = IndexMap(
-                    [tuple(g.mapping.to_ids(v.cpu())) for v in m.layers[2].data.node_sequence]
+                m.layers[2] = HigherOrderGraph.aggregate(
+                    edge_index=edge_index,
+                    node_sequence=node_sequence,
+                    edge_weight=edge_weight,
+                    first_order_mapping=g.mapping,
+                    n_first_order=g.n,
                 )
 
             for k in range(3, max_order + 1):
@@ -186,6 +201,7 @@ def from_temporal_graph(
                     edge_weight=edge_weight,
                     aggr="src",
                     save=cached or k == max_order,
+                    n_first_order=g.n,
                 )
                 if cached or k == max_order:
                     m.layers[k] = gk  # type: ignore[assignment]
@@ -251,17 +267,24 @@ def from_path_data(
         elif mode == "propagation":
             aggr = "src"
 
-        m.layers[1] = aggregate_edge_index(edge_index=edge_index, node_sequence=node_sequence, edge_weight=edge_weight)
-        m.layers[1].mapping = path_data.mapping
+        g1 = aggregate_edge_index(edge_index=edge_index, node_sequence=node_sequence, edge_weight=edge_weight)
+        g1.mapping = path_data.mapping
+        # Nodes that are not traversed by any path are not part of the aggregated graph,
+        # so the first-order node set can be larger than the order-1 layer.
+        n_first_order = max(path_data.mapping.num_ids(), g1.n)
+        m.layers[1] = HigherOrderGraph.from_aggregated_graph(
+            g1, first_order_mapping=path_data.mapping, n_first_order=n_first_order
+        )
 
         for k in range(2, max_order + 1):
             edge_index, node_sequence, edge_weight, gk = MultiOrderModel.iterate_lift_order(
                 edge_index=edge_index,
                 node_sequence=node_sequence,
-                mapping=m.layers[1].mapping,
+                mapping=path_data.mapping,
                 edge_weight=edge_weight,
                 aggr=aggr,
                 save=cached or k == max_order,
+                n_first_order=n_first_order,
             )
             if cached or k == max_order:
                 m.layers[k] = gk  # type: ignore[assignment]
@@ -563,7 +586,7 @@ def to_dbgnn_data(self, max_order: int = 2, mapping: str = "last") -> Data:
         edge_index_max_order = g_max_order.data.edge_index
         edge_weight = g.data.edge_weight
         edge_weight_max_order = g_max_order.data.edge_weight
-        bipartite_edge_index = generate_bipartite_edge_index(g, g_max_order, mapping=mapping, device=edge_index.device)
+        bipartite_edge_index = g_max_order.bipartite_edge_index(g, mapping=mapping, device=edge_index.device)
 
         if g.data.y is not None:
             y = g.data.y
diff --git a/src/pathpyG/utils/dbgnn.py b/src/pathpyG/utils/dbgnn.py
deleted file mode 100644
index a70ad8a7..00000000
--- a/src/pathpyG/utils/dbgnn.py
+++ /dev/null
@@ -1,46 +0,0 @@
-"""Utils for DBGNN models."""
-
-from typing import Optional
-
-import torch
-
-from pathpyG.core.graph import Graph
-
-
-def generate_bipartite_edge_index(
-    g: Graph, g2: Graph, mapping: str = "last", device: Optional[torch.device] = None
-) -> torch.Tensor:
-    """Generate edge_index for bipartite graph connecting nodes of a second-order graph to first-order nodes.
-    
-    The mapping strategy determines to which first-order nodes the second-order nodes are connected:
-    - "last": Connects each second-order node to the last node in its sequence.
-    - "first": Connects each second-order node to the first node in its sequence.
-    - "both": Connects each second-order node to both the first and last nodes in its sequence.
-
-    !!! warning "Only for Second-Order Graphs"
-        This function is intended to be used with second-order graphs only. 
-        It does not support the use of higher-order graphs, such as third-order graphs or beyond.
-
-    Args:
-        g (Graph): The first-order graph.
-        g2 (Graph): The second-order graph.
-        mapping (str, optional): The mapping strategy to use. Options are "last", "first", or "both". Defaults to "last".
-        device (torch.device, optional): The device to place the tensor on. Defaults to None.
-
-    Returns:
-        torch.Tensor: The edge_index tensor for the bipartite graph.
-    """
-    if mapping == "last":
-        bipartide_edge_index = torch.tensor([list(range(g2.n)), [v[1] for v in g2.data.node_sequence]], device=device)
-    elif mapping == "first":
-        bipartide_edge_index = torch.tensor([list(range(g2.n)), [v[0] for v in g2.data.node_sequence]], device=device)
-    else:
-        bipartide_edge_index = torch.tensor(
-            [
-                list(range(g2.n)) + list(range(g2.n)),
-                [v[0] for v in g2.data.node_sequence] + [v[1] for v in g2.data.node_sequence],
-            ],
-            device=device,
-        )
-
-    return bipartide_edge_index
diff --git a/tests/core/test_event_graph.py b/tests/core/test_event_graph.py
index 94dbf523..039bc734 100644
--- a/tests/core/test_event_graph.py
+++ b/tests/core/test_event_graph.py
@@ -7,6 +7,7 @@
 from torch_geometric.data import Data
 
 from pathpyG.core.event_graph import EventGraph
+from pathpyG.core.higher_order_graph import HigherOrderGraph
 from pathpyG.core.index_map import IndexMap
 from pathpyG.core.multi_order_model import MultiOrderModel
 from pathpyG.core.temporal_graph import TemporalGraph
@@ -218,6 +219,35 @@ def test_multi_order_model_construction(event_graph, temporal_graph):
         )
 
 
+def test_higher_order_graph_from_weighted_event_graph(temporal_graph):
+    """Aggregating an EventGraph respects the edge weights of the temporal graph.
+
+    Regression test: order 2 used to count each continuation once instead of carrying
+    the weight of the event it starts from, so it disagreed with the temporal-graph
+    route for every order but 2.
+    """
+    temporal_graph.data.edge_weight = torch.tensor([2.0, 5.0, 11.0, 7.0])
+    event_graph = EventGraph.from_temporal_graph(temporal_graph, delta=DELTA)
+
+    for k in (1, 2, 3):
+        from_eg = HigherOrderGraph.from_event_graph(event_graph, order=k)
+        from_tg = MultiOrderModel.from_temporal_graph(temporal_graph, delta=DELTA, max_order=k).layers[k]
+
+        assert from_eg.order == k
+        assert from_eg.nodes == from_tg.nodes
+        assert torch.equal(
+            from_eg.data.edge_index.as_tensor(),
+            from_tg.data.edge_index.as_tensor(),
+        )
+        assert torch.equal(from_eg.data.edge_weight, from_tg.data.edge_weight)
+
+    # The weights must actually reflect the temporal graph, not just agree with each other.
+    order_2 = HigherOrderGraph.from_event_graph(event_graph, order=2)
+    assert order_2.nodes == [("a", "b"), ("b", "c"), ("b", "d"), ("c", "e")]
+    # (a,b)->(b,c) carries the weight of event (a->b)@1, (b,c)->(c,e) that of (b->c)@2
+    assert order_2.data.edge_weight.tolist() == [2.0, 5.0]
+
+
 def test_to_device(event_graph):
     """Moving an EventGraph moves its underlying TemporalGraph too."""
     moved = event_graph.to(torch.device("cpu"))
diff --git a/tests/nn/test_dbgnn.py b/tests/nn/test_dbgnn.py
index 0de6d8a4..418f8084 100644
--- a/tests/nn/test_dbgnn.py
+++ b/tests/nn/test_dbgnn.py
@@ -5,7 +5,6 @@
 
 from pathpyG.core.multi_order_model import MultiOrderModel
 from pathpyG.nn.dbgnn import DBGNN
-from pathpyG.utils.dbgnn import generate_bipartite_edge_index
 
 
 def test_bipartite_edge_index(simple_walks):
@@ -17,13 +16,13 @@ def test_bipartite_edge_index(simple_walks):
     print(g2.data.edge_index)
     print(g2.mapping)
 
-    bipartite_edge_index = generate_bipartite_edge_index(g, g2, mapping="last")
+    bipartite_edge_index = g2.bipartite_edge_index(g, mapping="last")
     print(bipartite_edge_index)
     # ensure that A,C and B,C are mapped to C, C,D is mapped to D and C,E is mapped to E
 
     assert equal(bipartite_edge_index, tensor([[0, 1, 2, 3], [2, 2, 3, 4]]))
 
-    bipartite_edge_index = generate_bipartite_edge_index(g, g2, mapping="first")
+    bipartite_edge_index = g2.bipartite_edge_index(g, mapping="first")
     print(bipartite_edge_index)
     # ensure that A,C is mapped A, B,C is mapped to B, and C,D and C,E are mapped to C