Enumeration method for tree-like chemical compounds with benzene rings and naphthalene rings by breadth-first search order

Background Drug discovery and design are important research fields in bioinformatics. Enumeration of chemical compounds is essential not only for the purpose, but also for analysis of chemical space and structure elucidation. In our previous study, we developed enumeration methods BfsSimEnum and BfsMulEnum for tree-like chemical compounds using a tree-structure to represent a chemical compound, which is limited to acyclic chemical compounds only. Results In this paper, we extend the methods, and develop BfsBenNaphEnum that can enumerate tree-like chemical compounds containing benzene rings and naphthalene rings, which include benzene isomers and naphthalene isomers such as ortho, meta, and para, by treating a benzene ring as an atom with valence six, instead of a ring of six carbon atoms, and treating a naphthalene ring as two benzene rings having a special bond. We compare our method with MOLGEN 5.0, which is a well-known general purpose structure generator, to enumerate chemical structures from a set of chemical formulas in terms of the number of enumerated structures and the computational time. The result suggests that our proposed method can reduce the computational time efficiently. Conclusions We propose the enumeration method BfsBenNaphEnum for tree-like chemical compounds containing benzene rings and naphthalene rings as cyclic structures. BfsBenNaphEnum was from 50 times to 5,000,000 times faster than MOLGEN 5.0 for instances with 8 to 14 carbon atoms in our experiments.


