Searching Protein 3-D Structures in Faster Than Linear Time
Tetsuo Shibuya
Human Genome Center, Institute of Medical Science,
University of Tokyo
4-6-1 Shirokanedai, Minato-ku, Tokyo 108-8639, JAPAN,
E-mail: tshibuya@hgc.jp

Abstract
Searching for similar structures from a 3-D structure database of proteins is one of
the most important problems in post-genomic computational biology. To compare two
structures, we ordinarily use a measure called the RMSD (root mean square deviation)
as the similarity measure. We consider a very fundamental problem of finding all thesubstructures whose RMSDs to the query are within some given threshold, from a 3-D
structure database. The problem also appears in many other fields, such as computer
vision and robotics. In this paper, we propose the first algorithm that runs in faster thanlinear time in average. Our new algorithm runs in average-case O(m + N/m1−ϵ ), where
N is the database size, m is the query length, and ϵ is an arbitrary small constant suchthat 0 < ϵ < 1. It is a significant improvement over previous algorithms on the problem,
considering that the best known worst-case time complexity of the problem is O(N log m),
and the best known average-case (expected) time complexity of the problem was O(N ).

1

Introduction

With the aid of various state-of-the-art laboratory methods such as NMR (Nuclear Magnetic Resonance), more and more 3-D structures of bio-molecules, especially the proteins, aredetermined. For example, the number of structures contained in the PDB (Protein Structure Data Bank) database [3] was only around 1, 000 in 1993, but it becomes more than58, 000 in June 2009, and it is still increasing very rapidly (about 20% per year). Moreover,tremendous number of protein structures are being predicted on super computers these days.
Database searching on such databases has many applications in various ﬁelds of molecular biology [2, 12, 14, 21]. For example, we can predict the functions of some protein if its structure
is similar to another protein with known functions, as proteins with similar structures tend
to have similar functions. Hence, 3-D structure database searching becomes more and more
important in computational biology, and faster searching techniques are seriously needed.
A protein is a chain of amino acid residues, and its structure is usually represented by asequence of 3-D coordinates of representative atoms of residues (atoms named Cα are usuallyused). There are also many other important chain molecules in living cells, such as DNAs,
RNAs, and glycans, and they can also be represented by a sequence of 3-D coordinates. No twomolecular structures have exactly the same structure. Therefore we need to search for similar
structures from the structure databases, instead of searching for exactly the same structure. In
1this paper, we deal with a problem of searching for similar structures from structure databasesof chain molecules, each of which is represented by a sequence of 3-D coordinates.
There have been developed tremendous number of algorithms for comparing/searchingprotein structures, which can be categorized roughly into two types. One is a group of algorithms that compare two structures geometrically in the 3-D space, using the coordinatesof structures as their inputs. They consider that the structures are rigid or near-rigid, andsuperimpose (substructures of) the two structures by rotating and translating one of them.
The other is a group of algorithms that uses more abstract information of the structures, suchas the secondary structure elements (SSEs) [23, 24]. In this paper, we focus on the ﬁrst typeof the algorithms, i.e., we compare the sequences of coordinates without any abstraction.
Biologists ordinarily use a similarity measure called the RMSD (root mean square deviation) [1, 11, 18, 19, 21, 25, 26] to compare two structures of bio-molecules. There are also manyother measures, but many of them are just variants of the RMSD [17]. It is the most fundamental measure to determine geometric similarity between two same-length sequences of 3-D
coordinates, which is also used in various other ﬁelds, such as computer vision and robotics.
(The problem minimizing the RMSD is called the least squares ﬁtting problem in computer
vision.) It is deﬁned as the square root of the minimum value of the average squared distance
between each pair of corresponding atoms, over all the possible rotations and translations.
(See section 2.2 for more details.) In this paper, we consider the following problem related to
the RMSD, which is one of the most fundamental computational problems in computational
molecular biology.
3-D Substructure Search Problem
We are given a text structure P and a query structure Q, both of which are represented by
the 3-D coordinates of the residues of the structures. The problem is to ﬁnd all the positionsof the substructures of P whose RMSDs to Q are at most a given ﬁxed threshold c, without
considering insertions or deletions.1
In case the database consists of more than one structure, we can reduce the problem into
the above single-text problem by concatenating all the database structures into a single text
structure and ignoring substructures that cross over the boundaries of concatenated structures.
In general, c is set to a ﬁxed constant proportional to the distance between two adjacent atomsof the chain molecules. In the case of protein structures, the distance between two adjacentCα atoms is always around 3.8Å, while structural biologists say that two protein structuresare similar to each other if their RMSD is at most around 2Å.
History
In 1976, Kabsch [18, 19] proposed a sophisticated mathematical RMSD computation algorithm
based on the singular value decomposition that computes the RMSD between two structuresof size n in O(n) time. If we directly apply the algorithm to the 3-D substructure search
problem, it can be solved in O(N m) time, where N is the length of the text structure andthe m is the length of the query structure. In 1987, Schwartz and Sharir [25] proposed an
1

The problem becomes much harder if we take deletions and insertions into account. Even the simplest
version of the problem is known to be NP-hard [15].

