PHYSICAL REVIEW E 84, 026704 (2011)

Perturbation analyses of intermolecular interactions
Yohei M. Koyama,1,* Tetsuya J. Kobayashi,2 and Hiroki R. Ueda1,3,4,5,6,†
1

Laboratory for Synthetic Biology, Quantitative Biology Center, RIKEN, 2-2-3 Minatojima-minamimachi,
Chuo-ku, Kobe, Hyogo 650-0047, Japan
2
Institute of Industrial Science, The University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 153-8505, Japan
3
Laboratory for Systems Biology, Center for Developmental Biology, RIKEN, 2-2-3 Minatojima-minamimachi,
Chuo-ku, Kobe, Hyogo 650-0047, 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
6
Department of Mathematics, Kyoto University, Kyoto 606-8502, Japan(Received 5 November 2010; revised manuscript received 19 June 2011; published 12 August 2011)
Conformational fluctuations of a protein molecule are important to its function, and it is known thatenvironmental molecules, such as water molecules, ions, and ligand molecules, significantly affect the function
by changing the conformational fluctuations. However, it is difficult to systematically understand the role ofenvironmental molecules because intermolecular interactions related to the conformational fluctuations arecomplicated. To identify important intermolecular interactions with regard to the conformational fluctuations, wedevelop herein (i) distance-independent and (ii) distance-dependent perturbation analyses of the intermolecularinteractions. We show that these perturbation analyses can be realized by performing (i) a principal component
analysis using conditional expectations of truncated and shifted intermolecular potential energy terms and(ii) a functional principal component analysis using products of intermolecular forces and conditional cumulative
densities. We refer to these analyses as intermolecular perturbation analysis (IPA) and distance-dependentintermolecular perturbation analysis (DIPA), respectively. For comparison of the IPA and the DIPA, we apply
them to the alanine dipeptide isomerization in explicit water. Although the first IPA principal components
discriminate two states (the α state and PPII (polyproline II) + β states) for larger cutoff length, the separation
between the PPII state and the β state is unclear in the second IPA principal components. On the other hand, in thelarge cutoff value, DIPA eigenvalues converge faster than that for IPA and the top two DIPA principal components
clearly identify the three states. By using the DIPA biplot, the contributions of the dipeptide-water interactions toeach state are analyzed systematically. Since the DIPA improves the state identification and the convergence ratewith retaining distance information, we conclude that the DIPA is a more practical method compared with the IPA.
To test the feasibility of the DIPA for larger molecules, we apply the DIPA to the ten-residue chignolin folding in
explicit water. The top three principal components identify the four states (native state, two misfolded states, andunfolded state) and their corresponding eigenfunctions identify important chignolin-water interactions to eachstate. Thus, the DIPA provides the practical method to identify conformational states and their corresponding
important intermolecular interactions with distance information.
DOI: 10.1103/PhysRevE.84.026704

PACS number(s): 02.70.Rr, 36.20.Ey, 87.15.ap, 89.70.Cf

I. INTRODUCTION

Conformational fluctuations of a protein molecule areimportant to its function. Environmental molecules such as
solvents [1,2] and ligands [3–6] significantly affect proteinfunction by modifying its conformational fluctuations. However, it has been difficult to systematically analyze suchenvironmental effects on the conformational fluctuation ofproteins. Several methods have been developed to analyze theconformational fluctuations of a protein. For example, correlation analysis using atomic coordinates [7,8] and inter-residue
interaction energy [9,10] reveal communication within theprotein. The principal component analysis (PCA) [11] using
atomic coordinates [12–14], atomic pair distances [15], andmapped dihedral angles [16,17] decomposes the fluctuations
into large uncorrelated fluctuations. However, these methods*†

ym.koyama@gmail.com
uedah-tky@umin.ac.jp1539-3755/2011/84(2)/026704(20)are based only on the analysis of protein coordinates, so
direct information regarding the environmental molecules islost. Therefore, although the contributions of environmental
molecules might be implicitly represented in conformational
fluctuations of a protein, it is difficult to explicitly evaluate
their direct contribution to the conformational fluctuations ofa protein.
To understand the molecular conformational fluctuations
based on their atomic interactions, we previously introduced
the potential energy PCA (PEPCA) [18] in which we canidentify molecular conformational fluctuations (or states) by
the principal components and the important interactions by
the corresponding eigenvectors. Compared to other methods,
PEPCA can potentially be applied to evaluate the direct contribution of environmental molecules because PEPCA is based on
the analyzing potential energies, which can represent not only
intramolecular but also intermolecular atomic interactions.
However, directly applying PEPCA to investigate the contribution of environmental molecules presents some difficulties.
This is because PEPCA using intermolecular potential energy

026704-1©2011 American Physical Society

KOYAMA, KOBAYASHI, AND UEDA

PHYSICAL REVIEW E 84, 026704 (2011)

terms identify the intermolecular interactions that induce thelargest conformational change of the “whole” system, namely,
a protein and its environmental molecules (see Appendix A).
Therefore, PEPCA will collect the large fluctuations of theenvironmental molecules that may be irrelevant to conformational fluctuations of the target protein. To understandenvironmental effects on conformational fluctuations of thetarget protein, we need to focus on the conformational changeof the target protein rather than the whole system. In the present
study, we generalize our previous results [18] to evaluate
environmental effects on the conformational change of thetarget molecule.
This article is organized as follows: With respect to arbitrary
intermolecular perturbations, we show that changes in theconformational distribution of the target molecule due tothe perturbation can be evaluated by the variance of theconditional expectation of the perturbation potential energy.
By using this result, we introduce (i) distance-independentand (ii) distance-dependent perturbations of intermolecularinteractions. Then, we search the perturbations that induce
the largest change in the conformational distribution of thetarget molecule. We show that these can be solved by(i) PCA using conditional expectations of truncated and shiftedintermolecular potential energy terms and (ii) functional
principal component analysis (FPCA) [11,19] using products
of intermolecular forces and conditional cumulative densities.
We refer to these analyses as (i) intermolecular perturbation
analysis (IPA) and (ii) distance-dependent intermolecular
perturbation analysis (DIPA), respectively. For comparisonof the IPA and the DIPA, we applied them to the alanine
dipeptide isomerization in explicit water. We see that the
DIPA is more practical compared with the IPA. To test
the feasibility of the DIPA for larger molecules, we apply
the DIPA to the ten-residue chignolin folding in explicit water.
Finally, we discuss the computational cost of the DIPA for
protein molecules.
II. THEORY
A. Change in conformational distribution of target moleculedue to perturbation

We consider the change in the conformational distributionof the target molecule induced by a general perturbation. The
molecular system is described by the potential energy V (q,q ),where q are the coordinates of the target molecule and q
are the coordinates of the environmental molecules. In thiscase, the canonical distribution of the whole system at inverse
temperature β is represented asρ(q,q ) =1 −βV (q,q )
,
e
Z

(1)

where we defined the expectation value as−βV ≡ e−βV ρ(q,q )dqdq .
e

By using Eq. (3), the perturbed conformational distribution of
the target molecule ise−βV |q
ρ(q),
(5)
ρ  (q) = ρ  (q,q )dq =e−βV 
where the conditional expectation is defined as−βV (q,q )e|q ≡ e−βV (q,q ) ρ(q |q)dq .

(2)

The perturbed canonical distribution ρ  (q,q ) is represented asρ  (q,q ) =1 −βV  (q,q )e−βV (q,q )ρ(q,q ),e=Ze−βV (3)(6)

Next, we measure the change in conformational distributionof the target molecule due to the perturbation. For this purpose,
we use the Kullback-Leibler divergence (or relative entropy)[20–24]ρ  (q)dq  0.
(7)
D(ρ (q)||ρ(q)) ≡ ρ  (q) ln
ρ(q)This can be considered to be the expectation of the log-ratio[ln ρ  (q)/ρ(q)] under the perturbed equilibrium state. By using
the identity ρ  = ρ exp(ln ρ  /ρ), Eq. (7) can be expressed by
expectations under the unperturbed equilibrium state as

 
∞

ρ  (q) k
1

.
(8)
ln
D(ρ (q)||ρ(q)) =(k − 1)!
ρ(q)k=1
By applying cumulant expansions to Eq. (5), the ratio changeby the perturbation can be expanded as follows:ρ  (q)= −β (V | q − V )ρ(q)1+ [var(βV |q) − var(βV )] + · · · , (9)2where var(βV ) and var(βV |q) are the variance and theconditional variance of βV , respectively. By using Eq. (9)and the equality [25]lnvar(βV ) = var(β V |q) + var(βV |q)(10)

(sometimes called the law of total variance), Eq. (8) can be
expanded asD(ρ  (q)||ρ(q)) = 12 var(β V |q) + · · · .

(11)Thus, within the second-order approximation, we can evaluate
the change in the conformational distribution of the target
molecule by the variance of the conditional expectation of
the perturbation potential energy under the unperturbed equilibrium state. In Appendix A, we summarize the relationshipbetween the change in the conformational distribution of thewhole system and the target molecule by the perturbation.

where Z is the partition function. Next, we perturb the system
by adding a perturbation potential energy V (q,q ):V  (q,q ) = V (q,q ) + V (q,q ).

(4)

B. Intermolecular perturbation analysis

We develop a perturbation analysis of intermolecularinteractions. One natural perturbation of an interaction is tomultiply the original force with 1 + λ by introducing a perturbation parameter λ. With this perturbation, the interaction isstrengthened when λ > 0 and it is weakened when λ < 0. If
a perturbation to a certain interaction significantly affects the026704-2

PERTURBATION ANALYSES OF INTERMOLECULAR . . .

PHYSICAL REVIEW E 84, 026704 (2011)conformational distribution of the target molecule, we consider
it to be an “important” interaction. In the current ordinaryclassical force fields, intermolecular interactions consist ofvan der Waals interactions or electrostatic interactions. Since
they only depend on the atomic distance r, we denote a
nonbonded potential energy between two atoms as φ(r). To
avoid dealing the infinite long-range interaction, we only
perturb the interactions within some cutoff length rc . By
introducing a truncated and shifted potential energy [26]
φ(r) − φ(rc ) for r  rc ,φ(r; rc ) ≡(12)0
for r > rc ,the above-mentioned perturbation can be realized by adding
potential energy λφ(r; rc ) to the original potential energy.
By increasing the cutoff, we check the convergence of theperturbation analysis.
Then, we perturb all intermolecular interactions by introducing perturbation parameters λ = (λ1 , . . . ,λM )T . If thesystem includes atoms that have permutation symmetry [27],the perturbed potential energy is not invariant to the permutation of the atoms. To keep permutation invariance in theperturbed potential energy, perturbation parameters must havethe identical values when the corresponding potential energy
terms include atoms with permutation symmetry. With these
considerations, we can perturb all intermolecular interactions
with permutation invariance by adding perturbation potential
energyV (q,q ) =Mλk Vk (q,q ),covariance matrix is a positive semidefinite matrix, it can2,
be diagonalized with non-negative eigenvalues σ12 , . . . ,σMwhich are sorted in descending order, and the corresponding
orthogonal matrix U = (u1 , . . . ,uM ) [11,28] ascov(βV|q)ui = σi2 ui .

