PHYSICAL REVIEW E 78, 046702 共2008兲

Perturbational formulation of principal component analysis in molecular dynamics simulation
Yohei M. Koyama,1,2,* Tetsuya J. Kobayashi,3,1 Shuji Tomoda,2 and Hiroki R. Ueda1,4,5,†
1

Laboratory for Systems Biology, Center for Developmental Biology, RIKEN,
2-2-3 Minatojima-minamimachi, Chuo-ku, Kobe, Hyogo 650-0047, Japan2
Department of Life Sciences, Graduate School of Arts and Sciences, University of Tokyo, Komaba, Meguro-ku, Tokyo 153-8902, Japan
3
Institute of Industrial Science, The University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 153-8505, Japan
4
Functional Genomics Unit, Center for Developmental Biology, RIKEN,
2-2-3 Minatojima-minamimachi, Chuo-ku, Kobe, Hyogo 650-0047, Japan
5
Department of Bioscience, Graduate School of Science, Osaka University, Toyonaka, Osaka 560-0043, Japan
共Received 31 May 2008; revised manuscript received 18 August 2008; published 7 October 2008兲
Conformational fluctuations of a molecule are important to its function since such intrinsic fluctuations
enable the molecule to respond to the external environmental perturbations. For extracting large conformational fluctuations, which predict the primary conformational change by the perturbation, principal component
analysis 共PCA兲 has been used in molecular dynamics simulations. However, several versions of PCA, such as
Cartesian coordinate PCA and dihedral angle PCA 共dPCA兲, are limited to use with molecules with a single
dominant state or proteins where the dihedral angle represents an important internal coordinate. Other PCAs
with general applicability, such as the PCA using pairwise atomic distances, do not represent the physicalmeaning clearly. Therefore, a formulation that provides general applicability and clearly represents the physicalmeaning is yet to be developed. For developing such a formulation, we consider the conformational distribution change by the perturbation with arbitrary linearly independent perturbation functions. Within the second
order approximation of the Kullback-Leibler divergence by the perturbation, the PCA can be naturally interpreted as a method for 共1兲 decomposing a given perturbation into perturbations that independently contributeto the conformational distribution change or 共2兲 successively finding the perturbation that induces the largest
conformational distribution change. In this perturbational formulation of PCA, 共i兲 the eigenvalue measures the
Kullback-Leibler divergence from the unperturbed to perturbed distributions, 共ii兲 the eigenvector identifies the
combination of the perturbation functions, and 共iii兲 the principal component determines the probability change
induced by the perturbation. Based on this formulation, we propose a PCA using potential energy terms, andwe designate it as potential energy PCA 共PEPCA兲. The PEPCA provides both general applicability and clear
physical meaning. For demonstrating its power, we apply the PEPCA to an alanine dipeptide molecule invacuum as a minimal model of a nonsingle dominant conformational biomolecule. The first and second
principal components clearly characterize two stable states and the transition state between them. Positive andnegative components with larger absolute values of the first and second eigenvectors identify the electrostaticinteractions, which stabilize or destabilize each stable state and the transition state. Our result therefore
indicates that PCA can be applied, by carefully selecting the perturbation functions, not only to identify themolecular conformational fluctuation but also to predict the conformational distribution change by the perturbation beyond the limitation of the previous methods.
DOI: 10.1103/PhysRevE.78.046702

PACS number共s兲: 02.70.Rr, 36.20.Ey, 87.15.ap, 89.70.Cf

I. INTRODUCTION

Conformational fluctuations of a molecule are important
to its function, since such intrinsic fluctuations enable themolecule to respond to perturbations in the external environment 关1,2兴. For extracting a large conformational fluctuation,
which is believed to relate to the primary conformationalchange by the perturbation, principal component analysis
共PCA兲 关3兴, also called essential dynamics, has been generally
used in molecular dynamics 共MD兲 simulations 关4–8兴. In general, PCA is performed using Cartesian coordinates. Since
Cartesian coordinates depend on the overall molecular motion, translational and rotational motion have to be eliminated before the calculation of the covariance matrix. Molecular translation can be eliminated uniquely by translating*ym.koyama@gmail.com†

uedah-tky@umin.ac.jp1539-3755/2008/78共4兲/046702共11兲

the molecular center of mass to the origin. However, molecular rotation cannot be removed uniquely by least square fitting to a certain reference structure, except the molecule fluctuating around the single dominant state. This is because thefitting result generally depends on the selected reference
structure 关9–11兴. We refer to this problem as “frame fitting
problem” in this study. To avoid the frame fitting problem,
several PCAs using internal coordinates, which is independent of the overall molecular motion, have been proposedsuch as the dihedral angle PCA 共dPCA兲 关12,13兴 and a PCA
using pairwise atomic distances 关11兴. The dihedral angle represents an important internal coordinate for proteins. However, it is difficult to apply the dPCA to other molecules suchas nucleic acids. PCA using pairwise atomic distances havegeneral applicability, but their physical meaning is unclear.
Furthermore, these methods do not directly relate the identified large conformational fluctuations with the conformational change induced by the environmental perturbationsuch as ligand binding. Therefore, a formulation for charac-

046702-1©2008 The American Physical Society

KOYAMA et al.

PHYSICAL REVIEW E 78, 046702 共2008兲

terizing conformational fluctuations, which 共a兲 achieves general applicability, 共b兲 describes its physical meaning, and 共c兲
predicts the conformational change by the perturbation isunavailable.
Some previous studies more directly tried to relate conformational fluctuations with the conformational change induced by the environmental perturbation. Ikeguchi et al. 关14兴
applied linear response theory for predicting the average conformational change induced by constant force perturbations
from the covariance matrix of atomic coordinates in an unperturbed equilibrium state. Because they used Cartesianatomic coordinates, their method also inherits the frame fitting problem as in the case of Cartesian PCA, which limits
its applicability to molecules with a single dominant state.
Ming et al. 关15,16兴 proposed quantifying the conformational
distribution change by the perturbation with the KullbackLeibler divergence. They used the harmonic approximationof the unperturbed and the perturbed energy function, i.e.,
normal mode approximation, in which each equilibrium distribution becomes a multivariate normal distribution. By thisapproximation, they estimate the Kullback-Leibler divergence by deriving its analytical expression. Although applications of the Kullback-Leibler divergence are not limited tothe multivariate normal distributions, the normal mode approximation limits their method to conformational change
from a single-dominant state to another single-dominant stateby the perturbation. Therefore, a formulation that satisfies theabove mentioned three requirements 共a兲–共c兲 is still not available.
For developing such a formulation, we first introduce a
perturbation with arbitrary linearly independent perturbation
functions 共Sec. II A兲. We then quantify the conformational
distribution change by the Kullback-Leibler divergence,
which is naturally derived from large deviation theory 共Sec.
II B兲. Within the second order approximation of the
Kullback-Leibler divergence, the PCA can be naturally interpreted as a method for 共1兲 decomposing a given perturbation
into perturbations that independently contribute to the conformational distribution change or 共2兲 successively finding
the perturbation that induces the largest conformational distribution change 共Sec. II C兲. We then propose a PCA using
potential energy terms, which is designated as potential energy PCA 共PEPCA兲. The PEPCA gives general applicability
and clear physical meaning 共Sec. II D兲. For demonstrating its
power, we apply the PEPCA to an alanine dipeptide molecule in vacuum as a minimal model of a nonsingle dominant conformational biomolecule 共Sec. III兲. Finally, we discuss the theoretical implications of our perturbational
formulation of PCA, and the applicability of the PEPCA tolarger molecules 共Sec. IV兲.

