Correction: An Exact Expression to Calculate the Derivatives of Position-Dependent Observables in Molecular Simulations with Flexible Constraints.

[This corrects the article DOI: 10.1371/journal.pone.0024563.].

Predicting stabilizing mutations in proteins using Poisson-Boltzmann based models: study of unfolded state ensemble models and development of a successful binary classifier based on residue interaction energies.

In many cases the stability of a protein has to be increased to permit its biotechnological use. Rational methods of protein stabilization based on optimizing electrostatic interactions have provided some fine successful predictions. However, the precise calculation of stabilization energies remains challenging, one reason being that the electrostatic effects on the unfolded state are often neglected. We have explored here the feasibility of incorporating Poisson-Boltzmann model electrostatic calculations performed on representations of the unfolded state as large ensembles of geometrically optimized conformations calculated using the ProtSA server. Using a data set of 80 electrostatic mutations experimentally tested in two-state proteins, the predictive performance of several such models has been compared to that of a simple one that considers an unfolded structure of non-interacting residues. The unfolded ensemble models, while showing correlation between the predicted stabilization values and the experimental ones, are worse than the simple model, suggesting that the ensembles do not capture well the energetics of the unfolded state. A more attainable goal is classifying potential mutations as either stabilizing or non-stabilizing, rather than accurately calculating their stabilization energies. To implement a fast classification method that can assist in selecting stabilizing mutations, we have used a much simpler electrostatic model based only on the native structure and have determined its precision using different stabilizing energy thresholds. The binary classifier developed finds 7 true stabilizing mutants out of every 10 proposed candidates and can be used as a robust tool to propose stabilizing mutations.

The Quixote project: Collaborative and Open Quantum Chemistry data management in the Internet age.

Computational Quantum Chemistry has developed into a powerful, efficient, reliable and increasingly routine tool for exploring the structure and properties of small to medium sized molecules. Many thousands of calculations are performed every day, some offering results which approach experimental accuracy. However, in contrast to other disciplines, such as crystallography, or bioinformatics, where standard formats and well-known, unified databases exist, this QC data is generally destined to remain locally held in files which are not designed to be machine-readable. Only a very small subset of these results will become accessible to the wider community through publication.In this paper we describe how the Quixote Project is developing the infrastructure required to convert output from a number of different molecular quantum chemistry packages to a common semantically rich, machine-readable format and to build respositories of QC results. Such an infrastructure offers benefits at many levels. The standardised representation of the results will facilitate software interoperability, for example making it easier for analysis tools to take data from different QC packages, and will also help with archival and deposition of results. The repository infrastructure, which is lightweight and built using Open software components, can be implemented at individual researcher, project, organisation or community level, offering the exciting possibility that in future many of these QC results can be made publically available, to be searched and interpreted just as crystallography and bioinformatics results are today.Although we believe that quantum chemists will appreciate the contribution the Quixote infrastructure can make to the organisation and and exchange of their results, we anticipate that greater rewards will come from enabling their results to be consumed by a wider community. As the respositories grow they will become a valuable source of chemical data for use by other disciplines in both research and education.The Quixote project is unconventional in that the infrastructure is being implemented in advance of a full definition of the data model which will eventually underpin it. We believe that a working system which offers real value to researchers based on tools and shared, searchable repositories will encourage early participation from a broader community, including both producers and consumers of data. In the early stages, searching and indexing can be performed on the chemical subject of the calculations, and well defined calculation meta-data. The process of defining more specific quantum chemical definitions, adding them to dictionaries and extracting them consistently from the results of the various software packages can then proceed in an incremental manner, adding additional value at each stage.Not only will these results help to change the data management model in the field of Quantum Chemistry, but the methodology can be applied to other pressing problems related to data in computational and experimental science.

An exact expression to calculate the derivatives of position-dependent observables in molecular simulations with flexible constraints.