2algorithm based on the fast Fourier transform technique that runs in O(N log m) time, whichis the best-known worst-case time complexity.2 More recently, Shibuya [27] proposed a linear
expected-time (average-case time) algorithm, which was the best-known time complexity before
this paper, though its worst-case time complexity is O(N m). To analyze the average-case timecomplexity, he used the freely-jointed chain model [4, 10, 13, 22] as a model of random protein
structures. The model is is widely used for analyzing properties of various chain molecules (seesection 2.3). The model was also used in [28] to analyze the average time complexity of analgorithm on structures of chain molecules. According to [27], there is high consistency between
the theoretical results deduced from the freely-jointed chain model and the experimental resultson the whole PDB database, and it is very reasonable to use the model for analyzing the averagecomputational complexity of algorithms for molecular structures.
Shibuya [27] also proposed several preprocessing algorithms that achieves even faster expected query time complexity. One is an O(N log2 N )-time and O(N log N )-space prepro√
cessing algorithm with expected query time complexity of O(m + N/ m). Another is anO(N log N )-time and O(N )-space preprocessing algorithm with expected query time complexity of O( √Nm + m log(N/m)). But the expected query time complexity of our algorithm in thispaper is theoretically better than any of these algorithms, even though our algorithm does not
require any preprocessing.
Our results
The expected time complexity of the algorithm in [27] seems to be optimal, as we apparently
require Ω(N ) time to sequentially load all the data into memory. But in case we can randomly
access the input data, there could be algorithms faster than O(N ). For example, the famousBoyer-Moore algorithm [5], which is one of the most widely-used searching algorithms forordinary texts, achieves O(m + N/ min(m, |Σ|)) average time complexity, where N is the text
size, m is the query size and |Σ| is the alphabet size, by permitting random access to thecharacters in the text.
In this paper, we make a breakthrough by proposing the ﬁrst algorithm whose expected
time complexity is better than linear for the above 3-D substructure search problem. It runs
in expected O(m + N/m1−ϵ ) time, where ϵ is an arbitrary small constant such that 0 < ϵ <1, under the assumption that we can randomly access any data in the database (i.e., thecoordinates of any atoms in the text structure). Note that this time complexity is even better
than the query time complexity of the previous algorithms with preprocessing [27], though
our new algorithm does not require any preprocessing. Our analysis of the expected timecomplexity is based on the freely-jointed chain model, which has also been used in previous
work [27, 28] for the analysis of expected time complexities of previous algorithms.
The organization of this paper
In section 2, we describe the notations used in this paper and previous related work as preliminaries. In section 3, we describe our algorithm. In section 4, we analyze the expected(average-case) time complexity of our algorithm. In section 5, we conclude our results.

2

The original algorithm runs in O(N log N ) time. See [27] for the technique to improve it into O(N log m).

3

2

Preliminaries

2.1

Notations and Definitions

A structure of a chain molecule is represented by a notation like S = {⃗s1 , ⃗s2 , . . . , ⃗sn }, where⃗si denotes the 3-D coordinates of the i-th atom. The length n of S is denoted by |S|. A
structure S[i..j] = {⃗si , ⃗si+1 , . . . , ⃗sj } (1 ≤ i ≤ j ≤ n) is called a substructure of S. For
two structures S = {⃗s1 , ⃗s2 , . . . , ⃗sn } and T = {⃗t1 , ⃗t2 , . . . , ⃗tn′ }, the concatenated structure{⃗s1 , ⃗s2 , . . . , ⃗sn , ⃗t1 , ⃗t2 , . . . , ⃗tn′ } is denoted by S + T. R · S denotes the structure S rotated by therotation matrix R, i.e., R · S = {R⃗s1 , R⃗s2 , . . . , R⃗sn }.
⃗v t denotes the transpose of the vector ⃗v and AT denotes the transpose of the matrix A.
trace(A) denotes the trace of the matrix A. |⃗v | denotes the norm of the vector ⃗v . ⃗0 denotes thezero vector. ⟨x⟩ denotes the expected value of x. var(x) denotes the variance of x. P rob(X )denotes the probability of the event X .
In the rest of this paper, P = {⃗p1 , p⃗2 , . . . , p⃗N } denotes the text structure and Q ={⃗q1 , ⃗q2 , . . . , ⃗qm } denotes the query structure. Note that the size m of the query could be
an arbitrary number such that 1 < m ≤ N . The 3-D substructure search problem is to ﬁnd all
the positions i such that the RMSD (see section 2.2 for its deﬁnition) between P[i..i + m − 1]
and Q is at most a given ﬁxed threshold c.

2.2

RMSD: Root Mean Square Deviations

