Free energy from simulations

Why Does the E1219V Mutation Expand T-Rich PAM Recognition in Cas9 from Streptococcus pyogenes?

Popular RNA-guided DNA endonuclease Cas9 from Streptococcus pyogenes (SpCas9) recognizes the canonical 5'-NGG-3' protospacer adjacent motif (PAM) and triggers double-stranded DNA cleavage activity. Mutations in SpCas9 were demonstrated to expand the PAM readability and hold promise for therapeutic and genome editing applications. However, the energetics of the PAM recognition and its relation to the atomic structure remain unknown. Using the X-ray structure (precatalytic SpCas9:sgRNA:dsDNA) as a template, we calculated the change in the PAM binding affinity in response to SpCas9 mutations using computer simulations. The E1219V mutation in SpCas9 fine-tunes the water accessibility in the PAM binding pocket and promotes new interactions in the SpCas9:noncanonical T-rich PAM, thus weakening the PAM stringency. The nucleotide-specific interaction of two arginine residues (i.e., R1333 and R1335 of SpCas9) ensured stringent 5'-NGG-3' PAM recognition. R1335A substitution (SpCas9R1335A) completely disrupts the direct interaction between SpCas9 and PAM sequences (canonical or noncanonical), accounting for the loss of editing activity. Interestingly, the double mutant (SpCas9R1335A,E1219V) boosts DNA binding affinity by favoring protein:PAM electrostatic contact in a desolvated pocket. The underlying thermodynamics explain the varied DNA cleavage activity of SpCas9 variants. A direct link between the energetics, structures, and activity is highlighted, which can aid in the rational design of improved SpCas9-based genome editing tools.

Effect of Conformational Equilibrium on Solvation Properties of 1,2-DCE in Water: A Solvation Thermodynamics and 3D-RISM Study

The contributions of the enthalpy and entropy of solvation for the study of chemical and biological systems are important in the prediction, interpretation, and manipulation of these processes. The relation between solvation Gibbs energies, enthalpies, and entropies of solvation, and their rigorous relation with the conformational equilibrium, are derived for the first time and applied with a computational method, in accordance with the Solvation Thermodynamics previous results, to 1,2-dichloroethane solvation in water. The rigid conformer calculations in solution were performed by using PC+/3D-RISM approach, with the conformational averaged results for enthalpy and solvation Gibbs energy reproducing the experimental results quite successfully. A qualitative agreement in the entropy of solvation predictions was also observed.

Adapting free energy perturbation simulations for large macrocyclic ligands: how to dissect contributions from direct binding and free ligand flexibility

Large and flexible ligands gain increasing interest in the development of bioactive agents. They challenge the applicability of computational ligand optimization strategies originally developed for small molecules. Free energy perturbation (FEP) is often used for predicting binding affinities of small molecule ligands, however, its use for more complex ligands remains limited. Herein, we report the structure-based design of peptide macrocycles targeting the protein binding site of human adaptor protein 14-3-3. We observe a surprisingly strong dependency of binding affinities on relatively small variations in substituent size. FEP was performed to rationalize observed trends. To account for insufficient convergence of FEP, restrained calculations were performed and complemented with extensive REST MD simulations of the free ligands. These calculations revealed that changes in affinity originate both from altered direct interactions and conformational changes of the free ligand. In addition, MD simulations provided the basis to rationalize unexpected trends in ligand lipophilicity. We also verified the anticipated interaction site and binding mode for one of the high affinity ligands by X-ray crystallography. The introduced fully-atomistic simulation protocol can be used to rationalize the development of structurally complex ligands which will support future ligand maturation efforts.

Contribution of Electrostatics to the Kinetics of Electron Transfer from NADH-Cytochrome b5 Reductase to Fe(III)-Cytochrome b5

Brownian dynamics (BD) simulations provide here a theoretical atomic-level treatment of the reduction of human ferric cytochrome b5 (cyt b5) by NADH-cytochrome b5 reductaste (cyt b5r) and several of its mutants. BD is used to calculate the second-order rate constant of electron transfer (ET) between the proteins for direct correlation with experiments. Interestingly, the inclusion of electrostatic forces dramatically increases the reaction rate of the native proteins despite the overall negative charge of both proteins. The role played by electrostatic charge distribution in stabilizing the ET complexes and the role of mutations of several amino acid residues in stabilizing or destabilizing the complexes are analyzed. The complex with the shortest ET reaction distance (d = 6.58 Å) from rigid body BD is further subjected to 1 ns of molecular dynamics (MD) in a periodic box of TIP3P water to produce a more stable complex allowed by flexibility and with a shorter average reaction distance d = 6.02 Å. We predict a docking model in which the following ion-ion interactions are dominant (cyt b5r/cyt b5): Lys162-Heme O1D/Lys163-Asp64/Arg91-Heme O1A/Lys125-Asp70.