In this work, we introduce an algorithm to compute the derivatives of physical observables along the constrained subspace when flexible constraints are imposed on the system (i.e., constraints in which the constrained coordinates are fixed to configuration-dependent values). The presented scheme is exact, it does not contain any tunable parameter, and it only requires the calculation and inversion of a sub-block of the Hessian matrix of second derivatives of the function through which the constraints are defined. We also present a practical application to the case in which the sought observables are the Euclidean coordinates of complex molecular systems, and the function whose minimization defines the flexible constraints is the potential energy. Finally, and in order to validate the method, which, as far as we are aware, is the first of its kind in the literature, we compare it to the natural and straightforward finite-differences approach in a toy system and in three molecules of biological relevance: methanol, N-methyl-acetamide and a tri-glycine peptide.

Exact and efficient calculation of Lagrange multipliers in biological polymers with constrained bond lengths and bond angles: proteins and nucleic acids as example cases.

To accelerate molecular dynamics simulations, it is common to impose holonomic constraints on the hardest degrees of freedom. In this way, the time step used to integrate the equations of motion can be increased, thereby allowing longer total simulation times. The imposition of such constraints results in an aditional set of N(c) equations (the equations of constraint) and unknowns (their associated Lagrange multipliers), whose solution is closely related to any algorithm implementing the constraints in Euclidean coordinates. In this work, it is shown that, due to the essentially linear structure of typical biological polymers the algebraic equations that need to be solved involve a matrix which is not only sparse, but also banded if the constraints are indexed in a skilful way. This allows the Lagrange multipliers to be obtained through a noniterative procedure, which can be considered exact up to machine precision, and which takes O(N(c)) operations, instead of the usual O(N c3) for generic molecular systems. We develop the formalism, and describe the appropriate indexing for a number of model molecules. Finally, we provide a numerical example of the technique in a series of polyalanine peptides of different lengths. Although a use of the Lagrange multipliers without any modification in the solution of the underlying ordinary differential equations yields unstable integration algorithms, the central role of these quantities makes their efficient calculation useful for the improvement of methods that correctly enforce the exact satisfaction of the constraints at each time step. We provide several examples of this.

Exploring the free energy landscape: from dynamics to networks and back.

Knowledge of the Free Energy Landscape topology is the essential key to understanding many biochemical processes. The determination of the conformers of a protein and their basins of attraction takes a central role for studying molecular isomerization reactions. In this work, we present a novel framework to unveil the features of a Free Energy Landscape answering questions such as how many meta-stable conformers there are, what the hierarchical relationship among them is, or what the structure and kinetics of the transition paths are. Exploring the landscape by molecular dynamics simulations, the microscopic data of the trajectory are encoded into a Conformational Markov Network. The structure of this graph reveals the regions of the conformational space corresponding to the basins of attraction. In addition, handling the Conformational Markov Network, relevant kinetic magnitudes as dwell times and rate constants, or hierarchical relationships among basins, completes the global picture of the landscape. We show the power of the analysis studying a toy model of a funnel-like potential and computing efficiently the conformers of a short peptide, dialanine, paving the way to a systematic study of the Free Energy Landscape in large peptides.

Modified Ehrenfest Formalism for Efficient Large-Scale ab initio Molecular Dynamics.

We present in detail the recently derived ab initio molecular dynamics (AIMD) formalism [Alonso et al. Phys. Rev. Lett. 2008, 101, 096403], which due to its numerical properties, is ideal for simulating the dynamics of systems containing thousands of atoms. A major drawback of traditional AIMD methods is the necessity to enforce the orthogonalization of the wave functions, which can become the bottleneck for very large systems. Alternatively, one can handle the electron-ion dynamics within the Ehrenfest scheme where no explicit orthogonalization is necessary, however the time step is too small for practical applications. Here we preserve the desirable properties of Ehrenfest in a new scheme that allows for a considerable increase of the time step while keeping the system close to the Born-Oppenheimer surface. We show that the automatically enforced orthogonalization is of fundamental importance for large systems because not only it improves the scaling of the approach with the system size but it also allows for an additional very efficient parallelization level. In this work, we provide the formal details of the new method, describe its implementation, and present some applications to some test systems. Comparisons with the widely used Car-Parrinello molecular dynamics method are made, showing that the new approach is advantageous above a certain number of atoms in the system. The method is not tied to a particular wave function representation, making it suitable for inclusion in any AIMD software package.