1
␳共q兲 = e−V共q兲/kBT ,Z

II. THEORY

共1兲

where kB and Z are the Boltzmann constant and the partition
function 共normalizing constant兲, respectively. We perturb themolecular system using perturbation parameters ␭= 共␭1 , . . . , ␭ M 兲T and linearly independent perturbation functions f共q兲 = 关f 1共q兲 , . . . , f M 共q兲兴T as
M

V␭共q兲 = V共q兲 − kBT 兺 ␭i f i共q兲.

共2兲

i=1

Then, the perturbed canonical distribution ␳␭共q兲 is␳␭共q兲 =

1 −V 共q兲/k T ␭·f共q兲−␺共␭兲e ␭ B =e
␳共q兲,
Z␭

共3兲

Z␭
= ln具e␭·f典,
Z共4兲

where we define␺共␭兲 ⬅ ln

and 具¯典 indicates the average under ␳共q兲. We consider the ␭region where ␺共␭兲 exists. This representation enables us toconsider the perturbed distribution 关Eq. 共3兲兴 as an exponential family, which is broadly studied in statistics 关17–21兴. In
the theory of the exponential family, ␭, f共q兲, and ␺共␭兲 arecalled natural parameters, sufficient statistics, and cumulant
generating function, respectively. The nth derivative of thencumulant generating function 兩␺i1¯in共␭0兲 ⬅ ⳵␭⳵i ␺¯共␭兲⳵␭i1 兩␭=␭0 isnknown as the nth cumulant. In particular, the first cumulant is
the average ␺i共␭兲 = 具f i典␭, and the second cumulant is the covariance ␺ij共␭兲 = 具共f i − 具f i典␭兲共f j − 具f j典␭兲典␭, where 具¯典␭ indicates the average under ␳␭共q兲. 具¯典0 is equivalent to 具¯典.
B. Quantification of the conformational distributionchange by the perturbation

For quantifying the conformational distribution change induced by the perturbation, we use the Kullback-Leibler divergence 关15–17,21–24兴, which naturally appears in Sanov’s
theorem in large deviation theory 关21,23,25,26兴, and is defined as
D共␳␭ 储 ␳兲 ⬅

冕␳␭共q兲ln␳␭共q兲dq 艌 0.
␳共q兲

共5兲

The Kullback-Leibler divergence is non-negative, equals
zero if and only if two distributions are identical, and isgenerally asymmetric with respect to the exchange of the two
distributions. Using Eq. 共3兲, Eq. 共5兲 is represented as

A. Perturbation of a molecular system

M

We consider a molecular system that is described by thepotential energy V共q兲, where q = 共q1 , . . . , q3N兲T are the atomic
Cartesian coordinates and N is the number of atoms. Then,the canonical distribution ␳共q兲 at temperature T is

D共␳␭储␳兲 = ␭ · 具f典␭ − ␺共␭兲 = 兺 ␭i␺i共␭兲 − ␺共␭兲.

共6兲

i=1

By using the Taylor expansions of ␺i共␭兲 and ␺共␭兲 around theorigin, Eq. 共6兲 is expanded as

046702-2

PERTURBATIONAL FORMULATION OF PRINCIPAL …
M

D共␳␭储␳兲 =PHYSICAL REVIEW E 78, 046702 共2008兲M

1
1
␺ij共0兲␭i␭ j + 兺 ␺ijk共0兲␭i␭ j␭k + ¯ ,兺
2 i,j=1
3 i,j,k=1

D共␳␭储␳兲 ⬇

max
兩␭兩=␦

u j·␭=0,j=1,. . .,i−1

共7兲

1
1
= 共␦ui兲TC共␦ui兲 = ␦2␴2i . 共11兲2
2
From Eq. 共3兲, the probability change of q induced by the
perturbation ␭ can be represented as

C. Perturbational formulation of principal component analysis

Here, we consider a PCA with respect to the conformational distribution change induced by the perturbations formulated in the preceding sections. We denote the covariancematrix of f共q兲 in an unperturbed distribution ␳共q兲 as C,
where Cij = ␺ij共0兲 = 具共f i − 具f i典兲共f j − 具f j典兲典. Since the covariancematrix is a positive definite matrix 关we assume each varianceof f i共q兲 is nonzero, i.e., not constant兴, we can diagonalize thematrix using the orthogonal matrix U = 共u1 , . . . , u M 兲 and positive eigenvalues ␴2i . We assume the eigenvalues are sorted indescending order. On the basis of the eigenvectors, theKullback-Leibler divergence in the second order approximation 关Eq. 共7兲兴 is represented asln

共8兲
where we define␰ ⬅ UT␭,

共9兲

which are the perturbation parameters corresponding to theeigenvectors. Thus, within the second order approximationof the Kullback-Leibler divergence, the given perturbation ␭
is decomposed into the perturbations ␰i, which independently
contribute to the Kullback-Leibler divergence by 21 ␴2i ␰2i .
Therefore, the eigenvalue ␴2i measures the conformational
distribution change induced by the decomposed perturbation
␰i. Components of the eigenvector ui identifies the combinations of the perturbation functions that induces the decomposed perturbation ␰i.
Since the covariance matrix is symmetric, the ith eigenvector ui can be also characterized variationally as ␭, which
maximizes ␭TC␭ with constraints 兩␭兩 = const and u j · ␭ = 0, j
= 1 , . . . , i − 1 关3,27兴. The second order approximation of the
Kullback-Leibler divergence is in the quadratic form 21 ␭TC␭.
Therefore, we can interpret the diagonalization of the covariance matrix by successively finding the perturbation ␭
which maximizes the Kullback-Leibler divergence D共␳␭ 储 ␳兲under the constraints兩␭兩 = ␦,␦ such that the second order,
approximation of D共␳␭储␳兲 is valid,

共10a兲

u j · ␭ = 0,

共10b兲

j = 1, . . . ,i − 1.

Equivalently, it can be expressed as

␳␭共q兲= ␭ · f共q兲 − ␺共␭兲
␳共q兲

冉

1
= ␭ · f共q兲 − 具f典 · ␭ + ␭TC␭ + ¯
2冊M= g共q兲 · ␰ −M1
兺 ␴ 2␰ 2 − ¯
2 i=1 i i

冉冊1
= 兺 gi共q兲␰i − ␴2i ␰2i − ¯ ,
2
i=1

共12兲

共13兲

where we define g共q兲 = 关g1共q兲 , . . . , g M 共q兲兴T as
g共q兲 ⬅ UT关f共q兲 − 具f典兴.

M

1
1
1
D共␳␭储␳兲 ⬇ ␭TC␭ = 共UT␭兲T共UTCU兲共UT␭兲 = 兺 ␴2i ␰2i ,
2
2
2 i=1

兩␭兩=␦

1 T
␭ C␭
2

u j·␭=0,j=1,. . .,i−1

1. The Kullback-Leibler dithe nth order coefficient is 共n−2兲!nvergence can also be represented by 具f典␭ instead of ␭ asshown in Appendix A.

max

共14兲

These functions g共q兲 can be considered as the principal components of f共q兲 under ␳共q兲. Thus, in the second order approximation of the Kullback-Leibler divergence, the probability change induced by the perturbation ␰i is determined
only by the ith principal component gi共q兲 and the ith eigenvalue ␴2i . Moreover, for the perturbation where the first orderapproximation of the cumulant generating function is valid,
the probability change is determined only by the ith principal
component. Perturbation ␰i ⬎ 0 increases the probability of q
where gi共q兲 ⬎ 0, and decreases q where gi共q兲 ⬍ 0. Equation
共14兲 indicates that these principal components g共q兲 are
具gi典 = 0,

