Improved Algorithms for
Enumerating Tree-like Chemical
Graphs with Given Path Frequency
Yusuke Ishida* Liang Zhao Hiroshi Nagamochi
Graduate School of Informatics, Kyoto University, Japan

Tatsuya Akutsu
Institute for Chemical Research, Kyoto University, Japan

Outline1.

Problem Formulation

2.

Canonical Representation and Family tree

3.

Algorithm (Branch and Bound)–Detachment-cut

4.

Experimental Results for the first formulation

5.

H-less Single-bond Formulation–Hydrogen-cut

6.

Experimental Results for the second formulation

7.

Conclusions

1Background

inferencechemical compounds
applications:• structure determination
using mass-spectrum• drug design

given partial structures

Definition
A feature vector fK(G) :#occurrences of eachvertex-labeled pathof length 0,1,..., KHHOCCOH

Ha chemical compound G
f3(G) (K = 3)
length = 0
HOC

4
2
2

length = 1

length = 2

length = 3

OH
CHCOCC

COH 1
HCH 3
CCH 3
OCO 1
CCO 2

OCOH 1
CCOH 1OCCH 613212

Enumerating Chemical Multitree Problem
OutputInputlabel set
Σ = { H, O, C }valenceval (H) = 1val (O) = 2val (C) = 4

level of g
K=1featurevector g
HOCOHCHCOCC

4221321

HHOCCOHHHCCHOO

HHall Σ-labeled multitrees satisfying
the feature vector constraint andthe valence constraint

Previous Work• Aringhieri et al. [4OR, 2003]– designed two algorithms to generate all alkane
isomers.

• Fujiwara et al. [J. Chem. Inf. Model., 2008]– proposed a branch and bound algorithm forchemical multitree problem.
– gave H-less single-bond formulation of chemical
multitree problem and their algorithm can be also
applied to it.

3Canonical Representation for multitrees
T1,T2,T3 : isomorphic trees

canonical
depth representation
C012OC

O

H

H

H

C

H

C
O

H

H

CO

H

H

T2

T1

O

H

C

H

O

H

H

T3

Define the canonical representation as
the embedding with largest depth-label sequence.

Family treeφ• Add a vertex to T = φ.
• Repeatedly attach a newvertex to T based on
canonical representation.

family tree（#labels=2, #vertices=5 ）4Parent of multitrees
canonical rightmost canonicalrepresentation path representation
CCOCO

H

HT

canonical
representation

HHO

CO

HH

branching

P(T)

CHO

branchingCOH

HP(P(T))

Define the parent P(T) of T as the tree obtainedby removing the rightmost leaf from T.

If T is canonical, P(T) remains canonical.
Branching operation : O(1) time / node

Algorithm (Branch and Bound)φ

