Gyration- and Inertia-Tensor-Based Collective Coordinates for

ARTICLE pubs.acs.org/JPCA

Gyration- and Inertia-Tensor-Based Collective Coordinates for Metadynamics. Application on the Conformational Behavior of Polyalanine Peptides and Trp-Cage Folding Jirí Vymetal†,‡ and Jirí Vondrasek*,† †

Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic, Flemingovo nam. 2, 166 10 Prague 6, Czech Republic ‡ Faculty of Natural Sciences, Physical Chemistry, Charles University in Prague, Hlavova 2030 Prague 2, 128 40 Prague 2, Czech Republic ABSTRACT: Effective simulations of proteins, their complexes, and other aminoacid polymers such as peptides or peptoids are critically dependent on the performance of the simulation methods and their ability to map the conformational space of the molecule in question. The most important step in this process is the choice of the coordinates in which the conformational sampling will be executed and their uniqueness regarding the capability to unambiguously determine the energy minimum on the free-energy hypersurface. In the presented study, we show that metadynamics and chosen collective coordinates—the principal moments of the tensors of gyration and inertia, the principal radii of gyration around the principal axes, asphericity, acylindricity, and anisotropy—can be used as a powerful combination to map the conformational space of peptides and proteins. We show that the combination of these coordinates with metadynamics produces a powerful tool for the study of biologically relevant molecules.

’ INTRODUCTION Metadynamics—a method for the calculation of the freeenergy profiles from nonequilibrium molecular simulations— has been introduced1 and intensively developed over the past decade. A comprehensive review of the historical development of metadynamics including all of its variants and numerous examples of its application has been provided by Laio and Gervasio.2 One of the most prominent applications of metadynamics is connected with studies of the dynamical behavior of proteins and their complexes spanning a large repertoire of their features. The conformational sampling and free-energy profiling along with a study of the protein-folding process seem to be the most frequent ones. The situation in the computation of the protein characteristics and protein dynamics is complicated by three separate factors: time, the quality of the force field, and the choice of the collective coordinates. The algorithm assumes that the system can be described by a few collective coordinates and the location of the system in the conformational space can be calculated. Apparently the choice of the collective coordinates is of utmost importance. Some authors have suggested the use of essential coordinates3 (coordinates determined by essential dynamics), but their attempt, others have suggested a relative. Other authors have suggested a relatively free choice of the collective coordinates depending only on the nature of the problem studied,4 but this attempt seems to be useless for any comparison of the general protein features for a large set of proteins. To simulate proteins or their complexes adequately, we have to find, sometimes tediously, the most appropriate coordinates in which the freeenergy profile is evaluated unambiguously. r 2011 American Chemical Society

The description of general polymer molecules by the radius of the gyration or principal components of the gyration tensor has a long history in macromolecular chemistry.5 Later, more sophisticated descriptors, e.g., asphericity, acylindricity, and relative shape anisotropy, were suggested and explored.6 The components of the gyration tensors and derived indexes allow an evaluation of the overall shape of a molecule (spherical, prolate, or oblate) and reveal its symmetry. Although the main field for their utilization lies in polymer science, their application on the characterization of RNA structures,7 proteins and DNA,8 protein-nucleic acid complexes,9 and even the denatured state of proteins10 have recently been reported. Last but not least, the molecular-shape descriptors together with the related WHIM (weighted holistic invariant molecular) descriptors11 are often employed in QSAR modeling. The second principal problem of metadynamics and its application for proteins and peptides is connected with the accuracy and physical soundness of the force field utilized. It is known that, in general, force fields perform well12 regarding noncovalent interactions, although it is strongly dependent on the way the partial charges are calculated. The principal problem of some force fields is probably in the parametrization of their covalent intramolecular parameters, particularly in the parametrization of the dihedral terms. As has been demonstrated recently, Special Issue: Pavel Hobza Festschrift Received: July 11, 2011 Revised: August 20, 2011 Published: September 30, 2011 11455

dx.doi.org/10.1021/jp2065612 | J. Phys. Chem. A 2011, 115, 11455–11465

The Journal of Physical Chemistry A

ARTICLE

the current force fields are too helical,11 and the only way to improve such inconsistency is a complete reparameterization of the dihedral term according to the benchmark ab initio calculations.

manner, we prefer to use these modified CCs: pffiffiffiffi Si 0 ¼ Si

ð3Þ

’ MATERIALS AND METHODS

b0 ¼

pffiffiffi b

ð4Þ

