Overview
In drug discovery, detailed knowledge of protein-ligand interactions are key to understand the recognition processes between the two molecules, and techniques such as Molecular Dynamics (MD) are valuable when simulating the behavior of the complexes using atomistic detail. Several aspects can affect a protein dynamic: mutations, binding of different ligands, protonation states (either titratable residues or ligands), composition of the lipidic bilayer (in the case of membrane proteins), and temperature, to name just a handful. These same aspects can also affect the stability of a complex.
Analyzing and presenting MD results to help guide project outcomes can be challenging. This is particularly true when several MD should be compared, as within large Free Energy Perturbation (FEP) experiments. A common practice is to use Principal Component Analysis (PCA) to reduce the dimensionality of the data while keeping the most significant information. This analysis is frequently used to evaluate the conformational space that a protein can access, providing additional details when compared to other commonly used analysis tools like RMSD, RMSF, and radius of gyration (among others). As exemplified below, extracting useful information from PCA data strongly relies on the availability of adequate visual tools.
Visualizing results with PCA plots
We can use a visual analysis of the PCA to gain better understanding of the conformations accessed by a MD simulation. With PCA, we can transform the 3D coordinates from all frames in a trajectory file to a set of linear orthogonal vectors, called Principal Components (PC). To picture this in a simple system, in Figure 1 every blue dot describes the x and y coordinates of one atom for every step in a trajectory file, and the red arrows will be representative of the first two PC for this system.
Figure 1. Representation of a x and y coordinates for one atom system (blue dots), and a qualitative representation of the corresponding first 2 eigenvectors obtained from a PCA.
Now, let’s imagine a tridimensional system (e.g. MD simulation of a protein) with a hypothetical outcome for the analysis where we obtained just one PC that explains 100% of the variance: this would imply that the collective motion of all atoms can be expressed as a linear transformation of one component (i.e., translation of a rigid protein). In a more realistic example, if the variances obtained for the first two PCs are 52% and 29%, then 81% of the protein dynamics can be explained by following the evolution of these two PCs in a 2D plot (Figure 2B).
Rigorously speaking, to accomplish this analysis, a 3N x 3N covariance matrix is generated using the 3D coordinates (where N is the number of atoms), and from the diagonalization of this matrix we obtain the eigenvalues (variance for each coordinate of every atom), with the PCs being the corresponding eigenvectors (arranged in descending order of variance). By projecting the original 3D coordinates over the first two PC, we generate a 2D map that is representative of the protein dynamics, as described by the two largest uncorrelated motions found on the sampled conformations. This means we could, for example, decide if the protein achieved conformational convergence during the simulated time, compare different conditions, or confirm if the equilibration stage was correctly set.
Revealing hidden patterns in protein dynamics with PCA plots
To illustrate this analysis in a practical example, a 50 ns MD simulation is analyzed using RMSD and PCA (Figure 2). According to the RMSD analysis (Figure 2 A), conformations at 10, 30 and 45 ns of simulated time would be equivalent (orange stars). On the PCA, a 2D plot for the first two PC is presented, where we can follow the simulated time using the color bar (Figure 2 B). Analogously, conformations sampled at 10, 30 and 45 ns of simulated time are marked with white stars. As we can see, the PCA indicates that these conformations are not equivalent, achieving structural convergence in the conformational space sampled in the range from 40 to 50 ns. In other words, in this example the evidence obtained from the PCA plots points to the protein exploring 3 macrostates, which is not captured with RMSD.
Figure 2. A) RMSD analysis: From this plot we have no strong evidence indicative that the protein presents a different conformation between 10 and 50 ns of simulation time. B) PCA results for the same system, obtained by projecting the coordinates of the protein trajectory over the first and second eigenvectors (PC1 and PC2). Data points are colored according to the simulation time (in frames), which allows tracking the evolution of the conformations sampled during the simulation. In this plot we obtain more comprehensive information regarding the conformational space sampled by the protein (i.e., in this example the evidence points to the protein exploring 3 macrostates). White stars denote the same simulated times as the orange stars in panel A.
From knowing our system to selection of new candidates
With the increase of computational power, it has become common practice to assess the stability of bound drug candidate molecules using MD, and in this next scenario, we will show with PCA how different situations can have a significant impact in the protein conformations being sampled. In Figure 3, the projection of eight trajectories used in the full PCA analysis assessing the same protein in different conditions is depicted as grey dots, while individual results for two of these trajectories are colored following the simulated time (same as in Figure 2B). In this example, we analyze the difference between the results obtained for the same dimeric protein in complex with the co-crystallized compound in both active sites (Figure 3A), or with the second binding site unoccupied (Figure 3B). The PCA was performed using the positions of protein residues belonging to the active site occupied in both situations. The PCA plots indicates that the conformational space explored in Figure 3A is narrower and closer to the initial conformation, while in the case displayed in Figure 3B, the protein is experiencing a noteworthy restructuring. We can thus conclude that an allosteric effect connecting the active site from both subunits is revealed, evidenced by the significantly different set of values obtained for PC1 and PC2.
Figure 3. Tracking of conformational sampling in two different 200 ns MD runs for a dimeric protein, using the 3D coordinates for all heavy atoms from one binding site in the PCA. A) Trajectory including the co-crystallized compound in the active site of both protomers. B) Trajectory obtained from an MD simulation where only one binding site is occupied. In both images data points for each situation are colored according to the simulation time (in frames), and the space explored with all eight trajectories used in the PCA calculation (including different conditions for the dimer) is depicted with grey dots.
PCA can support FEP experiments in the identification of outliers in a congeneric series. In Figure 4, we have analyzed the data from an FEP study involving 28 ligand structures. For each complex, the last frame of each trajectory was projected onto the PC map (magenta circles, Figure 4B) that is defined by the initial 200 ns MD simulation of the complex between the target and the reference lead compound (Figure 4A). The scientist undertaking the work could use this analysis to identify outlier structures described by PC1 and PC2 and evaluate the source of the difference observed (e.g., induced-fit effect, non-converged results, conformational noise due to an initial bad alignment, etc.). In this example, the four ligands displaying a larger deviation (Figure 4 B, green oval) present either more than one substitution point, or some of the bulkiest R-group replacements.
Figure 4. A) Heatmap generated with the first two PCs obtained with the MD run used as a reference (preceding step to the FEP calculation), and the 41 FEP trajectory structures obtained from the last frame of each transformation. B) All the conformational space covered is depicted in gray, and the projections for the structures obtained from FEP are colored in magenta. The area corresponding to the most populated region in the heatmap (A) is delimited with a white oval (B).
Visualize PCA analysis within Flare™
Within Flare, we provide a wide variety of tools for the analysis of MD trajectory results, and its analytical power can be significantly extended by utilizing pyflare. Pyflare enables editing and running python scripts within the Flare GUI. For the analysis presented here, we have combined the object-oriented python toolkit MDAnalysis1-2 with pyflare scripts written by Cresset. If you are interested in using these scripts for your research, please contact us and our support team will be happy to provide them.
PCA of dynamics trajectories will shortly be available in Flare V9. Stay tuned for the release announcement!
References
- N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. J. Comput. Chem. 32 (2011), 2319–2327. doi:10.1002/jcc.21787
- R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. MDAnalysis: A Python package for the rapid analysis of molecular dynamics simulations. In S. Benthall and S. Rostrup, editors, Proceedings of the 15th Python in Science Conference, pages 98-105, Austin, TX, 2016. SciPy. doi:10.25080/Majora-629e541a-00e