Background
Enumeration of chemical compounds is important in bioinformatics, and has been adapted to several applications such as drug discovery and design [1][2][3], structure elucidation [4][5][6], and analyses of chemical spaces [7][8][9][10][11][12][13]. It is defined as a problem of generating all non-redundant chemical structures satisfying some constraints. For example, a chemical formula, which consists of the number of each atom included in the compound, is given as an input. There are several algorithms for enumerating chemical compounds from a chemical formula and most of them use a molecular graph to represent a chemical compound, where the nodes and edges of the graph refer to atoms and bonds of the chemical compound, respectively. Some of those algorithms are claimed to be able to enumerate various chemical structures without restriction of the structure, such as MOLGEN [14] and Open Molecule Generator (OMG) [15]. It was reported that OMG is able to deal with different valences for a kind of atom, and was not efficient for several instances compared with MOLGEN. While the remaining ones, such as EnuMol [16,17] as well as BfsSimEnum and BfsMulEnum [18], have a limitation of the structure of enumerated compounds, such as acyclic compounds for BfsSimEnum and BfsMulEnum and compounds with no cycle except for benzene rings for EnuMol, the methods consume significantly less computational time. There are also related application softwares, e.g. SmiLib [19] and CLEVER [20], that generate chemical compounds from given fragments. The limitation of these tools is that they require a library of desired chemical fragments, which can be generated by the enumeration tool.
Our previous methods, BfsSimEnum and BfsMulEnum, use a tree structure, instead of a general graph, to represent a chemical compound and call it a molecular tree so they can generate only tree-like chemical compounds. In this work, we develop BfsBenNaphEnum, which aims to reduce the limitation of previous methods by extending them such that they can enumerate chemical compounds containing only benzene rings and naphthalene rings as cyclic structures, which are six carbon atoms cyclic structures and ten carbon atoms bicyclic structures, respectively. Pólya proposed a group-theoretic method for isomer counting of single cyclic structures such as a benzene ring, a naphthalene ring, and an anthracene ring using the cycle index, from which many studies followed [21]. However, structures enumerated by these methods are restricted to certain types. Indeed, Meringer wrote that up to now the only way to calculate the number of isomers belonging to an arbitrary molecular formula is to use structure generators [22]. Suzuki et al. considered the problem of enumerating structures having monocyclic graph structures, each of which has exactly one cycle [23]. An enumeration method for tree-like chemical compounds containing only benzene rings as cyclic structures has been implemented on EnuMol web server (http://sunflower.kuicr.kyoto-u.ac.jp/tools/enumol/). On the other hand, our method can enumerate compounds containing naphthalene rings in addition to benzene rings. Moreover, the proposed algorithm can calculate the number of benzene rings and naphthalene rings from chemical formula, while users have to specify the number of benzene rings in EnuMol.
Chemical structures considered in this study can be represented by a molecular tree, where a benzene ring is converted to a node with valence six and a naphthalene ring is considered as two benzene nodes having a special bond. We name that special bond as a merge bond. Since a merge bond merges two carbon atoms of two benzene rings together, it reduces the number of carbon atoms with free valence electron of two benzene rings by two so we represent a merge bond by a double-edge. Moreover, benzene nodes cannot have double bonds with other nodes because they bond with other non-benzene atoms by a single bond [24]. This means that a double-edge represents a double bond if it connects two non-benzene nodes, while it represents a merge bond if it connects two benzene nodes. Therefore, bonds in a benzene ring and a naphthalene ring are considered as the same bond and Kekulé representation is not included in this work. Besides, this work uses a two-dimensional molecular tree to represent a chemical structure so it cannot deal with stereoisomers. For tautomeric, this work considers two structures in a pair of tautomeric as non-redundant compounds and generates both of them.
BfsSimEnum and BfsMulEnum are modified to return a set of molecular trees as the output, given a chemical formula, the number of benzene rings, and the number of naphthalene rings. After that, an attribute called carbon position list is added into benzene nodes in a molecular tree to represent the way that benzene nodes bond with their adjacent nodes. This attribute is important because bonding with different carbon atoms in a benzene ring may result in different chemical structures. Finally, for each molecular tree from BfsSimEnum and BfsMulEnum, we generate a set of molecular trees whose nodes adjacent to benzene nodes are labeled with a carbon position such that all chemical structures are enumerated without redundancy based on normal form rule.
For evaluating our proposed method, we perform computational experiments for several instances, and compare the execution time by our method with that by MOLGEN. We show that our proposed method is efficient for enumerating chemical compounds containing benzene rings and naphthalene rings, and is from 50 times to 5,000,000 times faster than MOLGEN for several instances in our experiments.

Enumeration problem
Let be a finite set of labels of atoms, for example, = {C,N,O,H}, where 'C' , 'N' , 'O' , and 'H' denote carbon, nitrogen, oxygen, and hydrogen atoms, respectively. A molecular graph is defined as a multi-graph G(V, E), where V is a set of nodes and E is a set of multi-edges, also denoted by V (G) and E(G), respectively. Each node is labeled with an atom-label in , while each edge represents the bond between two atoms and the multiplicity of edge represents the bond type. The degree of each node is equal to the valence of its atom. Let deg(v) and l(v) be the degree and the label of node v, respectively. Let val(l i ) be the valence of the atom represented by label l i in . It should be noted that there exist different valences for a kind of atom, for example, carbon atoms of CO 2 and CO. For this case, it is sufficient to put two distinct labels C and C (2) in , and to define val(C) = 4 and val C (2) = 2. Let num (G,l i ) be the total number of nodes labeled with label l i in molecular graph G. Then, the enumeration problem is defined as follows. Problem 1. Given the numbers n l i of atoms for all labels l i ∈ , the number n b of benzene rings, and the number n n of naphthalene rings, enumerate all non-redundant molecular graphs G such that num(G,l i ) = n l i for all l i ∈ , deg(v) = val(l(v)) for all nodes v ∈ V (G), and G includes exactly n b benzene rings, n n naphthalene rings, and no other cyclic structures. It must be noted that n b and n n can be zero.
In the case that the input chemical formula contains five or less carbon atoms, BfsStructEnum can enumerate only tree-like chemical compounds by specifying the number of benzene rings and the number of naphthalene rings to be zero. Because we enumerate molecular trees such that degree of each node equals to valence of atom label of that node, charged molecules cannot be enumerated automatically. However, they can still be enumerated by specifying a charged atom as a new kind of atom type with appropriate valence value.
Since our enumeration methods deal with a chemical compound as a node-labeled rooted ordered tree for efficient enumeration, we contract cyclic structures appearing in a molecular graph to single nodes. Concretely, we contract a benzene ring to a node, called benzene node, labeled with a special label 'b' , and contract a naphthalene ring to two benzene nodes connected by a special bond, called merge bond, represented by a double edge (see Fig. 1). Since six carbon atoms contained in a benzene ring are contracted into a benzene node, we need to remember which carbon atom in the benzene ring connects to its adjacent node in a molecular graph. Hence, we add an attribute called carbon position list to each benzene node. Figure 1b shows examples of carbon position lists using numbers assigned to carbon atoms in benzene rings in Fig. 1a. We call such a node-labeled rooted ordered tree whose benzene nodes are attributed with carbon position lists a carbon position-assigned molecular tree. We enumerate carbon position-assigned molecular trees instead of molecular graphs.

Center-rooted and left-heavy
In our previous work, we defined the normal form for molecular trees without any cyclic structures using centerrooted and left-heavy to avoid its redundant generation.
In this work, we also utilize center-rooted and left-heavy for carbon position-assigned molecular trees, of which properties do not depend on carbon position lists.
A molecular tree T is called center-rooted if its root is the center node (see Fig. 2a) or one endpoint of the center edge of the longest path in T (see Fig. 2b). The center can be either a node or an edge depending on the length of the longest path.
In order to define a left-heavy tree, atom-labels must be ordered so that they can be compared with each other, for example, b> C > N > O >H for = {b,C,N,O,H}, where 'b' denotes a special atom representing a benzene ring. Let T(u) be the ordered subtree rooted at u in T. Let u and v be two nodes in a molecular tree T, (u 1 , u 2 , . . . , u h ) and (v 1 , v 2 , . . . , v k ) be lists of child nodes of u and v, respectively. It is defined that Let mul(e) and mul(u,v) be the multiplicity of edge e = (u,v). Let (e 1 , e 2 , . . . , e m ) and (e 1 , e 2 , . . . , e m ) be two lists of edges in T(u) and T(v) in breadth-first search (BFS) order (see Fig. 4 It should be noted that center-rooted and left-heavy are different from centroid-rooted and left-heavy defined by Fujiwara et al. [16], for example, the molecular tree in Fig. 1b is center-rooted and is not centroid-rooted because the number of nodes in the left subtree by removing the root, 4, is more than (total number of nodes −1)/2 = (7 − 1)/2 = 3. In addition, their left-heavy is defined Fig. 1 Example of a molecular graph including benzene rings and naphthalene rings. a A molecular graph including one benzene ring and one naphthalene ring. b A rooted tree contracted from the left graph. It is noted that hydrogen atoms are omitted

Carbon position list
Let s = (v 1 , v 2 , . . . , v n ) be a list of nodes, |s| and s[i] denote the size and the i-th element of s, respectively. Let T sub (v 1 , v 2 ) be the left-heavy tree rooted at v 1 that consists of the connected component including v 1 when the by traversing a center-rooted left-heavy molecular tree T with BFS order, which is also denoted by index(v) if T is clear.

Proposition 1.
For a node v that has the parent node v p and a child node v c in a center-rooted molecular tree T, We define an equality T 1 = C T 2 for two rooted carbonposition assigned trees T 1 and T 2 if T 1 = m T 2 , and v is a list of lists, called a carbon position list explained later, for a benzene node v in T. For convenience, we define another equality T 1 = C T 2 by removing the condition that C T 1 r 1 = C T 2 r 2 for the roots r 1 and r 2 of T 1 and T 2 , respectively, from the conditions of T 1 = C T 2 , if r 1 and r 2 are benzene nodes.
For a node v having the parent v p and a child v c , Hence, only carbon position lists of descendent benzene nodes are needed to determine whether or not   [1] , v . Figure 6 shows examples of carbon positionassigned molecular trees, where benzene node v 1 in each tree has adjacent nodes v 2

Proposition 2. For a benzene node v that has the parent node v p in a center-rooted molecular tree T, A
Proof. If v has no child, it is clear because the adjacent node of v is only v p . We assume that v has a child v c . From . . , 6} and two carbon positions are assigned for a merge bond because a naphthalene ring shares two carbon atoms between two benzene rings.
In the examples of Fig. 6,

Definition 2. An adjacent node list A T
For a benzene node v 2 that is connected by a merge bond with the parent node v 1 , we suppose that the carbon atoms having positions 1,2 in v 2 are connected with the carbon atoms having positions x + 1, x in v 1 , respectively, where x takes an integer between 1 and 6, and x = (x mod 6) + 1 (see Fig. 7a). Here, consider the case that v 1 has the parent node v p . If T is in normal form (Definition 6), position 1 is assigned to the carbon atom connected with v p (Proposition 5). Then, from Proposi- for any child node v c of v 1 except v 2 , and the naphthalene ring is not symmetric. Consider the case that v 1 does not have a parent node, that is, v 1 is the root. If T sub (v 1 , v 2 ) = C T sub (v 2 , v 1 ), the naphthalene ring can be symmetric only with respect to the axis denoted by the dashed red line in Fig. 7a. Then, it is not needed to consider the other symmetry for the naphthalene ring.
Consider the case that Proposition 4). Then, a carbon position list C T (v 1 ,v 2 ) of a naphthalene ring consisting of two benzene nodes v 1 , v 2 is a list of lists determined from C T v 1 and C T v 2 according to the following rule, where C T in ascending order.  Fig. 7b).