Model Systems. The applicability of the collective coordinates (CCs) was investigated on a set of short neutral alanine peptides blocked by acetyl (Ace) and N-methyl-amine (Nme) terminal groups, namely Ace-Ala3-Nme, Ace-Ala6-Nme, AceAla9-Nme, and Ace-Ala12-Nme. The Trp-cage mini-protein was chosen as a model molecule for the evaluation of the folding abilities of the newly implemented CCs. In the case of the Trpcage with the sequence Asn-Leu-Tyr-Ile-Gln-Trp-Leu-Lys-AspGly-Gly-Pro-Ser-Ser-Gly-Arg-Pro-Pro-Pro-Ser, the zwitter-ionic forms were simulated with a protonated N-terminal amino group and deprotonated C-terminal carboxyl group. For Asp, Lys, and Arg, the natural protonation in pH 7 was chosen, resulting in the net charge of +1 for the whole molecule. All of the molecules were built de novo using the tLeap program from the AmberTools package14 in the fully extended conformation. Collective Coordinates. For the purpose of this study, we used the following descriptors as the CCs in metadynamics: the principal moments of the tensors of gyration and inertia, the principal radii of the gyration around the principal axes, asphericity, acylindricity, and anisotropy. Both the tensors of inertia and gyration are related to the distribution of mass in the object. In the case of (macro)molecules and molecular clusters, the distribution of mass is represented by the respective positions of atoms. Therefore, the discrete definition of tensors is usual. The gyration tensor is defined as 2 3 x2i x i yi xi zi 16 7 y2i yi zi 7 ð1Þ S¼ 6 4 xi yi 5 N xi zi yi zi z2i

c0 ¼

pffiffi c

ð5Þ

∑ ∑ ∑

∑ ∑ ∑

∑ ∑ ∑

where the summation is performed over N atoms and the coordinates x, y, and z are related to the geometrical center of the molecule/cluster. Similarly, the tensor of inertia is defined as 2 3 mi ðy2i þ z2i Þ mi xi yi mi xi zi 6 7 mi xi yi mi ðx2i þ z2i Þ mi xi zi 7 I ¼6 4 5 mi xi zi mi yi zi mi ðx2i þ y2i Þ

∑

∑ ∑

∑

∑

∑

∑

∑ ∑

ð2Þ where the coordinates are related to the center of total mass and weighted by the mass of the individual atoms mi. Both tensors are symmetric of rank 3 and nonsingular for a nonlinear arrangement of atoms. The eigenvectors of the inertia tensor define the principal axes of inertia, i.e., special rotational axes, which are mutually orthogonal and absent of the products of inertia (nondiagonal elements). The corresponding eigenvalues I1, I2 and I3, called the principal moments of inertia, reflect the distribution of mass along the principal axes. (For details, see the Appendix.) All of the CCs meet the requirements and conditions of continuity and differentiability. From a practical point of view, we propose to use slightly adapted forms of some descriptors. We suggest applying such a transformation that converts the physical dimension of the descriptor to the dimension of length. The relation A3 can be considered as a case of this transformation. In the same

The reason for this modification is the fact that we are primarily interested in a compactness of the structures studied. Using the principal moments directly, the extended or diffusive forms would be sampled excessively just because of the choice of the quadratic form of the coordinates (physical unit of the moments is proportional to m2). Similar arguments may be applied also for asphericity and acylindricity. Implementation. We introduced the principal radii of gyration rgi moments of gyration tensor S0 i, asphericity b0 , acylindricity c0 , and relative shape anisotropy k2 as the CCs for metadynamics into PLUMED, a portable plugin for free-energy calculations with molecular dynamics.15 PLUMED provides an abstract interface for coding new extensions such as CCs; it is therefore very suitable as a starting point for developing new simulation methods. Because of PLUMED’s broad spectrum of supported computational codes, our implementation can be accessed in various MD packages including Gromacs, Amber, NAMD, DL-POLY, and ACEMD. Moreover, what is available thanks to PLUMED’s capabilities are not only metadynamics and its variants but also umbrella sampling and steered dynamics methods. From a technical point of view, the value of a collective coordinate CC(X) must be evaluated as a function of the Cartesian coordinates X as well as its Cartesian derivatives δCC(X)/δxi. All of the collective variables in question can be calculated from both tensors, S and I . We have chosen an evaluation of the gyration tensor, because some relations become more straightforward. The key point—the diagonalization of the tensor—employs a specialized routine devoted to the diagonalization of the real symmetric matrices of rank 3.16 The evaluation of the other collective variables from S0 1, S0 2, and S0 3 exploits the previously outlined relations. Computation of derivatives is facilitated by utilizing the chain rules. The diagonalization of the tensors may be interpreted as a rotational transformation, and this fact is also used in projecting the derivatives in the Cartesian coordinate system. In order to provide a general and more flexible code, our implementation supports two methods of weighting the atomic contributions. Mass-weighted as well as nonweighted versions are available in PLUMED for all of the newly introduced CCs. Simulation Details. All of the presented results were obtained using Gromacs (version 4.5.3)17 enhanced by PLUMED15 (version 1.2) with newly implemented CCs. The ff03 force field18 was used in all of the simulations provided in the Gromacs package. In order to include the solvation effects, the implicit generalized Born model (with εr = 80) was utilized with the Still method of the calculation of Born radii19 as implemented in Gromacs. No cutoff was employed for the computation of the electrostatic and other nonbonded calculations; all of the pair interactions were included explicitly. The integration step of 2 fs was used along with the constraint of all of the bond lengths by means of the LINCS algorithm. In order to avoid the “flying ice cube” artifact, the total translational and rotational moments were regularly removed every 10 steps. All of the simulations were performed at a temperature of 300 K controlled 11456

