Skip to main content
  • Research Article
  • Open access
  • Published:

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

Abstract

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–3], structure elucidation [4–6], and analyses of chemical spaces [7–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.

Preliminaries

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 d e g(v) and l(v) be the degree and the label of node v, respectively. Let v a l(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 CO2 and CO. For this case, it is sufficient to put two distinct labels C and C (2) in Σ, and to define v a l(C)=4 and v a l(C (2))=2. Let n u m(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}}\phantom {\dot {i}\!}\) for all l i ∈Σ, d e g(v)=v a l(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 1 b shows examples of carbon position lists using numbers assigned to carbon atoms in benzene rings in Fig. 1 a. 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.

Fig. 1
figure 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

Center-rooted and left-heavy

In our previous work, we defined the normal form for molecular trees without any cyclic structures using center-rooted 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. 2 a) or one endpoint of the center edge of the longest path in T (see Fig. 2 b). The center can be either a node or an edge depending on the length of the longest path.

Fig. 2
figure 2

Illustration of center-rooted molecular trees. a Center of the longest path is a node. b Center of the longest path is an edge. The thick lines indicate one of the longest paths and the center node/edge is shown in red

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 T(u)> s T(v) if l(u)>l(v) (Fig. 3 a) or there exists an integer i such that T(u j )= s T(v j ) for all j<i and (T(u i )> s T(v i ) (Fig. 3 b) or i=k+1≤h (Fig. 3 c)). If T(u)> s T(v) or T(v)> s T(u) does not hold, it is said that T(u)= s T(v).

Fig. 3
figure 3

Illustration of three molecular trees such that T(u)> s T(v) or T(u)> m T(v). a l(u)>l(v). b l(u)=l(v), T(u 1)> s T(v 1). c l(u)=l(v), T(u 1)= s T(v 1), h=2>1=k. d T(u)= s T(v), m u l(e 1)>m u l(e1′)

Let m u l(e) and m u l(u,v) be the multiplicity of edge e=(u,v). Let (e 1,e 2,…,e m ) and (\(e^{\prime }_{1},e'_{2},\ldots,e'_{m}\)) be two lists of edges in T(u) and T(v) in breadth-first search (BFS) order (see Fig. 4), respectively. T(u)> m T(v) if T(u)> s T(v), or if T(u)= s T(v) and there exists an integer i such that m u l(e j )=m u l(e j′) for all j<i, and m u l(e i )>m u l(e i′) (Fig. 3 d). If T(u)> m T(v) or T(v)> m T(u) does not hold, it is said that T(u)= m T(v).

Fig. 4
figure 4

Illustration of breadth-first search (BFS) order. Numbers indicate BFS order for this example

Let c h i l d(v)=(v 1,v 2,…) be a list of all child nodes of node v in BFS order. It is defined that a molecular tree T is left-heavy if T(v i )≥ m T(v i+1) holds for all nodes v in T and all i=1,…,|c h i l d(v)|−1.

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. 1 b 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 using depth-first search order, not our breadth-first search order.

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 edge (v 1,v 2) is deleted from T (see Fig. 5). T sub(v 1,v 2)= m T(v 1) if v 1 is a child of v 2 in T. Let i n d e x(v,T) be the order of v∈V(T) by traversing a center-rooted left-heavy molecular tree T with BFS order, which is also denoted by i n d e x(v) if T is clear.

Fig. 5
figure 5

Illustration of subtree T sub(v 1,v 2). a A molecular tree T and T sub(v 1,v 2), which is surrounded by a red rectangle. b T sub(v 2,v 1)

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, T sub(v p ,v)≠ m T sub(v c ,v).

Proof.

The height of T sub(v p ,v) is larger than that of T sub(v c ,v) because T is center-rooted. Hence, T sub(v p ,v) is always different from T sub(v c ,v).