exist two integers i and j such that
This definition is applied to comparison of C T 1 (v 1 ,v 2 ) and C T 2 (v 1 ,v 2 ) for a naphthalene ring with v 1 and v 2 in the same way.
In the example of Fig. 6, T 1 and T 2 have the same tree structure, and C T 2 v 1 = ((1), (4), (2, 3)) < ( (3), (4), (1,2) Let Aut b and Aut n be the automorphism groups of a benzene ring and a naphthalene ring, respectively (see Fig. 9). Aut b is generated from rotation of π/3 radians and reflection. For in a benzene ring. Aut n is generated from rotation of π radians and reflection. We suppose that a list φ(C T v [i] ) of carbon positions for a map φ and i = 1, . . . , |C T v | is in ascending order by sorting elements of the list because all nodes (4), (1, 2)) and the reflection map φ b by the perpendicular bisector between carbon atoms of 1 and 2.

Normal form of a carbon position-assigned molecular tree
In order to prevent generating redundant molecular trees in enumeration, we define a normal form of a carbon position-assigned molecular tree.

Definition 5. Let P be a path in T consisting of n nodes
. P is called a symmetric path if the following conditions are satisfied.
for all i = 1, · · · , n 2 , where x is the largest integer less than or equal to x, , and v ∈ V 1 \V 2 means that v ∈ V 1 and v / ∈ V 2 . Proof. For a path (v 1 , · · · , v n ), v i+1 and v n−i must be the parent nodes of v i and v n−i+1 , respectively, for i = 1, · · · , n−1 2 if n is odd and for i = 1, · · · , n 2 − 1 if n is even due to the center rooted property. Therefore, if the length of path is odd, v n+1 2 is the parent node of both v n+1 2 −1 and v n+1 2 +1 , which means that the depth of v n+1 2 is less than that of any node in the path.
In the case that n is even, either v n 2 or v n 2 +1 has the least depth among all nodes in the path and another node is the child node of that node. Assume that between these two nodes the parent node is v a and the child node is v b . v a cannot have a parent node because the height of T sub (v p , v a ), where v p is the parent node of v a , cannot be equal to the height of T sub (v c , v b ) for any nodes v c that are adjacent to v b due the center-rooted condition, which means that T sub (v a , v b ) = m T sub (v b , v a ) cannot be hold and the first condition of symmetric path is violated. In other words, v a , which is either v n 2 or v n 2 +1 , is the root node of the tree if n is even.
We say that v 1 is left of v n for a symmetric path Figure 10 shows examples of symmetric paths, , and C T 6 v 4 = C T 6 v 6 . We define an inequality T 1 > C T 2 for carbon positionassigned molecular trees T 1 and T 2 if T 1 > m T 2 , or T 1 = m T 2 , and there exists an integer i such that v i is a benzene node, Fig. 7a