dx.doi.org/10.1021/jp2065612 |J. Phys. Chem. A 2011, 115, 11455–11465

The Journal of Physical Chemistry A Table 1. Protocols Used in BIASXMD Simulations Protocol

CCs

P1

Rg, S10 , S20 , S30 , b0 , c0 , k2, 1 neutral replica

P2

Rg, S10 , S20 , S30 , rg 1, rg 2, rg 3, 1 neutral replica

P3

Rg, S10 , S20 , S30 , 4 neutral replicas

P4

Rg,rg 1, rg 2, rg 3, 4 neutral replicas

by a Langevin stochastic integrator as implemented in Gromacs with a friction constant of 0.1 p s1. The simulation of all of the systems in this study was conducted in a way parallel to multiple-replica bias-exchange metadynamics (BIASXMD)20 including eight replicas. Exchanges between replicas were allowed according to the rule taken from the paper of Piana and Laio.20 The pairs of replicas, between which the exchanges were attempted, were selected randomly. Therefore, two basic protocols were devised to comprise all of the newly implemented CCs. Table 1 summarizes the CCs included in all of the protocols and the differences between them. The protocols P1 and P2 were designed for simulations of both peptides and Trp-cage. The additional protocols P3 and P4 were tested on folding simulations for Trp-cage. Each replica was biased by one one-dimensional collective variable, except for neutral replicas, where only Langevin dynamics without any bias were performed. The exchanges between the replicas were attempted by PLUMED each 2 ps (1000 steps) of simulated time. Each BIASXMD simulation of the alanine peptides was 100 ns per replica long, and the folding simulations of Trp-cage were 200 ns long. Metadynamics in our simulation were accomplished by the PLUMED plugin extended by our implementation of gyration and inertia tensor-based CCs. Specifically, a well-tempering algorithm was chosen with a “BIASFACTOR” parameter of 17.66, a starting height of the bias potential of 0.25 kJ/molm and a deposit period of 1 ps (500 steps). The width of 0.005 nm of the bias potential (the “SIGMA” parameter) was used for all of the CCs with a physical dimension of length (Rg, S10 , S20 , S30 , rg1, rg2, rg3, b0 , c0 ). For the dimensionless relative shape anisotropy k2, the value of 0.05 was selected. Only the mass-weighted variants of the newly implemented CCs were utilized in this study. The ability of PLUMED to evaluate the bias potential on a grid was exploited to accelerate the force calculation. Construction of Free-Energy Profiles. The free-energy profiles were calculated by reversing the bias potential from metadynamics. Time-domain averaging was used in the evaluation of the free energy from each biased replica according to the following simple algorithm: After the deposit of the bias potential with a height of 50 kJ/mol, the averaging started. The next 1000 gaussians were not summed directly with the previous bias potential but used as a correction. Although well-tempered metadynamics guarantee the final convergence of the bias potential, they must be tuned in the opposite way relative to the quick convergence to overcome high barriers. The proposed algorithm of the time averaging does not interfere with the self-convergence of well-tempered metadynamics, because both the averaged and nonaveraged free-energy profiles must provide the same surface in a long enough time. Furthermore, the free-energy profiles from the independent simulations were compared by means of a basis statistical analysis. The average profiles are provided as final results with their reliability estimates based on the variance between the individual profiles.

ARTICLE

Structural Analysis. The assignment of the secondary structures was performed by the DSSP program.21 Only the second half (50 ns) of each neutral replica was used for the analysis. The basic statistical treatments followed the above-mentioned procedure. The root-mean-square devistion (rmsd) was calculated by the g_rms utility from the Gromacs package. A visualization of the molecules was rendered by Pymol,22 and the plots and other graphics were elaborated using the matplotlib23 library.