共15兲

具gig j典 = ␴2i ␦i,j .

共16兲

Thus, the principal components make orthogonal functions
with respect to the inner product
具gi,g j典 ⬅

冕

gi共q兲g j共q兲␳共q兲dq.

共17兲

In summary, in our perturbational formulation, within thesecond order approximation of the Kullback-Leibler divergence, the PCA can be interpreted as a method for 共1兲 decomposing a given perturbation ␭ into perturbations that independently contribute to the conformational distributionchange measured by the Kullback-Leibler divergence or 共2兲
successively finding the perturbation that induces the largest
conformational distribution change under the constraints
关Eqs. 共10兲兴. Then, 共i兲 the eigenvalue ␴2i measures the
Kullback-Leibler divergence by the perturbation ␰i as 21 ␴2i ␰2i关Eq. 共8兲兴. 共ii兲 Eigenvector ui identifies the combination of
perturbation functions that induces the ith largest conformational distribution change. 共iii兲 The principal component

046702-3

PHYSICAL REVIEW E 78, 046702 共2008兲

KOYAMA et al.

gi共q兲 determines the probability change of q induced by the
perturbation ␰i 关Eq. 共13兲兴. In Appendix B, we give perturbational formulations of Cartesian PCA and dPCA.
D. Perturbation of potential energy terms

To develop a PCA that satisfies the three requirements
共a兲–共c兲 mentioned in the Introduction, we consider a perturbation of potential energy terms. We perturb some potential
energy terms
V共q兲 = 关V1共q兲, . . . ,VM 共q兲兴T

共18兲

in potential energy V共q兲 with perturbation term in Eq. 共2兲
M

M

i=1

i=1

− kBT 兺 ␭i f i共q兲 = − 兺 ␭iVi共q兲.

共19兲This perturbation changes the weight of the potential energy
term Vi by 共1 − ␭i兲. This gives perturbation functions
f共q兲 =

冉冊

V1共q兲
V M 共q兲 T, ... ,
.
k BT
k BT

共20兲

Thus, a PCA using these perturbation functions 共PEPCA兲 canbe interpreted as a method for 共1兲 decomposing a given potential energy terms perturbation 关Eq. 共19兲兴 into perturbations
that independently contribute to the conformational distribution change or 共2兲 successively finding the potential energy
terms perturbation 关Eq. 共19兲兴 that induces the largest conformational distribution change under the constraints Eqs. 共10兲.
Since the potential energy terms are always defined in MD
simulations, and they are invariant with respect to the overallmolecular motion, the requirement 共a兲 in the Introduction issatisfied. A potential energy term represents the molecularinteraction, therefore, the requirement 共b兲 is satisfied. The
requirement 共c兲 is automatically satisfied in our perturbational formulation.
Perturbation ␭i = 1 eliminates the potential energy term Vi.
␰ corresponding to ␭ = 共1 , . . . , 1兲T is冉兺
M

␰ = UT␭ = UT共1, . . . ,1兲T =i=1

M

Ui1, . . . , 兺 UiM
i=1

冊

TABLE I. Effects of the perturbation ␰ j in the potential energy
terms perturbation. Change in the weight of the potential energy
term Vi by ␰ j is determined by Eq. 共22兲. The probability change ofq by ␰ j is determined by Eq. 共13兲. The term “increase” or “decrease” indicates the probability change in the first order approximation of the cumulant generating function ␺共␭兲.

␰j ⬎ 0␰j ⬍ 0V共q兲 =

M

i=1

i=1

Uij ⬍ 0g j共q兲 ⬎ 0

g j共q兲 ⬍ 0

weakenstrengthen

strengthen
weakenincreasedecreasedecreaseincrease

Kb共rb − r̄b兲2 + 兺 Ka共␪a − ¯␪a兲2
兺b苸bondsa苸angles兺+

Kd关cos共nd␸d − ␥d兲 + 1兴d苸dihedrals

+

冉

A

B

ij

ij

冊

QQ

ij
ij
i j
+
.
兺
兺
12 −
r 6
共i,j兲苸nonbonds r
共i,j兲苸nonbonds ␧rij

共23兲Each term corresponds to the bond, angle, dihedral, van der
Waals, and electrostatic potential energy terms, respectively.
We can select all or some 共e.g., nonbonded interactions only兲
of the potential energy terms to perturb as冉

V共q兲 = 兵Kb共rb − r̄b兲2兩b 苸 some bonds其,T.

兵Ka共␪a − ¯␪a兲2兩a 苸 some angles其,兵Kd共cos共nd␸d − ␥d兲 + 1兲兩d 苸 some dihedrals其,再
再We standardize the direction of the jth eigenvector u j as
M
Uij 艌 0 for convenience, because the direction of the ei兺i=1
genvector does not change the result. If we perturb ␰ j direction, i.e., ␭ = ␰ ju j, then the perturbation term Eq. 共19兲 is represented asM

Uij ⬎ 0

⬎ 0, and g j共q兲 ⬎ 0 by Uij ⬍ 0. Therefore, the magnitude andthe sign of the eigenvectors provide important information in
the PEPCA. Moreover, we can investigate the detailed information of potential energy term, such as attractive or repulsive electrostatic interaction.
Next, we consider some concrete potential energy terms
to perturb. Currently, the most commonly used biomolecularpotential energy 关28兴 has a functional form

共21兲

− 兺 ␭iVi共q兲 = − 兺 ␰ jUijVi共q兲.

Probability of q

potential energy term Vi

共22兲

Therefore the perturbation ␰ j ⬎ 0 strengthens the ith potential
energy term Vi if Uij ⬍ 0, and it weakens the term if Uij ⬎ 0.
As discussed in Sec. II C, perturbation ␰ j ⬎ 0 increases theprobability of q if g j共q兲 ⬎ 0, and decreases it if g j共q兲 ⬍ 0, in
the first order approximation of ␺共␭兲. These results are summarized in Table I. Thus, we can identify the combination ofpotential energy terms stabilizing q, where g j共q兲 ⬍ 0 by Uij

Aij
rij12

−

Bij
rij6

冎

兩共i, j兲 苸 some nonbonds ,冎冊Q iQ j兩共i, j兲 苸 some nonbonds
␧rij

T

.

共24兲

III. NUMERICAL RESULTS

For demonstrating the power of the theory, we apply the
PEPCA to an alanine dipeptide molecule in vacuum. We performed 10 ns molecular dynamics simulation by integration
with a 2 fs time step. For generating a canonical ensemble at
300 K, Langevin dynamics with the collision frequency ␥= 1.0 ps−1 were used. We used molecular dynamics package
AMBER8 关29兴 to perform molecular dynamics simulation and046702-4

PERTURBATIONAL FORMULATION OF PRINCIPAL …

PHYSICAL REVIEW E 78, 046702 共2008兲

TABLE II. Index of the potential energy terms of the alanine
dipeptide. 89, 90, 101, and 102 are the improper torsion potential
energy terms. 236–276 van der Waals and 410–450 electrostatic
potential energy terms are 1–4 interactions, and they are divided by
2.0 and 1.2 共these are default values in AMBER8 关29兴兲, respectively.
Bond

Angle

Dihedralvan der Waals

Electrostatic

1–21

22–57

58–102