Efficient model chemistries for peptides. I. General framework and a study of the heterolevel approximation in RHF and MP2 with Pople split-valence basis sets.

We present an exhaustive study of more than 250 ab initio potential energy surfaces (PESs) of the model dipeptide HCO-L-Ala-NH(2). The model chemistries (MCs) investigated are constructed as homo- and heterolevels involving possibly different RHF and MP2 calculations for the geometry and the energy. The basis sets used belong to a sample of 39 representants from Pople's split-valence families, ranging from the small 3-21G to the large 6-311++G(2df,2pd). The reference PES to which the rest are compared is the MP2/6-311++G(2df,2pd) homolevel, which, as far as we are aware, is the most accurate PES in the literature. All data sets have been analyzed according to a general framework, which can be extended to other complex problems and which captures the nearness concept in the space of MCs. The great number of MCs evaluated has allowed us to significantly explore this space and show that the correlation between accuracy and computational cost of the methods is imperfect, thus justifying a systematic search for the combination of features in a MC that is optimal to deal with peptides. Regarding the particular MCs studied, the most important conclusion is that the potentially very cost-saving heterolevel approximation is a very efficient one to describe the whole PES of HCO-L-Ala-NH(2). Finally, we show that, although RHF may be used to calculate the geometry if a MP2 single-point energy calculation follows, pure RHF//RHF homolevels are not recommendable for this problem.

Quantum mechanical calculation of the effects of stiff and rigid constraints in the conformational equilibrium of the alanine dipeptide.

If constraints are imposed on a macromolecule, two inequivalent classical models may be used: the stiff and the rigid one. This work studies the effects of such constraints on the conformational equilibrium distribution (CED) of the model dipeptide HCO-L-Ala-NH(2)without any simplifying assumption. We use ab initio quantum mechanics calculations including electron correlation at the MP2 level to describe the system, and we measure the conformational dependence of all the correcting terms to the naive CED based in the potential energy surface that appear when the constraints are considered. These terms are related to mass-metric tensors determinants and also occur in the Fixman's compensating potential. We show that some of the corrections are non-negligible if one is interested in the whole Ramachandran space. On the other hand, if only the energetically lower region, containing the principal secondary structure elements, is assumed to be relevant, then, all correcting terms may be neglected up to peptides of considerable length. This is the first time, as far as we know, that the analysis of the conformational dependence of these correcting terms is performed in a relevant biomolecule with a realistic potential energy function.

Explicit factorization of external coordinates in constrained statistical mechanics models.

If a macromolecule is described by curvilinear coordinates or rigid constraints are imposed, the equilibrium probability density that must be sampled in Monte Carlo simulations includes the determinants of different mass-metric tensors. In this work, the authors explicitly write the determinant of the mass-metric tensor G and of the reduced mass-metric tensor g, for any molecule, general internal coordinates and arbitrary constraints, as a product of two functions; one depending only on the external coordinates that describe the overall translation and rotation of the system, and the other only on the internal coordinates. This work extends previous results in the literature, proving with full generality that one may integrate out the external coordinates and perform Monte Carlo simulations in the internal conformational space of macromolecules.

Definition of Systematic, Approximately Separable, and Modular Internal Coordinates (SASMIC) for macromolecular simulation.