’ RESULTS Simulation of Short Alanine Peptides. Short neutral alanine peptides with blocked terminal groups, namely, Ace-Ala3-Nme (Ala3), Ace-Ala6-Nme (Ala6), Ace-Ala9-Nme (Ala9), and AceAla12-Nme (Ala12), were investigated by means of multiplereplica BIASXMD. Two protocols, P1 and P2 (described in the Materials and Methods section), distinctive in the CCs used, were utilized simultaneously in order to involve all of the new implemented CCs. In total, 2 6 simulations were performed for each peptide, which allowed the statistical and critical evaluation of the results. The convergence was monitored, and the differences in the profiles from the individual independent runs served for the estimation of the statistical uncertainty of our results. The freeenergy profiles from all of the independent simulations are plotted in Figure 1, where the curves are compared to each other to demonstrate the robustness of the presented results. Naturally, the smallest deviations between individual runs were observed for shorter peptides (Ala3 and Ala6) while growing gradually for the longer ones (Ala9 and Ala12). The numerical estimates of standard deviation for individual free-energy surfaces on Figure 1 are 0.2 kJ.mol1 for Ala3, 0.2 kJ.mol1 for Ala6, 0.5 kJ.mol1 for Ala9, and 0.7 kJ.mol1 for Ala12. There was no difference in performance between the simulation protocols P1 and P2. Both P1 and P2 sampled the conformation space of alanine peptides comparably and provided equivalent free-energy profiles as well as the other investigated properties mentioned further. The one-dimensional free-energy profiles as a function of the individual employed CCs were obtained by postprocessing the bias potential of the corresponding replicas. The resulting profiles are plotted in Figure 2. The fact that almost all of the free-energy profiles in Figure 2 are very smooth without distinctive minima clearly manifests the intrinsic flexibility of the peptides or CCs. The only exception is the shortest peptide in our study Ala3, for which two or three minima on all of the free-energy profiles could be localized. We interpret this observation as the existence of well-recognizable conformational families. In the case of the longer peptides (Ala6, Ala9, Ala12), the continuum of the conformational states (in the sense of the collective variables used in Figure 2ag) was sampled. For their further characterization, the asphericity, acylindricity, and relative shape anisotropy (Figure 2hj) were suitable tools. The systematic shifts of the main minima in Figure 2h for Ala6, Ala9, and Ala12 suggest a deviation from spherical symmetry. This fact is also confirmed by the systematic shifts in the relative shape anisotropy (Figure 1j) toward the helices. Finally, the free-energy profiles in acylindricity (Figure 2j) showed common minima for its relatively small value, hence the structures possess high cylindrical symmetry. Such an interpretation can also be directly deduced from the profiles in the moments of gyration tensor and the principle radii of gyration (Figure 2b2g). These analyses clearly confirm the 11457



ARTICLE

Figure 1. Accuracy of the obtained free-energy profiles demonstrated on the radius of gyration. The red curve denotes the average free-energy profiles, and the black curves were derived from independent runs of BIASXMD. The variance between the profiles grows toward the longer sequences.

main global minimum for the Ala6, Ala9 and Ala12 as the helical conformer. However, the considerable structural deviations are energetically feasible according to the presented free-energy profiles. The indication of a free-energy barrier between the helical and collapsed conformation is elucidated by relative-shape anisotropy as a reaction coordinate. For Ala6, the collapsed structure seems to be iso-energetical while the helical state is strongly preferred by longer peptides. We have to also add that the systematic shifts in the allowed relative shape anisotropy parameter toward helices can be partially caused by the fact that α-helices are preferred conformers for alanine peptides in ff03 force field combined with the GBSA (Still radii) implicit solvation model. We should also add that using of implicit solvent may shift the bias known from explicit simulations. The prevalence of helical structures was also confirmed by visual inspections of the trajectories and the analysis of the secondary structure elements. The analyses of the secondary structures in peptides Ala6, Ala9, and Ala12 were performed on neutral replica in each independent BIASXMD run, which approximates the classical canonical distribution of energies. The comprehensive results are plotted in Figure 3. The conformation preferences for secondary structures distinguished by DSSP were calculated for every residue in the peptide (Figure 3A,C,E). Moreover, the cooperativity of the residues to create secondarystructure elements (helix and turn) was explored as shown in Figure 3B,D,E. The growing preferences for forming an α helix are noticeable for all of the residues in the peptide on its elongation. The statistically most helix-favoring residues lie in the middle of the sequence, and the propensities for helicity drop toward the boundaries. The opposite trends were observed for the preferences of residues to form hydrogen-bonded turns. Figure 3B,D,F shows the probability of a simultaneous occurrence of the given number of residues in a helix or turn. Panels D and F clearly show that the full-length helix is formed only with a probability of 30% or 20% for ala9 and ala12, respectively. Furthermore, this trend is not linear, but a significant preference for a four-residue helix was