103–276

277–450

Eigenvalue

the ff03 force field 关30兴 which has the functional form given
in Eq. 共23兲. We saved the coordinates at every 1 ps, and used
a total of 10,000 coordinates for the analysis. This molecule
has two stable states in our 10 ns simulation. The transition
between the two states occurred many times, and it equilibrated well within the two states. For performing the
PEPCA, we developed a program that outputs the potential
energy terms from the topology and coordinates files. We use
all the potential energy terms for Eq. 共18兲, and index each ofthem as in Table II. As indicated in Eq. 共20兲, all potential
energy terms are divided by kBT for PCA. PCA was performed with the statistical computing environment R 关31兴.
Figure 1 shows the eigenvalues of the covariance matrixof the potential energy terms divided by kBT. The largest
eigenvalue ␴21 is remarkably larger than the other eigenvalues. The second order approximation of the Kullback-Leibler
divergence is represented as Eq. 共8兲. Therefore, the perturbation corresponding to the largest eigenvalue changes the conformational distribution remarkably, and the other perturbations have less effect on the conformational distribution.
Therefore, we concentrate on analyzing the largest perturbation ␰1. Figure 2共a兲 shows the sign of the first principal
component g1共q兲 on a Ramachandran plot, which clearly
characterizes the two stable states. As summarized in Table I,
the perturbation ␰1 ⬎ 0 increases the probability of q where

Index of the eigenvalue
FIG. 1. Eigenvalues of the covariance matrix of the potential
energy terms divided by kBT. Eigenvalues are sorted in descendingorder. The inset figure shows the eigenvalues in log scale. First fiveeigenvalues are 461.5, 16.3, 9.3, 5.7, and 5.3.

g1共q兲 ⬎ 0 关red in Fig. 2共a兲兴, and decreases the probability ofq, where g1共q兲 ⬍ 0 关blue in Fig. 2共a兲兴, in the first order approximation. Therefore, the perturbation ␰1 controls the ratioof the two stable states. As summarized in Table I, we canidentify the combination of interactions stabilizing the conformation q, where g1共q兲 ⬎ 0 关red in Fig. 2共a兲兴 as the negative components of the first eigenvector 共Ui1 ⬍ 0兲. In details,
according to Table I, we first notify that the probability of theconformation q, where g1共q兲 ⬎ 0 increases if ␰1 ⬎ 0. We thennotify that the potential energy term Vi is strengthened under␰1 ⬎ 0 if Ui1 ⬍ 0. We thus identify the interactions 共Ui1 ⬍ 0兲
which stabilize the conformation q where g1共q兲 ⬎ 0. Figure
2共b兲 shows that smaller negative components of the first eigenvector are el-6-17, el-5-18, and el-8-16. According toTable III, el-8-16 is an attractive interaction while el-6-17
and el-5-18 are repulsive interactions. Conformation A 关Figs.
2共a兲 and 2共b兲兴 ascertains that the attractive interaction el-8-16
is an important stabilizing interaction. Interestingly, the conformation A also indicates that the extended ␤-like structureis enhanced by repulsive interactions el-6-17 and el-5-18.
Similarly, the positive components of the first eigenvector
共Ui1 ⬎ 0兲 stabilizes the conformation q, where g1共q兲 ⬍ 0
关blue in Fig. 2共a兲兴. The larger positive component el-6-18 is
an attractive interaction 关Figure 2共b兲 and Table III兴. Conformation B 关Figs. 2共a兲 and 2共b兲兴 ascertains that the attractive
interaction el-6-18 is an important stabilizing interaction.
Thus, the perturbation ␰1 controls the ratio of the two stablestates, i.e., the equilibrium constant of the two stable states.
The negative and the positive components of the first eigenvector clearly identifies the stabilizing interactions of the twostable states.
Next, we analyze the second largest perturbation ␰2, although its effect is smaller compared with that of the largest
perturbation ␰1. Figure 2共c兲 shows the sign of the second
principal component g2共q兲. Interestingly, g2共q兲 clearly characterizes the transition state region of the two stable states.
Therefore, the perturbation ␰2 controls the ratio between thetransition state and the two stable states. The smaller negative components of the second eigenvector 共Ui2 ⬍ 0兲, which
stabilize the conformation q where g2共q兲 ⬎ 0 关red in Fig.
2共c兲兴, are attractive interactions el-5-16 and el-6-15 关Figure
2共d兲 and Table III兴. Conformation TS 关Figs. 2共c兲 and 2共d兲兴
ascertains that attractive interactions el-5-16 and el-6-15 areimportant stabilizing interactions. On the other hand, thelarge positive component of the second eigenvector 共Ui2
⬎ 0兲, which stabilizes the conformation q where g2共q兲 ⬍ 0
关blue in Fig. 2共c兲兴, is a repulsive interaction el-6-16 关Figure
2共d兲 and Table III兴. The conformation TS indicates that therepulsive interaction el-6-16 destabilizes the transition state.
Therefore, el-6-16 relatively stabilizes the two stable states关blue in Fig. 2共c兲兴 by destabilizing the transition state 关red inFig. 2共c兲兴. Thus, the second largest perturbation ␰2 controls
the stability of the transition state of the two stable states,i.e., the rate constant of the transition between the two stablestates. The negative and the positive components of the second eigenvector identify the stabilizing interactions of thetransition state and the two stable states.
IV. DISCUSSION AND CONCLUSIONSIn this study, we clarified that within the second orderapproximation of the Kullback-Leibler divergence, a PCA

046702-5

PHYSICAL REVIEW E 78, 046702 共2008兲

KOYAMA et al.

(a)

(b)
Component of the 1st eigenvector

A

B
9
7 15
5 17

el-6-18

B

6 18

A
6
5

18
17
8 16

el-5-18

el-8-16el-6-17

Index of the potential energy term(c)Component of the 2nd eigenvector(d)TS

el-6-16

TS
15
16
5
6

el-5-16 el-6-15

Index of the potential energy term
FIG. 2. 共Color兲 The sign of 共a兲 the first 关g1共q兲兴 and 共c兲 the second 关g2共q兲兴 principal component on a Ramachandran plot. Red and blue
points show gi共q兲 艌 0 and gi共q兲 ⬍ 0, respectively. The dihedral angles ␾ and ␺ are defined by atom indices 5-7-9-15 and 7-9-15-17,
respectively. The 3D structures labeled A, B, and TS show snapshots of the MD simulation indicated by the black lines. The 3D structures
were rendered by VMD 关32兴. Cyan, white, red, and blue spheres indicate the C, H, O, and N atoms, respectively. The components of 共b兲 thefirst and 共d兲 the second eigenvectors. Label el-6-18 indicates the component of the electrostatic potential energy term between atom 6 andatom 18. Such labels are captioned where the absolute value of the component is larger than 0.3.

can be interpreted as a method for 共1兲 decomposing a given
perturbation into perturbations that independently contributeto the conformational distribution change or 共2兲 successively
finding the perturbation that induces the largest conformaTABLE III. Charge 共electron charge units兲 of the numbered atoms in Figs. 2共b兲 and 2共d兲.
5共C兲

6共O兲

8共H兲

15共C兲

16共O兲

17共N兲

18共H兲

0.512

−0.550

0.294

0.570

−0.555−0.4240.290