Predicting RNA Duplex Dimerization Free-Energy Changes upon Mutations Using Molecular Dynamics Simulations

The dimerization free energies of RNA-RNA duplexes are fundamental values that represent the structural stability of RNA complexes. We report a comparative analysis of RNA-RNA duplex dimerization free-energy changes upon mutations, estimated from a molecular dynamics simulation and experiments. A linear regression for nine pairs of double-stranded RNA sequences, six base pairs each, yielded a mean absolute deviation of 0.55 kcal/mol and an R(2) value of 0.97, indicating quantitative agreement between simulations and experimental data. The observed accuracy indicates that the molecular dynamics simulation with the current molecular force field is capable of estimating the thermodynamic properties of RNA molecules.

Trends in template/fragment-free protein structure prediction

Predicting the structure of a protein from its amino acid sequence is a long-standing unsolved problem in computational biology. Its solution would be of both fundamental and practical importance as the gap between the number of known sequences and the number of experimentally solved structures widens rapidly. Currently, the most successful approaches are based on fragment/template reassembly. Lacking progress in template-free structure prediction calls for novel ideas and approaches. This article reviews trends in the development of physical and specific knowledge-based energy functions as well as sampling techniques for fragment-free structure prediction. Recent physical- and knowledge-based studies demonstrated that it is possible to sample and predict highly accurate protein structures without borrowing native fragments from known protein structures. These emerging approaches with fully flexible sampling have the potential to move the field forward.

CHARMM: the biomolecular simulation program

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983.

Protein Structure Evaluation using an All-Atom Energy Based Empirical Scoring Function

Arriving at the native conformation of a polypeptide chain characterized by minimum most free energy is a problem of long standing interest in protein structure prediction endeavors. Owing to the computational requirements in developing free energy estimates, scoring functions--energy based or statistical--have received considerable renewed attention in recent years for distinguishing native structures of proteins from non-native like structures. Several cleverly designed decoy sets, CASP (Critical Assessment of Techniques for Protein Structure Prediction) structures and homology based internet accessible three dimensional model builders are now available for validating the scoring functions. We describe here an all-atom energy based empirical scoring function and examine its performance on a wide series of publicly available decoys. Barring two protein sequences where native structure is ranked second and seventh, native is identified as the lowest energy structure in 67 protein sequences from among 61,659 decoys belonging to 12 different decoy sets. We further illustrate a potential application of the scoring function in bracketing native-like structures of two small mixed alpha/beta globular proteins starting from sequence and secondary structural information. The scoring function has been web enabled at www.scfbio-iitd.res.in/utility/proteomics/energy.jsp.

Molecular mechanics methods for predicting protein–ligand binding

Ligand binding affinity prediction is one of the most important applications of computational chemistry. However, accurately ranking compounds with respect to their estimated binding affinities to a biomolecular target remains highly challenging. We provide an overview of recent work using molecular mechanics energy functions to address this challenge. We briefly review methods that use molecular dynamics and Monte Carlo simulations to predict absolute and relative ligand binding free energies, as well as our own work in which we have developed a physics-based scoring method that can be applied to hundreds of thousands of compounds by invoking a number of simplifying approximations. In our previous studies, we have demonstrated that our scoring method is a promising approach for improving the discrimination between ligands that are known to bind and those that are presumed not to, in virtual screening of large compound databases. In new results presented here, we explore several improvements to our computational method including modifying the dielectric constant used for the protein and ligand interiors, and empirically scaling energy terms to compensate for deficiencies in the energy model. Future directions for further improving our physics-based scoring method are also discussed.

A molecular thermodynamic view of DNA–drug interactions: a case study of 25 minor-groove binders

Developing a molecular view of the thermodynamics of DNA recognition is essential to the design of ligands for regulating gene expression. In a first comprehensive attempt at sketching an atlas of DNA-drug energetics, we present here a detailed thermodynamic view of minor-groove recognition by small molecules via a computational study on 25 DNA-drug complexes. The studies are configured in the MMGBSA (Molecular Mechanics-Generalized Born-Solvent Accessibility) framework at the current state of the art and facilitate a structure-energy component correlation. Analyses were conducted on both energy minimized structures of DNA-drug complexes and molecular dynamics trajectories developed for the purpose of this study. While highlighting the favorable role of packing, shape complementarity, and van der Waals and hydrophobic interactions of the drugs in the minor groove in conformity with experiment, the studies reveal an interesting annihilation of favorable electrostatics by desolvation. Structural modifications attempted on the ligands point to the requisite physico-chemical factors for obtaining improved binding energies. Hydrogen bonds predicted to be important for specificity based on structural considerations do not always turn out to be significant to binding in post facto analyses of molecular dynamics trajectories, which treat thermal averaging, solvent, and counterion effects rigorously. The strength of the hydrogen bonds retained between the DNA and drug during the molecular dynamics simulations is approximately 1kcal/mol. Overall, the study reveals the compensatory nature of the diverse binding free energy components, possible threshold limits for some of these properties, and the availability of a computationally viable free energy methodology which could be of value in drug-design endeavors.