The ith eigenvector perturbation λ = δui maximizes thevariance Eq. (15) under |λ| = δ and λ · u1 = 0, . . . ,λ · ui−1 =0. By using Eqs. (11), (16), and (17), the change in theconformational distribution of the target molecule due to
λ = δui becomesD(ρδui (q)||ρ(q)) = 12 δ 2 σi2 + · · · .

V (q,q ) = δ

Labels Ik and Jk are the set of atoms with permutation
symmetry in the target molecule and the environmental
molecules, respectively. Although a perturbation within theenvironmental molecule (such as water-water interactions) canchange the conformational distribution of the target molecule,for simplicity we only perturb interactions between the target
molecule and environment molecules (such as protein-waterinteractions) in this study.
To identify important intermolecular interactions, we finda perturbation λ that induces a large change in the conformational distribution. For this purpose, we search for λthat maximizes D (ρλ (q)||ρ(q)) under the constraint |λ| = δ,where ρλ (q) is the conformational distribution of the target
molecule under the perturbation described by Eq. (13). Within
the second-order approximation, this is achieved by finding λthat maximizes the variance in Eq. (11). By using Eq. (13), thevariance is represented asvar(βV |q) = var βM

λk Vk |q

(15)

k=1

= λT cov(βV|q)λ,(16)

where cov (β  V| q) is the covariance matrix whose (i,j ) entry
is the covariance between β  Vi | q and βVj |q. Since theMUki Vk (q,q ).

(19)

k=1

By using Eqs. (9) and (19), the change in the ratio (orpopulation shift) of the conformation q induced by the
perturbation λ = δui isln

k=1i∈Ik j ∈Jk(18)

Thus, we can find the important combinations of the interactions to the target molecule in the top eigenvectors. This
procedure can be considered to be the PCA [11] using β  V| qor −β  V| q.
Next, we consider the contribution of the ith eigenvector
to the conformational distribution of the target molecule.
The perturbation λ = δui corresponds to the potential energy
perturbation [Eq. (13)] as(13)

where V(q,q ) = [V1 (q,q ), . . . ,VM (q,q )]T are defined asφk (rij ; rc ).
(14)
Vk (q,q ) ≡(17)

ρδui (q)= δgi (q) + · · · ,
ρ(q)

(20)

where we have introducedgi (q) ≡ −βui · ( V| q − V) .

(21)

By definition, gi (q) is the ith principal component of the
PCA using −β  V| q. Note that these principal components
have the opposite sign with respect to the definition in ourprevious article [18], which uses βV(q). Performing PCA
using −β  V| q instead of β  V| q is useful to visualize theresults by the biplot [29,30], as we will see in the numericalresults. By using Eqs. (19) and (20), we can analyze in detailthe effects of the ith eigenvector on the target molecularconformational distribution. We consider the case δ > 0. From
Eq. (19), the perturbation δui strengthens the kth potential
energy term Vk if Uki > 0, and it weakens this term if Uki < 0.
From Eq. (20), the perturbation δui increases the ratio of q
if gi (q) > 0, and decreases it if gi (q) < 0 in the first-orderapproximation. These results are summarized in Table I.
In summary, the PCA using −β  V| q identifies the
combinations of the intermolecular interactions that areimportant for the conformational distribution of the target
molecule. We refer to the PCA using −β  V| q as the IPA.
The ith eigenvector ui is the ith important combination of
the interactions [Eq. (19)], the ith eigenvalue σi2 measures thechange in the conformational distribution [Eq. (18)], and theith principal component gi (q) represents the change in the ratioof the target molecular conformation q due to the perturbation[Eq. (20)]. The detailed contributions of the interactions to thetarget molecular conformation are given in Table I. If we only
perturb interactions within the target molecule, the equality

026704-3

KOYAMA, KOBAYASHI, AND UEDA

PHYSICAL REVIEW E 84, 026704 (2011)

TABLE I. Effects of the ith eigenvector perturbation. A change inthe weight of a truncated and shifted intermolecular potential energy
term Vk is determined by Eq. (19). The ratio change in q is determinedby Eq. (20). The term “increase” or “decrease” indicates the ratiochange for the first-order approximation.

invariance by adding the potential energy V (q,q ) with
perturbation parameters λ(r) = [λ1 (r), . . . ,λM (r)]T as
V (q,q ) =MVk (q,q ),(27)

φk (rij ; rc ,λk ).

(28)

k=1

The kth interaction
Uki > 0δ>0δ<0Uki < 0Strengthen
Weakenwhere

Ratio of q

WeakenStrengthen

gi (q) > 0
Increase
Decrease

Decrease
Increase V| q = V(q) holds. Therefore, IPA includes the PEPCA that
corresponds to the intramolecular perturbation analysis [18].
In this sense, IPA is a natural generalization of PEPCA.

i∈Ik ,j ∈Jk

The meaning of labels Ik and Jk are identical in Sec. II B.
Next, we quantify the perturbation with Eq. (11). By
introducing the conditional density niJk (r|q) of Jk atoms at
a distance r from the ith atom for a given q, the conditional
expectation of Eq. (28) can be expressed as  rc  rcVk | q =λk (r )Fk (r )dr niJk (r|q) dr.

C. Distance-dependent intermolecular perturbation analysis

We develop a perturbation analysis of intermolecularinteractions with distance dependence by generalizing the
IPA. We consider the distance-dependent perturbation of anintermolecular interaction between the ith atom in the target
molecule and the j th atom in the environmental molecules. We
represent the force acting on the ith atom from the j th atom
through van der Waals or electrostatic interaction as Fj i , andthe distance between the ith atom and the j th atom as rij . With
this notation, one natural distance-dependent perturbation isrepresented as[1 + λ(rij )]Fj i ,

(22)



Vk (q,q ) ≡gi (q) < 0

i∈Ik

0

r

(29)
By defining the conditional cumulative density
 r
NiJk (r|q) ≡niJk (r  |q) dr 

(30)

0

and using the equality
dNiJk (r|q)= niJk (r|q),dr
the integration by parts of Eq. (29) leads to rc
Vk | q =λk (r)fk (r|q) dr,

(31)

(32)

0

where the distance-dependent perturbation function λ(rij ) is
introduced. To derive the potential energy that realizes theperturbation given in Eq. (22), we first note that the force Fj i
is the gradient of the corresponding potential energy φ(rij ).
Therefore, if we introduce coordinates of the ith atom as(xi ,yi ,zi ), the force becomes


∂ ∂ ∂
φ(rij ) = F (rij ) ej i ,Fj i = −,,(23)∂xi ∂yi ∂zi

where we have defined f(r|q) = [f1 (r|q), . . . ,fM (r|q)]T asfk (r|q) ≡ Fk (r)NiJk (r|q).
(33)
i∈Ik

As shown in Appendix B, fk (r|q) is also characterized bya derivative of Vk | q used in the IPA. By using Eqs. (27)and (32), the variance in Eq. (11) is expressed asvar(β V | q)

where we have introduced
F (r) ≡ −= var β

dφ(r)
dr

(24)

and ej i is a three-dimensional unit vector directed from atom j
to atom i. By using Eqs. (23) and (24), the perturbation given
in Eq. (22) can be realized by adding a potential energyφ(r; rc ,λ) ≡

rc



r λ(r )F (r )drforr  rc ,0

forr > rc ,(25)

where we have introduced the cutoff distance rc for convenience in the following analyses. If we ignore the dependence
on distance by using λ(r) = λ, Eq. (25) reduces toφ(r; rc ,λ) = λφ(r; rc ),=k=1
 rc  rc

λk (r)fk (r|q) dr

(34)

0

λ(r)T cov(βf(r|q),βf(r  |q))λ(r  )drdr  , (35)00where cov(βf(r|q),βf(r  |q)) is an M × M covariance matrix
whose i,j element is the covariance between βfi (r|q) andβfj (r  |q) for given r and r  .
By introducing a covariance operator Ĉ [19] that is defined
as
 rc
Ĉλ ≡
cov(βf(r|q),βf(r  |q))λ(r  ) dr 
(36)0

and an inner product(26)

where φ(r; rc ) is defined in Eq. (12). The perturbation is
identical to that of the IPA. Therefore, Eq. (25) is a naturalgeneralization of the perturbation used in the IPA. Then,we perturb all intermolecular interactions with permutation

M  rcλ,λ  ≡ rc

λ(r)T λ (r) dr,

(37)

0

Eq. (35) can be represented as

026704-4var(βV |q) = λ,Ĉλ  0.

(38)

PERTURBATION ANALYSES OF INTERMOLECULAR . . .

PHYSICAL REVIEW E 84, 026704 (2011)

By using the equality(a) Perturbation analysis
cov(βf(r|q),βf(r |q)) = cov(βf(r |q),βf(r|q)),
T

(39)





Ĉλ,λ  = λ,Ĉλ .

ΔV

(40)

Therefore, Ĉ is a self-adjoint operator. Furthermore, the nonnegativity of the variance [Eq. (38)] indicates that Ĉ is the
positive-semidefinite operator. Thus, Ĉ can be diagonalized by
the orthonormal eigenfunctions ui (r) = [U1i (r), . . . ,UMi (r)]Twith the corresponding non-negative eigenvalues σi2 (which
are sorted in descending order) asĈui = σi2 ui .

qGenerally, the eigenfunctions ui (r) constitute an infinite seriesi = 1, . . . ,∞.
To identify important intermolecular interactions with theirdistance information, we determine λ(r) that maximizes
D(ρλ (q)||ρ(q)) under the√ constraint |λ| = δ, where the norm
|λ| is defined as |λ| ≡ λ,λ. Within the second-order approximation [Eq. (11)], the maximization is equivalent to maximizing var (β V | q) under the constraint |λ| = δ. Because
the variance can be represented as the quadratic form of the covariance operator Ĉ [Eq. (38)], the variance maximization canbe solved by the diagonalization [Eq. (41)] [19]. In particular,a perturbation λ(r) = δui (r) maximizes the variance under theconstraints |λ| = δ and λ,u1  = 0, . . . , λ,ui−1  = 0. By the
perturbation λ(r) = δui (r), Eq. (11) becomes Eq. (18). Thus,
eigenfunctions with larger eigenvalue represent the importantintermolecular interactions with their distance information.
We note that finding λ(r) that maximizes the functional
form of Eq. (34) under the constraints |λ| = δ and λ,u1  =0, . . . , λ,ui−1  = 0 can be considered as the FPCA [11,19]
using βf(r|q) or −βf(r|q). By the perturbation λ(r) = δui (r),Eq. (9) becomes Eq. (20), where gi (q) is defined asgi (q) ≡ ui (r),−β[f(r|q) − f(r|q)].