tional distribution change. To achieve general applicability
and provide clear physical meaning of the perturbation, we
introduced a perturbation of potential energy terms, whichleads to the PEPCA. The effectiveness of the PEPCA wasdemonstrated using alanine dipeptide as a minimal model ofa nonsingle dominant state behavior biomolecule.
To apply the PEPCA to larger molecules, we must investigate some practical aspects that were not addressed in thispaper. First, water molecules are not incorporated in our MD
simulation. For more practical simulation, it is important to
explicitly incorporate water molecules in MD simulation.

046702-6

PERTURBATIONAL FORMULATION OF PRINCIPAL …

PHYSICAL REVIEW E 78, 046702 共2008兲

Here we represent the biomolecular coordinates q and thewater coordinates q⬘. If we only perturb the biomolecules,
i.e., perturbation functions depend only on q, then the perturbed marginal distribution of q is␳␭共q兲 =

冕␳␭共q,q⬘兲dq⬘ =冕

e␭·f共q兲−␺共␭兲␳共q,q⬘兲dq⬘= e␭·f共q兲−␺共␭兲␳共q兲.

共25兲This is the same distribution as Eq. 共3兲. Therefore, our perturbational formulation of PCA is valid without any modification in explicit water simulation if we only perturb the
biomolecular coordinates.
Second, sufficient conformational sampling is required
because insufficient conformational sampling results in poor
estimations of the principal components. This problem has
been intensively challenged at least by three different approaches. The first is employment of more efficient samplingalgorithms such as extended ensemble algorithms 关33,34兴,the second is the development of the algorithm for improving
the parallel scalability of the cluster computer for MD simulation 关35,36兴, and the last is design of fast special-purpose
hardware for MD simulation 关37,38兴. By combining theseapproaches, the insufficient conformational sampling problem may be solved in the near future. At that time, analysis
method of the MD simulation data will become more important.
Third, computational feasibility of diagonalization of thecovariance matrix is a problem in PCA. The atom number Nof an ordinary protein is 103 – 104. If we perturb all theatomic potential energy terms in the protein molecule, thenthe size of the covariance matrix is N2 ⬃ 106 – 108. Such alarge matrix is difficult to diagonalize with the existing computational power. Instead of perturbing each atomic potential
energy term, we can perturb some potential energy terms at atime. For example, we can consider groups of atoms such asamino acids or bases, then perturb the intergroup nonbonded
potential energy terms simultaneously. If we consider aminoacid as one such group, then the group size, i.e., the aminoacid number Na is 102 – 103 and the size of the covariancematrix is Na2 ⬃ 104 – 106. Thus the diagonalization is feasible.
In this analysis, we can identify the important combinationsof the intergroup interactions contributing to conformational
distribution change. Note that Kong et al. 关39兴 used residuebased interaction-correlation matrix to study the signaltransduction mechanism of rhodopsin. In this respect, our
PEPCA method is considered to be a covariance matrix
counterpart of their interaction correlation matrix analysis.
As another way to reduce the matrix size, we can use theGram matrix diagonalization instead of the covariance matrix. This is because the formal equivalence between the
PCA using perturbation functions and the kernel PCA 关40兴
holds as shown in Appendix C. Since the Gram matrix depends only on the sample data size, the Gram matrix diagonalization may be useful when the number of the perturbation functions is larger than the sample data size. Thisequivalence between the PCA using perturbation functions
and the kernel PCA may also provide the physical meaning
to some dimension reduction methods. It is known that sev-

eral dimension reduction methods, such as metric multidimensional scaling 共metric MDS 关41兴, also known as principal coordinate analysis 关3,42兴兲, Isomap 关43兴, Laplacian
eigenmap 关44兴, and locally linear embedding 共LLE兲 关45兴 areinterpreted as kernel PCA 关46,47兴. Therefore, it may be interesting to interpret these dimensional reduction methods
with our perturbational formulation, because the perturbation
may have a clearer physical meaning than the kernel function. We also note that our perturbational formulation of
PCA provides the clearer limitation that is defined by thesecond order approximation of the Kullback-Leibler divergence in comparison with the kernel PCA.
Finally, we estimate the Kullback-Leibler divergence by
the perturbation from the unperturbed equilibrium state, i.e.,a perturbation estimator only reweights the samples from theunperturbed distribution. Therefore, if the perturbation
changes the conformational distribution into unsampled conformations in unperturbed equilibrium state, then the perturbation estimation is invalid. Recently, it has been observedthat some proteins “sample” the functional conformations
under the unperturbed conformational fluctuation, which is
known as a “pre-existing” behavior 关2,48兴. As demonstrated
in Sec. III, our perturbational formulation of PCA is well
suited for pre-existing behavior. If a molecular function follows the pre-existing behavior, then the principal components identify the functional substates, and their detailed
atomic interactions can be obtained from their eigenvectors.
From the components of the eigenvectors, we may be able toidentify the natural perturbation targets for molecular functions, such as the ligand binding sites, which induce the large
conformational distribution change 关16兴. Thus, while weneed more investigations for larger molecules, PCA can beapplied, by carefully selecting the perturbation functions, notonly to identify the molecular conformational fluctuation but
also to predict the conformational distribution change by the
perturbation beyond the limitation of the previous methods.
ACKNOWLEDGMENTS

We thank Mitsunori Takano, Hironori Nakamura, Makoto
Taiji, and Noriaki Okimoto for valuable discussions. We arealso thankful for the computational resources of the RIKENSuper Combined Cluster 共RSCC兲. This research was supported by the Center for Developmental Biology 共CDB兲
共H.R.U.兲 and by the Japan Society for the Promotion of Science for Young Scientists 17-4169 共T.J.K兲.
APPENDIX A: EXPECTATION PARAMETERREPRESENTATION OF THE KULLBACKLEIBLER DIVERGENCE

Equation 共6兲 indicates that D共␳␭ 储 ␳兲 is the Legendre transformation of ␺共␭兲. Therefore, 具f典␭ and ␭ have a one-to-one
correspondence. As a result, we can use 具f典␭ for specifying
the perturbed distribution ␳␭共q兲 instead of ␭. In this respect,
具f典␭ are called expectation parameters 关20兴. Because of
Cramér’s theorem in large deviation theory 关23,26兴, the Legendre transformation of the cumulant generating function
␺共␭兲 is also a rate function ␾共具f典␭兲. Therefore, the equality

046702-7

PHYSICAL REVIEW E 78, 046702 共2008兲

KOYAMA et al.













































































































FIG. 3. Hyperplanes determining the probability change of q by
the perturbation ␭ and contours of the Kullback-Leibler divergence.
The solid line shows the hyperplane ␭ · y − ␺共␭兲 = 0 which separates
the increases or decreases of the probability change of q by the
perturbation ␭. The solid and broken ellipses show the contours of
the Kullback-Leibler divergence of ␾共y兲 = C共␳ , ␳␭兲 and ␾共y兲
= D共␳␭ 储 ␳兲, respectively.

␾共具f典␭兲 = D共␳␭ 储 ␳兲

共A1兲

holds. This relation between the rate function and the
Kullback-Leibler divergence is a direct consequence of thecontraction principle in large deviation theory 关23,26兴, whichstates that the rate function is obtained by minimization of
the Kullback-Leibler divergence D共␳⬘ 储 ␳兲 under the probability distribution ␳⬘共q兲 which has given averages of statistics
f共q兲. Since the exponential family 关Eq. 共3兲兴 is the solution ofsuch minimization, the above equality holds. In this interpretation, ␭ indicate the Lagrange multipliers of the constraints
for the averages.
Equation 共3兲 indicates that the probability change of q by
the perturbation ␭ is determined by ␭ · f共q兲 − ␺共␭兲, and ␭ · y− ␺共␭兲 = 0 defines the hyperplane in which the probability ofq does not change by the perturbation ␭. We can show thatthe inequality
␭ · 具f典 艋 ␺共␭兲 艋 ␭ · 具f典␭tours of the Kullback-Leibler divergence is summarized inFig. 3.
Within the second order approximation of the cumulant
generating function ␺共␭兲, ␭ can be expressed by 具f典␭ as␭ ⬇ C−1共具f典␭ − 具f典兲.