Add a vertex toT = φ.
Attach a new vertex to T basedon canonical representation.
(Branching Operation)
Check each constraint anddecide whether to discard T.
(Bounding Operation)
If |T| = (#input vertices),
output T as a solution.

family tree（#labels=2, #vertices=5）5

Bounding Operations
For a current tree T,

the input

CC

test whether fK(T) ≦ g holds.
the feature vector of T

Input

C

①feature-vector-cut

OCH 1

O H H OT②bond-cut

test whether deg(v;T) ≦ val(ℓ (v)) for ∀v ∈ T holds.
#edges incident to v in T

the valence of the label of v

③detachment-cut

CCO H H

test whether T can be extended to multitrees
satisfying degree constraint.

TDetachment-cut
Test whether T can be extended to multitrees satisfying degree constraint.

input

residual

H 13O2N1C6NH 2OH 1CH 10NC 1CO 2CC 5

HON
CA
NH
CH
COCC
NACA

shrinkrightmostpath
CCCCO HCN

H H H H H H H Ha current tree T

41131131112O
C

A

H

N

G=(V, E)‫‏‬r(H) = 4r(O) = 1r(N) = 1r(C) = 3r(A) = 1

ρ(H) = (1,1,1,1)‫‏‬
ρ(O) = (2)‫‏‬ρ(N) = (2)‫‏‬
ρ(C) = (4,3,2)‫‏‬
ρ(A) = (3)‫‏‬ρ:4-3+1=2A6Detachment : A Reverse Operation of Contraction
The number r(v) of copies of v
and the degree ρ(v i)of each copy v i are specified.
r(b)=3
ρ(b)=(3,3,1)‫‏‬

baWa

3Wb

3

132d

3r(d)=2
ρ(d)=(2,3)‫‏‬r(a)=1
ρ(a)=(3)‫‏‬1r(c)=4
ρ(c)=(1,3,2,3)‫ ‏‬WccWd

332a ρ-detachment H of Ga multigraph G

Detachment-cutv

Test whether G has a connected
and loopless ρ-detachment.
4+3+2 = 98for v = C① ∑ρ(vi) ≧ deg(v;G) ∀v∈V1≦i≦r(v)#edges incident to v in G

condition for degree constraint75XH
ρ(H) = (1,1,1,1)‫ ‏‬C
ρ(O) = (2)‫‏‬ρ(N) = (2)‫‏‬N
ρ(C) = (4,3,2)‫ ‏‬A
ρ(A) = (3)‫‏‬
G=(V, E)‫‏‬

2② d(X,V;G) ≧ ∑r(v) + c(G-X) - 1 φ≠∀X ⊆Vv∈X#edges joining X and V #connected components in G - X

condition for connectivityOr(H) = 4r(O) = 1r(N) = 1r(C) = 3r(A) = 17

Experiment• Comparing the running time of our algorithmwith Fujiwara et al.’s [2008]• Instances from KEGG LIGAND database(Replacing each benzene ring by a virtual atom of valence 6)• Pentium4 3.00GHz• T.O. : time over 1800 (sec)Experimental Results for the first formulation
Formula#atomsFujiwara et al.’s algorithm
CPU time#solutions

N.F. : not found

our algorithm
CPU time#solutions

K(sec)

C16H22O4 1
37 2
3T.O.
3.11
3.25

N.F.
92

158.230.48
0.30

570,773
92

C17H28N2O 1
43 2
3T.O.
50.55
16.78

N.F.
55
1109.271.40
0.6173,711
55
1C21H28N2O5 146 2
3T.O.
51.72
4.26

N.F.
162500.78
3.51
0.32

70,170162

C24H38O4 1
61 2
3

T.O.
T.O.
T.O.

N.F.
N.F.
N.F.

T.O.
318.68188.13N.F.
1,198
8(sec)

8

H-less Single-bond Formulation [Fujiwara et al. 2008]

H

C

O

C

single-bondtransformationH – removaltransformation

H

C

H
O

O

H
O
C
OH
CH
CO
CC

H

422
1
321

new atomCC

O{(O,C),2}

CO

O
CRCO
RO
CC
RC

2211111O

R = {(O,C),2}

Hydrogen-cutinputOCRCOROCCRC

231
1
12
1

R = {(O,C),2}h*(ℓ ) : sum of #hydrogens adjacent
to each atom of label ℓ.
H0OHCH

C{(O,C),2}

CO

320H
Hh*(C) = 3 + 2 + 0 = 5h*(O) = 0 + 1 = 11H

9

Hydrogen-cuth(ℓ, T) : sum of #hydrogens that must be adjacent
to each atom of label ℓ in a current tree T
If h(ℓ ;T) > h*(ℓ ) holds for a label ℓ, discard T.
Assume h*(C) = 7.
2
1: a new vertex prightmostpath

CNCCCCC

21NC

2 CC

3 C

{(O,C),2}{(O,C),2}OO

The degrees ofthese vertices
won’t be changed.

Th(C ; T) = 3ppp

Tc discard
3 >78h(C ; Tc) = 6

Experimental Results for the second formulation
Formula#atomsFujiwara et al.’s algorithm
CPU time#solutions

our algorithm
CPU time#solutions

K(sec)

C21H28N2O5 119 23

222.290.11
0.09

70,170
162

9.03
0.02
0.01

70,170162

C24H38O4 125 2
3T.O.
23.36
15.87

5,305,243
1,198
8T.O.
8.10
5.66

60,257,365
1,198
8

C19H39O7P 2
29 3
4T.O.
184.5411.86

161111543.37
45.36
3.60

2,520
11

C21H39O7P 2
33 3
4

T.O.
T.O.
118.48

77
11
11T.O.
438.1925.65

1,7361311(sec)

10

Conclusions• We proposed a branch and bound algorithm with twonew bounding operations.
• For the first formulation, we can solve the problem
with about 25 non-hydrogen atoms for K ≧ 2.
• For the second formulation, we can solve the problem
with about 30 non-hydrogen atoms for K ≧ 2.

Future Work• Treat more general graphs (e.g., outerplanar graphs).
• Use other graph structures for representing partial
structures of input.

11

Definition
A feature vector fK(G) :#occurrences of each vertex-labeled pathof length 0,1,..., Kf1(G)HOCHHOHCH
HOOO
COHCOCCC

4220
13102322

H

H

O

CCO

H

Ha chemical compound G

12

Canonical Representation for multitrees
T1,T2,T3 : isomorphic trees
DL: depth label sequence

left-heavyＣＣＣＯ Ｈ Ｈ

Ｏ Ｈ Ｈ

ＣＣＨ Ｏ Ｈ

Ｏ

Ｈ Ｃ Ｈ

Ｏ Ｈ

ＨＨＯ Ｈ

T2

T1

DL(T1) = (0,C,1,C,2,O,2,H,1,O,2,H,1,H,1,H)‫‏‬

T3

DL(T2) = (0,C,1,C,2,O, DL(T3) = (0,C,1,O,2,H,2,H,1,H,1,O,2,H,1,H)‫ ‏‬1,H,1,C,2,O,2,H,1,H)‫‏‬The embedding with largest DL is left-heavy.
[Nakano, Uno (2003)]

We enumerate left-heavy trees.

Parent of multitrees
left-heavy

left-heavy

left-heavy

ＣＣＣＣＯ Ｈ Ｈ

Ｏ Ｈ Ｈ

parentTＣＯ Ｈ

Ｏ Ｈ Ｈ

P(T)

ＣＯＯ Ｈ Ｈ

P(P(T))

P(T) : tree obtained by removing the rightmost leaf from T
If T is left-heavy, P(T) remains left-heavy.
In branching operation,
we attach a new vertex to rightmost path.

13

Feature Vector (K=1) and ρ-detachmentΣ = { H, O, C }val (H) = 1
val (O) = 2
val (C) = 4HOCHHOHCHHOOOCOHC
OCCCr(H) = 4r(O) = 2r(C) = 242
20
13
102
322

ρ(H) = (1,1,1,1)‫‏‬
ρ(O) = (2,2)‫‏‬
ρ(C) = (4,4)‫‏‬

HOCCOHHHOHCGDetachment-cut
Test whether T will generate trees satisfying given constraint.
inputresidual
H 13O2N1C6NH 2OH 1CH 10
NC 1CO 2CC 5

Which vertices arecounted as “residual”?
New vertices canbe attached.
CCCCO H C

H H H H H

N
H

H
O
N
C

4+1+0+1+NHCH
COCC

1311NOHC

G=(V, E)‫‏‬r(H)= 4+r(O)= 1+r(N)= 0+r(C)= 1+H H

a current tree T14

Detachment-cut
Test whether T will generate trees satisfying given constraint.
inputresidual
H 13O2N1C6NH 2OH 1CH 10
NC 1
CO 2CC 5

Which vertices arecounted as “residual”?
New vertices canbe attached.
CCCC

O H C

H H H H H

CN
H

4
H 4+1
O 1+1
N 0+3
C 1+A
1
NH 1
CH 3
CO 1
CC 1
NA 1
CA 2NAHa current tree TC

G=(V, E)‫‏‬r(H)= 4+4r(O)= 1+1r(N)= 0+1r(C)= 1+3r(A)= 1

This edge will
be counted.

H HODetachment-cutX

Test whether ① and ② holds.
necessary conditions for“G has a connected and loopless ρ-detachment” N[Nagamochi 2006]

75H#edges joining X and V#connected components in G - X

condition for connectivity

4+3+2 = 9O

2①②A

8for v = C#edges incident to v in GC

G=(V, E)‫‏‬

ρ(H) = (1,1,1,1)‫‏‬r(H)=r(O)=ρ(O) = (2)‫‏‬
ρ(N) = (2)‫‏‬r(N)=ρ(C) = (4,3,2)‫ ‏‬r(C)=ρ(A) = (3)‫‏‬
r(A)=

4
1131

condition for degree specification15

Detachment-cutConsider the degree constraint ρ(v i)of each vertex v i for 1 ≦ i ≦ r(v).

O

For a vertex v i ∈T, ρ (v i ) = val (ℓ (v i )).
For a vertex v i ∈T,ρ (v i ) = val (ℓ (v i )) - deg (v i ;T) + 1.
Cρ:4-2+1=3CCCO H

H H H H H H H H

N

G=(V, E)‫‏‬

ρ:3-2+1=2N

HAA

ρ:4-3+1=2CC

ρ(H) = (1,1,1,1)‫ ‏‬r(H)=r(O)=ρ(O) = (2)‫‏‬
ρ(N) = (2)‫‏‬r(N)=ρ(C) = (4,3,2)‫ ‏‬r(C)=ρ(A) = (3)‫‏‬r(A)=41131

Test whether there exists multitreessatisfying the degree constraint
by considering a “detachment” of G.

a current tree T

Detachment-cut
Test whether T can be extended multitrees satisfying degree constraint.

input

residual

H 13O2
N1C6
NH 2OH 1
CH 10
NC 1
CO 2
CC 5

HONC
ANH
CH
COCC
NA
CA

rightmostpathCC

fixedCCO H

CN

H H H H H H H H

411
3
1
13
1
112OCA

H

N

G=(V, E)‫‏‬r(H) =4+4r(H)=r(O) =1+1r(O)=r(N)=1r(N)= 0+r(C) =1+3r(C)=r(A) = 1

a current tree T16

Detachment : A Reverse Operation of Contraction
WbWabWdadcWc
A multigraph GA detachment H of G17