(43)

By definition, gi (q) is the ith principal component of the FPCA
using −βf(r|q).
In summary, performing the FPCA using −βf(r|q) identifies the important combinations of the intermolecular interactions to the target molecule with their distance information.
We refer to the FPCA using −βf(r|q) as the DIPA. In amanner similar to the PEPCA and the IPA, we can interpretthe eigenfunctions, eigenvalues, and principal components of
DIPA in terms of the perturbation as follows: The ith eigenfunction perturbation λ(r) = δui (r) induces the largest change
in the conformational distribution of the target molecule under
the constraints |λ| = δ and λ,u1  = 0, . . . , λ,ui−1  = 0; theith eigenvalue σi2 measures the change in the conformational
distribution of the target molecule by 12 δ 2 σi2 [Eq. (18)]; the ith
principal component gi (q) represents the ratio change of theconformation q of the target molecule by δgi (q) [Eq. (20)].
Finally, we summarize perturbation analyses of atomic
interactions. The common idea is to find perturbations of

PCA using -β V(q) (PEPCA)(d) Distance-dependentintermolecular perturbation

PCA using -β<V|q> (IPA)

FPCA using -βf(r|q) (DIPA)

(41)

(42)

q

(c) Intermolecular perturbation

environmental molecules
Atom : target moleculeInteraction :bondednonbonded
Perturbation :strengthen
weakenThe orthonormality of the eigenfunctions is represented asui ,uj  = δi,j .

ρ'(q)

ρ(q)

we can show that(b) Intramolecular perturbation

FIG. 1. (Color online) (a) Perturbation analysis of atomic interactions. Atomic interactions are perturbed by adding a perturbation
potential energy V . The perturbation effect is measured by thechange in the conformation distribution of the target moleculeq. (b) Perturbation analysis of intramolecular interactions, which
can be solved by the PCA using −βV(q) (PEPCA) [18], whereV(q) are intramolecular potential energy terms. (c) Perturbation
analysis of intermolecular interactions, which can be solved by
performing the PCA using −β  V| q (IPA), where  V| q areconditional expectations of truncated and shifted intermolecularpotential energy terms. (d) Perturbation analysis of intermolecularinteractions with dependence on distance, which can be solved
by performing the FPCA using −βf(r|q) (DIPA), where f(r|q)
are products of intermolecular forces and conditional cumulative
densities as defined in Eq. (33).

atomic interactions that significantly change the conformational distribution of the target molecule [Fig. 1(a)]. The
difference is the three ways in which the atomic interactions are
perturbed: (i) If we consider perturbations of intramolecularinteractions [Fig. 1(b)], the perturbation analysis can be solvedby the PCA using −βV(q), where V(q) are intramolecularpotential energy terms (PEPCA) [18]. (ii) If we consider
perturbations of intermolecular interactions [Fig. 1(c)], the perturbation analysis can be solved by the PCA using −β  V| q,
where  V| q are conditional expectations of truncated andshifted intermolecular potential energy terms. (iii) If weconsider distance-dependent perturbations of intermolecularinteractions [Fig. 1(d)], the perturbation analysis can be solvedby the FPCA using −βf(r|q), where f(r|q) are products ofintermolecular forces and conditional cumulative densities
defined in Eq. (33) (DIPA).
D. IPA and DIPA procedure

For IPA and DIPA, we need to estimate the conditional
expectations and the conditional cumulative densities, respectively. These require additional molecular dynamics (MD)simulations of environmental molecules q with a fixed target
molecule q to sample q that follows ρ(q |q). For this purpose,first, we generate coordinates q and q that follow the canonicaldistribution ρ(q,q ). This can be done using an MD simulationthat can generate the canonical distribution such as a Langevindynamics [Fig. 2(a)]. Then, we can perform PEPCA [18] by

026704-5

KOYAMA, KOBAYASHI, AND UEDA

(b2)

× -β

PCA

V(q1)

V(q2)

V(qn)

q1,q'1q2,q'2qn,q'nq1,q'1,1

q2,q'2,1

qn,q'n,1

q1,q'1,n'q2,q'2,n'qn,q'n,n'

(d1)

V|q1

V|q2

V|qn

(e1)

f(r|q1)f(r|q2)f(r|qn)

× -β

PCA

IPA

× -β

FPCA

DIPA

(b1)

MD simulations
fixing target molecule(a) q0,q'0
(c)

(d2)

PHYSICAL REVIEW E 84, 026704 (2011)

(e2)

PEPCA

FIG. 2. Procedure to perform PEPCA, IPA, and DIPA. (a) An MD
simulation of the target molecule q and the environmental moleculesq . (b1) Calculation of intramolecular potential energy terms V(q).
(b2) PCA using −βV(q) (PEPCA [18]). (c) MD simulations of
the environmental molecules q with the target molecule q fixed.
(d1) Calculation of conditional expectations of truncated and shiftedintermolecular potential energy terms  V| q. (d2) PCA using−β V| q (IPA). (e1) Calculation of products of intermolecular forces
and conditional cumulative densities defined in Eq. (33). (e2) FPCA
using −βf(r|q) (DIPA).

roles regarding the states via intermolecular (peptide-water)interactions. We first apply PEPCA to understand the contribution of intrapeptide interactions to the peptide conformational
distribution. Next, we apply IPA and DIPA to understandthe contribution of peptide-water interactions to the peptide
conformational distribution.
MD simulations were performed by MD package AMBER
10 [31] with the ff03 force field [32] and TIP3P water [33].
Water molecules were included in a cubic box whose edge
length was determined to be the minimum distance between
peptide atoms and the faces of the box over 15 Å. As a result,2192 water molecules were contained within the unit cell.
Bonds involving hydrogen in the peptide were constrained by
SHAKE [34], and TIP3P water molecules were constrained
by SETTLE [35]. A direct-space cutoff distance of 10 Å
was used for the Lennard-Jones and electrostatic interactions.
Long-range electrostatic interactions were calculated using
particle mesh Ewald (PME) [36]. A 100 ps NpT (1 atm and
300 K) MD simulation was performed to equilibrate the water
density, which was stationary around 30 ps and resulted in a
(40.6 Å)3 cubic periodic box at 100 ps. Then a 10 ns NVT (300K) Langevin dynamics simulation with a collision frequencyof 1.0 ps−1 by integration with a 2 fs time step was performed.
The coordinates were saved every 1 ps, and 10 000 snapshots
were used for the analysis. From the Ramachandran plot
(Fig. 3), we see that the peptide has three stable states around(φ,ψ) = (−80◦ ,−20◦ ), (φ,ψ) = (−80◦ ,150◦ ), and (φ,ψ) =(−150◦ ,150◦ ). We call these states α, PPII, and β based on theircorresponding secondary structures of α helix, polyproline II
helix, and β sheet, respectively.

performing PCA using −βV(q) [Fig. 2(b)]. To get the concrete
expression for ρ(q |q), we divide the potential energy of thewhole system V (q,q ) into three terms: V (q,q ) = VT (q) +VTE (q,q ) + VE (q ). The terms VT (q), VTE (q,q ), and VE (q )are the target molecular term, target-environmental term, andenvironmental term, respectively. Then the distribution ρ(q |q)
is given asρ(q |q) =ρ(q,q )1 −β[VTE (q,q )+VE (q )]
e
=
.
Z(q)
ρ(q,q )dq

(44)This distribution can be considered to be a canonical distribution of q for a given q and the potential energy VTE (q,q ) +VE (q ) at inverse temperature β. Therefore, an MD simulation
with the given q and the potential energy VTE (q,q ) + VE (q )generates q following the distribution ρ(q |q). We then iteratesuch MD simulations for various different q, which can be
sampled from the original MD simulation [Fig. 2(c)]. After
this iteration, we can finally perform the IPA [Fig. 2(d)] andthe DIPA [Fig. 2(e)].
III. NUMERICAL RESULTS
A. PEPCA of the alanine dipeptide isomerizationin explicit water

For comparison of the IPA and the DIPA, we apply themto the alanine dipeptide isomerization in explicit water. In
this system, the target molecule (peptide) has multiple stablestates, and the environment molecules (water) have important

FIG. 3. (Color online) Ramachandran plot of the alanine dipeptide in explicit water. Angles φ and ψ are the dihedral angles indicatedby the curved arrows. Dots represent the scatter plot of φ and ψ(Ramachandran plot). The three-dimensional structures labeled α,
PPII, and β are snapshots at 120, 350, and 370 ps, respectively.
The structures were rendered by a molecular visualization program(VMD) [37].

026704-6

PERTURBATION ANALYSES OF INTERMOLECULAR . . .

PHYSICAL REVIEW E 84, 026704 (2011)

TABLE II. Number of the intrapeptide potential energy termswith permutation invariance. The permutation symmetry is relabelinghydrogens in a methyl group. The number of all potential energy terms
ignoring the permutation invariance is given in parentheses.
Potential type

No.

Bond
Angle
Dihedralvan der Waals
Electrostatic

15(21)
24(36)
31(45)
84(174)84(174)Total

238(450)After calculating V(q) with permutation invariance(Table II) for each of the 10 000 snapshots, PEPCA was performed [Fig. 2(b)]. The PCA was implemented by the singularvalue decomposition (SVD) after centering the data [11]. The
top three PEPCA eigenvalues overwhelm the other eigenvalues(Fig. 4). Therefore, we expect that the three corresponding
eigenvectors are sufficient to understand the intrapeptide
interactions that are important to the peptide conformational
distribution. In fact, we can clearly discriminate three stablestates by only the first two principal components that are
visualized when we plot the components of the first and thesecond eigenvectors (Uk1 ,Uk2 ) (squares in Fig. 5) and the first
and the second principal components [g1 (q),g2 (q)] (dots in
Fig. 5). This plot is known as the biplot [29,30]. As shown in
Fig. 5, the first principal components discriminate between the
α and the PPII + β states, and the second principal components
discriminate between the PPII and the β states. Therefore, we
can expect that the intrapeptide interactions that are importantto these three states are identified by the first and the second
eigenvectors.
The perturbation effects by the eigenvectors summarized inTable I are easily read from the biplot (Fig. 5) as follows.
First, we consider the first eigenvectors. If we strengthen