found. On the contrary, five and six helical residues are observed less frequently for both Ala9 and Ala12; in the case of Ala12, it is nine residues. A similar nonlinear dependence is found for hydrogen-bonded turns, which are usually populated in lesser counts. Ala9 and Ala12 share the same trend, but Ala6 exhibits different proportions between the doubles and triples of the residues involved in hydrogen-bonded turns. The presented results were obtained from all of the available neutral replicas and possess high statistical reliability because of the small deviations between individual runs. The error bars are not shown for panels A, C, and E for the sake of lucidity, but the standard deviation does not exceed 3% (typically below 1%). The cluster analysis of the conformers from neutral replica was performed with the aim of illustrating the representative structures and the discriminative ability of gyration tensor-based CCs. Because of the ambiguous criteria and poor performance of the common cluster algorithm on continuous sets of conformers, the choice of representative structures was more or less arbitrary. For the sake of visualization, the multidimensional descriptor of the conformation state composed from all of the individual CCs was reduced to a two-dimensional representation by means of a principal component analysis. Therefore, the resulting principal descriptors PC1 and PC2 cover the majority of the structural variance. Figures 47 depict this conformational space sampled in neutral replica for all of the studied peptides projected to PC1 and PC2. The selected structures of the peptides are drawn along with their position in the projected subspace. Also the corresponding values of CCs are listed in attached tables for each structure. The picture of the sampled conformational space for each peptide from Figures 47 completes the energetic picture in Figure 1. For the shortest peptide Ala3, Figure 4 shows three clusters that correspond to the minima on the free-energy profiles. In the case of the longer peptides, namely, Ala6, Ala9 and Ala12, the clustering in the projected space is less traceable and one dense region dominates. The analysis of the trajectories 11458



ARTICLE

Figure 2. The free-energy profiles of the investigated alanine peptides along the gyration radius-based CCs.

revealed that the prevalent regions contain helical conformers. The distance in the projected space from the ideal helical conformer (depicted as Structure 1 in Figures 5,6,7) reflects the deviation from a regular helix. The existence of one prevalent region is demonstrated on the free-energy profile as a single broad global minimum.

In the case of the short alanine peptides, we aimed not only for an extensive conformational sampling but also for obtaining precise free-energy profiles along the CCs and an unbiased distribution of the other molecular descriptors. We were able to obtain free-energy profiles with controlled precision. Moreover, using a neutral replica in the BIASXMD simulations made it 11459



ARTICLE

Figure 3. Analysis of the secondary-structure elements in alanine peptides. The propensity of each residue to a different secondary structure is plotted for Ala6 (A), Ala9 (C), and Ala12 (E) peptides. The cooperativity of the helices and turns is expressed as a probability of the simultaneous occurrence of N residues in the given conformational state.

possible to approximate the canonical equilibrium distribution with an enhanced sampling of the conformers because of its linking with metadynamics. For this reason, we succeeded in a reliable estimation of the distribution of the other structural descriptors that were not explicitly considered as collective variables in metadynamics. Folding Simulation of Trp-Cage. The series of the BIASXMD simulation were conducted in order to simulate the folding of the Trp-cage mini-protein. In addition to the P1 and P2 simulation protocols, we employed two extra protocols, P3 and P4, with the aim of enhancing the folding event of Trp-cage. The typical course of the performed folding simulation is shown in Figure 8. The rmsd between the native structure and the continuous trajectories obtained from the exemplary BIASXMD simulation executed in P1P4 protocols reveals the folding event. When the native-like fold with the lowest rmsd

(approximately 0.1 nm) was reached, it usually remained preserved to the end of the simulation. It is worth mentioning that an unfolding event is possible as well. However, some other structures seem to also be sufficiently stable: they possessed higher rmsd (typically 0.20.4 nm). It is questionable how far from a native state they are. The overall efficiency of all four of the folding protocols is depicted in Figures 9 and 10. Both panels show the distribution of the rmsd between structures sampled in all of the simulated neutral replicas (the native state should be trapped there) and the native structure from the PDB database (1L2Y). We compared both deviations for the main chain atoms (Figure 9) and in the residues involved in the protein core (Figure 10). This confrontation reveals that the native state was found almost in all of the simulations, because such structures possess low acceptable structural deviations for both conformations: the main chain and 11460



ARTICLE

Figure 4. Ala3 conformers.




the inner protein core. The only exceptions are two unproductive simulations conducted with the P1 and P2 protocols, which is a reason to improve the sampling strategy. All of the simulations performed under the P3 and P4 protocols were able to sample native-like conformers. However, the addition of three other neutral replicas in our BIASXMD scheme resulting in enhanced protocols 3 and 4 did not bring any substantial increase in the population of the native state (cf. the distributions of the rmsd on panels A, B, C, and D in Figures 9 and 10. The population of the

lower rmsd did not rise noticeably in the case of P3 and P4). Therefore, only the probability of the folding events seems to be increased by the presence of an extra amount of neutral replicas. The sampled native-like conformers are depicted in Figure 11. The experimental structure is superimposed with the most similar structures (in the sense of the rmsd between the main chain) obtained from all of our BIASXMD simulations. The rmsd of the main chain does not exceed 0.1 nm; moreover, the overwhelming majority of conformers possess the correct packing of 11461



ARTICLE

Figure 8. Typical rmsd from the native state in the course of a BIASXMD simulation.