We define an equality T 1= C T 2 for two rooted carbon-position assigned trees T 1 and T 2 if T 1= m T 2, and C v 1 T 1=C v 2 T 2 for all benzene nodes v 1∈V(T 1), where v 2∈V(T 2) satisfies i n d e x(v 1,T 1)=i n d e x(v 2,T 2), and \({C^{T}_{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}=_{\underline {C}}T_{2}\) by removing the condition that C r 1 T 1=C r 2 T 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 , T sub(v p ,v)≠ C T sub(v c ,v) if T sub(v p ,v)≠ m T sub(v c ,v). Hence, only carbon position lists of descendent benzene nodes are needed to determine whether or not \(T^{sub}(v_{c_{1}},v)=_{C} T^{sub}(v_{c_{2}},v)\) for child nodes \(v_{c_{1}}\) and \(v_{c_{2}}\) of v.

Definition 1.

An adjacent node list \({A^{T}_{v}}\) of a benzene node v in a carbon position-assigned molecular tree T is defined as a list of lists of nodes adjacent to v using carbon position lists of descendent benzene nodes such that

  • \(|{A^{T}_{v}}[\!i]|\leq |{A^{T}_{v}}[\!i+1]|\) for all i,

  • \(index \left ({A^{T}_{v}}[\!i][\!1]\right) < index \left ({A^{T}_{v}}[\!i+1][\!1]\right)\) if \(|{A^{T}_{v}}[\!i]|=|{A^{T}_{v}}[\!i+1]|\),

  • \(index \left ({A^{T}_{v}}[\!i][\!j]\right) < index \left ({A^{T}_{v}}[\!i][j+1]\right)\) for all i,j,

  • \({A^{T}_{v}}[\!i]=(v')\) if (v, v ′) is a merge bond for some i,

  • \(v'\in {A_{v}^{T}}[\!i]\) if (v, v ′) is not a merge bond, and \(T^{sub}(v',v)=_{C}T^{sub} \left ({A^{T}_{v}}[\!i][\!1],v\right)\).

Figure 6 shows examples of carbon position-assigned molecular trees, where benzene node v 1 in each tree has adjacent nodes v 2,v 3,v 4,v 5. Then, \(T_{1}^{sub}(v_{2},v_{1}) =_{C} T_{1}^{sub}(v_{3},v_{1}) \neq _{C} T_{1}^{sub}(v_{4},v_{1}) \neq _{C} T_{1}^{sub}(v_{5},v_{1})\) and i n d e x(v 4)<i n d e x(v 5), so we have \(A^{T_{1}}_{v_{1}}=((v_{4}),(v_{5}),(v_{2},v_{3}))\). Also for T 2, \(A^{T_{2}}_{v_{1}}=((v_{4}),(v_{5}),(v_{2},v_{3}))\). For T 3, \(A^{T_{3}}_{v_{1}}=((v_{2}),(v_{3}),(v_{4}),(v_{5}))\) because (v 2,v 1) is a merge bond. If (v 2,v 1) is not a merge bond and C v 2 T 3=C v 3 T 3, then \(A^{T_{3}}_{v_{1}}=((v_{4}),(v_{5}),(v_{2},v_{3}))\).

Fig. 6
figure 6

Examples of adjacent node lists and carbon position lists. a T 1. b T 2. c T 3. d Molecular graph of T 1. e Molecular graph of T 2. f Molecular graph of T 3. Red numbers represent carbon positions of node v 1

Proposition 2.

For a benzene node v that has the parent node v p in a center-rooted molecular tree T, \({A^{T}_{v}}[\!1]=(v_{p})\).

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 Proposition 1 and i n d e x(v p )<i n d e x(v c ), \({A^{T}_{v}}[\!1]=(v_{p})\) always holds.

A carbon position list \({C^{T}_{v}}\) of a benzene node v in T is a list of lists, where \({C^{T}_{v}}[\!i]\) is a list of carbon positions of the nodes in \({A^{T}_{v}}[\!i]\). It is sufficient to enumerate \({C^{T}_{v}}[\!i]\) in ascending order because each node in \({A^{T}_{v}}[\!i]\) has the same subtree. If \(\left ({A^{T}_{v}}[\!i][\!1],v\right)\) is a merge bond, \({C^{T}_{v}}[\!i]\) has two carbon positions instead of one as usual. It should be noted that \({C^{T}_{v}}[\!i] \subseteq \{1,\ldots,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, \(C^{T_{1}}_{v_{1}}=((3),(4),(1,2))\) for \(A^{T_{1}}_{v_{1}}=((v_{4}),(v_{5}),(v_{2},v_{3}))\), \(C^{T_{2}}_{v_{1}}=((1),(4),(2,3))\) for \(A^{T_{2}}_{v_{1}}=((v_{4}),(v_{5}),(v_{2},v_{3}))\), \(C^{T_{3}}_{v_{1}}=((1,2),(3),(5),(4))\) for \(A^{T_{3}}_{v_{1}}=((v_{2}),(v_{3}),(v_{4}),(v_{5}))\).

Definition 2.

An adjacent node list \(A^{T}_{(v_{1},v_{2})}\) for a naphthalene ring with two benzene nodes v 1, v 2, where (v 1,v 2) is a merge bond, is defined as a list of lists of nodes adjacent to v 1 or v 2 except v 1 and v 2 such that

  • \(|A^{T}_{(v_{1},v_{2})}[\!i]|\leq |A^{T}_{(v_{1},v_{2})}[\!i+1]|\) for all i,

  • \(index \left (A^{T}_{(v_{1},v_{2})}[\!i][\!1]\right) < index \left (A^{T}_{(v_{1},v_{2})}[\!i+1][\!1]\right)\) if \(|A^{T}_{(v_{1},v_{2})}[\!i]|=|A^{T}_{(v_{1},v_{2})}[\!i+1]|\),

  • \(index \left (A^{T}_{(v_{1},v_{2})}[\!i][\!j]\right) < index \left (A^{T}_{(v_{1},v_{2})}[\!i][j+1]\right)\) for all i,j,

  • \(v'\in A^{T}_{(v_{1},v_{2})}[\!i]\) if \(T^{sub}(v',bn(v'))=_{C}T^{sub} \left (A^{T}_{(v_{1},v_{2})}[\!i][\!1],bn(A^{T}_{(v_{1},v_{2})}[\!i][\!1])\right)\), where b n(v) is v 1 or v 2 that is adjacent to v.

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 \(\overline {x+1}, \overline {x}\) in v 1, respectively, where x takes an integer between 1 and 6, and \(\overline {x}=(x \mod 6)+1\) (see Fig. 7 a). 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 Proposition 1, T sub(v p ,v 1)≠ C T sub(v c ,v 2) for any child node v c of v 2, T sub(v p ,v 1)≠ C T sub(v c ,v 1) 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})\neq _{\underline {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. 7 a. Then, it is not needed to consider the other symmetry for the naphthalene ring.

Fig. 7
figure 7

Correspondence between carbon positions in a naphthalene ring. a Correspondence between carbon positions involved with a merge bond in two benzene rings. b Correspondence between carbon positions of a naphthalene ring and two benzene rings in the case of \(T^{sub}(v_{1},v_{2})=_{\underline {C}}T^{sub} (v_{2},v_{1})\). The upper benzene ring v 1 is the parent of the lower benzene ring v 2. \(\overline {x}\) denotes (x mod 6)+1. Blue, red, and green numbers are positions of \(C^{T}_{v_{1}}\), \(C^{T}_{v_{2}}\), and \(C^{T}_{(v_{1},v_{2})}\), respectively. The dashed red line denotes the symmetric axis of Ï• ref

Consider the case that \(T^{sub}(v_{1},v_{2})=_{\underline {C}}T^{sub}(v_{2},v_{1})\). We can prove that x=1 if T is in normal form (see 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}_{(v_{1}, v_{2})}[\!i]\) is a list of carbon positions of nodes in \(A^{T}_{(v_{1}, v_{2})}[\!i]\) in ascending order.

Definition 3.

Carbon positions in a naphthalene ring correspond to carbon positions in two benzene nodes v 1,v 2, where v 1 is the parent node of v 2, if \(T^{sub}(v_{1},v_{2})=_{\underline {C}}T^{sub}(v_{2},v_{1})\), as follows (see Fig. 7 b).

  • For the benzene ring of v 1, positions 1,2 are assigned to carbons of the merge bond in \(C^{T}_{v_{1}}\). Position i (i=3,…,6) in \(C^{T}_{v_{1}}\) corresponds to i−2 in \(C^{T}_{(v_{1}, v_{2})}\).

  • For the benzene ring of v 2, positions 1,2 are assigned to carbons of the merge bond in \(C^{T}_{v_{2}}\). Position i (i=3,…,6) in \(C^{T}_{v_{2}}\) corresponds to i+2 in \(C^{T}_{(v_{1}, v_{2})}\).

Figure 8 shows examples of carbon position lists for a naphthalene ring, where \(T^{\prime }_{4}\) is T 4 with \(C^{T'_{4}}_{v_{1}}=((1,2),(4),(3))\) and \(C^{T'_{4}}_{v_{2}}=((1,2),(4),(5))\), \(T^{\prime \prime }_{4}\) is T 4 with \(C^{T^{\prime \prime }_{4}}_{v_{1}}=((1,2),(4),(5))\) and \(C^{T^{\prime \prime }_{4}}_{v_{2}}=((1,2),(4),(3))\). Then, \(A^{T'_{4}}_{(v_{1},v_{2})}=A^{T^{\prime \prime }_{4}}_{(v_{1},v_{2})}=((v_{3},v_{5}),(v_{4},v_{6}))\), \(C^{T'_{4}}_{(v_{1},v_{2})}=((2,6),(1,7))\), and \(C^{T^{\prime \prime }_{4}}_{(v_{1},v_{2})}=((2,6),(3,5))\).

Fig. 8
figure 8

Example of carbon position lists for a naphthalene ring. a T 4. b Molecular graph of \(T^{\prime }_{4}\), which is T 4 with \(C^{T'_{4}}_{v_{1}}=((1,2),(4),(3))\), \(C^{T'_{4}}_{v_{2}}=((1,2), (4),(5))\). c Molecular graph of \(T^{\prime \prime }_{4}\), which is T 4 with \(C^{T^{\prime \prime }_{4}}_{v_{1}}=((1,2),(4),(5))\), and \(C^{T^{\prime \prime }_{4}}_{v_{2}}=((1,2),(4),(3))\)

Definition 4.

For carbon position lists \(C^{T_{1}}_{v}\), \(C^{T_{2}}_{v}\), where A v T 1=A v T 2, it is defined that C v T 1<C v T 2 if there exist two integers i and j such that

  • C v T 1[ i ′][j ′]=C v T 2[ i ′][j ′] for all i ′<i and all \(j' = 1,\ldots,|C^{T_{1}}_{v}[\!i']|\),

  • C v T 1[ i][j ′]=C v T 2[ i][j ′] for all j ′<j,

  • C v T 1[ i][ j]<C v T 2[ i][ j].

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 v 1 T 2=((1),(4),(2,3))<((3),(4),(1,2))=C v 1 T 1 because C v 1 T 2[ 1][ 1]=1<3=C v 1 T 1[ 1][ 1].

Let A u t b and A u t n be the automorphism groups of a benzene ring and a naphthalene ring, respectively (see Fig. 9). A u t b is generated from rotation of π/3 radians and reflection. For ϕ b ∈A u t b , v 1 is adjacent to v 2 in a benzene ring if and only if ϕ b (v 1) is adjacent to ϕ b (v 2) in a benzene ring. A u t n is generated from rotation of π radians and reflection. We suppose that a list \(\phi ({C^{T}_{v}}[\!i])\) of carbon positions for a map ϕ and \(i=1,\dots,|{C^{T}_{v}}|\) is in ascending order by sorting elements of the list because all nodes in \({A^{T}_{v}}[\!i]\) have the same subtree. For example, ϕ b (C v 1 T 1)=((6),(5),(1,2)) for \(C^{T_{1}}_{v_{1}}=((3),(4),(1,2))\) and the reflection map ϕ b by the perpendicular bisector between carbon atoms of 1 and 2.

Fig. 9
figure 9

Illustration of automorphism of a benzene ring and a naphthalene ring. a A benzene ring. b A naphthalene ring. Dashed lines indicate reflections, curves indicate rotations, where all automorphisms are not shown

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 (v 1,v 2,…,v n )(n≥2). P is called a symmetric path if the following conditions are satisfied.

  • \(T^{sub} \left (v_{\lfloor \frac {n}2\rfloor },v_{\lfloor \frac {n}2\rfloor +1}\right) {=\;}_{m}T^{sub} \left (v_{n-\lfloor \frac {n}2\rfloor +1},v_{n-\lfloor \frac {n}2\rfloor }\right)\),

  • \(index \left (v_{i},T^{sub}\left (v_{\lfloor \frac {n}2\rfloor },v_{\lfloor \frac {n}2\rfloor +1}\right)\right) = index \left (v_{n-i+1},T^{sub}\left (v_{n-\lfloor \frac {n}2\rfloor +1},v_{n-\lfloor \frac {n}2\rfloor }\right)\right)\) for all \(i=1,\cdots,\lfloor \frac {n}2\rfloor \), where ⌊x⌋ is the largest integer less than or equal to x,

  • \({C^{T}_{v}} = C^{T}_{v'}\) for all benzene nodes \(v\in V\left (T^{sub}\left (v_{\lfloor \frac {n}2\rfloor },v_{\lfloor \frac {n}2\rfloor +1}\right)\right) \backslash V \left (T^{sub}(v_{1},v_{2})\right)\), where \(v'\in V \left (T^{sub}\left (v_{n-\lfloor \frac {n}2\rfloor +1},v_{n-\lfloor \frac {n}2\rfloor }\right)\right)\) satisfies \(index \left (v', T^{sub}\left (v_{n-\lfloor \frac {n}2\rfloor +1},v_{n-\lfloor \frac {n}2\rfloor }\right)\right) = index \left (v, T^{sub}\left (v_{\lfloor \frac {n}2\rfloor },v_{\lfloor \frac {n}2\rfloor +1}\right)\right)\), and v∈V 1∖V 2 means that v∈V 1 and v∉V 2.

Proposition 3.

For a center-rooted molecular tree, either of \(v_{\frac {n}2}\) and \(v_{\frac {n}2+1}\) is the root if the length of a symmetric path (v 1,⋯,v n ) is even. Otherwise, the depth of \(v_{\frac {n+1}2}\) is less than that of any node in the path.

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,\cdots,\frac {n-1}{2}\) if n is odd and for \(i=1,\cdots,\frac {n}{2}-1\) if n is even due to the center rooted property. Therefore, if the length of path is odd, \(v_{\frac {n+1}{2}}\) is the parent node of both \(v_{\frac {n+1}{2}-1}\) and \(v_{\frac {n+1}{2}+1}\), which means that the depth of \(v_{\frac {n+1}2}\) is less than that of any node in the path.

In the case that n is even, either \(v_{\frac {n}{2}}\) or \(v_{\frac {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_{\frac {n}{2}}\) or \(v_{\frac {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 (v 1,…,v n ) when \(v_{n-\lfloor \frac {n}2\rfloor +1}\) is the root, or i n d e x(v 1)<i n d e x(v n ).

Figure 10 shows examples of symmetric paths, (v 2,v 1,v 3) in T 5 and (v 5,v 2,v 1,v 3) in T 6, where \(T_{5}^{sub}(v_{2},v_{1})=_{m}T_{5}^{sub}(v_{3},v_{1})\), \(T_{6}^{sub}(v_{2},v_{1})=_{m}T_{6}^{sub}(v_{1},v_{2})\), and C v 4 T 6=C v 6 T 6.

Fig. 10
figure 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_{5}^{sub}(v_{2},v_{1})=_{m}T_{5}^{sub} (v_{3},v_{1})\). b T 6, where (v 5,v 2,v 1,v 3) is a symmetric path, \(T_{6}^{sub}(v_{2},v_{1})=_{m}T_{6}^{sub}(v_{1},v_{2})\) and C v 4 T 6=C v 6 T 6

We define an inequality T 1> C T 2 for carbon position-assigned 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, C v i T 1>C v i′T 2, and C v j T 1=C v j′T 2 for all benzene nodes v j with j>i, where i n d e x(v k ,T 1)=i n d e x(v k′,T 2) for all k=1,…,|V(T 1)|.

Definition 6.

Let Ï• ref be the reflection map with the symmetric axis shown in Fig. 7 a. A carbon position-assigned 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. 1.

    T is center-rooted and left-heavy.

  2. 2.

    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. 3.

    Positions in each sublist of \({C^{T}_{v}}\) for each benzene node v are in ascending order.

  4. 4.

    \({C^{T}_{v}}\leq \phi _{b}\left ({C^{T}_{v}}\right)\) for all benzene nodes v that is not connected by a merge bond with the parent node and all ϕ b ∈A u t b .

  5. 5.

    For benzene nodes v 1,v 2 connected by a merge bond such that v 1 is the root of T,

    1. (a)

      \(C^{T}_{(v_{1},v_{2})}\leq \phi _{n} \left (C^{T}_{(v_{1},v_{2})}\right)\) for all ϕ n ∈A u t n if \(T^{sub}(v_{1},v_{2})=_{\underline {C}}T^{sub}(v_{2},v_{1})\), where \(C^{T}_{(v_{1},v_{2})}\) is related with \(C^{T}_{v_{1}}\) and \(C^{T}_{v_{2}}\) by Definition 3.

    2. (b)

      \(C^{T}_{v_{2}}\leq \phi_{ref} \left (C^{T}_{v_{2}}\right)\) if \(T^{sub}(v_{1},v_{2})\neq _{\underline {C}} T^{sub}(v_{2},v_{1})\) and \(C^{T}_{v_{1}}=\phi_{ref} \left (C^{T}_{v_{1}}\right)\).

  6. 6.

    T sub(v 1,v 2)≥ C T sub(v n ,v n−1) for all pairs v 1,v n of nodes such that the path (v 1,…,v n ) is a symmetric path, v 1 and v n (=v 2) are not connected by a merge bond, and v 1 is left of v n .

We call a tree in normal form a normal tree.

Figure 8 also shows molecular trees in normal form and not in normal form. For condition 4 of the definition, \(C^{T'_{4}}_{v_{1}}=((1,2),(4),(3))\leq \phi _{b} \left (C^{T'_{4}}_{v_{1}}\right)\), \(C^{T^{\prime \prime }_{4}}_{v_{1}} = ((1,2),(4),(5))\leq \phi _{b}\left (C^{T^{\prime \prime }_{4}}_{v_{1}}\right)\). \(T^{\prime }_{4}\) and \(T^{\prime \prime }_{4}\) satisfy conditions 1, 2, 3, and 4. For condition 5, \(C^{T'_{4}}_{(v_{1},v_{2})}=((2,6),(1,7))\leq \phi _{n}\left (C^{T'_{4}}_{(v_{1},v_{2})}\right)\), whereas \(C^{T^{\prime \prime }_{4}}_{(v_{1},v_{2})}=((2,6),(3,5))>((2,6),(1,7))=\phi_{rot}\left (C^{T^{\prime \prime }_{4}}_{(v_{1},v_{2})}\right)\) for rotation ϕ rot of π radians, and \(T^{\prime \prime }_{4}\) violates the condition. It is noted that \(T^{\prime \prime }_{4}\) is rotated by π radians from \(T^{\prime }_{4}\). For condition 6, v 1 and v 2 are connected by a merge bond. Thus, \(T^{\prime }_{4}\) is a normal tree, and \(T^{\prime \prime }_{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}(v_{1},v_{2})=_{\underline {C}}T^{sub}(v_{2},v_{1})\), positions 1,2 are assigned to the merge bond in the benzene ring of v 1. Furthermore, if \(C^{T}_{(v_{1},v_{2})}\leq \phi _{n}\left (C^{T}_{(v_{1},v_{2})}\right)\) for all ϕ n ∈A u t n , then \(C^{T}_{v_{1}}\leq \phi _{b}\left (C^{T}_{v_{1}}\right)\) for all ϕ b ∈A u t b .

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 is left-heavy, T(v l )≥ m T(v 2), and l(v l )=l(v 2)=‘b’ is needed. However, T(v l )= C T(v c ), where v c is the leftmost child of v 2, because \(T^{sub}(v_{1},v_{2})=_{\underline {C}}T^{sub}(v_{2},v_{1})=_{C}T(v_{2})\). Hence, T(v l )< m T(v 2). It contradicts the assumption, and v 2 is the leftmost child of v 1. Therefore, \(A^{T}_{v_{1}}[\!1]=(v_{2})\). From condition 4 of Definition 6, \(C^{T}_{v_{1}}[\!1]=(1,2)\), and positions 1,2 are assigned to the merge bond, that is x=1 in Fig. 7 a.

For a map ϕ b ∈A u t b other than the identity and reflection map ϕ ref for a benzene ring, \(C^{T}_{v_{1}}<\phi _{b}\left (C^{T}_{v_{1}}\right)\) because each of ϕ b (1) and ϕ b (2) is at least 2. From \(C^{T}_{(v_{1},v_{2})}\leq \phi_{ref}\left (C^{T}_{(v_{1},v_{2})}\right)\) and the correspondence between \(C^{T}_{v_{1}}\) and \(C^{T}_{(v_{1},v_{2})}\), \(C^{T}_{v_{1}}\leq \phi_{ref}\left (C^{T}_{v_{1}}\right)\). Therefore, \(C^{T}_{v_{1}}\leq \phi _{b}\left (C^{T}_{v_{1}}\right)\) for all ϕ b ∈A u t b .

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, \({C^{T}_{v}}[\!1][\!1]=1\). Otherwise, from Definition 3, \({C^{T}_{v}}[\!1][\!1]=1\).

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.

Proposition 6.

For a normal tree T with a path (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 and v n from T, where v 1 is left of v n . If there is a non-identity map ϕ of the automorphism group of G ′ satisfying ϕ(v i )=v n−i+1 for all i=1,…,n, then T sub(v 1,v 2)≥ C T sub(v n ,v n−1), where ϕ in G ′ is naturally extended to T.

Proof.

If \(T^{sub}\left (\!v_{\lfloor \frac {n}2\rfloor },v_{\lfloor \frac {n}2\rfloor +1}\right)\!\!>_{m} \!T^{sub}\left (v_{n-\lfloor \frac {n}2\rfloor +1},v_{n-\lfloor \frac {n}2\rfloor }\right)\), then T sub(v 1,v 2)> m T sub(v n ,v n−1), and T sub(v 1,v 2)> C T sub(v n ,v n−1). We assume \(T^{sub}\left (v_{\lfloor \frac {n}2\rfloor },v_{\lfloor \frac {n}2\rfloor +1}\right)=_{m} T^{sub}\left (v_{n-\lfloor \frac {n}2\rfloor +1},v_{n-\lfloor \frac {n}2\rfloor }\right)\). If the path (v 1,…,v n ) is a symmetric path, T sub(v 1,v 2)≥ C T sub(v n ,v n−1) from condition 6. We assume that (v i+1,…,v n−i ) is a symmetric path for some i, and i n d e x(v i ,T sub(v i+1,v i+2))>i n d e x(v n−i+1,T sub(v n−i ,v n−i−1)) (see Fig. 11). Then,

$$ \begin{aligned} T^{sub}(v_{i+1},v_{i+2})&=_{m}T^{sub}(v_{n-i},v_{n-i-1}),\\ T^{sub}(v_{i+1},v_{i+2})&\geq_{C}T^{sub}(v_{n-i},v_{n-i-1}). \end{aligned} $$
((1))
Fig. 11
figure 11

Illustration of an automorphism ϕ in the proof. The red path indicates (v 1,…,v n ), where ϕ(v i )=v n−i+1 for all i=1,…,n

Let u j and w j be child nodes of v i+1 and v n−i , respectively. Then, \(v_{i}=u_{j_{2}\phantom {\dot {i}\!}}\) and \(v_{n-i+1}=w_{j_{1}\phantom {\dot {i}\!}}\), where j 1=i n d e x(v n−i+1,T sub(v n−i ,v n−i−1)) and j 2=i n d e x(v i ,T sub(v i+1,v i+2)). If v i+1 and v n−i are benzene nodes, \(T(u_{j_{1}})=_{C}T(v_{i})\), \(T(v_{n-i+1})=_{C}T\left (w_{j_{2}}\right)\), and T(v i )= C T(v n−i+1) because \(C^{T}_{v_{i+1}}=C^{T}_{v_{n-i}}\) and ϕ(v i )=v n−i+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 ϕ, T(u j )= C T(ϕ(u j )) for all u j ≠v i . If i n d e x(ϕ(u j+l ))<i n d e x(ϕ(u j )) for u j ,u j+l ≠v i and l>0, T(u j )≥ C T(u j+l )= C T(ϕ(u j+l ))≥ C T(ϕ(u j ))= C T(u j ). It means T(u j )= C T(u j+l ). We assume that i n d e x(ϕ(u j ))<i n d e x(ϕ(u j+l )) for all u j ≠v i , that is, ϕ(u j )=w j+1 for all j=j 1,…,j 2−1. Then,

$$ \begin{aligned} T\left(u_{j}\right)\,&{=}_{C}T\left(w_{j+1}\right)\leq_{C} T\left(w_{j}\right),\\ \text{and}\ \ T(v_{i})\,&{\leq}_{C} T\left(u_{j_{2}-1}\right){=}_{C} T\left(w_{j_{2}}\right). \end{aligned} $$
((2))

If T sub(v i+1,v i+2)> C T sub(v n−i ,v n−i−1), 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, T sub(v i+1,v i+2)= C T sub(v n−i ,v n−i−1), and T(v i )= C T(v n−i+1). Also for the case that (v i+1,…,v n−i ) is a symmetric path for some i and i n d e x(v i ,T sub(v i+1,v i+2))<i n d e x(v n−i+1,T sub(v n−i ,v n−i−1)), then T(v i )= C T(v n−i+1). Thus, T sub(v 1,v 2)≥ C T sub(v n ,v n−1).

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)=l(ψ(v 1)) for all v 1∈V(G 1), (ψ(v 1),ψ(v 2))∈E(G 2) if and only if (v 1,v 2)∈E(G 1), and m u l(ψ(v 1),ψ(v 2))=m u l(v 1,v 2).

Consider the case that the automorphism group A u t(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 different from T 2, there is a benzene node v 1∈V(T 1) such that

$$\begin{array}{@{}rcl@{}} C^{T_{1}}_{v_{1}}\neq C^{T_{2}}_{\psi(v_{1})}. \end{array} $$
((3))

If v 1 is not connected by a merge bond with the parent node, there is a non-identity map ϕ b ∈A u t b such that \(C^{T_{1}}_{v_{1}}=\phi _{b}\left (C^{T_{2}}_{\psi (v_{1})}\right)\) 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 v p T 1=C ψ(v p )T 2. If \(T^{sub}\left (v_{p},v_{1}\right) {=}_{\underline {C}}T^{sub}\left (v_{1},v_{p}\right)\), then v p is the root, and there is a non-identity map ϕ n ∈A u t n such that \(C^{T_{1}}_{(v_{p},v_{1})}=\phi _{n}\left (C^{T_{2}}_{(\psi (v_{p}),\psi (v_{1}))}\right)\) because T 1 and T 2 represent the same molecular graph. It contradicts condition 5a. Otherwise, \(T^{sub}\left (v_{p},v_{1}\right)\neq _{\underline {C}}T^{sub}\left (v_{1},v_{p}\right)\). 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, \(C^{T_{1}}_{v_{p}}=\phi_{ref}\left (C^{T_{1}}_{v_{p}}\right)\) and \(C^{T_{1}}_{v_{1}}=\phi_{ref}\left (C^{T_{2}}_{\psi (v_{1})}\right)\) because T 1 and T 2 represent the same molecular graph. It contradicts condition 5b.

Consider the case that there is an element ϕ∈A u t(G 1) such that ϕ(v 1)=v 2 for v 1 and v 2 belonging to distinct benzene rings. Since T 1 is different from T 2, there is a benzene node v 1∈V(T 1) such that

$$\begin{array}{@{}rcl@{}} C^{T_{1}}_{v_{1}}\neq C^{T_{2}}_{\psi(v_{1})}. \end{array} $$
((4))

Here, we suppose that conditions 3, 4, 5 are satisfied for all benzene nodes in T 1 and T 2. Then, there is a path from v 1 to ϕ(v 1)=v n , (v 1,…,v n ), in T 1. Since T 1 and T 2 represent the same molecular graph,

$$\begin{array}{@{}rcl@{}} &&T_{1}^{sub}(v_{1},v_{2})=_{C}T_{2}^{sub}(\psi(v_{n}),\psi(v_{n-1}))~\text{and}\\ &&T_{1}^{sub}(v_{n},v_{n-1})=_{C}T_{2}^{sub}(\psi(v_{1}),\psi(v_{2})). \end{array} $$
((5))

Here, we can assume that v 1 is left of v n and ψ(v 1) is left of ψ(v n ) without loss of generality. Then, from Proposition 6, for paths of (v 1,…,v n ) and (ψ(v 1),…,ψ(v n )),

$$\begin{array}{@{}rcl@{}} &&T_{1}^{sub}(v_{1},v_{2})\geq_{C}T_{1}^{sub}(v_{n},v_{n-1})\,\,\text{and}\\ &&T_{2}^{sub}(\psi(v_{1}),\psi(v_{2}))\geq_{C}T_{2}^{sub}(\psi(v_{n}),\psi(v_{n-1})) \end{array} $$
((6))

because T 1 and T 2 are normal trees. There is no carbon position lists that satisfy Eqs. (4), (6) and (7).

Therefore, T 1 does not represent the same molecular graph as T 2.

Methods

We propose an algorithm BfsBenNaphEnum for enumerating chemical compounds containing benzene rings and naphthalene rings as cyclic structures. BfsBenNaphEnum 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 v a l(b)=6, and executes BfsSimEnum to generate all non-redundant molecular trees T such that \(num(T, l_{i})=n_{l_{i}\phantom {\dot {i}\!}}\) for l i ∈Σ except l i =b,C and n u m(T,b)=n b +2n n , n u m(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+\sum _{l_{i}\in \Sigma, l_{i}\neq 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 BfsBenNaphEnum 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]|\geq 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.

Fig. 12
figure 12

Illustration of benzene rings having each carbon position list in Table 1. a ((1,2,3)). b ((1,2,4)). c ((1,3,5)). d ((1,2,3),(4,5,6)). e ((1,2,4),(3,5,6)). f ((1,3,5),(2,4,6)). g ((1,2,3,4)). h ((1,2,3,5)). i ((1,2,4,5)). j ((1,2,3,4,5)). k ((1,2,3,4,5,6)). Solid and dashed lines correspond to \({A^{T}_{v}}[\!1]\) and \({A^{T}_{v}}[\!2]\), respectively

Table 1 Carbon position lists for \({A^{T}_{v}}\), where v is the root, and \(|{A^{T}_{v}}[\!1]|\geq 3\)

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 \(A^{T_{1}}_{v_{1}}=((v_{4}),(v_{5}),(v_{2},v_{3}))\). First, carbon position 1 is assigned to \(A^{T_{1}}_{v_{1}}[\!1][\!1]=v_{4}\), that is, \(C^{T_{1}}_{v_{1}}[\!1][\!1]=1\). Since v 1 is the root and \(|A^{T_{1}}_{v_{1}}[\!1]|=1<3\), 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 v 3, each position from 2 to 6 that is more than the position assigned to v 2 except the position assigned to v 5 is examined (line 27) because v 2 and v 3 have the same subtree and condition 3 must be satisfied. Thus, \(C^{T_{1}}_{v_{1}}\!=((1),(2),(3,4)),((1),(2),(3,5)),((1),(2),(3,6)),\dots, ((1), (3),(2,4)),((1),(3),(2,5)),((1),(3),(2,6)),\dots,((1), (6),(4,5))\) are examined, where ((1),(6),(2,3)),((1),(6),(2,4)),((1),(5),(2,3)) and so on are discarded in the next step.

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).

Since an input of this part, that is, an output of BfsMulEnum, satisfies conditions 1, 2 of the normal form, BfsBenNaphEnum always outputs normal trees. In ASSIGN_CHILD, a distinct carbon position list is always assigned, and all patterns are assigned (line 28). Hence, BfsBenNaphEnum outputs all distinct normal trees.

Theorem 1.

BfsBenNaphEnum outputs all non-redundant molecular graphs that are solutions of Problem 1.

Figure 13 shows another example T 7 of molecular trees. T 7 includes four benzene nodes v 5, v 4, v 3, v 2 in reverse BFS order, and edges (v 2,v 4), (v 3,v 5) are merge bonds. First, our algorithm assigns carbon position lists for \(A^{T_{7}}_{v_{5}}=((v_{3}),(v_{7}))\) as \(C^{T_{7}}_{v_{5}}=((1,2),(3)),((1,2),(4)), ((1,2),(5)),((1,2),(6))\). In a similar way, for \(A^{T_{7}}_{v_{4}}= ((v_{2}),(v_{6}))\), \(C^{T_{7}}_{v_{4}}=((1,2),(3)),((1,2),(4)),((1,2),(5)), ((1,2),(6))\). For \(A^{T_{7}}_{v_{3}}=((v_{1}),(v_{5}))\), \(C^{T_{7}}_{v_{3}}=((1), (2,3)),((1),(3,4)),((1),(4,5)), ((1),(5,6))\) are examined. In line 5 of ASSIGN_CHILD, ((1),(4,5)) and ((1),(5,6)) are discarded because ϕ b (((1),(4,5)))=((1),(3,4)), ϕ b (((1),(5,6)))=((1),(2,3)) for the reflection map ϕ b with respect to the axis through positions 1 and 4, and these violate condition 4. In a similar way, for \(A^{T_{7}}_{v_{2}}=((v_{1}),(v_{4}))\), \(C^{T_{7}}_{v_{2}}=((1),(2,3)),((1),(3,4))\) are assigned. After carbon position lists are assigned to all benzene nodes, condition 6 is confirmed in line 4 of ASSIGN. If C v 2 T 7≠C v 3 T 7, then there is one symmetric path, \({\mathcal {P}}=\{(v_{2},v_{3})\}\), and T 7(v 2)≥ C T 7(v 3) must be satisfied. It means that C v 4 T 7=C v 5 T 7=((1,2),(3)),((1,2),(4)),((1,2),(5)),((1,2),(6)) and C v 2 T 7=((1),(3,4))>C v 3 T 7=((1),(2,3)), or C v 4 T 7>C v 5 T 7 and C v 2 T 7≠C v 3 T 7. Hence, there are \(4+ {4 \choose 2}\cdot 2=16\) structures. If C v 2 T 7=C v 3 T 7=((1),(2,3)) (or C v 2 T 7=C v 3 T 7=((1),(3,4))), then \({\mathcal {P}}=\{(v_{2},v_{3}),(v_{4},v_{5})\}\), and both of T 7(v 2)≥ C T 7(v 3) and T 7(v 4)≥ C T 7(v 5), that is, C v 4 T 7≥C v 5 T 7, must be satisfied. Hence, there are 4+3+2+1=10 structures. In total, 16+10·2=36 structures are generated by BfsBenNaphEnum for T 7.

Fig. 13
figure 13

Example of a molecular tree T 7

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 well-known 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/.

Table 2 Results on execution time (sec), the number of enumerated structures by BfsBenNaphEnum and MOLGEN, and the number of chemical compounds exist in PubChem database for several instances

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 MOLGEN 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.

Fig. 14
figure 14

Relation between the number of enumerated structures and the computational time (sec)

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 non-redundant 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. BfsBenNaphEnum 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.

References

  1. Ward RA, Kettle JG. Systematic enumeration of heteroaromatic ring systems as reagents for use in medicinal chemistry. J Med Chem. 2011; 54(13):4670–7.

    Article  CAS  PubMed  Google Scholar 

  2. Blum LC, Reymond JL. 970 million druglike small molecules for virtual screening in the chemical universe database gdb-13. J Am Chem Soc. 2009; 131(25):8732–3.

    Article  CAS  PubMed  Google Scholar 

  3. Mishima K, Kaneko H, Funatsu K. Development of a new de novo design algorithm for exploring chemical space. Mol Inform. 2014; 33(11-12):779–89.

    CAS  Google Scholar 

  4. Funatsu K, Sasaki S. Recent advances in the automated structure elucidation system, chemics. utilization of two-dimensional NMR spectral information and development of peripheral functions for examination of candidates. J Chem Inform Comput Sci. 1996; 36(2):190–204.

    Article  CAS  Google Scholar 

  5. Meringer M, Schymanski EL. Small molecule identification with MOLGEN and mass spectrometry. Metabolites. 2013; 3:440–62.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Koichi S, Arisaka M, Koshino H, Aoki A, Iwata S, Uno T, Satoh H. Chemical structure elucidation from 13C NMR chemical shifts: Efficient data processing using bipartite matching and maximal clique algorithms. J Chem Inform Model. 2014; 54:1027–35.

    Article  CAS  Google Scholar 

  7. Bytautas L, Klein DJ, Schmalz TG. All acyclic hydrocarbons: Formula periodic table and property overlap plots via chemical combinatorics. New J Chem. 2000; 24(5):329–36.

    Article  CAS  Google Scholar 

  8. Faulon J, Visco DP, Roe D. Enumerating molecules. Rev Comput Chem. 2005; 21:209.

    CAS  Google Scholar 

  9. Koch MA, Schuffenhauer A, Scheck M, Wetzel S, Casaulta M, Odermatt A, Ertl P, Waldmann H. Charting biologically relevant chemical space: A structural classification of natural products (sconp). Proc Natl Acad Sci U S A. 2005; 102(48):17272–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Mauser H, Stahl M. Chemical fragment spaces for de novo design. J Chem Inf Model. 2007; 47(2):318–24.

    Article  CAS  PubMed  Google Scholar 

  11. Andricopulo AD, Guido RV, Oliva G. Virtual screening and its integration with modern drug design technologies. Curr Med Chem. 2008; 15(1):37–46.

    Article  PubMed  Google Scholar 

  12. Reymond JL, van Deursen R, Blum LC, Ruddigkeit L. Chemical space as a source for new drugs. MedChemComm. 2010; 1(1):30–8.

    Article  CAS  Google Scholar 

  13. Bürgi JJ, Awale M, Boss SD, Schaer T, Marger F, Viveros-Paredes JM, Bertrand S, Gertsch J, Bertrand D, Reymond JL. Discovery of potent positive allosteric modulators of the α3β2 nicotinic acetylcholine receptor by a chemical space walk in chembl. ACS Chem Neurosci. 2014; 5(5):346–59.

    Article  PubMed  PubMed Central  Google Scholar 

  14. Gugisch R, Kerber A, Kohnert A, Laue R, Meringer M, Rücker C, Wassermann A. MOLGEN 5.0, a molecular structure generator. Sharjah, United Arab Emirates: Bentham Science Publishers Ltd.; 2012.

    Google Scholar 

  15. Peironcely JE, Rojas-Chertó M, Fichera D, Reijmers T, Coulier L, Faulon JL, Hankemeier T. OMG: Open Molecule Generator. J Cheminformatics. 2012; 4:21.

    Article  CAS  Google Scholar 

  16. Fujiwara H, Wang J, Zhao L, Nagamochi H, Akutsu T. Enumerating treelike chemical graphs with given path frequency. J Chem Inf Model. 2008; 48(7):1345–57.

    Article  CAS  PubMed  Google Scholar 

  17. Shimizu M, Nagamochi H, Akutsu T. Enumerating tree-like chemical graphs with given upper and lower bounds on path frequencies. BMC Bioinformatics. 2011; 12:14–3.

    Article  Google Scholar 

  18. Zhao Y, Hayashida M, Jindalertudomdee J, Akutsu T. Breadth-first search approach to enumeration of tree-like chemical compounds. J Bioinformatics Comput Biol. 2013; 11:1343007.

    Article  Google Scholar 

  19. Schüller A, Hähnke V, Schneider G. SmiLib v2.0: A Java-based tool for rapid combinatorial library enumeration. QSAR Comb Sci. 2007; 26(3):407–10.

    Article  Google Scholar 

  20. Song CM, Bernardo PH, Chai CLL, Tong JC. CLEVER: Pipeline for designing in silico chemical libraries. J Mol Graph Model. 2009; 27(5):578–83.

    Article  CAS  PubMed  Google Scholar 

  21. Trinajstić N. Chemical Graph Theory, 2nd edn. Boca Raton, Florida: CRC Press; 1992, pp. 275–391. Chap. 11 Isomer Enumeration.

    Google Scholar 

  22. Meringer M. Handbook of Chemoinformatics Algorithms. Boca Raton, Florida: CRC Press; 2010, pp. 233–67. Chap. 8 Structure Enumeration and Sampling.

    Google Scholar 

  23. Suzuki M, Nagamochi H, Akutsu T. Efficient enumeration of monocyclic chemical graphs with given path frequencies. J Cheminformatics. 2014; 6:31.

    Article  Google Scholar 

  24. Hardinger SA, University of California LADoC. Biochemistry: Chemistry 14D: Organic Reactions and Pharmaceuticals : Course Thinkbook, Lecture Supplements, Concept Focus Questions, OWLS Problems, Practice Problems. Plymouth, MI 48170: Hayden-McNeil Pub; 2008.

    Google Scholar 

Download references

Acknowledgements

This work was partially supported by Grants-in-Aid #26240034, #24500361, and #25-2920 from MEXT, Japan.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Morihiro Hayashida.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

JJ and MH developed, implemented the methods, and drafted the manuscript. YZ and TA participated in the discussions during the development of the methods and helped draft the manuscript. All authors read and approved the final manuscript.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Jindalertudomdee, J., Hayashida, M., Zhao, Y. et al. Enumeration method for tree-like chemical compounds with benzene rings and naphthalene rings by breadth-first search order. BMC Bioinformatics 17, 113 (2016). https://0-doi-org.brum.beds.ac.uk/10.1186/s12859-016-0962-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12859-016-0962-4

Keywords