A set of rules is defined to systematically number the groups and the atoms of polypeptides in a modular manner. Supported by this numeration, a set of internal coordinates is defined. These coordinates (termed Systematic, Approximately Separable, and Modular Internal Coordinates--SASMIC) are straightforwardly written in Z-matrix form and may be directly implemented in typical Quantum Chemistry packages. A number of Perl scripts that automatically generate the Z-matrix files are provided as supplementary material. The main difference with most Z-matrix-like coordinates normally used in the literature is that normal dihedral angles ("principal dihedrals" in this work) are only used to fix the orientation of whole groups and a different type of dihedrals, termed "phase dihedrals," are used to describe the covalent structure inside the groups. This physical approach allows to approximately separate soft and hard movements of the molecule using only topological information and to directly implement constraints. As an application, we use the coordinates defined and ab initio quantum mechanical calculations to assess the commonly assumed approximation of the free energy, obtained from "integrating out" the side chain degree of freedom chi, by the Potential Energy Surface (PES) in the protected dipeptide HCO-L-Ala-NH2. We also present a subbox of the Hessian matrix in two different sets of coordinates to illustrate the approximate separation of soft and hard movements when the coordinates defined in this work are used. (PACS: 87.14.Ee, 87.15.-v, 87.15.Aa, 87.15.Cc)

A physically meaningful method for the comparison of potential energy functions.

In the study of the conformational behavior of complex systems, such as proteins, several related statistical measures are commonly used to compare two different potential energy functions. Among them, the Pearson's correlation coefficient r has no units and allows only semiquantitative statements to be made. Those that do have units of energy and whose value may be compared to a physically relevant scale, such as the root-mean-square deviation (RMSD), the mean error of the energies (ER), the standard deviation of the error (SDER) or the mean absolute error (AER), overestimate the distance between potentials. Moreover, their precise statistical meaning is far from clear. In this article, a new measure of the distance between potential energy functions is defined that overcomes the aforementioned difficulties. In addition, its precise physical meaning is discussed, the important issue of its additivity is investigated, and some possible applications are proposed. Finally, two of these applications are illustrated with practical examples: the study of the van der Waals energy, as implemented in CHARMM, in the Trp-Cage protein (PDB code 1L2Y) and the comparison of different levels of the theory in the ab initio study of the Ramachandran map of the model peptide HCO-L-Ala-NH2.

Distance-d covering problems in scale-free networks with degree correlations.

A number of problems in communication systems demand the distributed allocation of network resources in order to provide better services, sampling, and distribution methods. The solution to these issues is becoming more challenging due to the increasing size and complexity of communication networks. We report here on a heuristic method to find near-optimal solutions to the covering problem in real communication networks, demonstrating that whether a centralized or a distributed design is to be used relies upon the degree correlations between connected vertices. We also show that the general belief that by targeting the hubs one can efficiently solve most problems on networks with a power-law degree distribution is not valid for assortative networks.

Relevant distance between two different instances of the same potential energy in protein folding.

In the context of complex systems and, particularly, of protein folding, a physically meaningful distance is defined which allows to make useful statistical statements about the way in which energy differences are modified when two different instances of the same potential-energy function are used. When the two instances arise from the fact that different algorithms or different approximations are used, the distance herein defined may be used to evaluate the relative accuracy of the two methods. When the difference is due to a change in the free parameters of which the potential depends on, the distance can be used to quantify, in each region of parameter space, the robustness of the modeling to such a change and this, in turn, may be used to assess the significance of a parameters' fit. Both cases are illustrated with a practical example: the study of the Poisson-based solvation energy in the Trp-Cage protein (PDB code 1L2Y).

Improved routing strategies for Internet traffic delivery.

We analyze different strategies aimed at optimizing routing policies in the Internet. We first show that for a simple deterministic algorithm the local properties of the network deeply influence the time needed for packet delivery between two arbitrarily chosen nodes. We next rely on a real Internet map at the autonomous system level and introduce a score function that allows us to examine different routing protocols and their efficiency in traffic handling and packet delivery. Our results suggest that actual mechanisms are not the most efficient and that they can be integrated in a more general, though not too complex, scheme.