the interactions when Uk1 < 0 (left squares in Fig. 5) andweaken when Uk1 > 0 (right squares in Fig. 5), the ratio of thepeptide conformation q will increase when g1 (q) < 0 (leftdots in Fig. 5) or decrease when g1 (q) > 0 (right dots in
Fig. 5). Similarly, if we strengthen the interactions indicatedby the right squares and weaken the left squares, the peptide
conformation will increase in the right dots and decrease in theleft dots, respectively. We can also apply the same discussion
in the second eigenvector by considering the bottom and the
top direction in the biplot. Thus, the perturbation effects by the
first and the second eigenvectors are systematically understoodby the biplot.
By using these biplot properties, we first consider the
perturbation effects of the first eigenvector, whose principal
components discriminate between the α and the PPII + β
states. Smaller negative components of the first eigenvector(left squares in Fig. 5) are the electrostatic interactions 2N-3H,
1C-3N, 1O-3H, and 2N-2O. Among them, 2N-3H, 1C-3N,
and 1O-3H are attractive interactions and more favorable forthe α state (left dots in Fig. 5). The 2N-2O interaction isrepulsive, and more unfavorable to the PPII + β states (rightdots in Fig. 5). This can be confirmed in the structure because
2N and 2O atoms are separate in the α state. Larger positive
components of the 1st eigenvector (right squares in Fig. 5)
are the electrostatic interactions 2H-2O, 1C-3H, and 1O-3N.
The interaction 2H-2O is attractive and more favorable for
the PPII + β states. The interactions 1C-3H and 1O-3N are
repulsive and more unfavorable to the α state. Their atoms are
separate in the PPII + β states.
Next, we consider the perturbation effects of the second
eigenvector, whose principal components discriminate between the PPII and the β state. Smaller negative componentsof the second eigenvector (bottom squares in Fig. 5) are theattractive electrostatic interactions 1O-2C and 1C-2O, andmore favorable for the PPII state (bottom-right dots in Fig. 5).
The larger positive component of the second eigenvector (topsquares in Fig. 5) is the repulsive electrostatic interaction
1O-2O. This is more unfavorable to the PPII state, and its atoms
are separate in the β state (top-right dots in Fig. 5). In summary,by using the PEPCA biplot we can identify the three peptide
states in explicit water by the first and the second principal
components and their important intrapeptide interactions by
the components of the first and the second eigenvectors.

σ i2

B. IPA of the alanine dipeptide isomerization in explicit wateri
FIG. 4. PEPCA eigenvalues of the alanine dipeptide in explicit
water.

We apply the IPA to understand the contribution of
the peptide-water interactions to the peptide conformational
distribution. To perform MD simulations with fixed peptide
coordinates q [Fig. 2(c)], we used the belly option of the sander
program of AMBER 10 [31]. Although the original belly-option
code turns off the forces acting on the target molecule derived
from the potential energy, the random forces still act on thetarget molecule in the Langevin dynamics. Therefore, weslightly modified the code to turn off the random forces to
fix the target molecule completely. For initial coordinates,
we used 1000 conformations from the original 10 ns MD
simulation [Fig. 2(a)] with 10 ps intervals. Starting from these
initial conformations, 1000 NVT (300 K) MD simulations with
fixed peptide coordinates q were performed [Fig. 2(c)]. After

026704-7

KOYAMA, KOBAYASHI, AND UEDA

PHYSICAL REVIEW E 84, 026704 (2011)
Uk1

α

β

1O-2O

2H-3N
2N-2O
2N-3H

1O-3N

3H

2C

2H

2O

Uk2
1C-3H
2N-3N
2H-3H

1C
1O

2N

2H-2O

g2(q)

2N

3N

1C

1C-2C

2O
2C
3N

2H

3H
1O

1C-3N

PPII
3H
2H

1O-3H

1C-2O

1O-2C

2N
1C

2C

3N

2O

1O

g1(q)

FIG. 5. (Color online) The PEPCA biplot of the alanine dipeptide in explicit water. Dots show the scatter plot of the first and thesecond principal components [g1 (q),g2 (q)]. Left, bottom-right, and top-right structures represent the snapshot at 120 ps (in the α state),
350 ps (in the PPII state), and 370 ps (in the β state), respectively. Squares show the scatter plot of the top ten components of the first
and the second eigenvectors (Uk1 ,Uk2 ). The label “1O-3H” indicates the electrostatic interactions between 1O atom and 3H atom shown instructures.

calculating truncated and shifted intermolecular interactions(Table III), we performed the IPA [Fig. 2(d)].
Figure 6 shows the IPA eigenvalues with different cutoff
distances and sampling times for the conditional expectations. We can see the eigenvalues decrease and converge
by increasing the sampling time. Therefore, the IPA surely
gives the definite result with sufficient sampling time for theconditional expectations. Figure 7 shows the IPA eigenvalues
with different values for rc . The eigenvalues converge for
11 Å  rc  15 Å except the third eigenvalue. Figure 8
shows the IPA biplots using different values for rc . For
rc = 2 to 4 Å, states can not be identified from the principal
components. For rc = 5 to 6 Å, the first principal components
roughly discriminate the α state and PPII + β states. For rc 
7 Å, the separation of the α and PPII + β states becomes
clearer. However, the separation of the PPII state and theβ state is unclear. Therefore, we cannot discuss importantintermolecular interactions to discriminate the PPII state andthe β state.
We now identify the important intermolecular interactions
to the α state and PPII + β states by using the biplot with rc =15 Å. Smaller negative components of the first eigenvector [leftsquares in Fig. 8 (rc = 15 Å)] are 3H-H, 3N-O, 1C-H, 1O-O,
TABLE III. Number of peptide-water potential energy termswith permutation invariance. The permutation symmetries consistof relabeling hydrogens in a methyl group, hydrogens in a watermolecule, and water indexes. The symbols O and H in the left column
indicate interactions with oxygen or hydrogen atoms of water. Van
der Waals interactions about hydrogens in water are not incorporated
in the TIP3P water model [33].
Potential typeNo.

van der Waals (O)Electrostatic (H)Electrostatic (O)

16
16
16

Total

48

2O-H, and 2H-H electrostatic interactions, where -H and -O
indicate interactions with the hydrogen and oxygen atoms ofwater, respectively. Because 3H-H, 3N-O, 1C-H, 1O-O, and2H-H electrostatic interactions are repulsive, they represent
more unfavorable interactions for the PPII + β states than the
α state. The attractive interaction 2O-H prefers the α state to the
PPII + β states. Larger positive components [right squares inFig. 8 (rc = 15 Å)] are 3H-O, 1O-H, 3N-H, 2N-H, 1C-O, 2CO, and 2O-O electrostatic interactions. Attractive electrostaticinteractions 3H-O, 1O-H, 3N-H, 2N-H, 1C-O, and 2C-O prefer
the PPII + β states to the α state. The repulsive electrostaticinteraction 2O-O is a more unfavorable interaction for the αstate than the PPII + β states.
C. DIPA of the alanine dipeptide isomerization in explicit water

To incorporate the distance information systematically, weapply DIPA to the alanine dipeptide in explicit water. The
simple method to implement FPCA using −βf(r|q) (DIPA) is
to perform PCA with discretized functional data, as describedin Appendix C. We discretize the interval [0,rc ] with NB binsasrl ≡ lr, r ≡ rc /NB , l = 0, . . . ,NB .

(45)

To obtain f(rl |q) [Eq. (33)], we estimate the force Fk (rl ) andthe conditional cumulative density NiJk (rl |q). The force Fk (rl )can be calculated from the force field used. The conditional
cumulative density NiJk (rl |q) must be estimated from the MD
simulations that fix the coordinates q of the target molecule[Fig. 2(c)], which are identical to the MD simulations used toperform IPA. We can estimate NiJk (rl |q) by averaging thenumber of Jk atoms within rl of the ith atom. The PCA
was implemented by performing a SVD after centering thediscretized data [11]. By using the correspondence between
the PCA and the FPCA that is given in Table VII, we obtain
the eigenfunctions, eigenvalues, and principal components of
DIPA.

026704-8

PERTURBATION ANALYSES OF INTERMOLECULAR . . .

PHYSICAL REVIEW E 84, 026704 (2011)

FIG. 6. Dependence of the IPA eigenvalues of alanine dipeptide in explicit water on the cutoff distance rc and the sampling time. IPAs were
performed from 10 ps (10 sample) to 25 ns (25 000 sample) averages for the conditional expectations  V| q.

We investigate the convergence of DIPA by varying theconditions for the conditional cumulative densities. First,
we check the dependence of the DIPA eigenvalues on thesampling time to estimate the conditional cumulative densities(Fig. 9). Upon increasing the sampling time, the eigenvalues
decrease and converge. In particular, the top three eigenvalues
converge within 1 ns. This convergence rate is faster thanthat for the IPA using rc = 15 Å (Fig. 6). Second, weinvestigate the dependence of the DIPA eigenvalues on the
cutoff distance rc (Fig. 10). Upon increasing rc , the eigenvalues
increase and converge for rc  10 Å in DIPA. We attribute

this monotonically increasing property in DIPA to the factthat a distance-dependent perturbation with a longer cutoff
distance includes shorter distances, so that a larger change
in the conformational distribution can be induced. Since
DIPA converges for rc  10 Å, the cutoff value rc = 15 Åis sufficient for further analysis. Finally, we checked thedependence of the DIPA eigenvalues on the stride r of
the interval [0,rc ] for the conditional cumulative densities(Fig. 11). The top eigenvalues are robust against variation inr. In particular, for r  0.2 Å, the eigenvalues are almost
identical. In summary, allowing sufficient sampling time, alonger cutoff distance rc , and a small stride r for conditional
cumulative densities, DIPA gives definite results.
We analyze in detail the converged result of DIPA. Because
the top three eigenvalues converge within 1 ns (Fig. 9), we use
the results of DIPA with 1 ns MD simulations in the following
analysis. The first principal components discriminate between
the α and the PPII + β states, and the second principal
components discriminate between the PPII and the β states(dots in Fig. 12). Thus, the top two eigenfunctions suffice
to identify important peptide-water interactions for the threestates. Next, we consider the dependence on distance, whichcannot be directly interpreted from the IPA biplot. By using
Eqs. (37) and (42), the normality of the ith DIPA eigenfunction
is represented as
M  rc
k=1

Uki (r)2 dr = 1.

(46)

0

By using Eq. (46), we can quantify the contribution of the
distance r to the ith eigenfunction by introducing
FIG. 7. Dependence of the IPA eigenvalues of alanine dipeptide in
explicit water on the cutoff distance rc . The conditional expectations
were estimated from 25 ns (25 000 sample) averages.
026704-9