. A carbon positionassigned molecular tree T that contains a carbon position list C T v for each benzene node v is in normal form if the following conditions are satisfied.
1. T is center-rooted and left-heavy.

T(v) ≥ m T sub (r, v) if the center of the longest path in
T with the root r is the edge (r, v). 3. Positions in each sublist of C T v for each benzene node v are in ascending order. 4. C T v ≤ φ b C T v for all benzene nodes v that is not connected by a merge bond with the parent node and all φ b ∈ Aut b . 5. For benzene nodes v 1 , v 2 connected by a merge bond such that v 1 is the root of T, We call a tree in normal form a normal tree.
(v 1 ,v 2 ) for rotation φ rot of π radians, and T 4 violates the condition. It is noted that T 4 is rotated by π radians from T 4 . For condition 6, v 1 and v 2 are connected by a merge bond. Thus, T 4 is a normal tree, and T 4 is not a normal tree.

Proposition 4. For a normal tree T with a benzene node v 1 that is connected by a merge bond with its child node v 2 and satisfies T sub
Proof. We assume that there exists a node v l as a left sibling of v 2 , and v l is the leftmost child of v 1 . Since T Fig. 10 Examples of symmetric paths. The red lines denote symmetric paths. a T 5 , where (v 2 , v 1 , v 3 ) is a symmetric path, and T sub . It contradicts the assumption, and v 2 is the leftmost child of v 1 . Therefore, 2), and positions 1,2 are assigned to the merge bond, that is x = 1 in Fig. 7a.
For a map φ b ∈ Aut b other than the identity and reflection map φ ref for a benzene ring, (v 1 ,v 2 ) and the correspondence between C T v 1 and