Figure 9. Fraction of the folded structures. The rmsd of the main chain atoms between structures occurring in neutral replica and the native structure (1L2Y) was binned and plotted as histograms. The distributions from all of the simulations performed by means of four different protocols (A, B, C, D) are shown. The occurrence of structures with low rmsd proves that native-like conformers are sampled.

the inner core formed by Tyr3, Trp6, Pro12, Pro17, Pro18, and Pro19. The free-energy profiles along the utilized CCs derived from BIASXMD showed high variance between individual runs. Therefore, any free-energy profiles are presented. A much more extensive BIASXMD simulation would have had to be performed in order to reach reliable converged freeenergy profiles.

We expect that a substantial improvement of folding efficiency could be achieved by an extension of the BIASXMD protocol by additional protein-specific CCs and by a further optimization of the parameters involved in the simulation (the number of replicas, frequency of exchanges between the replicas, width and height of the gaussians, and their deposit rate). Such a complex and demanding procedure has only recently been addressed for BIASXMD.24 11462



ARTICLE

Figure 10. Fraction of folded structures. The rmsd of the residues in the hydrophobic core of Trp-cage between structures sampled in neutral replicas and the native structure (1L2Y). Four different protocols are compared (A, B, C, D). The non-negligible fractions of low rmsd indicate the sampling of a conformation with a native-like hydrophobic core.

Figure 11. Superposition of Trp-cage structures. The native state (ribbon) 1L2Y is compared with structures obtained from the simulations. From each BIASXMD run, the most similar structure was selected, aligned, and drawn along the native structure.

’ DISCUSSION The gyration- and inertia-tensor-based CCs are capable of an overall geometrical description of the molecules as follows from their definition. We cannot expect them to capture the local structural details, which do not influence the total three-dimensional shape of the molecule. We may ask how much information about the shape is stored in the gyration or inertia tensor. Both tensors are symmetrical and may be diagonalized in a suitably rotated coordinate system. Therefore, only the values on the diagonals bear information on the molecule independent of its orientation in space. This trio of values bears the maximal amount of information available. However, metadynamics utilizing the three-dimensional coordinate

begins to be rather ineffective, and the acquisition of properly converged free-energy profiles becomes troublesome. Therefore, the one- or two-dimensional reduction is required with a minimal loss of useful information. The Suter’s geometry descriptors—asphericity, acylindricity, and relative-shape anisotrophy—fulfill these requirements. All of them have well-defined geometrical meanings. A simultaneous use of all of the CCs may be questionable because of their mutual dependency and redundancy. However, metadynamics generally as a method does not suffer from the utilization of dependent CCs. Since the application of multidimensional CCs is prohibitive, the complete information from the tensor cannot be recovered without loss even by exploiting all of the one-dimensional CCs as we did in this study. A better and still feasible description of the system in question could be achieved by employing two-dimensional metadynamics. A suitable combination of CCs may reveal a large portion of the correlation between the individual components of the gyration or inertia tensor. Such a correlation would significantly improve the knowledge about the free-energy profile of the molecule and reveal more subtle details of the geometrical preferences. Concerning the efficiency of our simulation, no comprehensive tuning was intended. We used the same parameters for all of the simulations assuming their robustness. The only optimization in our study was introducing additional P3 and P4 protocols with more neutral replicas. We expected an improved folding rate because of the greater relaxation mediated by classical dynamics in a neutral replica, which is believed (by us) to facilitate the local arrangement of the backbone and core packing toward the native state. According to our results from P3 and P4 protocols, these prospects were not met or had only a minor effect on the overall performance of the folding simulation. However, we still expect that a significant improvement in performance will be achieved by the combination of global shape descriptors from this study 11463



ARTICLE

with other CCs designed for protein folding such as backbone helicity, number of contacts, dihedral correlation or the content of secondary structures. The other factors that inevitably affect the quality of the results in this study are the force field and solvation model. Most of the current force fields for proteins and peptides including ff03 were found to favor the helix excessively.25 In spite of the unbalanced propensity for secondary structures, successful applications of ff03 on Trp-cage folding were reported.26,27 Moreover, the justification for implicit solvents to describe the solvation of different peptide conformers, protein dynamics,28,29 and protein folding30 correctly remains no less controversial. None of these issues have been addressed in this study.