Si (r) ≡M
k=1

Uki (r)2  0,(47)

KOYAMA, KOBAYASHI, AND UEDA

PHYSICAL REVIEW E 84, 026704 (2011)

FIG. 8. (Color online) Dependence of the IPA biplot of the alanine dipeptide in explicit water on the cutoff distance rc . The conditional
expectations are estimated from 25 ns MD simulations (25 000 sample). Directions of the eigenvectors and corresponding signs of the principal
components are determined to be consistent with the PEPCA biplot (Fig. 5). The meaning of biplots is same as in Fig. 5. Labels α, PPII, and β
indicate principal components at identical snapshots in Fig. 5. The label “A-H” (“A-O”) indicates electrostatic interactions between dipeptide
atom A (shown on structures of Fig. 5) and hydrogen (oxygen) atoms of water molecules. Components of eigenvectors (squares) are selectedfrom the union of the top ten components of the first and the second eigenvectors.

which satisfy
 rc

Si (r) dr = 1.

(48)

0

FIG. 9. Dependence of the DIPA eigenvalues of the alanine
dipeptide in explicit water on the sampling time for conditional
cumulative densities (rc = 15 Å and r = 0.05 Å). DIPA was
performed every 10 ps from 10 ps to 25 ns MD simulations.

For r  10 Å, S1 (r) and S2 (r) are almost zero (Fig. 13).
This indicates that long-range intermolecular interactions
over 10 Å have no preference with respect to the threestates. These results also explain the convergence of the
top two eigenvalues for rc  10 Å in Fig. 10. Although

FIG. 10. Dependence of the DIPA eigenvalues of the alanine
dipeptide in explicit water on the cutoff distance rc . Conditional
cumulative densities are estimated from 1 ns MD simulations(r = 0.05 Å).

026704-10

PERTURBATION ANALYSES OF INTERMOLECULAR . . .

PHYSICAL REVIEW E 84, 026704 (2011)

FIG. 11. Dependence of the DIPA eigenvalues of the alanine
dipeptide in explicit water on the stride r of the interval [0,rc ].
Conditional cumulative densities are estimated from 1 ns MD
simulations (rc = 15 Å).

FIG. 13. Contribution of the distance r to the first and the second
DIPA eigenfunctions [Eq. (47)] of the alanine dipeptide in explicit
water. Conditional cumulative densities are estimated from 1 ns MD
simulations (rc = 15 Å and r = 0.05 Å).

long-range intermolecular interactions cause the slow convergence in IPA (Fig. 6), these smaller contributions oflong-range interactions are considered to be responsible for the
fast convergence in DIPA. The contributions S1 (r) and S2 (r)are maximized at r = 5.55 Å and r = 4.55 Å, respectively.
Therefore, hydration shell around these distances has animportant role to determine the preference among the threestates.
Since components of the PEPCA or IPA eigenvectors are
represented as points (squares in Figs. 5 and 8), eigenvector
components are compactly plotted in a biplot. However,components of DIPA eigenfunctions are represented as curves

FIG. 12. (Color online) DIPA biplots of the alanine dipeptide in explicit water. Conditional cumulative densities are estimates from 1 ns
MD simulations (rc = 15 Å and r = 0.05 Å). Dots show the scatter plot of the first and the second principal components [g1 (q),g2 (q)]. Labelsα, PPII, and β indicate principal components at identical snapshots in Fig. 5. The curves are components of the first and the second DIPA
eigenfunctions [Uk1 (r),Uk2 (r)] (0 Å  r  15 Å). Label “A-O” (“A-H”) indicates that the curve shows the electrostatic interactions between
A atom of the peptide (Fig. 5) and oxygen (hydrogen) atoms of water molecules. Numbers on the curves represent r in units of Å.
026704-11

KOYAMA, KOBAYASHI, AND UEDA

PHYSICAL REVIEW E 84, 026704 (2011)
D. PEPCA of the chignolin folding in explicit water

FIG. 14. Contribution of the kth component to the first and thesecond DIPA eigenfunctions [Eq. (49)] of the alanine dipeptide in
explicit water. Conditional cumulative densities are estimated from
1 ns MD simulations (rc = 15 Å and r = 0.05 Å). Components
are selected from the union of the top nine components of the first
and the second eigenfunctions. The significance of the labels (suchas “1C-O”) is identical to that used in Fig. 12.

because the component of DIPA eigenfunctions depends on
distance. Therefore, representing all components in a biplot iscomplicated. To circumvent this, we select some componentsthat are large contributors to the eigenfunction and plot themon different biplots. By using Eq. (46), we may quantify thecontribution of the kth component to the ith eigenfunction byintroducing
 rc
Ski ≡
Uki (r)2 dr  0,
(49)

To test the feasibility of the DIPA for larger molecules, we
use the chignolin folding in explicit water. The chignolin is aten-residue peptide (Gly1-Tyr2-Asp3-Pro4-Glu5-Thr6-Gly7Thr8-Trp9-Gly10) that folds into β hairpin (Protein Data Bank
(PDB) ID 1UAO) [38]. In previous studies [39–41], folding
simulations of the chignolin were successfully performed. The
protocol of the following MD simulation is identical to thatof the alanine dipeptide (Sec. III A). The first structure in
1UAO was used for the initial coordinates of the chignolin.
Water molecules were included in a cubic box whose edge
length was determined to be the minimum distance between
peptide atoms and the faces of the box over 20 Å. As a result,
6437 water molecules were contained within the unit cell. To
neutralize the whole system, 12 Na+ and 10 Cl− ions wereadded. Małolepsza et al. [27] pointed out that some improper
torsion potentials of AMBER force field are not invariant toexchange of symmetrical atoms and these could be avoided
by reordering the atoms for the torsion angle. We applied
their program [42] to the topology file. This reordered atomsof improper torsions of Tyr2, Asp3, Glu5, and C-terminalof Gly10. After energy minimization and 100 ps NpT MD
simulation at 1 atm and 312 K (melting temperature) [38], thebox size became (58.4 Å)3 . Then, a 1 μs NVT (312 K) Langevindynamics simulation with a collision frequency 1.0 ps−1 byintegration with a 2 fs time step was performed. The root
mean square deviation (RMSD) of the chignolin (Fig. 15)
shows native to non-native transitions three times.
After calculating the potential energy terms with permutation invariance (Table IV) with 1 ns interval (1000 structures),the PEPCA was performed. The PCA was implemented by
the SVD after centering the potential energy terms. Two large
gapped eigenvalues are observed in Fig. 16. Therefore, we
can expect that the top two principal components identify
conformational states and their corresponding eigenvectors
can identify important intramolecular interactions to each0which satisfyMSki = 1.

(50)

k=1

Figure 14 shows Ski with large values for the first and thesecond eigenfunctions.Then, we represent these top components by biplots (Fig. 12). First, the attractive interactionthat prefers the α state to the PPII + β states is 2O-H
electrostatic interaction. Second, the attractive interaction that
prefers the PPII state to the β state is 2H-O electrostaticinteraction. Finally, the attractive interactions that prefer theβ state to the PPII state are 1C-O, 1O-H, 2C-O, and 2O-H
electrostatic interactions. By comparing components of the
IPA eigenvectors for rc = 15 Å [Fig. 8 (rc = 15 Å)] withcomponents of the DIPA eigenfunctions (Fig. 12), we see thatthey correspond. For example, 1C-H components point to theleft and 1C-O components point to the right for both IPA and
DIPA. Thus, we can confirm the consistency between the IPA
and the DIPA.

FIG. 15. RMSD of Cα atoms from the first structure of thechignolin in explicit water.

026704-12

PERTURBATION ANALYSES OF INTERMOLECULAR . . .

PHYSICAL REVIEW E 84, 026704 (2011)

TABLE IV. Number of the chignolin potential energy terms withpermutation invariance. The number of all potential energy terms
ignoring the permutation invariance is given in parentheses.
Potential type

No.

Bond
Angle
Dihedralvan der Waals
Electrostatic

116 (141)
204 (249)
313 (403)
6147 (9063)
6147 (9063)

Total

12 927 (18 919)

state. Figure 17 shows the PEPCA principal components. By
comparing Fig. 15, we can see that the first principal components correspond to the native and non-native transitions.
The second principal components indicate that the non-native
structures around 0.1 μs are different from the non-native
structures around 0.6 and 0.75 μs.
We analyze the conformational states by using the PEPCA
biplot (Fig. 18). In the left side of the biplot [g1 (q) < 0],dense dots are observed. By comparing Figs. 15 and 17,these dots correspond to the native structures. A snapshot
from the state also confirms the native β-hairpin structure(dashed box labeled N in Fig. 18). We refer to the state as
the “N” state (native state). Dots in the right side of the
biplot [g1 (q) > 0] corresponds to the non-native state. In thebottom-right [g1 (q) > 0 and g2 (q) < 0] of the biplot, we cansee two dense dots. Structures from these states (dashed boxes
labeled M1 and M2 in Fig. 18) share the same backbone
structures. Although the M1 structure has no contact between
N- and C-terminal residues, the M2 structure has contact
between them. By comparing the backbone of the native
β hairpin, these structures can be recognized as misfolded
structures. Therefore, we refer to these states as the “M1”state and the “M2” state (misfolded state), respectively. In the
top-right [g1 (q) > 0 and g2 (q) > 0] of the biplot, we can see
scattered dots. A structure from these dots (dashed box labeled
U in Fig. 18) has unfolded backbone structure. We refer to this

FIG. 16. The PEPCA eigenvalues of the chignolin in explicit water.

FIG. 17. The top three PEPCA principal components of thechignolin in explicit water.

state as the “U” state (unfolded state). Thus, we can identify
the four states, that are the native (N) state [g1 (q) < 0], two
misfolded (M1 and M2) states [g1 (q) > 0 and g2 (q) < 0], andthe unfolded (U) state [g1 (q) > 0 and g2 (q) > 0] in the PEPCA
biplot.
Next we identify preferable intramolecular interactions toeach state by the PEPCA biplot (Fig. 18). First, we consider
the interactions that prefer the N state to the non-native state.
These interactions are located on the left side of the biplot(Uk1 < 0). We can see that the bottom-left (Uk1 < 0 and
Uk2 < 0) components are attractive electrostatic interactions
for the N- and the C-terminal residues except 1N-10O
repulsive electrostatic interaction. We can confirm 1H-10O
and 1O-10H contacts in the N structure. The top-left (Uk1 <0 and Uk2 > 0) components are 3OD-6HG1, 6OG1-8HG1,
3H-8O, and 3OD-8HG1 electrostatic interactions. These are
attractive electrostatic interactions among Asp3, Thr6, andThr8 residues. We can also confirm 3OD-6HG1, 6OG1-8HG1,
and 3H-8O atoms contact in the N structure.
Second, we identify important interactions to the M1 andthe M2 states. The N state and the M1 + M2 states are
discriminated by the bottom-right direction. As shown inAppendix D, we can understand a perturbation effect of thelinear combination of the first and the second eigenvectors
from the biplot. This indicates that the favorable interactions
for the M1 and the M2 states can be found in the bottom-right
direction. In the bottom-right direction, there are repulsive
electrostatic interactions among Asp3, Thr6, and Thr8 residues
that destabilize the native state. Attractive electrostatic interactions in the bottom-right directions are 3H-7O, 1O-9H, and2C-7O. From the M1 and the M2 structures in Fig. 18, we cansee 3H-7O and 1O-9H atoms contact and these contacts are notobserved in the N structure. The M1 state and the M2 state are
discriminated by the top-right and bottom-left directions. The
attractive electrostatic interaction 1H-10O prefers the M2 state026704-13