Bias and error in estimates of equilibrium free-energy differences from nonequilibrium measurements

In 1997, Jarzynski proved a remarkable equality that allows one to compute the equilibrium free-energy difference DeltaF between two states from the probability distribution of the nonequilibrium work W done on the system to switch between the states, e-DeltaF/kappaT =e-W/kappaT [Jarzynski, C. (1997) Phys. Rev. Lett. 87, 2690-2693]. The Jarzynski equality provides a powerful free-energy difference estimator from a set of N irreversible experiments and is closely related to free-energy perturbation, a common computational technique for estimating free-energy differences. Despite the many applications of the Jarzynski estimator, its behavior is only poorly understood. In this article we derive the large N limit for the Jarzynski estimator bias, variance, and mean square error that is correct for arbitrary perturbations. We then analyze the properties of the Jarzynski estimator for all N when the probability distribution of work values is Gaussian, as occurs, for example, in the near-equilibrium regime. This allows us to quantitatively compare it to two other free-energy difference estimators: the mean work estimator and the fluctuation-dissipation theorem estimator. We show that, for near-equilibrium switching, the Jarzynski estimator is always superior to the mean work estimator and is even superior to the fluctuation-dissipation estimator for small N. The Jarzynski-estimator bias is shown to be the dominant source of error in many cases. Our expression for the bias is used to develop a bias-corrected Jarzynski free-energy difference estimator in the near-equilibrium regime.

Biomolecular simulations: recent developments in force fields, simulations of enzyme catalysis, protein-ligand, protein-protein, and protein-nucleic acid noncovalent interactions

Computer modeling has been developed and widely applied in studying molecules of biological interest. The force field is the cornerstone of computer simulations, and many force fields have been developed and successfully applied in these simulations. Two interesting areas are (a) studying enzyme catalytic mechanisms using a combination of quantum mechanics and molecular mechanics, and (b) studying macromolecular dynamics and interactions using molecular dynamics (MD) and free energy (FE) calculation methods. Enzyme catalysis involves forming and breaking of covalent bonds and requires the use of quantum mechanics. Noncovalent interactions appear ubiquitously in biology, but here we confine ourselves to review only noncovalent interactions between protein and protein, protein and ligand, and protein and nucleic acids.

Structure-based drug design approaches for predicting binding affinities of HIV1 protease inhibitors

Computational assessment of the binding affinity of enzyme inhibitors prior to synthesis is an important component of computer-assisted drug design (CADD) paradigms. The free energy perturbation (FEP) methodology is the most accurate means of estimating relative binding affinities between two inhibitors. However, due to its complexity and computation-intensive nature, practical applications are restricted to analysis of structurally-related inhibitors. Accordingly, there is a need for methods that enable rapid assessment of large number of structurally-unrelated molecules in a suitably accurate manner. In this review, the FEP method is compared with regression-based methods that employ multivariate models to assess the advantages of each in the estimation of relative binding affinities of inhibitors to an enzyme. Semiquantitative predictions of relative binding free energies of human immunodeficiency virus 1 (HIV1) protease inhibitors are also presented and compared with the corresponding FEP results. The results indicate that the regression-based methods and the FEP method are useful in the semi-quantitative and quantitative assessment of relative binding affinities of enzyme inhibitors, respectively, prior to synthesis.

Protein strain in blue copper proteins studied by free energy perturbations