’ CONCLUSION The investigated CCs for metadynamics are able to unambiguously distinguish between the conformers of only small peptides such as Ala3. The flexibility of longer peptides and the number of possible conformers inevitably smudge the resolution of such global shape descriptors. However, the geometrically distinctive conformers such as helices are still easily recognizable, especially by means of CC acylindricity and relative-shape anisotropy. Our results suggest that CCs used in our study may facilitate successful folding of Trp-cage mini-protein, but not the complete folding process. Some key aspects of the folding pathway do not seem to be accelerated by global shape CCs, and such events are a bottleneck of our simulations. The global character of the newly implemented CCs in PLUMED makes them attractive for protein-folding simulation. They naturally facilitate the change of shape of flexible molecules. In the course of folding simulations, they provide collapsing to compact structures as well as loosening to extended forms. ’ APPENDIX: THEORY This is a summary of the well-known relations between tensors of gyration and inertia, radii of gyration, and shape descriptors developed by Suter. For additional details, see the original paper.6 If we assume the same mass of all of the atoms for inertia tensor or alternatively weigh all of the atomic contributions to the gyration tensor by the mass of the individual atom, then both tensors can be diagonalized in the same reference frame. 2 3 S1 0 0 6 7 7 ðA1Þ V 1 SV ¼ Sdiag ¼ 6 4 0 S2 0 5 0 0 S3 2

V 1 IV ¼ I diag

I1 6 ¼6 40 0

0 I2 0

3 0 7 07 5 I3

ðA2Þ

We are obeying the convention of indexing the eigenvalues according to their magnitudes. Since both tensors are positivesemidefinite, all eigenvalues must be larger than or equal to zero. Index 1 indicates the largest eigenvalue, and index 3 represents the smallest. The eigenvectors of the gyration tensor can be interpreted as the geometrical axes of a rotational ellipsoid that effectively approximate the average spatial distribution of atoms with respect to the center of a molecule/cluster. The eigenvalues—the

principal moments of the gyration tensor S1, S2, S3—then correspond to the squares of the length (moments) of the associated elliptic axes. An equivalent way of expressing the distribution of mass is to use the radius of gyration around the given axes pffiffiffiffiffi I ax ðA3Þ rgax ¼ m where Iax stands for the moment of inertia around a certain axis, and m is for the total mass of an object. Unfortunately, the similarly named scalar quantity is routinely used to measure size, i.e., the radius of the gyration: qffiffiffiffiffiffiffiffiffiffiffiffiffi mi ri2 ðA4Þ Rg ¼ m

∑

Here the atomic distances ri are related to the center of mass. We note that the nonweighted form of Rg is also commonly used in the statistical physics of polymers. qffiffiffiffiffiffiffiffi ri2 ðA5Þ Rg ¼ N Assuming the consistent weighting of the atomic contribution in the previously defined quantities, the relations between them can be easily demonstrated. Dealing with particles with the same mass mi or alternatively weighting the particles’ contributions to the gyration tensor by miN/m, the proportion between S and I become obvious

∑

I ¼ mðRg2 Eð3Þ SÞ

ðA6Þ

E(3) states the identity matrix of rank 3. The eigenvalues of both tensors are therefore linked I1 ¼ mðS2 þ S3 Þ, I2 ¼ mðS1 þ S3 Þ, I3 ¼ mðS1 þ S2 Þ

ðA7Þ

and so are their traces ðA8Þ

TrI ¼ 2mTrS ¼ 2mRg2

The principal radii of gyration around the principal axes follow pffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi I1 ¼ S2 þ S3 ðA9Þ rg1 ¼ m and finally obey Rg2

ðA10Þ 2 Suter proposed additional descriptors in order to emphasize and distinguish between certain shape properties of the objects. All are defined in terms of the principal components of the gyration tensor. Asphericity b measures the deviation from the spherical symmetry: 2 2 2 rg1 þ rg2 þ rg3 ¼

1 b ¼ S1 ðS2 þ S3 Þ 2

ðA11Þ

Similarly, acylindricity c expresses the divergence of the cylindrical symmetry: c ¼ S2 S3

ðA12Þ

The last descriptor, the relative shape anisotropy k , reflects both the symmetry and dimensionality of an object. Limited between 2

11464



ARTICLE

values of 0 and 1, k2 reaches 1 for an ideal linear arrangement and drops to zero in the case of highly symmetric configurations (at least tetrahedral symmetry). For planar symmetric objects, anisotropy converges to the value of 1/4 S1 S2 þ S1 S3 þ S2 S3 k2 ¼ 1 3 ðS1 þ S2 þ S3 Þ2