KOYAMA, KOBAYASHI, AND UEDAU

PHYSICAL REVIEW E 84, 026704 (2011)

1C
2N

1H

M1

7C
7O

8OG1
8O
10C

1O
3N

10O

1N

3OD

1O-9N

3OD

10O
10C
10H
7O
8O
8HG1

6HG1
6OG1

3OD
5OE

10C

2N-10O
1H-10C
1C-10C

10O 10H 9H 8O 7O

6OG1

1O-10O

6OG1-8HG1
3OD-6HG1
3OD-8HG1
3H-8O

M2
Uk2

1N 2N
1H 1C
1O
2C 3H

g2(q)

3N-7O

3N 3H

1H

Uk1
N

2C
1O 1C2N

2C-10O
1O-10H
1O-10C
1HA-10O
2C-7O
1N-10O
1O-9H
3OD-5OE
1C-10O
3H-7O

3OD-6OG1
3OD-10O
3N-8O
3OD-8O
3OD-8OG1
6OG1-8OG1

1C
2C 3N 3H 3OD
1N 1O 2N
5OE
1H
10O

10C

1H-10O

10H

9H 8O 8OG1

6OG1
7O

g1(q)
FIG. 18. (Color online) The PEPCA biplot of the chignolin in explicit water. Dots show the scatter plot of the first and the second
principal components [g1 (q),g2 (q)]. Squares show the scatter plot of the components of the first and the second eigenvectors (Uk1 ,Uk2 ). Top 20
components of the first and the second eigenvectors are shown. Labels such as 1H-10O indicate the electrostatic interaction between the 1H
atom and the 10O atom in the dashed boxes. Structures represent the snapshot at 0.68 μs (in the N state), 0.73 μs (in the M1 state), 0.75 μs (inthe M2 state), and 0.114 μs (in the U state), respectively. In M1 and M2 structures, the side-chain of Trp9 is not drawn.

to the M1 state. We can confirm that there is 1H-10O contact
in the M2 structure and there is not in the M1 structure.
Finally, we identify favorable interactions for the U state.
Components in the top-right direction are 1O-10O, 2N-10O,
1H-10C, and 1C-10C repulsive electrostatic interactions. Since
the U state does not have stable intramolecular interactions,these results are reasonable. In summary, we can identify thefour states (N, M1, M2, and U) by the top two principal
components and their important intramolecular interactions
by corresponding eigenvectors.

E. DIPA of the chignolin folding in explicit water

To identify important intermolecular interactions, we apply
the DIPA to the chignolin folding in explicit water. The
protocol to perform the DIPA is identical to that of thealanine dipeptide (Sec. III C). For initial coordinates, we used1000 conformations from the original 1 μs MD simulation[Fig. 2(a)] with 1 ns intervals. Starting from these initialconformations, 1000 NVT (312 K) MD simulations with
fixed chignolin coordinates q were performed [Fig. 2(c)].
Conditional cumulative densities were calculated using rc =25 Å and r = 0.05 Å. The number of the chignolin-solvent
interaction with permutation invariance is M = 798 (Table V)

and the number of bins for the conditional cumulative densities
is NB = rc /r = 500. The PCA variable number for the
FPCA (Appendix C) is p = MNB = 399 000 and the sample
size is n = 1000. Since the p × n data matrix is large, the
SVD implementation for the PCA is difficult to apply. Insteadof the SVD, we implemented the PCA by diagonalizing an
n × n centered Gram matrix of the data matrix [18,43].
There are three large DIPA eigenvalues and they are almost
converged within 100 ps MD simulations (Fig. 19). Therefore,
we can expect that the top three DIPA eigenfunctions identify
important chignolin-solvent interactions. Interestingly, theconvergence rate is faster than that of the alanine dipeptide(Fig. 9). We discuss this point in Sec. IV.

TABLE V. Number of chignolin-solvent potential energy termswith permutation invariance.
Potential typevan der Waals (Na+ , Cl− , O)
Electrostatic (Na+ , Cl− , H, O)
Total

026704-14

No.
114 × 3114 × 4798

PERTURBATION ANALYSES OF INTERMOLECULAR . . .

PHYSICAL REVIEW E 84, 026704 (2011)

FIG. 19. Dependence of the DIPA eigenvalues of the chignolinin explicit water on the sampling time for conditional cumulative
densities (rc = 25 Å and r = 0.05 Å). DIPA was performed every
10 ps from 10 to 500 ps sampling time.

The top two principal components (Fig. 20) exhibit similar
dynamics to that of the PEPCA principal components (Fig. 17).
This indicates that the four states (N, M1, M2, and U)
can be identified by the top two principal components andtheir important chignolin-solvent interactions can be identified
by the corresponding eigenfunctions. The third principal
components discriminate between the U state and the otherstates. As shown in Fig. 17, the third PEPCA principal
components do not show this discrimination. Since the U stateis exposed to the solvent compared with the other folded three

FIG. 20. Top three DIPA principal components of the chignolinin explicit water. Conditional cumulative densities are estimated from100 ps MD simulations (rc = 25 Å and r = 0.05 Å).

FIG. 21. Contribution of the distance r to DIPA eigenfunctions[Eq. (47)] of the chignolin in explicit water. Conditional cumulative
densities are estimated from 100 ps MD simulations (rc = 25 Å andr = 0.05 Å).

states (N, M1, and M2), the third eigenfunction is expected torepresent the chignolin-solvent interactions that stabilize these
exposed atoms. We confirm this later by using biplots.
The first and the second DIPA eigenfunctions are 0 forr  20 Å (Fig. 21). Therefore, the cutoff value rc = 25 Å issufficient to converge the DIPA. The eigenfunctions have large
value for 5 Å  r  8 Å. This suggests that the solvation shellaround these distances have different contribution to each state.
To identify important intermolecular interactions to eachstate, we select the large contributors to each eigenfunction(Fig. 22). Based on Ski , we show biplots using the first andthe second or the third principal components (Fig. 23). First,
the attractive interaction that prefers the N state to the M1 +M2 states is 7O-H electrostatic interaction. From structures inFig. 18, we can see that 7O atom is exposed to the solvent in
the N state and it contacts 3H atom in the M1 + M2 states.
Second, the attractive interaction that prefers the M1 + M2
states to the N states is 8O-H electrostatic interaction. The 8O
atom contacts 3H in the N state and it is exposed to the solventin the M1 + M2 states. Third, the attractive interactions that
prefer the M1 state to the M2 state are 10C-O and 10O-H
electrostatic interactions. The 10O atom contacts 1H in the M2state and it is exposed to the solvent in the M1 state. Finally,
the attractive interaction that prefers the U state to the folded
states (N, M1, and M2) is 3OD-H electrostatic interaction. The
3OD atoms contact the peptide atoms in the N, M1, and M2
states and they are exposed to the solvent in the U state.
In summary, the DIPA identify the four states (N, M1,
M2, and U) by the top three principal components and theirimportant intermolecular interactions by their corresponding
eigenfunctions. Although we only use the information about
intermolecular interactions, the DIPA can correctly identify thechignolin conformational states. This can be possible becausesome atoms in the chignolin expose to the solvent with state
dependence.

026704-15

KOYAMA, KOBAYASHI, AND UEDA

PHYSICAL REVIEW E 84, 026704 (2011)
IV. DISCUSSION AND CONCLUSIONS

N
M1
M2
1O-O

Uk1(r)

N

M1
8O-O

10C-O

g1(q)
Uk1(r)

U

N
M1
M2

8O-H
M2

7O-O

g1(q)

7C-H

U

10O-H

N

Uk3(r)
g2(q)

g3(q)

M1
M2

M1
M2

N

Uk1(r)

U

7O-H
U

M1
M2

g1(q)

Uk1(r)

U

3OD-O

g1(q)

g1(q)

7C-O

U
N

Uk1(r)

g1(q)

N
M1

10C-H

10O-O

g1(q)

g1(q)

Uk2(r)

M1
M2
1C-H

3OD-H
U

Uk2(r)g2(q)N

Uk1(r)

Uk3(r)
g3(q)1O-H

Uk2(r)g3(q)U

Uk2(r)
g2(q)

g2(q)

1C-O

Uk1(r)

Uk2(r)g2(q)Uk1(r)

Uk3(r)

FIG. 22. Contribution of components to DIPA eigenfunctions[Eq. (49)] of the chignolin in explicit water. Conditional cumulative
densities are estimated from 100 ps MD simulations (rc = 25 Å andr = 0.05 Å). Components are selected from the union of the top
nine components of the first, the second, and the third eigenfunctions.
Label A-O (A-H) represents the electrostatic interactions between
the atom A of the chignolin (dashed boxes in Fig. 18) and oxygen(hydrogen) atoms of water molecules.

To identify conformational states of the target moleculeand intermolecular interactions that are important for eachstate, we developed perturbation analyses of intermolecular interactions. We show (i) distance-independent and(ii) distance-dependent perturbation analyses can be realized
by performing (i) a PCA using conditional expectations oftruncated and shifted intermolecular potential energy termsand (ii) a FPCA using products of intermolecular forces andconditional cumulative densities. We refer to these analyses as
IPA and DIPA, respectively.
For comparison of the IPA and the DIPA, we appliedthem to the alanine dipeptide in explicit water. Although
the first IPA principal components identify the α state and
the PPII + β states for larger cutoff length, the separationof the PPII state and the β state is unclear in the second
IPA principal components. DIPA eigenvalues converge with asufficient sampling time, longer cutoff distance rc , and smaller
stride r of the interval [0,rc ] for the conditional cumulative
densities. With this convergence condition, DIPA identifies
three peptide states from the top two eigenfunctions. The
separation between the PPII and the β state is clearer than thatfor IPA. Although long-range peptide-water interactions cause
the slow convergence in IPA, they have small contributions to
the top two eigenfunctions in DIPA. This fact is considered tobe responsible for the fast convergence of DIPA. Thus, DIPA
improves IPA by introducing a dependence on distance.
To show the feasibility of the DIPA for larger molecules,
we applied the DIPA to the ten-residue chignolin folding
in explicit water. Interestingly, the DIPA converges faster