Proposition 5. For a benzene node v of a normal tree T, C T v [1] [1] is always equal to 1.
Proof. If v is not connected by a merge bond with the parent node, from condition 4, C T v must be the least possible carbon position list. Hence,

Lemma 1. Given a molecular graph G without cyclic structures except benzene rings and naphthalene rings, G can be represented by a normal tree.
Proof. We can assign numbers to carbons in benzene rings and naphthalene rings of G such that the conditions of Definition 6 are satisfied.

Lemma 2.
Given two different molecular graphs G 1 and G 2 , they cannot be represented by the same normal tree.
Proof. We can unambiguously obtain a molecular graph from a normal tree by replacing all benzene nodes with benzene rings according to its carbon position lists. (v 1 , . . . , v n ), G is the molecular graph obtained from the tree T by removing T sub (v 1 , v 2 ) and T sub (v n , v n−1 ) except v 1 Fig. 11). Then,

and v n from T, where v 1 is left of v n . If there is a nonidentity map φ of the automorphism group of G satisfying
Let u j and w j be child nodes of v i+1 and v n−i , respectively. Then, v i = u j 2 and v n−i+1 = w j 1 , We assume that v i+1 and v n−i are not benzene nodes. For child nodes u j of v i+1 , T u j ≥ C T u j+1 because u j , v i+1 , u j+1 is a symmetric path. Also for child nodes w j of v n−i , T w j ≥ C T w j+1 . From the definition of φ, . It means T(u j )= C T u j+l . We assume (2) , then there is an integer j (j 1 ≤ j ≤ j 2 ) such that T(u j ) > C T(w j ), and it contradicts Eq. (2). Therefore,