ðA13Þ

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ REFERENCES (1) Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99 (20), 12562–12566. (2) Laio, A.; Gervasio, F. L. Metadynamics: A method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science. Rep. Prog. Phys. 2008, 71 (12). (3) Spiwok, V.; Kralova, B.; Tvaroska, I. Continuous metadynamics in essential coordinates as a tool for free energy modelling of conformational changes. J. Mol. Model. 2008, 14 (11), 995–1002. (4) Henin, J.; Fiorin, G.; Chipot, C.; Klein, M. L. Exploring multidimensional free energy landscapes using time-dependent biases on collective variables. Journal of Chemical Theory and Computation. 2010, 6 (1), 35–47. (5) Solc, K. Shape of a random-flight chain. J. Chem. Phys. 1971, 55 (1), 335–&. (6) Theodorou, D. N.; Suter, U. W. Geometrical considerations in model systems with periodic boundaries. J. Chem. Phys. 1985, 82 (2), 955–966. (7) Hyeon, C.; Thirumalai, D. Unfolding RNA secondary structure by force. Biophys. J. 2004, 86 (1), 313A. (8) Rawat, N.; Biswas, P. Size, shape, and flexibility of proteins and DNA. J. Chem. Phys. 2009, 131, 16. (9) Rawat, N.; Biswas, P. Shape, flexibility and packing of proteins and nucleic acids in complexes. Phys. Chem. Chem. Phys. 2011, 13 (20), 9632–9643. (10) Dima, R. I.; Thirumalai, D. Asymmetry in the shapes of folded and denatured states of proteins. J. Phys. Chem. B 2004, 108 (21), 6564–6570. (11) Todeschini, R.; Lasagni, M. New molecular descriptors for 2D and 3D structures - Theory. J. Chemom. 1994, 8 (4), 263–272. (12) Berka, K.; Laskowski, R. A.; Hobza, P.; Vondrasek, J. Energy matrix of structurally important side-chain/side-chain interactions in proteins. J. Chem. Theory Comput. 2010, 6 (7), 2191–2203. (13) Best, R. B.; Buchete, N. V.; Hummer, G. Are current molecular dynamics force fields too helical?. Biophys. J. 2008, 95 (1), L7–L9. (14) Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. The Amber biomolecular simulation programs. J. Comput. Chem. 2005, 26 (16), 1668–1688. (15) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M. PLUMED: A portable plugin for free-energy calculations with molecular dynamics. Comput. Phys. Commun. 2009, 180 (10), 1961–1972. (16) Kopp, J. Efficient numerical diagonalization of hermitian 3 3 matrices (vol 19, pg 523, 2008). International Journal of Modern Physics C 2008, 19 (5), 845. (17) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 2008, 4 (3), 435–447. (18) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G. M.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang,

J. M.; Kollman, P. A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations. J. Comput. Chem. 2003, 24 (16), 1999–2012. (19) Qiu, D.; Shenkin, P. S.; Hollinger, F. P.; Still, W. C. The GB/SA continuum model for solvation. A fast analytical method for the calculation of approximate Born radii. J. Phys. Chem. A 1997, 101 (16), 3005–3014. (20) Piana, S.; Laio, A. A bias-exchange approach to protein folding. J. Phys. Chem. B 2007, 111 (17), 4553–4559. (21) Kabsch, W.; Sander, C. Dictionary of protein secondary structure - Pattern-recognition of hydrogen-bonded and geometrical features. Biopolymers 1983, 22 (12), 2577–2637. (22) Delano, W. L. PyMOL molecular viewer: Updates and refinements. In Abstracts of Papers of the American Chemical Society; American Chemical Society: Washington, DC, 2009; Vol. 238. (23) Hunter, J. D. Matplotlib: A 2D graphics environment. Comput. Sci. Eng. 2007, 9 (3), 90–95. (24) Cossio, P.; Marinelli, F.; Laio, A.; Pietrucci, F. Optimizing the performance of bias-exchange metadynamics: Folding a 48-residue LysM domain using a coarse-grained model. J. Phys. Chem. B 2010, 114 (9), 3259–3265. (25) Best, R. B.; Buchete, N. V.; Hummer, G. Are current molecular dynamics force fields too helical?. Biophys. J. 2008, 95 (1), L7–L9. (26) Beck, D. A. C.; White, G. W. N.; Daggett, V. Exploring the energy landscape of protein folding using replica-exchange and conventional molecular dynamics simulations. J. Struct. Biol. 2007, 157 (3), 514–523. (27) Paschek, D.; Nymeyer, H.; Garcia, A. E. Replica exchange simulation of reversible folding/unfolding of the Trp-cage miniprotein in explicit solvent: On the structure and possible role of internal water. J. Struct. Biol. 2007, 157 (3), 524–533. (28) Zhou, R. H.; Berne, B. J. Can a continuum solvent model reproduce the free energy landscape of a beta-hairpin folding in water?. Proc. Natl. Acad. Sci. U.S.A. 2002, 99 (20), 12777–12782. (29) Fan, H.; Mark, A. E.; Zhu, J.; Honig, B. Comparative study of generalized Born models: Protein dynamics. Proc. Natl. Acad. Sci. U. S. A. 2005, 102 (19), 6760–6764. (30) Rhee, Y. M.; Sorin, E. J.; Jayachandran, G.; Lindahl, E.; Pande, V. S. Simulations of the role of water in the protein-folding mechanism. Proc. Natl. Acad. Sci. U.S.A. 2004, 101 (17), 6456–6461.

11465


Gyration- and Inertia-Tensor-Based Collective Coordinates for

Recommend Documents