Free energy perturbations have been performed on two blue copper proteins, plastocyanin and nitrite reductase. By changing the copper coordination geometry, force constants, and charges, we have estimated the maximum energy with which the proteins may distort the copper coordination sphere. By comparing this energy with the quantum chemical energy cost for the same perturbation on the isolated copper complex, various hypotheses about protein strain have been tested. The calculations show that the protein can only modify the copper-methionine bond length by a modest amount of energy-<5 kJ/mol-and they lend no support to the suggestion that the quite appreciable difference in the copper coordination geometry encountered in the two proteins is a result of the proteins enforcing different Cu-methionine bond lengths. On the contrary, this bond is very flexible, and neither the geometry nor the electronic structure change appreciably when the bond length is changed. Moreover, the proteins are rather indifferent to the length of this bond. Instead, the Cu(II) coordination geometries in the two proteins represent two distinct minima on the potential surface of the copper ligand sphere, characterized by different electronic structures, a tetragonal, mainly sigma-bonded, structure in nitrite reductase and a trigonal, pi-bonded, structure in plastocyanin. In vacuum, the structures have almost the same energy, and they are stabilized in the proteins by a combination of geometric and electrostatic interactions. Plastocyanin favors the bond lengths and electrostatics of the trigonal structure, whereas in nitrite reductase, the angles are the main discriminating factor. Proteins 1999;36:157-174.

Helix-capping interaction in lambda Cro protein: a free energy simulation analysis

The stability mutant Tyr-26-->Asp was studied in the Cro protein from bacteriophage lambda using free energy molecular dynamics simulations. The mutant was calculated to be more stable than the wild type by 3.0 +/- 1.7 kcal/mol/monomer, in reasonable agreement with experiment (1.4 kcal/mol/monomer). Moreover, the aspartic acid in the mutant was found to form a capping interaction with the amino terminus of the third alpha-helix of Cro. The simulations were analyzed to understand better the source of the stability of this helix-capping interaction and to examine the results in light of previous explanations of stabilizing helix caps--namely, a model of local unsatisfied hydrogen bonds at the helix termini and the helix macrodipole model. Analysis of the simulations shows that the stabilizing effect of this charged helical cap is due both to favorable hydrogen bonds with backbone NH groups at the helix terminus and to favorable electrostatic interactions (but not hydrogen bonds) with their carbonyls (effectively the next row of local dipoles in the helix). However, electrostatic interactions are weak or negligible with backbone dipolar groups in the helix further away from the terminus. Moreover, the importance of other local electrostatic interactions with polar side chains near the helix terminus, which are neglected in most treatments of this effect, are shown to be important. Thus, the results support a model that is intermediate between the two previous explanations: both unsatisfied hydrogen bonds at the helix terminus and other, local preoriented dipolar groups stabilize the helix cap. These findings suggest that similar interactions with preoriented dipolar groups may be important for cooperativity in other charge-dipole interactions and may be employed to advantage for molecular design.

Molecular modeling methods. Basic techniques and challenging problems

An overview is presented of computer modeling and simulation methods that play an increasing role in drug design: quantum chemical methods, molecular mechanics, molecular dynamics and Brownian dynamics. The application of molecular dynamics for the prediction of thermodynamic properties like free energy differences and binding constants is discussed. The Brownian dynamics method is presented in connection with the calculation of effective electrostatic forces using the Poisson-Boltzmann equation, which allows one to sample ligand-binding geometries and to predict the kinetics of diffusion-limited enzyme reactions. New techniques that have recently been extensively developed, such as the global energy minimization and quantum-classical dynamics methods, are also introduced. The molecular modeling methods are illustrated with selected examples.

On the interpretation of biochemical data by molecular dynamics computer simulation

The application of computer simulation to molecular systems of biochemical interest is reviewed. It is shown that computer simulation is a tool complementary to experimental methods, which can be used to access atomic details inaccessible to experimental probes. Examples are given in which computer simulation augments the experimental information by providing an atomic picture of high resolution with respect to space, energy or time. The usefulness of a computer simulation largely depends on its quality. The most important factors that limit the accuracy of simulated results are discussed. The accuracy of different simulation studies can differ by orders of magnitude. The accuracy will depend on the type of biomolecular system and process studied. It will also depend on the choice of force field, the simulation set-up and the protocol that is used. A list of quality-determining factors is given, which may be useful when interpreting simulation studies appearing in the literature.

Accurate prediction of the stability and activity effects of site-directed mutagenesis on a protein core

Theoretical prediction of the structure, stability and activity of proteins, an important unsolved problem in molecular biology, would be of use for guiding site-directed mutagenesis and other protein-engineering techniques. X-ray diffraction studies have provided extensive structural information for many proteins, challenging theorists to develop reliable techniques able to use such knowledge as a base for prediction of mutants' characteristics. Here we report theoretical calculation of stabilization energies for 78 triple-site sequence variants of lambda repressor characterized experimentally by Lim and Sauer. The calculated energies correlate with the mutants' measured activities; active and inactive mutations are discriminated with 92% reliability. They correlate even more directly with the mutants' thermostabilities, correctly identifying two of the mutants to be more stable than the wild type.