M2

FIG. 23. (Color online) DIPA biplots of the chignolin in explicit water. Conditional cumulative densities are estimated from 100 ps MD
simulations (rc = 25 Å and r = 0.05 Å). Dots show the scatter plot of the first and the second or the third DIPA principal components. Labels
N, M1, M2, and U indicate principal components at identical snapshots in Fig. 18. The curves are components of the first and the second or thethird DIPA eigenfunctions. A label “A-O” (“A-H”) indicates that the corresponding curve shows the electrostatic interactions between the A atomof the chignolin (dashed boxes in Fig. 18) and oxygen (hydrogen) atoms of water molecules. Numbers on the curves represent r in units of Å.
026704-16

PERTURBATION ANALYSES OF INTERMOLECULAR . . .

PHYSICAL REVIEW E 84, 026704 (2011)

than that of the alanine dipeptide. The top three DIPA
principal components identify the four states (native state, two
misfolded states, and unfolded state) and the corresponding
eigenfunctions identify important chignolin-water interactions
to each state. Although the DIPA only use information about
intermolecular interactions, it can identify the conformational
states by recognizing atoms in the chignolin that expose to thesolvent with state dependence.
Top two DIPA eigenvalues of the chignolin (Fig. 19) exhibit
faster convergence than that of the alanine dipeptide (Fig. 9).
To explain this behavior, we analyze the relative error of
eigenvalues due to the finite sampling estimation. We assume
that the DIPA is performed according to the procedure shownin Appendix C. For convenience, we treat fk (rl |q) as a vector[f1 (r1 |q), . . . ,fM (rNB |q)]T in the following part. We denotethe n sample estimator of fk (rl |q) as fk,n (rl |q). When n islarge, the multivariate central limit theorem indicates1
fk,n (rl |q) = fk (rl |q) + √ k (rl |q),n(51)

where k (rl |q) follows the multivariate normal distribution
k (rl |q) ∼ N (0, cov(fk,1 (rl |q)|q)).

(52)

For simplicity, we assume fk (rl |q) and k (rl  |q) are uncorrelated and replace cov(fk,1 (rl |q)|q) with its averagecov (fk,1 (rl |q)|q). Under these assumptions, the covarianceof Eq. (51) becomes√cov(−βfk,n (rl |q) r)√√1= cov(−βfk (rl |q) r) + cov(−βfk,1 (rl |q) r|q).
n(53)
With Table VII, we note that the ith eigenvalue of the second
covariance matrix in Eq. (53) is σi2 . We denote the ith
eigenvalue of the first covariance matrix and the average of2
the third covariance matrix in Eq. (53) as σi,nand σi2 ,respectively. By applying Weyl’s inequality [28] to Eq. (53),the relative error of the ith eigenvalue is bounded as
2
2
σi,n
− σi2
1 σMN1 σ12
B


.
n σi2
n σi2
σi2

(54)

The first inequality in Eq. (54) indicates that finite sampling
eigenvalues are overestimated. The second inequality in
Eq. (54) shows the relative error is in inverse proportion to
the sampling time. These properties qualitatively explain themonotonically decreasing behavior of eigenvalues with respect
to the sampling time in Figs. 9 and 19. The second inequality
in Eq. (54) also indicates that the sampling time required toachieve some relative error is proportional to σ12 /σi2 . In thethird covariance matrix in Eq. (53), the variation is due to thenumber fluctuation of atoms within rl because the intermolecular force is constant at rl . Since the excluded volume of thechignolin is larger than that of the alanine dipeptide, the number of atoms within rl is smaller in the chignolin. Therefore, thesmaller number fluctuation can be expected in the chignolin.
This may induce the smaller value of σ12 relative to the conformational fluctuation of the target molecule (σi2 ). Thus, one
possibility of the faster convergence of the top two eigenvaluesof the chignolin is due to the larger excluded volume and thelarger conformational fluctuation of the chignolin.
Here we discuss the applicability of DIPA to protein
molecules. To perform DIPA, we need to perform (i) a longtime MD simulation of the whole system, (ii) MD simulations
that fix the target molecule, and (iii) FPCA using products ofintermolecular forces and conditional cumulative densities.
First, we consider the feasibility of (i) the long-time MD
simulation of the whole system. Of course, to identify the
conformational states by DIPA, the states have to be sampled
by the MD simulation. In many protein molecules, the timescale of functional motions is submicrosecond to millisecond,so the corresponding MD simulation is required.
Second, we consider the feasibility of (ii) MD simulations
that fix the target molecule. Before performing the MD simulations, we confirm that the functional states are sampled in the
MD simulation. For this purpose, we can use PEPCA [18] by
identifying molecular states and their important intramolecularinteractions. After confirmation of the molecular states, we
determine the number n of MD simulations to perform (for thisstudy, we performed n = 1000 MD simulations). If the numberof functional states is small, n = 1000 to 10 000 is sufficient.
We must also determine the sampling time to estimate
conditional cumulative densities. In the chignolin, the top
three eigenvalues converge within 100 ps MD simulations.
Because the number of water molecules around an atom of
the peptide and a protein is not significantly different, we
can expect that a few ns MD simulations also suffice for theprotein. In the current massive computing clusters, it is feasible
to perform 1000 to 10 000 independent MD simulations witha few ns.
Finally, we discuss the feasibility of (iii) the FPCA using
products of intermolecular forces and conditional cumulative
densities. Since we have assumed n = 1000 to 10 000,
as per the discussion above, we can perform PCA using
diagonalization of an n × n centered Gram matrix as used inSec. III E. The centered Gram matrix can be calculated fromthe Gram matrix. The (i,j ) component of the Gram matrix canbe calculated from the ith and the j th structures. Therefore,
the Gram matrix calculation can be easily decomposed andparallelized. In summary, the most time-consuming part forprotein molecules is (i) the long-time MD simulation of thewhole system, which is currently a common problem in thefield of MD simulations. Although micro- to millisecond MD
simulations are difficult, the recent breakthrough involving
the special purpose machine [44] is promising. The otherparts [(ii) and (iii)] are feasible with current computational
power.
We can now analyze the molecular conformational fluctuations based on their intramolecular interactions (PEPCA) andintermolecular interactions (DIPA). The PEPCA is useful forunderstanding the molecular functional states based on theiramino acid residues or bases interactions, whereas the DIPA
is useful for understanding the regulation mechanism of thestates by the environmental molecules. An important aspectof PEPCA and DIPA is that they can be applied to ordinaryclassical mechanical force fields in MD simulation. Furthermore, both results can be systematically visualized by biplots.
Therefore, we believe that PEPCA and DIPA provide general
analysis methods to understand the molecular conformational

026704-17

KOYAMA, KOBAYASHI, AND UEDA

PHYSICAL REVIEW E 84, 026704 (2011)fluctuations based on their intra- and intermolecular interactions.

TABLE VI. Second-order approximations of changes in theconformational distribution of the system [D(ρ  (q,q )||ρ(q,q ))] andtarget molecule [D(ρ  (q)||ρ(q))] due to intra- [V (q)] and inter[V (q,q )] molecular perturbations.

ACKNOWLEDGMENTS

We thank Itoshi Nikaido for his valuable discussions. We
are also thankful to the RIKEN Super Combined Cluster(RSCC) and the RIKEN Integrated Cluster of Clusters (RICC)
for the computational resources. This research was supportedby the Research Fund for Computational Science in RIKEN(H.R.U.), the Center for Developmental Biology (CDB)(H.R.U.), the Quantitative Biology Center (QBiC) (H.R.U.),a Grant-in-Aid for Scientific Research on Innovative Areas(T.J.K), and Aihara Project, the FIRST program from JSPS,
initiated by CSTP (T.J.K).
APPENDIX A: CHANGES IN THE CONFORMATIONAL
DISTRIBUTION OF THE SYSTEM AND TARGET
MOLECULE DUE TO A PERTURBATION

Changes in the conformational distribution of the system[D(ρ  (q,q )||ρ(q,q ))] and target molecule [D(ρ  (q)||ρ(q))]
due to a perturbation V (q,q ) can be connected by theequality [20,22]D(ρ (q,q )||ρ(q,q ))








= D(ρ (q)||ρ(q)) + D(ρ (q |q)||ρ(q |q))V ,(A1)

which is called the chain rule. The second term on the righthand side is called the conditional Kullback-Leibler divergence(or conditional relative entropy) and is explicitly defined asD(ρ (q |q)||ρ(q |q))V ρ  (q |q)  ρ (q) dq  0. (A2)dq≡ρ  (q |q) ln
ρ(q |q)
By the non-negative nature of the Kullback-Leibler divergence, Eq. (A2) is also non-negative and Eq. (A1) indicatesthe inequality
D(ρ  (q,q )||ρ(q,q ))  D(ρ  (q)||ρ(q)).

(A3)

Therefore, the change in the conformational distribution of
the whole system due to the perturbation is always larger thanthe change in the conformational distribution of the target
molecule. By using Eqs. (3) and (5), the perturbed conditional
distribution ρ  (q |q) is expressed as

D (ρ  (q,q )||ρ(q,q ))D (ρ  (q)||ρ(q))
D (ρ  (q |q)||ρ(q |q))V

e−βV (q,q )ρ(q |q).
e−βV |q







D(ρ (q,q )||ρ(q,q )) = D(ρ (q)||ρ(q)).

V (q,q )1var(βV )2
1
var(βV )2

