Perturbational formulation of principal component analysis in molecular dynamics simulation

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 physical meaning clearly. Therefore, a formulation that provides general applicability and clearly represents the physical meaning 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 contribute to 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, and we 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 in vacuum 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 and negative components with larger absolute values of the first and second eigenvectors identify the electrostatic interactions, 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 the molecular conformational fluctuation but also to predict the conformational distribution change by the perturbation beyond the limitation of the previous methods.

Figure 1

Eigenvalues of the covariance matrix of the potential energy terms divided by $k_{B}T$. Eigenvalues are sorted in descending order. The inset figure shows the eigenvalues in log scale. First five eigenvalues are 461.5, 16.3, 9.3, 5.7, and 5.3.

Figure 2

(Color) The sign of (a) the first $\left[g_1\left(\mathbf{q}\right)\right]$ and (c) the second $\left[g_2\left(\mathbf{q}\right)\right]$ principal component on a Ramachandran plot. Red and blue points show $g_i\left(\mathbf{q}\right)⩾0$ and $g_i\left(\mathbf{q}\right)<0$, respectively. The dihedral angles $\varphi$ and $\psi$ 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) the first and (d) the second eigenvectors. Label el-6-18 indicates the component of the electrostatic potential energy term between atom 6 and atom 18. Such labels are captioned where the absolute value of the component is larger than 0.3.

Figure 3

Hyperplanes determining the probability change of $\mathbf{q}$ by the perturbation $\mathbf{\lambda }$ and contours of the Kullback-Leibler divergence. The solid line shows the hyperplane $\mathbf{\lambda }\cdot \mathbf{y}-\psi \left(\mathbf{\lambda }\right)=0$ which separates the increases or decreases of the probability change of $\mathbf{q}$ by the perturbation $\mathbf{\lambda }$. The solid and broken ellipses show the contours of the Kullback-Leibler divergence of $\varphi \left(\mathbf{y}\right)=C(\rho ,{\rho }_{\mathbf{\lambda }})$ and $\varphi \left(\mathbf{y}\right)=D({\rho }_{\mathbf{\lambda }}\parallel \rho )$, respectively.

Figure 4

Contours of the Kullback-Leibler divergence and eigenvectors. Within the second order approximation of the Kullback-Leibler divergence, the contours of the Kullback-Leibler divergence make ellipsoids and eigenvectors corresponding to the principal axes of ellipsoids [Eqs. (8, A6, A8)]. (a) $\mathbf{\lambda }$ (natural parameters in the exponential family) representation of the Kullback-Leibler divergence. The Kullback-Leibler divergence in the second order approximation with $\mathbf{\lambda }$ representation is Eq. (8). Axes ${\xi }_1$ and ${\xi }_2$ correspond to the first and second eigenvectors, respectively. $\mathbf{\lambda }=0$ corresponds to the unperturbed distribution. (b) ${⟨\mathbf{f}⟩}_{\mathbf{\lambda }}$ (expectation parameters in the exponential family) representation of the Kullback-Leibler divergence, i.e., rate function $\varphi \left({⟨\mathbf{f}⟩}_{\mathbf{\lambda }}\right)$ [Eq. (A1)]. The second order approximation of the Kullback-Leibler divergence with ${⟨\mathbf{f}⟩}_{\mathbf{\lambda }}$ representation is Eq. (A6). Axes $z_1$ and $z_2$ correspond to the first and second principal components, respectively [Eq. (A8)]. $\mathbf{y}=⟨\mathbf{f}⟩$ corresponds to the unperturbed distribution.