The RMSD (root mean square deviation) [1, 11, 18, 19, 25, 26] is the most widely-used geometric similarity measure between two sequences of 3-D coordinates that represent molecular
structures. The RMSD between two 3-D coordinates sequences S √= {⃗s1 , ⃗s2 , . . . , ⃗sn } and T =∑⃗⃗⃗{t1 , t2 , . . . , tn } is deﬁned as the minimum value of ER,⃗v (S, T) = 1 ni=1 |⃗si − (R · ⃗ti + ⃗v )|2nover all the possible rotation matrices R and translation vectors ⃗v . Let RM SD(S, T) denotethe minimum value, and let R̂(S, T) and ⃗vˆ(S, T) denote the rotation matrix and the translation
vector that minimizes ER,⃗v (S, T).
Kabsch [18, 19] proposed an eﬃcient linear-time algorithm to compute RM SD(S, T),R̂(S, T) and ⃗vˆ(S, T) as follows. If the rotation matrix R is ﬁxed, ER,⃗v (S, T) is known tobe minimized when the centroid (center of mass) of R · T is translated to the centroid of S by
the translation vector ⃗v , regardless of what the rotation matrix R is. It means that ⃗vˆ(S, T)
can be computed in linear time if we are given R̂(S, T). Moreover, it also means that the
problem of computing the RMSD can be reduced to a problem of ﬁnding R (i.e., R̂(S, T))′ (S, T) = ∑n |⃗⃗ 2that minimizes ER
i=1 si − R · ti | , by translating both S and T so that both oftheir centroids are moved to the origin of the coordinates, which can be done in linear time.
If both structures have been already translated so that both centroids are moved to the ori∑
gin, we can compute R̂(S, T) in linear time as follows [1, 18, 19]. Let J = ni=1 ⃗si · ⃗ti t . Clearly,∑n′ (S, T) can be described assi t⃗si + ⃗ti t⃗ti ) −J can be computed in O(n) time. Then ER
i=1 (⃗
2 · trace(R · J), and trace(R · J) is maximized when R = V U T , where U ΛV is the singularvalue decomposition (SVD) of J. Thus R̂(S, T) can be obtained from J in constant time, as
J is a 3 × 3 matrix and the SVD can be computed in O(d3 ) time for a d × d matrix [16]. Notethat there are degenerate cases where det(V U T ) = −1, which means that V U T is a reﬂection
matrix. See [1, 11] for the details of the degenerate cases. Finally, we can compute the RMSD

4

in linear time once we have obtained R̂(S, T). In total, we can compute the RMSD in O(n)time.
There exists an easily-computable lower bound of the RMSD [27], as in the following.
Let Ulef t denote {⃗u1 , ⃗u2 , . . . , ⃗u⌊k/2⌋ } and Uright denote {⃗u⌊k/2⌋+1 , ⃗u⌊k/2⌋+2 , . . . , ⃗u2·⌊k/2⌋ } for astructure U = {⃗u1 , ⃗u2 , . . . , ⃗uk }. Let G(U) denote the centroid of the structure U,√ i.e., G(U) =1 ∑kui . Let F (U) denote |G(Ulef t ) − G(Uright )|/2, and let D(S, T) denote 2 · |Slef t |/|S| ·i=1 ⃗k|F (S) − F (T)| for two structures such that |S| = |T|. Note that |Slef t |/|S| ≈ 1/2. It is known

that D(S, T) is always smaller than or equal to RM SD(S, T). Our algorithm in section 3
utilizes this lower bound.

2.3

Freely-Jointed Chain Molecular Model

To analyze the average-case (expected) time complexity of any algorithm, we need some model
for random inputs. For example, we usually use just random character strings as a model ofrandom inputs to the algorithms on textual strings in theoretical computer science.
The freely-jointed chain model [4, 10, 13, 22] is a simple and well-known fundamental model
for average chain molecular structures, and is often used in analyses of the chain molecularbehaviors in theoretical molecular physics. The model is also called the random-walk chain
model, or just the ideal chain model. In the model, we assume that the chain molecules can beconsidered as random walks. It ignores many physical/chemical constraints, but it is known
to reﬂect the behavior of real molecules very well. All the previous expected time complexity
analyses of the algorithms on protein 3-D structures are based the model [27, 28]. According
to [27], the experimental results on the PDB database consist with their theoretical result
based on the freely-jointed chain model. Following them, we also assume in this paper thatthe structures in the database (i.e., text structures) follow the freely-jointed chain model.3
Consider a chain molecule S = {⃗s0 , ⃗s2 , . . . , ⃗sn } of length n+1, in which the distance between
any two adjacent atoms is ﬁxed to some constant r.4 In the freely-jointed chain model, a bond
between two adjacent atoms, i.e., ⃗bi = ⃗si+1 −⃗si , is considered as a random vector that satisﬁes|⃗bi | = r, and ⃗bi is considered to be independent from any other bond ⃗bj (j ̸= i). In this paper,
we let r = 1 by considering the distance between two adjacent atoms as the unit of distance. If
n is large enough, the distribution of the end-to-end vector ⃗sn − ⃗s0 is known to converge to theGaussian distribution in 3-D space, in which ⟨⃗sn − ⃗s0 ⟩ = ⃗0 and ⟨|⃗sn − ⃗s0 |2 ⟩ = n [4, 10]. In the
distribution, the probability (or probability density) that ⃗sn − ⃗s0 is located at some position(x, y, z) is given as follows:
Wn (x, y, z)dxdydz = (3

3 3 −3(x2 +y2 +z 2 )/2n
)2 edxdydz.
2πn(1)

The Algorithm

We describe our new algorithm in this section. We ﬁrst introduce in section 3.1 an eﬃcientlycomputable lower bound for the RMSD between each text substructure and the query structure. In section 3.2, we describe a new algorithm that utilizes the lower bound. The algorithm
3
4

On the other hand, we give no assumption on the query structures.
In the case of proteins, the distance between two adjacent Cα atoms is fixed to 3.8Å.

5

P[i..i+m-1]i'= ( -1) mod
i

i,

P

i,1

p1 p2

Text Structure P

w

P

w

i,

P

i,2

i, k

p +1
i,h

P[i+w-1-i']

Query Structure Q

Q[w-i']

Q

i,1

Q

Q

i,2

i, k

Figure 1: Division of the text and query structures for computing the lower bound Bi of
RM SD(P[i..i + m − 1], Q). The black ellipses represent the atoms whose coordinates are used
for computing the lower bounds, while the white ellipses represent the dismissed atoms.
ﬁrst searches for all the substructures with lower bounds that are at most the given threshold c, as candidates of the similar substructures. Then the algorithm checks each candidate
whether its RMSD value is actually at most the threshold c. There are two keys for achieving
the better expected time. One is that the lower bounds can be checked without accessing allof the coordinates of the text structure, as is described in section 3.1. The other is that not
all the lower bounds are actually computed, as is described in section 3.2.

3.1

A Lower Bound of the RMSD

To achieve the time complexity better than linear time, we cannot access all the coordinatesof the text structure P.5 Let k be some positive integer such that k ≤ m/2, and w be somem+2
positive integer such that w ≤ 2·(k+1). Note that we will let in section 4 w be some integer in
1−ϵΘ(m 1 ) and k be some constant such that k > 4/ϵ1 , where ϵ1 is an arbitrary small constantsuch that 0 < ϵ1 < 1, to achieve the expected time complexity of O(m + N/m1+ϵ ) (ϵ is anarbitrary small constant such that ϵ1 < ϵ < 1). In our algorithm, we use only coordinates ofatoms whose indices are multiples of w, i.e., P[j · w] for some j, to obtain lower bounds of theactual RMSD values.
From now on, we present a lower bound (Bi in the following) for RM SD(P[i..i+m−1], Q),
which satisﬁes the above requirement, i.e., it must be computable with only the text coordim⌋. Let p⃗i,j denotenates whose indices are multiples of w. Let i′ = (i − 1) mod w, and h = ⌊ k·w′′
P[i−1+j ·w −i ] and ⃗qi,j denote Q[j ·w −i ]. Let Pi,j = {⃗pi,(j−1)h+1 , p⃗i,(j−1)h+2 , . . . , p⃗i,j·h }, andQi,j = {⃗qi,(j−1)h+1 , ⃗qi,(j−1)h+2 , . . . , ⃗qi,j·h }. Figure 1 shows the locations of these structures on
P and Q. In the following, we describe a lower bound for RM SD(P[i..i + m − 1], Q) that canbe computed byusing only the coordinates in Pi,1 +Pi,2 +. . .+Pi,k and Qi,1 +Qi,2 +. . .+Qi,k .
√Let Bi beh ∑k2
j=1 {D(Pi,j , Qi,j )} , where D(S, T) is the lower bound of RM SD(S, T)mintroduced in section 2.2. Let R̂i,j and ⃗vˆi,j denote the rotation matrix and the translation
vector that optimize RM SD(Pi,j , Qi,j ). Similarly, let R̂i and ⃗vˆi denote the rotation matrix
5

Moreover, we cannot compute all the lower bounds for all the positions of the text structure. We will handle
this problem in section 3.2.

6

and the translation vector that optimize RM SD(P[i..i + m − 1], Q). Then, we can prove thatBi ≤ RM SD(P[i..i + m − 1], Q) as follows:
Bi

v
u
k
uh ∑
u
= t
{D(Pi,j , Qi,j )}2

m j=1

≤

v
u
k
uh ∑
u
t
{RM SD(Pi,j , Qi,j )}2

m j=1

v
u
k
h
uh ∑
1∑
u
= t
{
{⃗
pi,(j−1)h+ℓ − R̂i,j · (⃗qi,(j−1)h+ℓ − ⃗vˆi,j )}2 }m j=1 h ℓ=1

v
u
k
h
uh ∑
1∑
u
≤ t
{
{⃗
pi,(j−1)h+ℓ − R̂i · (⃗qi,(j−1)h+ℓ − ⃗vˆi )}2 }m j=1 h ℓ=1v
uu j=h·ku1 ∑= t{⃗pi,j − R̂i · (⃗qi,j − ⃗vˆi )}2

m j=1

≤

v
u ∑
u1 m
t
{⃗
pi+j−1 − R̂i · (⃗qj − ⃗vˆi )}2

m j=1= RM SD(P[i..i + m − 1], Q).

3.2

(2)

The Search Algorithm

The strategy of our algorithm is very simple and is as follows: We ﬁrst enumerate all the
positions i such that Bi ≤ c as candidate positions of the substructures similar to the query,where c is the given threshold. Then we check whether the RMSD between the query andeach candidate substructure, i.e., RM SD(P[i..i + m − 1], Q) for a candidate at position i, is
actually less than or equal to c. But we cannot compute all the Bi values for all the positions of
the text structure if we want to achieve faster algorithm than the linear time, as it apparently
requires at least Ω(N ) time. From now on, we discuss how to enumerate all the positions suchthat Bi ≤ c for the given threshold c without computing all the Bi values.
Let P▽
pi·w , p⃗(i+1)·w , . . . , p⃗(i+h−1)·w } for integers i such that 1 ≤ i ≤ N
i denote {⃗w −h+1, and△let Qj denote {⃗qj , ⃗qj+w , . . . , ⃗qj+(h−1)w } for integers j such that 1 ≤ j ≤ m−h·w+1.6 Let Ld be
the set of integers i such that Pi,1 = P▽
d . Explicitly, Ld = {d · w − w + 1, d · w − w + 2, . . . , d · w},
which means |Ld | = w. Let Md denote a set of integers i such that i ∈ Ld and Bi ≤ c. Let∑⌊N/w⌋−h+1
M all = M1 ∪ M2 ∪ . . . ∪ M⌊N/w⌋−h+1 . Notice that |M all | = i|Mi |, as Md and Md′′
allare disjoint sets if d ̸= d . M is the set of the candidate positions that we want to obtain,i.e., Bi ≤ c √iﬀ i ∈ M all . From now on, we consider how to eﬃciently obtain the√set M all .
h ∑km2
If Bi = m
j=1 {D(Pi,j , Qi,j )} ≤ c, it is easy to see that D(Pi,j , Qi,j ) ≤h · c for eachj (1 ≤ j ≤ k). Let Kd,j be a set of integers i such that i ∈ Ld and D(Pi,j , Qi,j ) ≤√m
h · c.

Let Kd = Kd,1 ∩ Kd,2 ∩ . . . ∩ Kd,k . Notice that Md ⊆ Kd and consequently M all ⊆ K all , where
K all = K1 ∪ K2 ∪ . . . ∪ K⌊N/w⌋−h+1 .7 Thus, we consider how to compute the set K all at ﬁrst.
In our algorithm, we ﬁrst obtain the set Kd for each d (1 ≤ d ≤ ⌊ N
w ⌋ − h + 1) as follows.
√

Recall that D(S, T) =672 · |Slef t |/|S|·|F (S)−F (T|. (See section 2.2 for the deﬁnition of F (·).)

Note that Pi,j = P▽
, and Qi,j = Q△
, where i′ = (i − 1) mod w.
⌊(i−1)/w⌋+(j−1)h+1(j−1)h·w+w−i′′
If d ̸= d , Kd and Kd′ are disjoint sets as Ld and Ld′ are disjoint sets.

7

Thus Kd,j canbe obtained by ﬁnding all i ∈ Ld such that F (Pi,j )−c′ ≤ F (Qi,j ) ≤ F (Pi,j )+c′ ,√mNwhere c′ = 2·⌊h/2⌋· c. We can compute the values F (P▽d ) for all d (1 ≤ d ≤ w − h + 1)◦
in O( Npi·w , p⃗(i+1)·w , . . . , p⃗(i+h′ −1)·w }, where h′ = ⌊h/2⌋.
w ) time as follows. Let Pi denote {⃗▽▽ lef t▽ right
As F (Pd ) = |G(Pd) − G(Pd)|/2, we can compute the value of F (P▽d ) in constant′ + 1). (See section 2.2
time if we are given all the vectors G(P◦i ) for all i (1 ≤ i ≤ N−hw
for the deﬁnition of G(·).) Moreover, all these G(P◦i ) values can be computed in O( Nw ) time,△◦◦′as G(Pi+1 ) = G(Pi ) + (⃗p(i+h′ )·w − p⃗i·w )/h . On the other hand, F (Qℓ ) can be computed△in O(|Qℓ |) = O(h) time for each ℓ. Thus, we can obtain the values F (Q△ℓ ) for all ℓ =(j − 1)h · w + w − i′ such that 0 ≤ i′ < w and 1 ≤ j ≤ k, in O(h · w · k) = O(m) time.
We can enumerate all the integers in Kd = Kd,1 ∩ Kd,2 ∩ . . . ∩ Kd,k by utilizing the kdimensional range search [6, 7] from a set of w points in the k-dimensional space. We can ﬁnd all
the i ∈ Kd in O(|Kd | + logk−1 w) time after an O(w logk−1 w)-time preprocessing on the query∑ N −h+1

k−1w
structure [6, 7]. Thus K all can be obtained in O(m+w logk−1 w+ N
w+ d=1|Kd |) =
w log
k−1
N
all
O(m + (w + w ) log
w + |K |) time in total.
Then we can obtain M all , i.e., the set of all the positions i such that i ∈ K all and Bi ≤ c,by checking whether the Bi value is actually at most the given threshold c for all i ∈ K all .
We can check each Bi value only in O(k) time, as we have already computed all the necessary△
allvalues of F (P▽d ) and F (Qj ). It takes O(k · |K |) time in total.
Finally, we can obtain the positions with the RMSDs actually smaller than or equal toc, by checking whether RM SD(P[i..i + m − 1], Q) is actually at most the threshold c for alli ∈ M all , which takes O(m) time for each with Kabsch’s algorithm described in section 2.2. It
takes O(m · |M all |) time in total.
k−1
w + m · |K all |),Thus, the total time complexity of our algorithm is O(m + (w + Nw ) logallall
all
as |M | ≤ |K |. It could be O(N m) at worst, as |K | could be at most O(N ). But we willshow in the next section that such bad cases must be very rare under the assumption that the
text structure P follows the freely-jointed chain model.8 We will show that the expected timecomplexity of our algorithm is only O(m + N/m1−ϵ ), where ϵ is an arbitrary small positive
constant such that 0 < ϵ < 1, by setting w and k appropriately.

4

The Expected Time Complexity of the Algorithm

Consider a structure S = {⃗s1 , ⃗s2 , . . . , ⃗s1+(h−1)·w } of length 1 + (h − 1) · w, which follows therandom walk model. Let ⃗bi = ⃗si+1 − ⃗si . bi is a random vector that satisﬁes |⃗bi | = 1, and⃗bi and ⃗bj are independent if i ̸= j. ⃗si can be represented as ⃗s1 + ∑i−1 ⃗bj . Let S▽ denotej=1{⃗s1 , ⃗s1+w , ⃗s1+2w , . . . , ⃗s1+(h−1)·w }, and h′ = ⌊h/2⌋. Let H(U) = (G(Ulef t ) − G(Uright ))/2 forsome structure U.9 Then the following equation holds:
H(S▽ ) =1{G(S▽lef t ) − G(S▽right )}2

8

The query structure Q does not have to follow any model. But the same can be said if the query structure
Q follows the freely-jointed chain model, instead of the text structure P.
9
Notice that F (U) = |H(U)|.

8

′

′

(i−1)·w(i−1)·w
h
2h
∑
∑
∑
1 ∑
⃗
⃗bj )}
{
(⃗
s
+
b
)
−
(⃗
s
+
1
j
1
2 · h′ i=1
′
j=1
j=1
i=h +1

=

′

′

2h
−1
i·w
h
i·w
∑
∑
∑
1 ∑
′
⃗
⃗bj )}.
= −{ i·(bj ) +(2h − i) · (2 · h′ i=1′
i=h +1j=1+(i−1)w
j=1+(i−1)w(3)⃗bi }/2h′ if i ≤ h′ ·w, and {(2h′ −1−⌊ i−1 ⌋)·⃗bi }/2h′ if i > h′ ·w. Then
Let ⃗b′i denote {(1+⌊ i−1
w ⌋)·
w
∑
2h′ w ⃗ ′H(S▽ ) can be described as i=1
bi . Let zi denote the z-coordinate of ⃗bi , and zi′ denote the zcoordinate of ⃗b′i . It is easy to see that ⟨zi ⟩ = 0 and var(zi ) = 1/3, as ⃗bi is a random vector that
∑

√∑

′

′

2+δ

·w
2h ·w
′
′
2+δ /
′
satisﬁes |⃗bi | = 1. Then, consider the value Ah,w = 2h,i=1 ⟨|zi − ⟨zi ⟩|⟩i=1 var(zi )where δ is some positive constant. According to Lyapunov’s central limit theorem [20], the∑ ′w ⃗′
distribution of H(S▽ ) = 2hi=1 bi converges to the Gaussian distribution if the above Ah,w′
converges to 0 as 2h · w grows up to inﬁnity, for some δ (δ > 0). It can be proved as follows:

Ah,w =

′ ·w
2h
∑

′ ·w
2h
∑
1
′ 2+δ
⟨|zi |⟩⟨|zi′ |⟩2
≤√
√∑
2+δ
2+δ
∑
′
′
2h ·w
2h ·w
′ 2
′ 2
i=1
i=1
i=1 ⟨|zi |⟩i=1 ⟨|zi |⟩

1

′ ·w
2h
∑

= {

′

⟨|zi′ |⟩2 }−δ/2 = {

i=1

= {

′

h
2h
−1
∑
∑
w
2
·
{
i
+
(2h′ − i)2 }}−δ/212h′ 2 i=1
i=h′ +1

w −δ/2
1 ′
h ·w+
}
→ 0
18
36h′

(2h′ · w → ∞).

(4)

The same discussion can be done for the other two axes (x and y), which are independent
to each other. Thus we conclude that H(S▽ ) converges to a Gaussian distribution in 3-D
space. It is easy to see that ⟨H(S▽ )⟩ = 0. The variance of H(S▽ ) can be computed as follows:var(H(S▽ )) = ⟨|H(S▽ )|2 ⟩ − ⟨|H(S▽ )|⟩2 = ⟨|H(S▽ )|2 ⟩
′ −1
2h
∑

′

2h
∑
i

′

2h − i ⃗ 2· ⃗bi +· bi | ⟩= w · ⟨|
′
2h
h′
i=1
i=h′ +1
′

=

′

h
2h
−1
∑
w
1
1
w ∑
2
{
i
+
(2h′ − i)2 } = h′ · w +≈ h′ · w,′2′
4h i=1612h
6i=h′ +1(5)

It means that the distribution of H(S▽ ) is the same as that of the end-to end vectors of
random walks of length h′ · w/6. Thus, the probability distribution of H(S▽ ) is:
Zh,w (x, y, z)dxdydz = (9 3 −9(x2 +y2 +z 2 )/h′ w)2 edxdydz.
πh′ w(6)Consequently, the probability (probability density) that |H(S▽ )| = r is:
Zh,w (r)dr = 4πr2 (9 3 −9r2 /h′ w
)2 e
dr.
πh′ w

(7)

Integrating Zh,w (r)dr, we obtain P rob(x ≤ F (S▽ ) ≤ y) = P rob(x≤ |H(S▽ )| ≤ y) =√∫y1′
r=x Zh,w (r)dr. Zh,w (r) takes the maximum value at rmax = 3 h w, and Zh,w (rmax ) =9√−1 / πh′ w. Thus P rob(x ≤ F (S▽ ) ≤ y) is at most (y − x) · Z−112e√h,w (rmax ) = 12e (y −x)/ πh′ w for any x and y such that x < y.
To evaluate the expected time complexity of the algorithm proposed in the previous section,
we have to investigate the expected size of Kd under the assumption that the text structure P
follows the freely-jointed
chain model.10 Recall that Kd,j is a set of integers i such that i ∈ Ld√
and D(Pi,j , Qi,j ) ≤ m
h · c. Thus ⟨|Kd,j |⟩ = p · |Ld | = p · w, where p = P rob(D(Pi,j , Qi,j ) ≤√√m· c). As D(S, T) = |F (S) − F (T)|, p equals P rob(F (Qi,j ) − m
h · c ≤ F (Pi,j ) ≤ F (Qi,j ) +
√h
√
m
m
·
c).
According
to
the
above
discussion,
p
is
at
most
p
=
Z
(r
)
·
2
·
max
max
h,w
h
h ·c =√√−1′
24c · e · m/(π · h · h · w), for any i, j and query Q. The value pmax is in O(k w/m), as
h = ⌊m/(k · w)⌋, h′ = ⌊h/2⌋ and c is a ﬁxed constant. As the text structure P follows thefreely-jointed chain model, Pi,j and Pi,j ′ (j ̸= j ′ ) are built by independent random walks. It
means that Kd,j and Kd,j ′ (j ̸= j ′ ) are independently chosen from Ld . Thus the expected sizeof Kd = Kd,1 ∩ Kd,2 ∩ . . . ∩ Kd,k is at most w · pkmax , which is in O((c′ · k)k · w · (w/m)k/2 ),where c′ is some constant. Thus ⟨|K all |⟩ is in O((c′ · k)k · N · (w/m)k/2 ).
Thus, the expected time complexity of the algorithm described in the last section is O(m +k−1(w + Nw +(c′ ·k)k ·N ·wk/2 ·m1−k/2 ). Let w be some integer in Θ(m1−ϵ1 ) and k be somew ) log
constant such that k > 4/ϵ1 , where ϵ1 is an arbitrary small constant such that 0 < ϵ1 < 1.
k−1mThen the above expected time complexity becomes O(m + N · log). The ‘logk−1 m’ termm1−ϵ1
ϵcan be replaced by m 2 where ϵ2 is an arbitrary small constant such that 0 < ϵ2 < 1 − ϵ1 . Letϵ = ϵ1 + ϵ2 . Then, the above time complexity can be simpliﬁed to O(m + N/m1−ϵ ). Note that
we can let ϵ be an arbitrary small constant such that 0 < ϵ < 1.

5

Concluding Remarks

We proposed a new algorithm for the 3-D substructure search problem, which is the problem
to ﬁnd all the similar substructures from a 3-D structure database of chain molecules. The
algorithm runs in average-case O(m + N/m1−ϵ ) time where N is the database size, m is thequery size and ϵ is an arbitrary small constant such that 0 < ϵ < 1, while the best-known
expected time complexity of the problem was O(N ). Our algorithm is the ﬁrst algorithm thatruns in faster than linear expected time.
We improved the expected time complexity, but the coeﬃcient before the complexity islarge in case ϵ is small. It is an open problem whether we can reduce it or not. It is also an
open problem whether we can improve the worst-case time complexity. Our algorithm can also
be applied to similar problems in diﬀerent dimensions. It would be very interesting to apply
it to many other problems in higher dimensions.

Acknowledgement
The author would like to thank Jesper Jansson, Gregory Kucherov, and Kunihiko Sadakane
for invaluable discussions. This work was partially supported by the Grant-in-Aid from the
Ministry of Education, Culture, Sports, Science and Technology of Japan.
10

The same discussion can be done if the query structure Q follows the freely-jointed chain model instead of
the text structure P.

10

References[1] K. S. Arun, T. S. Huang, and S. D. Blostein. Least-squares ﬁtting of two 3-D point sets.
IEEE Trans Pattern Anal. Machine Intell., Vol. 9, pp. 698–700, 1987.
[2] Z. Aung and K.-L. Tan. Rapid retrieval of protein structures from databases. Drug Discovery Today, Vol. 12, pp. 732–739, 2007.
[3] H. M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T. N. Bhat, H. Weissig, I. N.
Shindyalov, and P. E. Bourne. The protein data bank. Nucl. Acids Res., Vol. 28, pp.
235–242, 2000.
[4] R. H. Boyd and P. J. Phillips. The Science of Polymer Molecules: An Introduction Concerning the Synthesis, Structure and Properties of the Individual Molecules That Constitute Polymeric Materials, Cambridge University Press, 1996.
[5] R. S. Boyer and J. S. Moore. A fast string searching algorithm, Commun. ACM, Vol. 20,
pp. 762–772, 1977.
[6] B. Chazelle. Filtering search: a new approach to query-answering, SIAM J. Comput., Vol.
15, No. 3, pp.703–724, 1986.
[7] B. Chazelle. A functional approach to data structures and its use in multidimensional
searching, SIAM J. Comput., Vol. 17, No. 3, pp. 427–453, 1988.
[8] J. W. Cooley and J. W. Tukey. An algorithm for the machine calculation of complex
Fourier series, Math. Comput., Vol. 19, pp. 297–301, 1965.
[9] J. Dayantis and J.-F. Palierne. Monte Carlo precise determination of the end-to-end distribution function of self-avoiding walks on the simple-cubic lattice. J. Chem. Phys., Vol.
95, pp. 6088–6099, 1991.
[10] P.-G. de Gennes. Scaling Concepts in Polymer Physics, Cornell University Press, 1979.
[11] D. W. Eggert, A. Lorusso, and R. B. Fisher. Estimating 3-D rigid body transformations:a comparison of four major algorithms. Machine Vision and Applications, Vol. 9, pp.
272–290, 1997.
[12] I. Eidhammer, I. Jonassen, and W. R. Taylor. Structure comparison and structure patterns. J. Computational Biology, Vol. 7, No. 5, pp. 685–716, 2000.
[13] P. J. Flory. Statistical Mechanics of Chain Molecules, Interscience, New York, 1969.
[14] M. Gerstein. Integrative database analysis in structural genomics. Nat. Struct. Biol.,Suppl., pp. 960–963, 2000.
[15] D. Goldman, S. Istrail and C. H. Papadimitriou. Algorithmic aspects of protein structure
similarity. Proc. 40th Annual Symposium on Foundations of Computer Science, pp. 512–
522, 1999.

11

[16] G. H. Golub and C. F. Van Loan. Matrix Computation. 3rd eds., John Hopkins University
Press, 1996.
[17] H. Hasegawa and L. Holm. Advances and pitfalls of protein structural alignment. Current
Opinion in Structural Biology, Vol. 19, pp. 341–348, 2009.
[18] W. Kabsch. A solution for the best rotation to relate two sets of vectors. Acta Cryst., Vol.
A32, pp. 922–923, 1976.
[19] W. Kabsch. A discussion of the solution for the best rotation to relate two sets of vectors.
Acta Cryst., Vol. A34, pp. 827–828, 1978.
[20] O. Kallenberg. Foundations of Modern Probability, Springer-Verlag, 1997.
[21] P. Koehl. Protein structure similarities. Current Opinion in Structural Biology, Vol. 11,
pp. 348–353, 2001.
[22] H. A. Kramers. The behavior of macromolecules in inhomogeneous ﬂow. J. Chem. Phys.,Vol. 14, No. 7, pp. 415–424, 1946.
[23] E. Krissinel and K. Henrick. Secondary-structure matching (SSM), a new tool for fast
protein structure alignment in three dimensions. Acta Cryst. Sect., Vol. D60, pp. 2256–
2268, 2004.
[24] A. C. R. Martin. The ups and downs of protein topology: rapid comparison of protein
structure. Protein Engineering, Vol. 13, pp. 829–837, 2000.
[25] J. T. Schwartz and M. Sharir. Identiﬁcation of partially obscured objects in two and three
dimensions by matching noisy characteristic curves. Intl. J. of Robotics Res., Vol. 6, pp.
29–44, 1987.
[26] T. Shibuya. Eﬃcient substructure RMSD query algorithms. J. Comput. Biol., Vol. 14,
No. 9, pp. 1201–1207, 2007.
[27] T. Shibuya, Searching protein 3-D structures in linear time. Proc. Conference on Researchin Computational Molecular Biology (RECOMB ’09), LNBI 5541, pp. 1–15, 2009.
[28] T. Shibuya, J. Jansson, and K. Sadakane, Linear-Time Protein 3-D Structure Searchingwith Insertions and Deletions, Proc. 9th Workshop on Algorithms in Bioinformatics, 2009,
to appear.

12