Lemma 3.
Given two different normal trees T 1 and T 2 , T 1 does not represent the same molecular graph as T 2 .
Proof. We assume that T 1 represents the same molecular graph as T 2 . Let G 1 and G 2 be molecular graphs transformed from T 1 and T 2 , respectively, where each carbon in benzene rings and naphthalene rings is connected with adjacent atoms according to carbon position lists of T 1 and T 2 . From the assumption, there is an isomorphism ψ from G 1 to G 2 . It means that l(v 1 Consider the case that the automorphism group Aut(G 1 ) of G 1 has only elements φ such that φ(v 1 ) = v 2 for v 1 and v 2 belonging to distinct benzene rings. Let T(G) be the molecular tree without carbon position lists, obtained from G by contracting benzene rings and naphthalene rings to benzene nodes, and satisfying conditions 1, 2 of Definition 6. We suppose that maps ψ and φ in G 1 are naturally extended to T(G 1 ). Since T 1 is dif- If v 1 is not connected by a merge bond with the parent node, there is a non-identity map because T 1 and T 2 represent the same molecular graph. It contradicts condition 4 of Definition 6. Suppose that v 1 is connected by a merge bond with the parent node v p and C T 1 v p = C T 2 ψ(v p ) . If T sub v p , v 1 = C T sub v 1 , v p , then v p is the root, and there is a non-identity map φ n ∈ Aut n such that C T 1 (v p ,v 1 ) = φ n C T 2 (ψ(v p ),ψ(v 1 )) because T 1 and T 2 represent the same molecular graph. It contradicts condition 5a. Otherwise, T sub v p , v 1 = C T sub v 1 , v p . If v p is not the root, then T 1 does not represent the same molecular graph as T 2 because T sub v a , v p , where v a is the parent of v p , is different from other subtrees connected to the naphthalene ring. It contradicts the assumption. If v p is the root, because T 1 and T 2 represent the same molecular graph. It contradicts condition 5b.

Methods
We propose an algorithm BfsBenNaphEnum for enumerating chemical compounds containing benzene rings and naphthalene rings as cyclic structures. BfsBen-NaphEnum utilizes our previously developed algorithms BfsSimEnum, BfsMulEnum [18], and assigns carbon position lists.

Modification of BfsSimEnum and BfsMulEnum
Suppose that the numbers n l i of atoms with label l i for all l i ∈ , the numbers n b , n n of benzene rings and naphthalene rings are given. BfsBenNaphEnum introduces a special label 'b' representing a benzene node to with b > l i ∈ and val(b) = 6, and executes BfsSimEnum to generate all non-redundant molecular trees T such that num(T, l i ) = n l i for l i ∈ except l i = b, C and num(T, b) = n b + 2n n , num(T, C) = n C − 6n b − 10n n . At this time, all edges of enumerated trees are single because BfsSimEnum generates only simple trees. Then, we modify BfsMulEnum to assign n n merge bonds to edges between benzene nodes in each tree enumerated by BfsSimEnum in addition to adding 1 + l i ∈ ,l i =b num(T, l i )(val(l i ) − 2)/2 bonds to edges between usual nodes. It should be noted that multiple bonds cannot be assigned to edges connected to benzene nodes since a carbon atom in benzene rings and naphthalene rings is connected with another adjacent atom by a single bond.

Assignment of carbon positions for molecular trees
In this algorithm, we traverse along the tree T from the rightmost deepest benzene node to the root in reverse BFS order because an adjacent node list depends on carbon position lists of descendant nodes. For each benzene node v we found, we assign a carbon position list not to violate the conditions of normal form.
The pseudocode of assignment part in BfsBen-NaphEnum is given in Algorithms 1 and 2. We always assign carbon position 1 to the first node in A T v (line 20 in ASSIGN function) due to Proposition 5, which is the parent node of v if v is not the root (Proposition 2). If v is the root and |A T v [1] | ≥ 3, we assign carbon position lists in Table 1 (see also Fig. 12) to v immediately for the sake of efficiency. Carbon position lists in Table 1 satisfy condition 4 of the normal form, and all the cases are included in the table.
For other carbon positions from 2 to 6, we use ASSIGN_CHILD to assign such positions to the remaining adjacent nodes. For example, let T 1 in Fig. 6 be output without any carbon position list by BfsMulEnum. T 1 has a benzene node v 1 , and Since v 1 is the root and Table 1 is not used, and the other nodes v 5 , v 2 , v 3 are assigned by ASSIGN_CHILD. For v 5 , each carbon position from 2 to 6 is examined (line 26 in ASSIGN_CHILD). For v 2 , each position from 2 to 6 except the position assigned to v 5 is examined (line 27).
For each benzene node v, after assignment of a carbon position list to A T v , whether or not C T v violates conditions 4, 5 of the normal form is confirmed (lines 5, 11, 14 in ASSIGN_CHILD). After carbon position lists are assigned to all benzene nodes, condition 6 is confirmed (line 4 in ASSIGN).

Results
In this section, we show that our proposed method can enumerate chemical compounds with benzene rings and naphthalene rings correctly and efficiently. For the evaluation, although MOLGEN 3.5 is more suitable than MOLGEN 5.0 to enumerate tree-like compounds because MOLGEN 3.5 offered the possibility to define substructures like benzene or naphthalene as macro atoms but MOLGEN 3.5 cannot handle all the cases provided in Table 2, we compared proposed tool with MOLGEN 5.0. Thereby, we implemented it and installed another wellknown general purpose structure generator, MOLGEN 5.0, on a computer with 3.47 GHz intel Xeon CPU and 23.5 GiB memory, and compared their computational time. The implementation of BfsBenNaphEnum is available on our supplementary web site, http://sunflower. kuicr.kyoto-u.ac.jp/jira/bfsenum/.
Since MOLGEN can enumerate chemical compounds without restriction on the structure, we must specify a benzene ring and a naphthalene ring as a substructure so that the enumerated structures contain only benzene rings and naphthalene rings as cyclic structures. As can be seen from Table 2, where 'n' and 'b' denote a naphthalene ring and a benzene ring, respectively, BfsBenNaphEnum enumerated chemical compounds much faster than MOL-GEN while giving the same number of enumerated structures. BfsBenNaphEnum was from 50 times to 5,000,000 times faster than MOLGEN for instances with 8 to 14 carbon atoms. Table 2 also compares the number of discovered compounds in PubChem, which are not limited to tree-like chemical compounds, with the number of compounds enumerated by the proposed algorithm for several chemical formulas. When the number of carbon atoms is large (greater than 8 in this case), the number of discovered compounds is much less than the number of enumerated compounds. This implies that there are still a numerous number of unknown compounds to be discovered, which possibly include some essential compounds. In this study, we examined chemical formulas including up to two benzene rings and one naphthalene ring because MOLGEN was not able to output results in practical time for chemical formulas including more benzene rings and naphthalene rings.
We plotted the relation between the number of enumerated structures and the computational time for both methods in Fig. 14, where both x-axis and y-axis are in a log scale. It is seen from the figure that the execution time of BfsBenNaphEnum is much smaller than that of MOLGEN.

Discussion
Our algorithm is limited to tree-like chemical structures without any cyclic structures except benzene rings and naphthalene rings while MOLGEN does not have such limitation. Therefore, in the future, we would like to extend the algorithm such that it can enumerate more complex cyclic structures, such as polycyclic aromatic compounds and nucleotides. Besides, in order to make enumeration tools practical, we need to rank enumerated structures because a large number of structures are usually enumerated. For that purpose, it might be useful to  employ drug likeness filters such as Lipinski RO5, and QED score. Incorporation of such filters into our system is also important future work.

Conclusions
We proposed a way to represent a benzene ring in a molecular tree by regarding it as a new defined atom with valence six and introducing a new attribute named carbon position list to benzene nodes. Carbon position of an atom specifies which carbon in a benzene ring that the corresponding atom bonds with. We also proposed a new kind of bond called merge bond that merges two benzene rings together to form a naphthalene ring. With merge bond a molecular tree can represent a structure containing naphthalene rings without defining new kind of atom. Moreover, since a benzene ring and a naphthalene ring are symmetric structures, we defined a rule to assign carbon position lists such that no redundant structures due to the symmetry of a benzene ring and a naphthalene ring are enumerated. The algorithm of this work consists of two main steps. Given the number of benzene rings, the number of naphthalene rings as well as a chemical formula, BfsSimEnum and BfsMulEnum are applied such that they can enumerate molecular trees with benzene nodes. Next, the new extension BfsBenNaphEnum assigns carbon position lists to benzene nodes in normal molecular trees.
To show the performance of our algorithm, all nonredundant chemical structures were enumerated for several chemical formulas by BfsBenNaphEnum and MOLGEN 5.0, a well-known general purpose structure generator. It is shown that our algorithm is reliable since it generated the same number of structures as MOLGEN, while expended much less computational time. BfsBen-NaphEnum was from 50 times to 5,000,000 times faster than MOLGEN for instances with 8 to 14 carbon atoms in our experiments. This is mainly because the number of nodes decreases from six to one for each benzene ring and from ten to two for each naphthalene ring in a chemical structure and because we enumerate chemical structures in the form of trees instead of graphs.