共A5兲

With this, Eq. 共7兲 can be represented as
1
␾共具f典␭兲 = D共␳␭ 储 ␳兲 ⬇ 共具f典␭ − 具f典兲TC−1共具f典␭ − 具f典兲.
2

共A2兲

holds. The left inequality is from Eq. 共4兲 with Jensen’s inequality. The right inequality follows Eq. 共6兲 with the nonnegativity of the Kullback-Leibler divergence 关Eq. 共5兲兴. Notethat this inequality is considered to be a multidimensional
perturbation version of the Gibbs-Bogoliubov inequality
关49兴. The inequality 共A2兲 shows that the hyperplane ␭ · y− ␺共␭兲 = 0 is located between ␭ · 共y − 具f典兲 = 0 and ␭ · 共y − 具f典␭兲= 0. Because ␭ is a gradient of ␾共y兲 at y = 具f典␭, ␭ · 共y − 具f典␭兲= 0 is the tangent hyperplane of the contour of ␾共y兲 at y
= 具f典␭. The inequality 共A2兲 also indicates that for a certain0 艋 ␣共␭兲 艋 1, the equality␺共␭兲 = ␭ · 具f典␣共␭兲␭FIG. 4. Contours of the Kullback-Leibler divergence and eigenvectors. Within the second order approximation of the KullbackLeibler divergence, the contours of the Kullback-Leibler divergence
make ellipsoids and eigenvectors corresponding to the principal
axes of ellipsoids 关Eqs. 共8兲, 共A6兲, and 共A8兲兴. 共a兲 ␭ 共natural parameters in the exponential family兲 representation of the KullbackLeibler divergence. The Kullback-Leibler divergence in the second
order approximation with ␭ representation is Eq. 共8兲. Axes ␰1 and ␰2
correspond to the first and second eigenvectors, respectively. ␭ = 0
corresponds to the unperturbed distribution. 共b兲 具f典␭ 共expectation
parameters in the exponential family兲 representation of the
Kullback-Leibler divergence, i.e., rate function ␾共具f典␭兲 关Eq. 共A1兲兴.
The second order approximation of the Kullback-Leibler divergencewith 具f典␭ representation is Eq. 共A6兲. Axes z1 and z2 correspond to
the first and second principal components, respectively 关Eq. 共A8兲兴.
y = 具f典 corresponds to the unperturbed distribution.

共A6兲
Note that冑共具f典␭ − 具f典兲TC−1共具f典␭ − 具f典兲共A7兲

is known as the Mahalanobis distance 关3,41兴. With Eq. 共14兲,
Eq. 共A6兲 can be represented asM

具gi典 2
1
1
␾共具f典␭兲 ⬇ 具g典␰TUTC−1U具g典␰ ⬇ 兺 2␰ .
2
2 i=1 ␴i

共A8兲

共A3兲

With Eqs. 共8兲, 共A6兲, and 共A8兲, contours of the KullbackLeibler divergence with ␭ and 具f典␭ representation are shownin Fig. 4.

共A4兲

APPENDIX B: PERTURBATIONAL FORMULATIONOF CARTESIAN PCA AND DPCA

holds. We can show that this equality leads to
D共␳␣共␭兲␭ 储 ␳兲 = D共␳␣共␭兲␭ 储 ␳␭兲 ⬅ C共␳, ␳␭兲,

where C共␳ , ␳␭兲 is known as the Chernoff information
关24,50,51兴. Therefore, ␭ · y − ␺共␭兲 = ␭ · 共y − 具f典␣共␭兲␭兲 = 0 is thetangent hyperplane of the contour of ␾共y兲 = C共␳ , ␳␭兲 at y
= 具f典␣共␭兲␭. The relationship of the hyperplanes and the con-

We consider perturbational formulations of Cartesian
PCA and dPCA. In Cartesian PCA, PCA is performed in
Cartesian coordinate space such as the NC C␣ atoms of theprotein molecule046702-8

PERTURBATIONAL FORMULATION OF PRINCIPAL …

f共q兲 = 共q1, . . . ,q3NC兲T .

PHYSICAL REVIEW E 78, 046702 共2008兲

共B1兲

Nd

共j兲
k共j兲
兺
l kl cos共␥l − ␥l兲 = 0,
l=1

Then, the perturbation term in Eq. 共2兲 is represented as
M

3NC− kBT 兺 ␭i f i共q兲 = − kBT 兺 ␭iqi .
i=1

共B2兲

i=1

This means that the coordinate qi is perturbed by the constantforce kBT␭i. Therefore, the Cartesian PCA can be interpreted
as a method for 共1兲 decomposing a given constant forces
perturbation 关Eq. 共B2兲兴 into perturbations that independently
contribute to the conformational distribution change or 共2兲
successively finding the constant forces perturbation 关Eq.
共B2兲兴 that induces the largest conformational distributionchange under the constraints Eqs. 共10兲. The former interpretation is already known 关14兴. However, it is always necessary
to pay attention to the frame fitting problem when using theCartesian PCA.
To avoid the frame fitting problem, dihedral angles
␸1 , . . . , ␸Nd were used in PCA 关52兴. However, the direct useof the dihedral angle is problematic because of the periodicity of the dihedral angles. dPCA 关12,13兴 solves this periodicity problem by performing PCA in the mapped dihedral
angle spacef共q兲 = 共cos ␸1,sin ␸1, . . . ,cos ␸Nd,sin ␸Nd兲T

共B3兲

instead of the dihedral angle space. Then, the perturbation
term in Eq. 共2兲 can be represented as
M

Nd

i=1

i=1

− kBT 兺 ␭i f i共q兲 = − kBT 兺 共␭2i−1 cos ␸i + ␭2i sin ␸i兲j = 1, . . . ,i − 1,

共B7b兲

in terms of k and ␥. Thus, the dPCA can be interpreted as amethod for 共1兲 decomposing a given dihedral angles perturbation 关Eq. 共B4兲兴 into perturbations that independently contribute to the conformational distribution change or 共2兲 successively finding the dihedral angles perturbation 关Eq. 共B4兲兴
that induces the largest conformational distribution change
under the constraints Eqs. 共B7兲. Above derivation of the perturbational formulation of dPCA also demonstrates the general procedure for inventing a PCA which is free from theframe fitting problem as follows. First, introduce a physicallymeaningful perturbation which perturbing some internal coordinates. Second, divide the perturbation as the perturbation
term in Eq. 共2兲. Finally, perform PCA using f共q兲. Introduction of the PEPCA in Sec. II D also follows this procedure.
APPENDIX C: RELATIONSHIPTO THE KERNEL PCA

We noticed a formal equivalence between the PCA using
the perturbation functions developed in Sec. II C and thekernel PCA 关40兴, because each methods can be interpreted as
the PCA in a mapped space f共q兲 共this space is called feature
space in the kernel method兲 instead of q. To show thisequivalence, we briefly introduce the kernel PCA in our notation. We consider n point data 共q1 , . . . , qn兲, and their centered mapped data matrixD̃ ⬅ 关f̃共q1兲, . . . ,f̃共qn兲兴,

