Computational Generation of Substrate-Specific Molecular Cages
Abstract
In this paper, we propose a method to build molecular cages designed to capture a specific substrate. We model a cage as a graph of atoms with coordinates in space, and several constraints on their edges (degree, length and angle). We use a simple method to place binding patterns which are able to interact with certain parts of the substrate. We then propose an algorithm which considers all possible ways of connecting these binding patterns and try to construct the smallest possible molecular paths realizing these connections. We investigate many variants of our method in order to obtain the most efficient algorithm, able to build cages of more than a hundred atoms.
Keywords and phrases:
Enumeration, Molecular Cage, Cheminformatics, Geometric Algorithms, Experimental AlgorithmsCopyright and License:
2012 ACM Subject Classification:
Theory of computation Design and analysis of algorithms ; Applied computing ChemistrySupplementary Material:
Software (Source code): https://github.com/NoeDemange/MolecularCagesGeneration [11]archived at
swh:1:dir:4724ee746285e388bfbb126e84babf8a21e1aecc
Acknowledgements:
Authors want to thank Olivier David for his helpful comments.Editors:
Martin Aumüller and Irene FinocchiSeries and Publisher:
Leibniz International Proceedings in Informatics, Schloss Dagstuhl – Leibniz-Zentrum für Informatik
1 Introduction
Molecular cages are a class of discrete, three-dimensional molecular structures formed by the self-assembly of building blocks into closed frameworks that define an internal cavity capable of hosting guest molecules. These architectures are characterized by well-defined size, shape, and connectivity determined by their constituent components and bonding topology. Since pioneering works of Nobel-prize winning chemists D.J. Cram, J.-M. Lehn and C.J. Pedersen [10, 19, 23], recent years have experienced spectacular growth of interest for molecular architectures possessing a defined inner-space [20]. The design and synthesis of such molecular cages remain major challenges, as they involve exploring a vast combinatorial space of possible building blocks, connectivities, and three-dimensional arrangements. Anticipating which combinations will lead to a target architecture that can be reliably realized through self-assembly requires navigating a complex design space, where small variations in topology or geometry may result in fundamentally different outcomes. Despite this complexity, successfully synthesized molecular cages exhibit a broad spectrum of applications [32], ranging from the absorption of gases, such as carbon dioxide () and methane [15], to their roles in medicine [27, 29] and the containment of hazardous substances, exemplified by the storage of white phosphorus [21].
In this context, cheminformatics and computational modeling provide essential tools to address the complexity of molecular cage design. Although such models offer only coarse abstractions of chemical reality, they enable the systematic exploration of large design spaces that are inaccessible to intuition or experimentation alone. In particular, molecular structures are commonly represented as graphs in which vertices represent atoms and edges represent chemical bonds. This graph structure can be further enriched with additional information, such as atom and bond types as well as spatial coordinates. Such graph-based representations, capturing both connectivity and geometry, have long provided a foundation for reasoning about molecular structure, similarity and assembly [24, 28].
More recently, advances in algorithmic approaches and data-driven methods, exemplified by breakthroughs such as AlphaFold in protein structure prediction [17], have demonstrated the potential of computational models to provide guidance to experimental chemists, even in the presence of significant physical simplifications. Within the community of experimental algorithms, several recent contributions have similarly combined heuristic search, combinatorial modeling, and empirical evaluation to study complex chemical problems [1, 12, 2].
Two complementary paradigms exist for constructing molecular cages for substrates. A substrate denotes the molecule that is intended to be captured by a molecular cage which is a connected molecular structure that encloses an internal cavity capable of accommodating the substrate. The first paradigm, traditionally used in supramolecular chemistry, follows a host-first strategy: a cage is designed and synthesized a priori, and its ability to encapsulate different substrates is evaluated afterward. This approach has been widely studied, especially for highly symmetric cages that are easier to synthesize and characterize [4, 32, 31]. Although it reliably produces stable and experimentally accessible structures, it often results in limited substrate specificity, since the guest is not considered during the design process. Consequently, selecting an appropriate host for a given target substrate remains largely empirical.
The second paradigm adopts a guest-driven strategy. Starting from a target substrate, the objective is to design a molecular cage that is both compatible with and specific to that substrate. In this context, a molecule is considered as a cage for a given substrate if it is shape-complementary and forms favorable non-covalent interactions with it. Compatibility requires chemical interactions between host and guest, while specificity implies preferential binding, mainly driven by geometric complementarity between the cage interior and the shape of the substrate.
1.1 Contributions
We introduce a guest-driven method for generating molecular cages that targets a given substrate. Our approach guarantees substrate specificity by positioning binding patterns around the substrate, following a workflow introduced in [7], as explained in Section 2.
We propose several variants of an algorithm to construct molecular paths that connect these binding patterns in Section 3. These methods aim to minimize path length, which naturally enforces geometric complementarity between the constructed cage and the shape of the substrate, thereby increasing substrate recognition specificity. Moreover, shorter paths facilitate synthesis in the laboratory, by limiting the structural complexity.
To structure this construction process, we introduce an intermediate combinatorial object in Section 4, the interconnection tree, which specifies which binding patterns must be connected. We provide an efficient enumeration algorithm, both in terms of theoretical complexity and practical performance, to generate such trees and to use them as a backbone for cage assembly.
Beyond the proposed modeling framework and algorithms, a major contribution of this work lies in the extensive experimental evaluation of our methods. In Section 5, we systematically explore a wide range of parameters, including discretization levels, distance-based construction heuristics, branching factor in exhaustive search, and the use of interconnection tree weights for guidance. These experiments are conducted on real molecular substrates extracted from the Cambridge Structural Database [9], allowing us to identify parameter settings that are effective in practice. As a result, our approach can automatically generate molecular cages for a wide diversity of substrates within reasonable computation times.
2 Modeling and Definitions
This section introduces the core objects used throughout the paper. Since we operate on molecular structures, we adopt a graph-based model that explicitly incorporates chemical constraints. Atoms and covalent bonds are represented as vertices and edges, respectively, and the embedding of the graph in three-dimensional space is used to capture geometric and chemical feasibility.
In general, substrates may contain a wide variety of atom types. In this work, we restrict ourselves to organic molecules. For the construction of molecular cages, we consider a limited set of atom types, namely carbon, hydrogen, oxygen, and nitrogen, with a predominance of carbon and hydrogen. This restriction simplifies the model while remaining expressive enough for the targeted applications. The algorithms themselves are generic and can be extended to additional atom types.
We introduce the concept of a molecular graph, which encodes all the relevant information about the molecules under study (see an example in Figure 1).
Definition 1.
A molecular graph is a triple , where is a graph with vertex set (atoms) and edge set (covalent bonds). The function assigns coordinates to each vertex, and the function assigns to each vertex an atom type from a finite set .
A covalent bond corresponds to a pair of atoms sharing electrons, and its length depends on the types of the atoms involved. We denote by the length of the covalent bond of atoms of types and . In this article, we use a simplified model with two bond lengths: for bonds involving hydrogen atoms, and for bonds between non-hydrogen atoms.
The collision distance models steric repulsion, which prevents atoms from being placed too close to each other due to overlap of their electron clouds. We denote by the collision distance of atoms of types and . In this article, we use a uniform collision distance of .
The degree of a vertex corresponds to the number of covalent bonds it can form and is bounded by the valence of its atom type. Since we consider organic molecules, this degree is at most 4, but may be smaller depending on the atom type (see Table 2). Formally, the degree of a vertex is bounded by .
To lighten the presentation, we denote the distance between two vertices by instead of .
We define the notion of a chemically realistic molecular graph to capture the set of chemical constraints that such a graph must satisfy.
Definition 2.
| Atom | #bonds |
|---|---|
| Carbon | 4 |
| Nitrogen | 3 |
| Oxygen | 2 |
| Hydrogen | 1 |
| Geometry | 2D | 3D | Angle | Margin |
| Tetrahedral | ![]() |
![]() |
109.5° | 3° |
| Triangular | ![]() |
120° | 2° | |
| Linear | 180° | – |
Binding Patterns
Efficient substrate recognition requires a molecular cage to establish multiple intermolecular interactions in order to ensure both binding strength and selectivity. We focus on organic substrates and consider two key classes of non-covalent interactions: hydrogen bonds and – stacking interactions.
A hydrogen bond [3] involves a hydrogen atom bound to an electronegative donor and a second electronegative acceptor atom, typically oxygen or nitrogen. Its geometry is constrained in both distance and orientation (Figure 3). In our model, when the substrate provides a donor (resp. acceptor), we introduce a complementary acceptor (resp. donor) as a hydrogen binding pattern. An atom set can act as donor or acceptor, but not both, and can participate in at most one hydrogen bond.
A – stacking interaction [22] occurs between approximately planar and parallel aromatic systems at a characteristic distance. We model such interactions by planar cyclic subgraphs representing aromatic rings (Figure 3).
Binding pattern placement follows the approach of Bricage [7]. For a given substrate, multiple valid placements may exist. We first identify all substrate sites eligible for hydrogen bonding or – stacking based on atom types and local geometry. For each site, we generate candidate placements satisfying the geometric constraints of the interaction and discard those that collide with the substrate.
– stacking patterns are placed first, as they are spatially extended, strongly constrain the geometry, and rarely interfere with one another. Hydrogen patterns are then considered. Their placement is subject to geometric and chemical incompatibilities, which are encoded in a conflict graph whose vertices represent patterns and edges represent conflicts. Selecting a compatible set of patterns corresponds to computing an independent set in this graph. We enumerate all maximal independent sets using the Bron–Kerbosch algorithm [8]. Maximal independent sets correspond to sets of binding patterns that cannot be extended by adding another compatible pattern, and therefore represent configurations that maximize the number of interactions. In practice, we generate all maximal independent sets and then select representative ones to construct candidate cages. In future work, one could systematically evaluate all maximal sets or also consider non-maximal sets when a smaller number of interactions is desired.
This process yields a disconnected molecular graph of mutually compatible binding patterns, which are in interaction with the substrate. The graph is typically disconnected, since the binding patterns do not form covalent bonds between each other at this stage. The next step is therefore to connect these patterns using molecular paths in order to obtain a connected molecular cage. Section 3 describes how these patterns can be connected by molecular paths to form a complete cage while satisfying chemical and geometric constraints.
Connection Patterns and Paths
Cage construction relies on basic building blocks, called connection patterns. A connection pattern is a molecular graph with two distinguished vertices, denoted and , corresponding to atoms that can form covalent bonds with other patterns. These vertices define entry and exit points, allowing patterns to be assembled into molecular paths that connect binding patterns. The vertices and may coincide, as in the case of a single central atom, such as the tetrahedral carbon extensively used in this work.
Definition 3.
Let be a sequence of connection patterns, and let and denote the distinguished vertices of . We define as the union of the graphs , augmented with edges for all . If is a molecular graph, the sequence is called a molecular path.
3 Construction of a Molecular Path
In this section, we introduce the concept of constructing a molecular path and the underlying Molecular Path Construction problem.
We consider two molecular graphs: representing a substrate and representing binding patterns and a partially constructed cage, together with two vertices and of corresponding to atoms of two binding patterns to be connected. A molecular path connects to if it can be attached to at these vertices and satisfies the following structural constraints:
-
1.
is a chemically realistic molecular graph;
-
2.
let be the distinguished vertex of , ;
-
3.
let be the distinguished vertex of , ;
-
4.
for all and all , .
The first constraint ensures chemical plausibility of the constructed structure, while constraints (2) and (3) enforce the attachment of the path to the binding patterns at atoms and . A weak bond is an attractive, non-covalent interaction either within a molecule or between molecules. While different types of weak interactions exist and may correspond to different distances, we adopt a simplified model in which such interactions occur at a fixed distance of . Accordingly, in Constraint (4), denotes the minimum distance required for weak interactions. Since is chosen such that , this constraint is stronger than the non-collision requirement implied by Constraint (1) and ensures that the only weak interactions present are those explicitly modeled by the binding patterns.
In practice, constructing a molecular path that exactly satisfies all geometric constraints is generally impossible, due to the flexibility of molecular conformations. We therefore allow small deviations at the final attachment to , bounded by for bond angles and for bond lengths. These tolerances are fixed for simplicity and because the considered paths are short (fewer than 15 atoms). For longer paths, it would be natural to scale them with path length. Rather than distributing these deviations along the entire structure, as in real molecules, we approximate this behavior by concentrating them at a single atom, which simplifies both construction and analysis.
We quantify the resulting distortion using the normalized root mean square deviation (NRMSD) of bond angles relative to VSEPR theory and of bond lengths relative to ideal covalent distances, providing a global measure of deviation from an ideal configuration. The terminal vertex must satisfy the above angular and bond length tolerances with respect to the target, yielding a worst-case NRMSD of , computed as
We now define the Molecular Path Construction problem as follows: given the substrate , a partial cage , and two vertices , find a molecular path connecting to that minimizes in lexicographic order, where denotes the length of .
For clarity of presentation, we restrict Molecular Path Construction to a single type of connection pattern: a tetrahedral carbon atom (see Figure 4(b)). This choice is sufficient to generate realistic cages, and all proposed methods naturally extend to multiple pattern types.
3.1 Valid Positions of the Next Connection Pattern
To solve the Molecular Path Construction problem, we propose a depth-first exploration of the solution space, guided by distance-based heuristics. The molecular path is constructed incrementally by successively adding connection patterns: At each step, a new pattern is appended to the end of a partial path , producing a longer partial solution.
A key challenge arises from the continuous nature of 3D space: even when geometric constraints are fixed, there is an infinite number of admissible positions for the next atom. To make the search tractable, we first discretize the space of possible positions around the relevant edge rotation. In this sense, the algorithm can be interpreted as a branch-and-bound search over the discretized rotation space, where partial paths are extended only along positions that are geometrically feasible and heuristically promising.
The main problem addressed in the remainder of this section is how to generate well-distributed, legal positions for the next connection pattern and how to select among them those that are most likely to yield shorter, chemically plausible paths.
First, we explain how to determine the legal positions to add an atom, according to VSEPR theory and fixed bond lengths. Once two consecutive atoms are known, the next atom must lie on a circle in 3D space (Figure 4(a)).
When the position of this next atom is fixed, the two missing hydrogen atoms required to complete the tetrahedral carbon pattern, centered on , the central atom of the pattern , are added automatically (Figure 4(b)). We must also check whether these hydrogens satisfy collision constraints. Thus, placing a single atom induces the placement of two additional atoms, increasing the complexity of collision detection.
A first, naive approach consists in discretizing the placement circle and independently testing each candidate position for collisions with existing atoms. This results in a large number of queries to the Spherical Range Emptiness (or Threshold Nearest Neighbor) problem [5], defined as follows: given a finite set , a query point , and a threshold , decide whether .
When the two additional hydrogen atoms are taken into account, the number of such queries becomes too large. Instead, we exploit the simple geometric structure of the feasible positions on the placement circle. Each existing atom induces a forbidden angular interval on the circle corresponding to positions that violate the collision constraint. Moreover, the number of atoms that can generate such intervals is bounded by a constant, since only atoms lying within collision distance of the circle need to be considered. As a result, the set of valid positions can be represented by a bounded set of angular intervals.
For each atom, we compute and exclude the forbidden interval, yielding a set of admissible angular regions. The same construction is applied to the two hydrogen atoms, whose positions depend deterministically on the chosen angle via fixed angular offsets derived from the tetrahedral angle (). The final set of valid placements is obtained by intersecting the admissible angular regions of all three atoms we place.
By storing the atoms of in an appropriate spatial data structure, we can identify the atoms close to the placement circle in constant time, since the number of atoms within a fixed radius is bounded by a constant. Since all subsequent operations involve only a bounded number of atoms and angular intervals, the overall complexity of this procedure is constant. Further details of this construction are provided in Appendix A.
3.2 Selecting the Best Positions
The admissible angular intervals correspond to infinitely many feasible placements. To obtain a finite set of candidates, we discretize these intervals and guide the selection of angles using distance-based heuristics that estimate how promising a placement is. We consider three such heuristics, each inducing a different discretization strategy. The first relies on the Euclidean distance to the target, which is inexpensive to compute but ignores collision constraints and, therefore, provides only a coarse estimate. The second uses a distance computed on a discretized collision-aware grid, yielding a more accurate evaluation of the remaining path length but at a significantly higher computational cost. The third heuristic combines these two distances in order to balance efficiency and accuracy, leveraging the speed of the Euclidean metric only when relevant.
Euclidean distance.
The Euclidean distance from a candidate position to the target vertex is inexpensive to compute, but ignores obstacles. It performs well when is in direct line of sight, but may fail in cluttered environments where collisions prevent straight-line connections as in Figure 5.
To generate candidate positions, we exploit the fact that , the angle minimizing the Euclidean distance to , can be computed analytically. Let denote the total angular width of the valid intervals and the desired number of angular samples. The angular step is , bounded below by some minimum angular spacing. The first candidate angle is chosen as closely as possible to within the intervals of valid positions. The following positions are computed from the first by incrementing the angular step symmetrically on both sides, up to the limit imposed by the branching factor. If a sampled angle lies in a forbidden interval, it is moved to the closest valid angle in the sampling direction, without reversing direction.
Discretized Distance.
To better account for obstacles, we use a collision-aware distance computed on a discretized representation of space. Candidate positions are first selected within the valid angular intervals as follows: each interval initially contains a single uniformly placed sample, and additional samples are iteratively added to the interval whose current subdivision contains the longest segment. This adaptive refinement yields a well-distributed set of candidate positions.
We then evaluate each candidate by estimating its distance to the target using a voxel-based grid that excludes colliding cells. The voxel grid is constructed from the bounding box of the molecule, defined by the minimum and maximum coordinates along each axis and expanded by a margin of . The grid resolution is set to , and grid vertices that lie closer than a given threshold to any atom (modeled as spheres) are marked as unusable. Each free voxel corresponds to a node in a graph, and edges connect voxels in the 26-neighborhood, with weights equal to the Euclidean distance between voxel centers. The graph is augmented with the candidate positions and the target, each connected to their nearest grid neighbors. Shortest-path distances are computed using the algorithm [14]. Since all candidates share the same target, we also evaluate the single-source multi-target variant () [16].
Hybrid Distance.
We also consider a hybrid heuristic that either uses the Euclidean distance when no obstacle blocks visibility from the current connection pattern to the target, and switches to the discretized distance otherwise. Visibility is determined by checking whether the line segment between the current connection pattern and the target intersects any atom, where atoms are modeled as spheres.
3.3 Bounding the Search Space
Branching Factor.
To limit the combinatorial explosion of the visited partial molecular path, we limit the number of positions retained at each step by a branching factor. Only this number of positions with the best heuristic scores are kept and explored further. This balances completeness and efficiency, allowing us to find valid molecular paths while keeping computation times reasonable.
Size Pruning.
To prevent unbounded exploration during molecular path construction, we introduce cut-offs that limit the number of patterns in a path. In its simplest form, this limit is a fixed maximum number of patterns defined a priori. Although this guarantees termination, it does not favor the shortest paths. To do so, we refine this limit by keeping track of the length (in number of patterns) of the shortest path found so far between the start and end vertices. Any newly constructed partial path whose length exceeds this bound is discarded. We refer to this rule as the Min Length Cut.
We further strengthen this criterion by estimating the minimum number of additional patterns required to complete the current partial path. This estimate is obtained from the Euclidean distance to the destination and the minimal number of patterns needed to span that distance. If the resulting projected path length exceeds the length of the shortest path found so far, the branch is pruned. We call this rule the Projected Length Cut. Together, these cuts effectively restrict the exploration on the shortest molecular paths.
4 Interconnection Trees
We have shown how to generate a molecular path between two vertices. To construct a complete cage interacting with the substrate, all binding patterns must be interconnected to form a single connected structure. In this section, we introduce the notion of an interconnection tree, which formalizes how binding patterns can be linked. We then present an efficient enumeration algorithm for interconnection trees. Combined with the molecular path construction described in the previous section, this approach enables the generation of complete molecular cages.
4.1 Interconnection Tree
We model the selection of binding patterns to connect as a graph-theoretic problem, abstracting away geometric and chemical constraints. We consider a complete multipartite graph , where the vertex set is partitioned into and . Each part represents a binding pattern, and each vertex in corresponds to an atom that can serve as an endpoint of a molecular path.
We define the function that maps each vertex to the part to which it belongs. This definition naturally extends to edges by setting . The quotient graph of (see [6]) is then defined as
Our objective is to connect all parts in while ensuring that each vertex of is used at most once, since an atom can participate in only one molecular path. Moreover, we aim to minimize the number of connections, which leads to simpler structures and facilitates chemical synthesis. This leads to the following definition, where the symbol denotes disjoint union.
Definition 4.
An interconnection tree of a multipartite graph is a set such that is a matching in and is a spanning tree of .
Examples of interconnection trees are shown in Figure 6.
In general, deciding whether a multipartite graph admits an interconnection tree is -hard, by reduction from the Hamiltonian path problem (see Appendix B). However, for complete multipartite graphs, existence admits a simple characterization.
Lemma 5 (Proof in Appendix B).
Let be a complete multipartite graph. Then admits an interconnection tree if and only if .
4.2 Enumeration of Interconnection Trees
Let denote the set of interconnection trees of . We propose an efficient algorithm to enumerate all elements of . By efficient, we mean that the delay between the output of two consecutive solutions is bounded by the size of a solution, a standard complexity measure in enumeration problems [26].
A classical approach to enumeration relies on binary partitioning: solutions are divided according to whether they contain a given edge. Let denote the graph obtained from by removing an edge , and let denote the contraction of , where vertices and are identified and their corresponding parts are merged.
Lemma 6.
Let be a multipartite graph. A set is an interconnection tree of containing an edge if and only if is an interconnection tree of .
Let be a set of edge sets and let be an edge , we write . By Lemma 6, we can decompose into two disjoint sets as follows
A similar decomposition is used for enumerating spanning trees [30], yielding an algorithm with linear delay and constant amortized delay. However, such algorithm requires efficiently deciding whether the recursive subproblems are non-empty in order to avoid useless recursion (a technique known as flashlight search). In our setting, even if is complete multipartite, the graph is not, and deciding whether is empty becomes -complete.
To address this issue, we refine the decomposition so that all recursive subproblems remain complete multipartite graphs, allowing us to apply Lemma 5. We fix an arbitrary ordering of the vertices of that respects the partition: if , , and , then . For a vertex , we let denote the graph obtained from by removing all vertices . This leads to the following decomposition:
| (1) |
Equation 1 can be proved inductively. By induction hypothesis, we have the equation with an union over some initial segment of the edges and an additional term equal to . We let be the first edge in the order after and we apply Lemma 6 to and the edge to prove the induction.
Equation 1 directly yields a recursive enumeration algorithm. With appropriate data structures, we can maintain in constant time whether each part is non-empty and whether the condition of Lemma 5 is satisfied.
Theorem 7 (Proof in Appendix B).
Let be a complete multipartite graph with parts. Then can be enumerated with worst-case delay and amortized delay .
In practice, the number of interconnection trees may be very large. To prioritize trees that are more likely to yield short molecular paths, we associate a weight to each edge, corresponding to the spatial distance between the two atoms. The weight of an interconnection tree is defined as the sum of the weights of its edges.
Ideally, one would enumerate interconnection trees in non-decreasing order of weight. However, computing a minimum-weight interconnection tree is already -hard. When the number of interconnection trees remains moderate, a practical strategy, called ordered, is therefore to enumerate all trees, sort them by weight, and subsequently construct molecular paths. This approach is justified because we show in Section 5 that the cost of generating interconnection trees is several orders of magnitude smaller than that of computing molecular paths.
5 Experimental Evaluation on Real Substrates
All experiments were conducted on an AMD Ryzen 7 PRO 250 CPU with 30 GB of RAM. The code was compiled with GCC 13.3.0 using the -O2 optimization level. Unless stated otherwise, all experiments use the same parameter settings. Pruning is set to Projected Length Cut, the branching factor to 3, the number of angular samples to 12, the angular separation to 15∘, the distance function to Hybrid, and the grid step to . The maximum path length is fixed to 15, and at most solutions are generated. The source code and datasets used in this work are available in our GitHub repository111https://github.com/NoeDemange/MolecularCagesGeneration.
The initial substrates are extracted from the Cambridge Structural Database (CSD) [9], with the help of a chemist, we have selected 20 small organic molecules (fewer than 100 atoms) with varying characteristics in terms of size and shape. In what follows, we use only a representative subset of the molecules tested. Each experiment takes as input a substrate together with its associated binding patterns. Substrates are identified by their CSD Database Identifier in uppercase; when a lowercase letter is appended, it indicates that additional paths have been introduced so that only a single molecular path remains to be constructed.
To the best of our knowledge, there is currently no existing algorithmic approach for generating molecular cages from a given substrate, which prevents direct comparison with prior work. Our evaluation is therefore conducted on instances proposed by a chemist. While molecular cages have been experimentally constructed for some substrates, these structures are not generated computationally, and no standardized metrics or computational baselines are available to quantitatively compare their quality with our approach.
5.1 Molecular Path Generation
In this section, we evaluate the influence of the main parameters involved in the generation of a single molecular path.
Pruning.
We first compare the baseline algorithm with the two pruning strategies introduced in Section 3.3, namely the Min Length Cut and the Projected Length Cut, in order to assess their impact on computational efficiency.
The results are reported in Table 3. As expected, both cuts significantly reduce the number of explored solutions by favoring shorter paths. The Projected Length Cut explores the fewest partial paths and consistently achieves the lowest running time, making it the most effective pruning strategy in practice.
| Instance | Cut | #Paths | #Partial paths | Time (ms) |
|---|---|---|---|---|
| ABABELa | Baseline | 45 | 30 839 | 145 |
| Min Length | 10 | 7 009 | 35 | |
| Projected Length | 10 | 6 225 | 22 | |
| BAHSUYa | Baseline | 5 | 138 | 1.8 |
| Min Length | 5 | 120 | 1.8 | |
| Projected Length | 5 | 120 | 1.8 | |
| YILLAGc | Baseline | 41 | 5 662 | 47 |
| Min Length | 10 | 1 276 | 12 | |
| Projected Length | 10 | 1 134 | 8 |
| Instance | #Pos. | Path length | #Paths | NRMSD | Time (ms) |
|---|---|---|---|---|---|
| ABABELa | 1 | N/A | N/A | N/A | 2.5 |
| 2 | 10 | 1 | 0.87 | 6.3 | |
| 3 | 10 | 10 | 0.53 | 22 | |
| 4 | 10 | 67 | 0.33 | 99 | |
| BAHSUYa | 1 | N/A | N/A | N/A | 2.3 |
| 2 | N/A | N/A | N/A | 1.7 | |
| 3 | 5 | 2 | 0.82 | 1.8 | |
| 4 | 5 | 3 | 0.82 | 10 | |
| YILLAGc | 1 | N/A | N/A | N/A | 1.8 |
| 2 | N/A | N/A | N/A | 2.1 | |
| 3 | 8 | 5 | 0.76 | 6.9 | |
| 4 | 7 | 32 | 0.56 | 15 |
Branching factor.
We then study the impact of the branching factor of the path generation algorithm by varying the number of retained candidate positions from 1 to 4. The results are summarized in Table 3.
When only one or two positions are retained, no valid solution is found, as the construction of paths lacks sufficient degrees of freedom to avoid obstacles and to respect the margins of the final attachment to . Increasing the branching factor to four significantly increases both the number of solutions and the computation time. A branching factor of offers the best compromise, yielding a reasonable number of solutions, moderate computational cost, and good path quality in terms of NRMSD. Consequently, we retain three positions in the remainder of this work.
Number of angular samples.
We next evaluate the influence of the number of angular samples used to discretize the angular intervals when selecting candidate positions. We test values of 24, 12, 8, and 6 samples; these values have been chosen by dividing by , , and respectively. The results are shown in Table 4.
Increasing the number of samples does not systematically improve solution quality, while reducing it does not necessarily reduce the computation time. Using too many samples often yields candidate positions that are geometrically close and offer limited diversity. Conversely, too few samples reduce the accuracy of position selection and may lead to longer paths. Moreover, when the number of samples is less than , it may be smaller than the number of valid angular intervals, and some intervals may contain no candidate position, further limiting exploration. In practice, we use 12 angular samples, although the choice between 8 and 12 remains inconclusive based on our experiments. A more extensive empirical study would be required to determine the optimal number of samples or to adapt this parameter dynamically to the input instance.
Minimum angular spacing.
We then analyze the impact of the minimum angular spacing between candidate positions. The tested values are 5, 10, 15, and 20 degrees, with results reported in Table 4. The choice of angular spacing has a limited influence on performance. However, a separation of 15 degrees consistently provides a good compromise between path quality and computational cost across all tested instances.
| Instance | #Pos. | Path length | #Paths | NRMSD | Time (ms) |
|---|---|---|---|---|---|
| ABABELa | 24 | 10 | 3 | 0.92 | 25 |
| 12 | 10 | 10 | 0.53 | 24 | |
| 8 | 9 | 2 | 1.21 | 17 | |
| 6 | 11 | 6 | 0.41 | 44 | |
| BAHSUYa | 24 | 5 | 1 | 0.82 | 2.6 |
| 12 | 5 | 2 | 0.82 | 2.1 | |
| 8 | 6 | 1 | 0.80 | 3.0 | |
| 6 | 7 | 6 | 0.74 | 5.1 | |
| YILLAGc | 24 | 8 | 5 | 0.24 | 5.8 |
| 12 | 8 | 5 | 0.76 | 8.2 | |
| 8 | 7 | 6 | 0.45 | 4.9 | |
| 6 | 7 | 1 | 1.16 | 5.2 |
| Instance | Angular spacing | Path length | #Paths | NRMSD | Time (ms) |
|---|---|---|---|---|---|
| ABABELa | 5 | N/A | N/A | N/A | 2.4 |
| 10 | 10 | 12 | 0.44 | 26 | |
| 15 | 10 | 10 | 0.53 | 22 | |
| 20 | 10 | 9 | 0.34 | 21 | |
| BAHSUYa | 5 | 6 | 4 | 0.51 | 2.1 |
| 10 | 6 | 5 | 0.51 | 1.6 | |
| 15 | 5 | 2 | 0.81 | 1.4 | |
| 20 | 5 | 1 | 1.29 | 1.7 | |
| YILLAGc | 5 | 8 | 2 | 0.43 | 3.2 |
| 10 | 7 | 3 | 0.76 | 2.6 | |
| 15 | 8 | 5 | 0.76 | 8.0 | |
| 20 | 7 | 1 | 1.03 | 4.9 |
Distance types.
Finally, we evaluate the impact of the distance function used to guide path generation. We compare the three proposed strategies: the Euclidean distance, the discretized distance, and the Hybrid distance which combines both. As an alternative implementation of to compute the discretized distance, we also evaluate , which computes the same distance but to all candidate positions at the same time. Since both methods return identical distances, we report their results in a single row, listing the execution times of both algorithms.
| Instance | Distance | #Results | Path length | #Paths | NRMSD | Time (ms) |
|---|---|---|---|---|---|---|
| ABABELa | A* / SSMTA* | 8 | 10 | 6 | 0.68 | 100 / 107 |
| Euclidean | 12 | 10 | 12 | 0.64 | 11 | |
| HYBRID | 10 | 10 | 10 | 0.53 | 23 | |
| BAHSUYa | A* / SSMTA* | 7 | 5 | 2 | 0.90 | 10 / 9 |
| Euclidean | 2 | 5 | 2 | 0.79 | 0.32 | |
| HYBRID | 5 | 5 | 2 | 0.82 | 1.9 | |
| YARZUN03a | A* / SSMTA* | 107 | 11 | 29 | 0.39 | 78 / 76 |
| Euclidean | 0 | NA | NA | NA | 0.24 | |
| HYBRID | 93 | 11 | 26 | 0.15 | 29 | |
| YILLAGc | A* / SSMTA* | 8 | 7 | 2 | 0.89 | 13 / 14 |
| Euclidean | 1 | 5 | 1 | 0.93 | 0.18 | |
| HYBRID | 10 | 8 | 5 | 0.76 | 8.3 |
As shown in Table 5, the Euclidean distance is by far the fastest to compute, but it fails to consistently produce valid paths, particularly in the presence of obstacles, as illustrated by the absence of results for instance YARZUN03a (shown in Figure 5). In contrast, both and are significantly more robust, at the cost of substantially higher computation times. Despite its theoretical advantages, does not yield a noticeable speed-up over in our experiments.
The Hybrid distance offers the best compromise between robustness and efficiency. It successfully generates valid paths in most instances while remaining considerably faster than purely graph-based distances. For this reason, we adopt the Hybrid distance in the remainder of this work.
5.2 Interconnection Tree Enumeration
We evaluate the practical performance of our interconnection tree enumeration algorithm on several instances corresponding to realistic molecular cages. Each instance is modeled as a complete multipartite graph with parts and vertices per part, and is denoted by in the experiments. The chosen values of and reflect the sizes encountered in our cage construction process.
The goal of this experiment is to assess the efficiency of the enumeration algorithm in terms of total runtime, amortized delay between successive solutions, and the overhead induced by storing interconnection trees and sorting them by weight. In our implementation, we adopt simpler data structures, which do not ensure constant amortized delay, but lead to improved practical delays on the small instances considered here.
| Instance | #Trees | Time (ms) | Delay (ns) | Sorting overhead (%) | Storage overhead (%) |
|---|---|---|---|---|---|
| (3,3) | 162 | 0.006 | 36 | 116 | 66 |
| (3,6) | 3 240 | 0.07 | 22 | 294 | 190 |
| (5,3) | 174 960 | 5.8 | 34 | 255 | 164 |
| (5,6) | 107 308 800 | 2 685 | 25 | 402 | 188 |
| (7,3) | 525 404 880 | 16 742 | 32 | N/A | N/A |
| (7,6) | 4 119 428 086 | TO | 29 | N/A | N/A |
The results reported in Table 6 confirm that the amortized delay is very small and effectively constant in practice. As expected, storing and sorting interconnection trees is more expensive than generating them, and the relative cost of sorting increases with the number of trees. Nevertheless, the time spent per interconnection tree remains negligible compared to the time required to compute the molecular paths associated with a given tree. Therefore, when the number of interconnection trees is reasonable, typically below , we can afford to enumerate all trees and sort them by weight.
5.3 Full Molecular Cage Generation
In this section, we evaluate the full pipeline on real substrates, where multiple molecular paths must be generated to connect the selected binding patterns. We evaluate the early edge removal heuristic and two generation strategies and, independently of these choices, we analyze the resulting generated cages. Figure 7 illustrates cages produced by our approach, visualized using PyMOL [25]. The molecular cages behave as expected: binding sites are positioned to interact with the substrate, and the connecting paths define a shape complementary to that of the substrate. No path obstructs the internal cavity intended to host the substrate.
Ethyl propionate.
Early edge removal heuristic.
We introduce the Early Edge Removal heuristic to accelerate the exploration of interconnection trees. When the algorithm fails to construct a molecular path for a given edge, all remaining interconnection trees containing this edge are skipped. This pruning strategy eliminates large portions of the search space that are unlikely to yield valid or compact cages.
| Instance | Mode | Edge Removal | #Cages | #Trees | Average NRMSD | MNoA | Total (s) |
|---|---|---|---|---|---|---|---|
| ABCLUA10 | OTF | no | 10 | 1 | 0.87 | 203 | 0.04 |
| OTF | yes | 10 | 1 | 0.87 | 203 | 0.04 | |
| ORD | no | 10 | 1 | 0.95 | 101 | 11.6 | |
| ORD | yes | 10 | 1 | 0.95 | 101 | 11.5 | |
| ABINOS | OTF | no | 10 | 166 | 0.63 | 151 | 15 |
| OTF | yes | 10 | 14 | 0.59 | 133 | 18 | |
| ORD | no | 10 | 4 | 0.99 | 91 | 5.9 | |
| ORD | yes | 10 | 3 | 0.99 | 91 | 5.8 | |
| ACANIL01 | OTF | no | 0 | 33 | N/A | N/A | TO |
| OTF | yes | 10 | 11 | 0.87 | 99 | 5.4 | |
| ORD | no | 10 | 500 | 0.70 | 87 | 14 | |
| ORD | yes | 10 | 14 | 0.70 | 87 | 3 |
| Instance | Mode | #Trees | #Cages | Average NRMSD | MNoA | Total (s) |
|---|---|---|---|---|---|---|
| ABCLUA10 | OTF | 1 | 17 424 | 0.67 | 191 | TO |
| ORD | 34 | 15 964 | 0.74 | 95 | TO | |
| ABINOS | OTF | 66 | 8 219 | 0.43 | 85 | TO |
| ORD | 174 | 6 834 | 0.61 | 61 | TO | |
| ACANIL01 | OTF | 50 | 564 | 0.57 | 90 | TO |
| ORD | 264 | 2 253 | 0.52 | 69 | TO | |
| ALFUCO | OTF | 26 | 280 | 0.64 | 100 | TO |
| ORD | 63 | 2 528 | 0.43 | 61 | TO | |
| BAHSUY | OTF | 6 | 156 | 0.17 | 19 | 0.06 |
| ORD | 6 | 156 | 0.17 | 19 | 0.05 |
Table 7 compares generation modes with and without this heuristic on real instances, using a time limit of 120 seconds and a maximum of 10 returned solutions. The heuristic proves highly effective: it prunes interconnection trees that would either fail entirely or lead to long paths, enabling the discovery of significantly smaller cages in substantially less time.
Generation strategies.
Finally, we compare two strategies for interconnection tree generation: on-the-fly and ordered. As shown in Table 7, when the time budget does not allow to try to construct paths for all trees enumerated, processing interconnection trees in increasing weight order yields smaller cages than the on-the-fly approach. When sufficient time is available to enumerate all trees and try to construct paths on each, both strategies find the same solutions, but the cost of sorting is negligible compared to the cost of molecular path construction.
6 Conclusion
In this paper, we introduced a computational framework to generate molecular cages tailored to a target substrate. By combining the placement of binding patterns, the construction of short molecular paths, and the enumeration of interconnection trees, our approach enables the automatic generation of chemically realistic cages while promoting geometric complementarity and substrate specificity. Extensive experimental results demonstrate that the proposed heuristics, in particular hybrid distance estimations and pruning strategies on interconnection trees, allow the construction of compact cages within reasonable computation times.
Several directions for future work emerge from this study. First, distance computation could be further improved by adapting algorithms and the spatial data structures to our specific problem, in order to obtain distance estimates that are both more accurate and faster to compute. Second, methods to rigidify the cage could be explored, for instance by completing the structure to reduce the length of molecular paths, leading to more stable and synthetically accessible architectures. Finally, allowing multiple binding patterns to be connected to a single attachment point could help further reduce molecular path lengths and increase the compactness of the generated cages.
References
- [1] Ylène Aboulfath, Dominique Barth, Thierry Mautor, Dimitri Watel, and Marc-Antoine Weisser. Polymorphic Cycle Basis in a Sequence of Graphs to Analyze the Structural Evolution of a Molecular Dynamic Trajectory. In 23rd International Symposium on Experimental Algorithms (SEA 2025), volume 338 of Leibniz International Proceedings in Informatics (LIPIcs), pages 1:1–1:14. Schloss Dagstuhl – Leibniz-Zentrum für Informatik, 2025. doi:10.4230/LIPIcs.SEA.2025.1.
- [2] Dmytro Antypov, Argyrios Deligkas, Vladimir Gusev, Matthew J. Rosseinsky, Paul G. Spirakis, and Michail Theofilatos. Crystal Structure Prediction via Oblivious Local Search. In 18th International Symposium on Experimental Algorithms (SEA 2020), pages 21:1–21:14. Schloss Dagstuhl – Leibniz-Zentrum für Informatik, 2020. doi:10.4230/LIPIcs.SEA.2020.21.
- [3] Elangannan Arunan, Gautam R. Desiraju, Roger A. Klein, Joanna Sadlej, Steve Scheiner, Ibon Alkorta, David C. Clary, Robert H. Crabtree, Joseph J. Dannenberg, Pavel Hobza, Henrik G. Kjaergaard, Anthony C. Legon, Benedetta Mennucci, and David J. Nesbitt. Definition of the hydrogen bond (IUPAC Recommendations 2011). Pure and Applied Chemistry, 83(8):1637–1641, 2011.
- [4] Dominique Barth, Olivier David, Franck Quessette, Vincent Reinhard, Yann Strozecki, and Sandrine Vial. Efficient generation of stable planar cages for chemistry. In International Symposium on Experimental Algorithms (SEA), pages 235–246. Springer, 2015. doi:10.1007/978-3-319-20086-6_18.
- [5] Jon L Bentley. A survey of techniques for fixed radius near neighbor searching. Technical report, Stanford University, 1975.
- [6] Alain Bretto, Alain Faisant, and François Hennecart. Elements of Graph Theory: From Basic Concepts to Modern Developments, volume 24 of EMS Textbooks in Mathematics. EMS Press, 1 edition, 2022. doi:10.4171/etb/24.
- [7] Marie Bricage. Modélisation et Algorithmique de graphes pour la construction de structures moléculaires. PhD thesis, Université Paris Saclay, 2018.
- [8] Coen Bron and Joep Kerbosch. Algorithm 457: finding all cliques of an undirected graph. Communications of the ACM, 16(9):575–577, 1973.
- [9] The Cambridge Crystallographic Data Centre. Cambridge structural database. URL: https://www.ccdc.cam.ac.uk/.
- [10] Donald J. Cram. The Design of Molecular Hosts, Guests, and Their Complexes (Nobel Lecture). Angewandte Chemie International Edition in English, 27(8):1009–1020, 1987.
- [11] Noé Demange. MolecularCagesGeneration. Software, swhId: swh:1:dir:4724ee746285e388bfbb126e84babf8a21e1aecc (visited on 2026-05-27). URL: https://github.com/NoeDemange/MolecularCagesGeneration, doi:10.4230/artifacts.25788.
- [12] Coline Gianfrotta, Vladimir Reinharz, Dominique Barth, and Alain Denise. A Graph-Based Similarity Approach to Classify Recurrent Complex Motifs from Their Context in RNA Structures. In 19th International Symposium on Experimental Algorithms (SEA 2021), pages 19:1–19:18. Schloss Dagstuhl – Leibniz-Zentrum für Informatik, 2021. doi:10.4230/LIPIcs.SEA.2021.19.
- [13] Ronald J. Gillespie. The valence-shell electron-pair repulsion (VSEPR) theory of directed valency. Journal of Chemical Education, 40(6):295, 1963.
- [14] Peter E. Hart, Nils J. Nilsson, and Bertram Raphael. A Formal Basis for the Heuristic Determination of Minimum Cost Paths. IEEE Transactions on Systems Science and Cybernetics, 4(2):100–107, 1968. doi:10.1109/TSSC.1968.300136.
- [15] James R. Holst, Abbie Trewin, and Andrew I. Cooper. Porous organic molecules. Nature Chemistry, 2(11):915–920, 2010.
- [16] Htoo Htoo, Yutaka Ohsawa, Noboru Sonehara, and Masao Sakauchi. Incremental Single-Source Multi-Target A* Algorithm for LBS Based on Road Network Distance. IEICE Transactions on Information and Systems, E96.D(5):1043–1052, 2013. doi:10.1587/TRANSINF.E96.D.1043.
- [17] John Jumper, Richard Evans, Alexander Pritzel, Tim Green, Michael Figurnov, Olaf Ronneberger, Kathryn Tunyasuvunakool, Russ Bates, Augustin Žídek, Anna Potapenko, et al. Highly accurate protein structure prediction with alphafold. nature, 596(7873):583–589, 2021.
- [18] Richard M. Karp. Reducibility among combinatorial problems. In Complexity of Computer Computations: Proceedings of a symposium on the Complexity of Computer Computations, pages 85–103. Springer, 1972. doi:10.1007/978-1-4684-2001-2_9.
- [19] Jean-Marie Lehn. Supramolecular Chemistry – Scope and Perspectives Molecules, Supermolecules, and Molecular Devices (Nobel Lecture). Angewandte Chemie International Edition in English, 27(1):89–112, 1987.
- [20] Jean-Marie Lehn. Dynamic combinatorial chemistry and virtual combinatorial libraries. In Essays in Contemporary Chemistry, pages 307–326. Wiley, 2001.
- [21] Prasenjit Mal, Boris Breiner, Kari Rissanen, and Jonathan R. Nitschke. White phosphorus is air-stable within a self-assembled tetrahedral capsule. Science, 324(5935):1697–1699, 2009.
- [22] Emmanuel A Meyer, Ronald K Castellano, and François Diederich. Interactions with aromatic rings in chemical and biological recognition. Angewandte Chemie International Edition, 42(11):1210–1250, 2003.
- [23] Charles J. Pedersen. The Discovery of Crown Ethers (Noble Lecture). Angewandte Chemie International Edition in English, 27(8):1021–1027, 1987.
- [24] John W Raymond and Peter Willett. Maximum common subgraph isomorphism algorithms for the matching of chemical structures. Journal of computer-aided molecular design, 16(7):521–533, 2002. doi:10.1023/A:1021271615909.
- [25] Schrödinger, LLC. The PyMOL molecular graphics system, version 2.5.0. Open-source build, 2021.
- [26] Yann Strozecki. Enumeration complexity: Incremental time, delay and space. Habilitation à diriger des recherches, 2021.
- [27] Bruno Therrien. Drug delivery by water-soluble organometallic cages. Chemistry of Nanocontainers, pages 35–55, 2012.
- [28] Nenad Trinajstic. Chemical graph theory. CRC press, 2018.
- [29] Robert A. Tromans, Tom S. Carter, Laurent Chabanne, Matthew P. Crump, Hongyu Li, Johnathan V. Matlock, Michael G. Orchard, and Anthony P. Davis. A biomimetic receptor for glucose. Nature Chemistry, 11(1):52–56, 2019.
- [30] Takeaki Uno. Constant time enumeration by amortization. In Workshop on Algorithms and Data Structures, pages 593–605. Springer, 2015. doi:10.1007/978-3-319-21840-3_49.
- [31] Xinchun Yang, Zakir Ullah, J. Fraser Stoddart, and Cafer T. Yavuz. Porous Organic Cages. Chemical Reviews, 123(8), 2023.
- [32] Gang Zhang and Michael Mastalerz. Organic cage compounds – from shape-persistency to function. Chemical Society Reviews, 43:1934–1947, 2014.
Appendix A Collision-Free Angular Intervals on a Circle problem
To construct the circle, we define a center point from the positions and . The next vertex () is constrained to lie on a circle of radius in a plane spanned by two orthonormal vectors and . Thus, any candidate position can be written as
For each point in the set of positions of atoms, the collision constraint defines one or more forbidden angular intervals. The valid placement domain for the next vertex is obtained by subtracting each forbidden interval to the full domain.
As previously discussed, once the next atom is positioned, the two missing hydrogen atoms can be added to complete the local bonding pattern. However, the position of the next vertex is not known in advance, but we know that an hydrogen atom satisfies the following constraint: the angle formed by the next vertex, the starting vertex , and the hydrogen atom is .
Following the same geometric construction used for the next atom, we determine a circle on which each hydrogen atom must lie. This circle has center , radius , and lies in a plane spanned by two orthonormal vectors and . Any feasible position of a hydrogen atom can therefore be parameterized by an angle as
Thus, we can also compute the valid angular intervals for the hydrogen atoms. Since , cannot be directly used. Instead, the angular shift must be determined using the formula for the angle between two vectors: .
From this relation, we compute a value giving the geometric correspondence between the angular parameters of the two circles and we write . Thus, we obtain the valid angular intervals for all three positions (the next atom and the two hydrogens), as shown in Figure 8, and we intersect them to find the final set of which are valid placement region.
The number of intervals created by possible collisions can easily be bounded as follows. All atoms in collision with an atom whose center lies on the circle are included in a region at distance at most of this circle. Hence, the region containing the potential interfering atoms is of volume and each atom is of volume . Hence, the number of such atoms is bounded by the quotient of these values, . Then, taking the intersection of the three set of intervals, multiplies this number by at most three. This is a crude bound and in practice, there are less than intervals.
Appendix B Interconnection tree
Proof of NP-Completeness.
The problem of deciding whether there is an interconnection tree can be generalized to arbitrary multipartite graphs.
Problem 1 (Interconnection Tree Problem (IT)).
Given a multipartite graph , decide whether there exists an interconnection tree of .
In contrast to complete multipartite graphs, the problem is -complete for general multipartite graphs.
Theorem 8.
The interconnection tree problem is -complete.
Proof.
Let us prove that IT NP. A witness is any set of edges of . We can check that is a matching and that is a spanning tree of in linear time in . Hence, IT NP.
We now prove that IT is NP-hard by a reduction from the Hamiltonian path problem [18]. We consider an instance of the Hamiltonian path problem.
We build a multipartite graph such that every vertex of is represented by a part containing two vertices. More formally, for each , we define and with
The graph is a multipartite graph, where the parts are the ’s. This graph can be built in time linear in .
() We prove that if has a Hamiltonian path then has an interconnection tree. Let be a Hamiltonian path of . We let . Since is a Hamiltonian path, no vertex is repeated in the path. Hence, by construction of , its edges have all distinct vertices, which means it is a matching of .
By construction, each edge of is of the form . Hence, is a path of the quotient graph. Because is Hamiltonian it contains all vertices of , therefore connect all vertices of the quotient graph . Therefore, is a Hamiltonian path of and thus is also a spanning tree.
() We prove that if has an interconnection tree then has a Hamiltonian path. Let be a matching of G and is a spanning tree of .
Since each part contains exactly two vertices, there are at most two edges of using a vertex of the same part. Thus, is a connected subgraph of and each of its vertices is of degree at most two. Therefore, is a path connecting all parts of .
For each consecutive pair in the path , the matching must contain an edge between a vertex of and a vertex of . By construction of , either all edges between and are of the form or . These edges exist in only if . Hence, for all , the pair is an edge in . Therefore, is a Hamiltonian path in .
Proof of Lemma 5.
Here, we give a proof of Lemma 5, which gives a simple characterization of the complete multipartite graphs admitting an interconnection tree.
Lemma (Lemma 5).
Let be a complete multipartite graph. Then admits an interconnection tree if and only if .
Proof.
Assume that is an interconnection tree of . It has edges, because is a spanning tree of , which has for vertices the parts of . Since is a matching, no vertex is used twice in , hence there are exactly vertices in . Thus, contains at least vertices.
We now prove by induction on that a complete multipartite graph with has an interconnection tree. If there is a single part, we always have then is an interconnection tree. If has exactly two parts and , we have and because the parts are non-empty and is an interconnection tree of .
Let us assume the property for parts and let us consider with parts such that . W.l.o.g., let the largest part of be . Because and , we have . We let and (which exists because the parts are non-empty). The multipartite graph has non-empty parts, and it satisfies . By the induction hypothesis used on , we have an interconnection tree of and by Lemma 6, is an interconnection tree of , which proves the induction.
Algorithm.
We give here a description of the algorithm described in Section 4. The variable represents the input, it encodes all vertices with their parts stored in order in an array. Function Parts(V) returns the number of parts of and the number of vertices.
We use a function Update, which modifies such that and are merged as well as their parts. We also use the function Remove, which removes vertex from . All these operations can be implemented in constant time.
Complexity of Algorithm 1 is given in the following theorem, along with the proof omitted from the main part of the text.
Theorem.
Let be a complete multipartite graph with parts. Then can be enumerated with worst-case delay and amortized delay .
Proof.
In Equation 1, all graphs are complete, hence we can use the algorithm derived from Lemma 5 to test whether they have an interconnection tree. Moreover, the depth is only , since the number of parts is decremented at each recursive call. Hence, we have a delay in where is the number of vertices (we do not need to represent the edges explicitly since the graph is complete multipartite).
In fact, the delay can be made better by maintaining the proper data structure and doing the union in order of increasing . First, the criterion of Lemma 5 can be maintained in constant time. We maintain the number of partitions which always decrease by one in a subgraph and the number of vertices which decreases by two because of the contraction and by one each time we consider a contraction from a new vertex.
Second, we can compute from the previous in constant time. To do so, we store the part as a list of the original parts which constitutes it and each original part is stored as a linked list of its element to allow for efficient removal of elements. We also maintain a table of the vertices pointing to the list each vertex belongs to and its position in this list.
When we compute from , we have to add back the vertex to its part and to remove the vertex in its part. Then, we remove the part of to the list representing and insert the part to which belongs. We also need to compute from , where is the element after in the order. The operations are the same as in the previous case, and we also remove from its part. The operations are in constant number and in constant time (insertion in a list), hence the complexity for a single recursive call is in . Therefore, the delay of the algorithm is in .
We now analyze the amortized delay of this algorithm. Assume that part is initially chosen to be the largest and will thus always be the largest in the recursive calls. We charge the computation time associated to computing the subgraph in the enumeration algorithm and to detect the first empty set to all its children divided equally. This cost is constant because of the previous analysis and denoted by .
When then , for any . Moreover, when , there are one or two parts and the interconnection tree is produced in one step and constant time . Therefore, if we consider an interconnection tree, it has received a cost of from its ancestor level upper in the tree. Hence, the cost by interconnection tree is bounded by .