1
var(βV )2
1
V | q)var(β21var(βV |q)2

0

D(ρ  (q,q )||ρ(q,q )) = 12 var(βV ) + · · · .

(A6)

By using Eqs. (5) and (A4), Eq. (A2) can be expanded asD(ρ  (q |q)||ρ(q |q))V = 12 var(βV |q) + · · · .

(A7)Therefore, using Eqs. (A6), (11), and (A7), the second-orderapproximation of the chain rule [Eq. (A1)] gives
1
var(βV ) = 12 var(βV |q) + 12 var(βV |q).
2(A8)

We note that this equality is the law of total variance [Eq. (10)].
The results are summarized in Table VI.
APPENDIX B: RELATIONSHIP OF FUNCTIONS USEDIN THE IPA AND THE DIPA

By using the conditional density niJk (r|q) and Eq. (14),Vk | q is represented as  rcVk |q =[φk (r) − φk (rc )]niJk (r|q) dr.
(B1)
i∈Ik

0

The integration by parts of Eq. (B1) leads to rc
Vk |q =fk (r|q) dr.

(B2)

0

Thus, Vk | q (used in the IPA) is related by integration of
fk (r|q) (used in the DIPA). This is also expressed as
fk (r|q) =d Vk | q
,
drc
rc =r

(B3)

which means fk (r|q) is a derivative of Vk | q.
APPENDIX C: PROCEDURE TO PERFORM FPCA

(A4)

If we consider an intramolecular perturbation V (q), Eq. (A4)can be expressed as ρ  (q |q) = ρ(q |q) and Eq. (A2) becomes
0. Therefore, the chain rule [Eq. (A1)] leads to the equalityV (q)

Next we consider the second-order approximation ofeach term in the chain rule [Eq. (A1)]. By using Eq. (3),D(ρ  (q,q )||ρ(q,q )) can be expanded asρ  (q |q) =Perturbation potential energy

Distribution changeby the perturbation(A5)

Thus, the changes in the conformational distribution of
the system and target molecule due to the intramolecular
perturbation V (q) are equal.

The simple method to perform FPCA is to apply PCA withdiscretized functional data [11,19]. We describe this methodin our setting. By discretization Eq. (45), Eq. (34) can beapproximated asvar(βV |q) ≈ var β

NB
M 


λk (rl )fk (rl |q)r .

(C1)

k=1 l=1

The right-hand side of Eq. (C1) can be considered the variance√of the linear combination of√the variable −βfk (rl |q) rwith the coefficient λk (rl ) r for k = 1, . . . ,M and026704-18

PERTURBATION ANALYSES OF INTERMOLECULAR . . .

PHYSICAL REVIEW E 84, 026704 (2011)

l = 1, . . . ,NB . Therefore, we can consider the corresponding
PCA. We denote the ith eigenvector, eigenvalue,and prin√
cipal component of the PCA as Uki (rl ) r, σi2 , and gi (q),respectively. If we use n sample target molecular coordinates
q to perform the PCA, then i = 1, . . . , min(n,MNB ). Then
we can show that the PCA estimates FPCA
√ using −βf(r|q)
as follows: First, the eigenvectors Uki (rl ) r of the PCA are
orthonormal:
NB
M 


√
√
(Uki (rl ) r)(Ukj (rl ) r) = δi,j .

TABLE VII. Estimate of FPCA using√ −βfk (r|q)(k = 1, . . . ,M)by performing PCA using −βfk (rl |q) r(k = 1, . . . ,M and l =1, . . . ,NB ). The ith√ eigenfunction Uki (r) is√estimated from the ith
eigenvector Uki (rl ) r of PCA dividing by r. The ith eigenvalue
σi2 and principal component gi (q) of FPCA are estimated by the ith
eigenvalue σi2 and principal component gi (q) of PCA.
FPCA using−βfk (r|q)
k = 1, . . . ,M(C2)

k=1 l=1

These equalities converge to the orthonormal conditions of theeigenfunctions Uki (r) according toM  rc
k=1

Uki (r)Ukj (r) dr = δi,j

ith eigenfunction (vector)
ith eigenvalue
ith principal component

(C3)

0

for n,NB → ∞. Second, the eigenvalue σi2 of the PCA isexpressed as

NB
M 

√
√
σi2 = var
(Uki (rl ) r)(−βfk (rl |q) r) ,APPENDIX D: PERTURBATION EFFECTS BY A LINEARCOMBINATION OF TWO EIGENVECTORS

We consider a unit vector u = (u1 , . . . ,uM )T , which is alinear combination of the two eigenvectors ui and uj :u ≡ cos θ ui + sin θ uj .

k=1 l=1(C4)
which converges to the eigenvalue of the FPCA
σi2 = var β rc M
0Uki (r)fk (r|q) dr

(C5)

k=1

for n,NB → ∞. Finally, the principal component gi (q) of the
PCA is expressed as
NB
M 

√
gi (q) =(Uki (rl ) r)(−β(fk (rl |q)√−fk (rl |q)) r),

k=1

Uki (r)(fk (r|q) − fk (r|q)) dr

By using Eqs. (11), (16), and (17), the change in theconformational distribution of the target molecule induced byλ = δu becomesD(ρδu (q)||ρ(q)) = 12 δ 2 cos2 θ σi2 + sin2 θ σj2 + · · · . (D2)
By using Eqs. (9), (13), and (21), the ratio change of the target
molecular conformation q by the perturbation can be expandedasρδu (q)ln= δg(q) + · · · ,
(D3)
ρ(q)
g(q) ≡ cos θgi (q) + sin θgj (q).

(C6)

(D4)

By introducing n ≡ (cos θ, sin θ ), Eqs. (D1) and (D4) arerepresented by the projection to n as

which converges to the principal component of the FPCA
M  rc(D1)where we define

k=1 l=1

gi (q) = −βUki (r)σi2
gi (q)

PCA using√−βfk (rl |q) r
k = 1, . . . ,M
l = 1, . . . ,NB√
Uki (rl ) r2σigi (q)(C7)0

uk = n · (Uki ,Ukj )(D5)g(q) = n · [gi (q),gj (q)].

(D6)

for n,NB → ∞. Thus, we can estimate FPCA us= 1, . . . ,M) by performing PCA usinging −βfk (r|q)(k√−βfk (rl |q) r(k = 1, . . . ,M and l = 1, . . . ,NB ), which issummarized in Table VII.

Thus, uk and g(q) are obtained on the biplot by projection
onto n. In summary, the perturbation effects of the linear
combination of eigenvectors can be understood from the biplot
as well as each eigenvector.

[1] P. W. Fenimore, H. Frauenfelder, B. H. McMahon, and F. G.
Parak, Proc. Natl. Acad. Sci. USA 99, 16047 (2002).
[2] P. W. Fenimore, H. Frauenfelder, B. H. McMahon, and R. D.
Young, Proc. Natl. Acad. Sci. USA 101, 14408 (2004).
[3] K. Gunasekaran, B. Ma, and R. Nussinov, Proteins 57, 433(2004).
[4] K. Henzler-Wildman and D. Kern, Nature (London) 450, 964(2007).

[5] R. G. Smock and L. M. Gierasch, Science 324, 198 (2009).
[6] D. D. Boehr, R. Nussinov, and P. E. Wright, Nat. Chem. Biol.
5, 789 (2009).
[7] T. Ichiye and M. Karplus, Proteins 11, 205 (1991).
[8] P. H. Hünenberger, A. E. Mark, and W. F. van Gunsteren, J. Mol.
Biol. 252, 492 (1995).
[9] Y. Kong and M. Karplus, Structure 15, 611(2007).

026704-19

KOYAMA, KOBAYASHI, AND UEDA

PHYSICAL REVIEW E 84, 026704 (2011)[10] Y. Kong and M. Karplus, Proteins 74, 145 (2009).
[11] I. T. Jolliffe, Principal Component Analysis, 2nd ed. (Springer,
New York, 2002).
[12] A. Kitao, F. Hirata, and N. Go, Chem. Phys. 158, 447 (1991).
[13] A. E. Garcı́a, Phys. Rev. Lett. 68, 2696 (1992).
[14] A. Amadei, A. B. M. Linssen, and H. J. C. Berendsen, Proteins
17, 412 (1993).
[15] R. Abseher and M. Nilges, J. Mol. Biol. 279, 911 (1998).
[16] Y. Mu, P. H. Nguyen, and G. Stock, Proteins 58, 45 (2005).
[17] A. Altis, P. H. Nguyen, R. Hegger, and G. Stock, J. Chem. Phys.
126, 244111 (2007).
[18] Y. M. Koyama, T. J. Kobayashi, S. Tomoda, and H. R. Ueda,
Phys. Rev. E 78, 046702 (2008).
[19] J. O. Ramsay and B. W. Silverman, Functional Data Analysis,
2nd ed. (Springer, New York, 2006).
[20] S. Kullback, Information Theory and Statistics (Wiley & Sons,
New York, 1959).
[21] S. Amari and H. Nagaoka, Methods of Information Geometry(AMS and Oxford University Press, New York, 2000).
[22] T. M. Cover and J. A. Thomas, Elements of Information Theory,
2nd ed. (Wiley Interscience, New York, 2006).
[23] D. Ming and M. E. Wall, Proteins 59, 697 (2005).
[24] D. Ming and M. E. Wall, J. Mol. Biol. 358, 213 (2006).
[25] D. R. Brillinger, Ann. Inst. Stat. Math. 21, 215 (1969).
[26] D. Frenkel and B. Smit, Understanding Molecular Simulation,
2nd ed. (Academic Press, San Diego, 2002), Chap. 3.2.2.
[27] E. Małolepsza, B. Strodel, M. Khalili, S. Trygubenko, S. N.
Fejer, and D. J. Wales, J. Comput. Chem. 31, 1402 (2010).
[28] R. A. Horn and C. R. Johnson, Matrix Analysis (Cambridge
University Press, Cambridge, 1985).

[29] K. R. Gabriel, Biometrika 58, 453 (1971).
[30] J. C. Gower and D. J. Hand, Biplots (Chapman & Hall, London,
1996).
[31] D. A. Case et al., AMBER10 (University of California, San
Francisco, 2008).
[32] 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, andP. Kollman, J. Comput. Chem. 24, 1999 (2003).
[33] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey,and M. L. Klein, J. Chem. Phys. 79, 926 (1983).
[34] J.-P. Ryckaert, G. Ciccotti, and H. J. C. Berendsen, J. Comput.
Phys. 23, 327 (1977).
[35] S. Miyamoto and P. A. Kollman, J. Comput. Chem. 13, 952(1992).
[36] T. Darden, D. York, and L. Pedersen, J. Chem. Phys. 98, 10089(1993).
[37] W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graphics 14,
33 (1996).
[38] S. Honda, K. Yamasaki, Y. Sawada, and H. Morii, Structure 12,1507 (2004).
[39] M. M. Seibert, A. Patriksson, B. Hess, and D. van der Spoel,
J. Mol. Biol. 354, 173 (2005).
[40] D. Satoh, K. Shimizu, S. Nakamura, and T. Terada, FEBS Lett.
580, 3422 (2006).
[41] A. Suenaga, T. Narumi, N. Futatsugi, R. Yanai, Y. Ohno, N.
Okimoto, and M. Taiji, Chem. Asian J. 2, 591 (2007).
[42] [http://www-wales.ch.cam.ac.uk/∼wales/perm-prmtop.ff03.py].
[43] C. M. Bishop, Pattern Recognition and Machine Learning(Springer, New York, 2006), Chap. 12.1.4.
[44] D. E. Shaw et al., Commun. ACM 51, 91 (2008).

026704-20