共C1兲

Nd

= − kBT 兺 ki cos共␸i − ␥i兲.

共B4兲where centered mapped data is

To clarify the physical meaning of the perturbation, we introduced k = 共k1 , . . . , kNd兲T and ␥ = 共␥1 , . . . , ␥Nd兲T as polar coordinates of ␭ such that 共␭2i−1 , ␭2i兲T = 共ki cos ␥i , ki sin ␥i兲T.
Since ␭ and 共k , ␥兲 show one-to-one correspondence, theconstraints Eqs. 共10兲 can also be represented by k and ␥
instead of ␭. For this purpose, we calculatef̃共qi兲 ⬅ f共qi兲 −

i=1

n

r

D̃ = 兺 冑n␴iuivTi ,

共B5兲

i=1

and if we represent the jth eigenvector u j by k共j兲 and ␥共j兲,
then
Nd

共j兲
共j兲
共j兲
u j · ␭ = 兺 关共k共j兲
l cos ␥l 兲共kl cos ␥l兲 + 共kl sin ␥l 兲共kl sin ␥l兲兴

=兺

共j兲
k共j兲
l kl cos共␥l − ␥l兲.

where r, 冑n␴i, U = 共u1 , . . . , ur兲, and V = 共v1 , . . . , vr兲 are therank of D̃, singular values, left singular vectors, and rightsingular vectors, respectively. In this representation, the covariance matrix of the mapped data isrWith these expressions, the constraints Eqs. 共10兲 are represented as兩k兩 = ␦,1C ⬅ D̃D̃T = 兺 ␴2i uiuTi .
n
i=1

共B6兲

l=1

␦ such that the second orderapproximation of D共␳k,␥储␳兲 is valid, 共B7a兲

共C3兲

i=1

l=1
Nd

共C2兲

We represent the centered mapped data matrix D̃ by singularvalue decomposition 关27兴 asNd兩␭兩2 = 兺 关共ki cos ␥i兲2 + 共ki sin ␥i兲2兴 = 兩k兩2

1
兺 f共ql兲.
n l=1

共C4兲

C is an M ⫻ M matrix, whose size depends on the dimension
M of the mapped space. The above equation indicates thatthe eigenvalues and the eigenvectors of C give ␴2i and ui,
respectively. In our PCA using perturbation functions, thisdiagonalization is performed. On the other hand, in the kernel PCA, the centered Gram matrix G̃

046702-9

PHYSICAL REVIEW E 78, 046702 共2008兲

KOYAMA et al.
rnG̃ ⬅ D̃TD̃ = 兺 n␴2i vivTi

共C5兲

i=1

uj =

is diagonalized instead of the covariance matrix C. G̃ is an ⫻ n matrix, whose size depends on the data size n. With
Eqs. 共C1兲, 共C2兲, and 共C5兲, the element of the centered Grammatrix G̃ij can be represented as

g j共q兲 ⬅ u j · f̃共q兲 =

n

n

共C6兲

where we introduced the kernel function k共q , q⬘兲
k共q,q⬘兲 ⬅ f共q兲 · f共q⬘兲

共C7兲

and the centered kernel function k̃共q , q⬘兲k̃共q,q⬘兲 ⬅ f̃共q兲 · f̃共q⬘兲.

共C11兲

.

n
Vijf̃共qi兲 · f̃共q兲
兺i=1

冑n ␴ j=n
Vijk̃共qi,q兲
兺i=1

冑n ␴ j.

Thus, in the kernel PCA, the eigenvalues ␴2i and the principal

components g共q兲 can be calculated by the diagonalization of
the centered Gram matrix G̃, which can be evaluated by the
kernel function. On the other hand, Eq. 共C11兲 indicates thatthe eigenvectors ui require the explicit functional form off共q兲.
For example, a homogeneous polynomial kernel of degree
two, i.e., c = 0 and d = 2 in Eq. 共C9兲, is represented as
3NC

k共q,q⬘兲 = 共q · q⬘兲 = 兺 q2i qi⬘2 + 2

共C8兲

Thus, the centered Gram matrix G̃ 关and also centered kernel
function k̃共q , q⬘兲兴 can be calculated using the kernel function
k共q , q⬘兲, and the explicit functional form of f共q兲 is not required. As a first selection of the kernel function, the polynomial kernel
k共q,q⬘兲 = 共q · q⬘ + c兲d冑n ␴ j

共C12兲

n

1
1
k共ql,q j兲 + 2 兺 k共ql,qm兲,兺n l=1
n l,m=1

冑n ␴ j =With this expression, the jth principal component g j共q兲 isrepresented as

1G̃ij = f̃共qi兲 · f̃共q j兲 = k̃共qi,q j兲 = k共qi,q j兲 − 兺 k共qi,qm兲
n m=1
−

D̃v j

Vijf̃共qi兲
兺
i=1

共C9兲

2

i=1

= f共q兲 · f共q⬘兲,兺

qiq jqi⬘q⬘j

1艋i⬍j艋3NC

共C13兲

where
f共q兲 = 共兵q2i 兩i = 1, . . . ,3NC其,兵冑2qiq j兩1 艋 i ⬍ j 艋 3NC其兲T .
共C14兲

are often used in the kernel method. In particular, the Gaussian kernel has infinite dimensional feature space. We can see
from Eq. 共C5兲, the eigenvalues and the eigenvectors of thecentered Gram matrix G̃ are n␴2i and vi, respectively. Therefore, by dividing eigenvalues n␴2i of G̃ by data size n, we
can obtain the eigenvalue ␴2i of C. The eigenvectors to be
determined are not vi but ui. From Eq. 共C3兲, u j can be represented as

In our perturbational formulation of PCA, this perturbation
functions f共q兲 give the harmonic perturbation around the origin. The polynomial kernel 关Eq. 共C9兲兴 and the Gaussian kernel 关Eq. 共C10兲兴 in Cartesian coordinates were applied in the
MD simulation 关53兴. However, using Cartesian coordinates
leads to the frame fitting problem, as in the case of the Cartesian PCA. It is desirable to use the invariant kernel function with respect to the overall molecular motion. For thispurpose, it may be useful to perform the kernel PCA with akernel function such as the polynomial kernel 关Eq. 共C9兲兴 andthe Gaussian kernel 关Eq. 共C10兲兴 with respect to the mappeddihedral angle space 关Eq. 共B3兲兴 or potential energy term
space 关Eq. 共20兲兴.

关1兴 S. Hammes-Schiffer and S. J. Benkovic, Annu. Rev. Biochem.
75, 519 共2006兲.
关2兴 K. Henzler-Wildman and D. Kern, Nature 共London兲 450, 964
共2007兲.
关3兴 I. T. Jolliffe, Principal Component Analysis, 2nd ed. 共Springer,
New York, 2002兲.
关4兴 A. Kitao, F. Hirata, and N. Go, Chem. Phys. 158, 447 共1991兲.
关5兴 A. E. García, Phys. Rev. Lett. 68, 2696 共1992兲.
关6兴 A. Amadei, A. B. M. Linssen, and H. J. C. Berendsen, Proteins
17, 412 共1993兲.
关7兴 A. Kitao and N. Go, Curr. Opin. Struct. Biol. 9, 164 共1999兲.
关8兴 H. J. C. Berendsen and S. Hayward, Curr. Opin. Struct. Biol.

10, 165 共2000兲.
关9兴 P. H. Hünenberger, A. E. Mark, and W. F. van Gunsteren, J.
Mol. Biol. 252, 492 共1995兲.
关10兴 M. Karplus and T. Ichiye, J. Mol. Biol. 263, 120 共1996兲.
关11兴 R. Abseher and M. Nilges, J. Mol. Biol. 279, 911 共1998兲.
关12兴 Y. Mu, P. H. Nguyen, and G. Stock, Proteins 58, 45 共2005兲.
关13兴 A. Altis, P. H. Nguyen, R. Hegger, and G. Stock, J. Chem.
Phys. 126, 244111 共2007兲.
关14兴 M. Ikeguchi, J. Ueno, M. Sato, and A. Kidera, Phys. Rev. Lett.
94, 078102 共2005兲.
关15兴 D. Ming and M. E. Wall, Proteins 59, 697 共2005兲.
关16兴 D. Ming and M. E. Wall, J. Mol. Biol. 358, 213 共2006兲.

and the Gaussian kernel
k共q,q⬘兲 = e−兩q − q⬘兩 /2␴ ,
2

2

共C10兲

046702-10

PERTURBATIONAL FORMULATION OF PRINCIPAL …

PHYSICAL REVIEW E 78, 046702 共2008兲

关17兴 S. Kullback, Information Theory and Statistics 共John Wiley &Sons, New York, 1959兲.
关18兴 O. Barndorff-Nielsen, Information and Exponential Families
In Statistical Theory 共John Wiley & Sons, Chichester, 1978兲.
关19兴 L. D. Brown, Fundamentals of Statistical Exponential Families: With Applications in Statistical Decision Theory 共Instituteof Mathematical Statistics, Hayward, CA, 1986兲.
关20兴 S. Amari and H. Nagaoka, Methods of Information Geometry
共AMS and Oxford University Press, New York, 2000兲.
关21兴 H.-O. Georgii, in Entropy, edited by A. Greven, G. Keller, andG. Warnecke 共Princeton University Press, Princeton, 2003兲, pp.
37–54.
关22兴 S. Kullback and R. A. Leibler, Ann. Math. Stat. 22, 79 共1951兲.
关23兴 S. R. S. Varadhan, in Entropy 关21兴, pp. 199–214.
关24兴 T. M. Cover and J. A. Thomas, Elements of Information
Theory, 2nd ed. 共Wiley-Interscience, New York, 2006兲.
关25兴 I. N. Sanov, Mat. Sb. 42, 11 共1957兲.
关26兴 A. Dembo and O. Zeitouni, Large Deviations Techniques andApplications, 2nd ed. 共Springer, New York, 1998兲.
关27兴 R. A. Horn and C. R. Johnson, Matrix Analysis 共Cambridge
University Press, Cambridge, 1985兲.
关28兴 J. W. Ponder and D. A. Case, Adv. Protein Chem. 66, 27
共2003兲.
关29兴 D. A. Case, T. A. Darden, T. E. Cheatham, III, C. L. Simmerling, J. Wang, R. E. Duke, R. Luo, K. M. Merz, B. Wang, D. A.
Pearlman, M. Crowley, S. Brozell, V. Tsui, H. Gohlke, J. Mongan, V. Hornak, G. Cui, P. Beroza, C. Schafmeister, J. W. Caldwell, W. S. Ross, and P. A. Kollman, AMBER8, University ofCalifornia, San Francisco, 2004.
关30兴 Y. Duan, C. Wu, S. Chowdhury, M. C. Lee, G. Xiong, W.
Zhang, R. Yang, P. Cieplak, R. Luo, T. Lee, J. Caldwell, J.
Wang, and P. Kollman, J. Comput. Chem. 24, 1999 共2003兲.
关31兴 http://www.r-project.org/
关32兴 W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graphics
14, 33 共1996兲.
关33兴 Y. Iba, Int. J. Mod. Phys. C 12, 623 共2001兲.
关34兴 Y. Okamoto, J. Mol. Graphics Modell. 22, 425 共2004兲.
关35兴 J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid,
E. Villa, C. Chipot, R. D. Skeel, L. Kale, and K. Schulten, J.
Comput. Chem. 26, 1781 共2005兲.
关36兴 K. J. Bowers, E. Chow, H. Xu, R. O. Dror, M. P. Eastwood, B.
A. Gregersen, J. L. Klepeis, I. Kolossváry, M. A. Moraes, F. D.

Sacerdoti, J. K. Salmon, Y. Shan, and D. E. Shaw, in Proceedings of the ACM/IEEE Conference on Supercomputing 共SC06兲
共ACM Press, New York, 2006兲.
关37兴 M. Taiji, T. Narumi, Y. Ohno, N. Futatsugi, A. Suenaga, N.
Takada, and A. Konagaya, in Proceedings of the ACM/IEEE
Conference on Supercomputing 共SC03兲 共ACM Press, New
York, 2003兲.
关38兴 D. E. Shaw, M. M. Deneroff, R. O. Dror, J. S. Kuskin, R. H.
Larson, J. K. Salmon, C. Young, B. Batson, K. J. Bowers, J. C.
Chao, M. P. Eastwood, J. Gagliardo, J. P. Grossman, C. R. Ho,
D. J. Ierardi, I. Kolossváry, J. L. Klepeis, T. Layman, C.
McLeavey, M. A. Moraes, R. Mueller, E. C. Priest, Y. Shan, J.
Spengler, M. Theobald, B. Towles, and S. C. Wang, Commun.
ACM 51共7兲, 91 共2008兲.
关39兴 Y. Kong and M. Karplus, Structure 共London兲 15, 611 共2007兲.
关40兴 B. Schölkopf, A. Smola, and K. Müller, Neural Comput. 10,1299 共1998兲.
关41兴 K. V. Mardia, J. T. Kent, and J. M. Bibby, Multivariate Analysis 共Academic Press, London, 1979兲.
关42兴 J. C. Gower, Biometrika 53, 325 共1966兲.
关43兴 J. B. Tenenbaum, V. de Silva, and J. C. Langford, Science 290,
2319 共2000兲.
关44兴 M. Belkin and P. Niyogi, Neural Comput. 15, 1373 共2003兲.
关45兴 S. T. Roweis and L. K. Saul, Science 290, 2323 共2000兲.
关46兴 C. K. I. Williams, Mach. Learn. 46, 11 共2002兲.
关47兴 J. Ham, D. D. Lee, S. Mika, and B. Schölkopf, in Proceedingsof the Twenty-First International Conference on Machine
Learning 共ACM Press, New York, 2004兲, pp. 369–376.
关48兴 K. Gunasekaran, B. Ma, and R. Nussinov, Proteins 57, 433
共2004兲.
关49兴 C. Chipot and A. Pohorille, in Free Energy Calculations:
Theory and Applications in Chemistry and Biology, edited by
C. Chipot and A. Pohorille 共Springer, Berlin, 2007兲 pp. 33–75.
关50兴 H. Chernoff, Ann. Math. Stat. 23, 493 共1952兲.
关51兴 S. Sinanović and D. H. Johnson, Signal Process. 87, 1326
共2007兲.
关52兴 D. M. F. van Aalten, B. L. de Groot, J. B. C. Findlay, H. J. C.
Berendsen, and A. Amadei, J. Comput. Chem. 18, 169 共1997兲.
关53兴 X.-w. Chen and J. Chen, in Proceedings of the IEEE International Joint Conference on Neural Networks 共IJCNN’05兲
共IEEE Press, New York, 2005兲, Vol. 1, pp. 521–526.

046702-11