Force field development phase II: Relaxation of physics-based criteria or inclusion of more rigorous physics into the representation of molecular energetics.

In the previous paper, we reviewed the origins of energy based calculations, and the early science of FF development. The initial efforts spanning the period from roughly the early 1970s to the mid to late 1990s saw the development of methodologies and philosophies of the derivation of FFs. The use of Cartesian coordinates, derivation of the H-bond potential, different functional forms including diagonal quadratic expressions, coupled valence FFs, functional form of combination rules, and out of plane angles, were all investigated in this period. The use of conformational energetics, vibrational frequencies, crystal structure and energetics, liquid properties, and ab initio data were all described to one degree or another in deriving and validating both the FF functional forms and force constants. Here we discuss the advances made since in improving the rigor and robustness of these initial FFs. The inability of the simple quadratic diagonal FF to accurately describe biomolecular energetics over a large domain of molecular structure, and intermolecular configurations, was pointed out in numerous studies. Two main approaches have been taken to overcome this problem. The first involves the introduction of error functions, either exploiting torsion terms or introducing explicit 2-D error correction grids. The results and remaining challenges of these functional forms is examined. The second approach has been to improve the representation of the physics of intra and intermolecular interactions. The latter involves including descriptions of polarizability, charge flux aka geometry dependent charges, more accurate representations of spatial electron density such as multipole moments, anisotropic nonbond potentials, nonbond and polarization flux, among others. These effects, though not as extensively studied, likely hold the key to achieving the rigorous reproduction of structural and energetic properties long sought in biomolecular simulations, and are surveyed here. In addition, the quality of training and validation observables are evaluated. The necessity of including an ample set of energetic and crystal observables is emphasized, and the inadequacy of free energy as a criterion for FF reliability discussed. Finally, in light of the results of applications of the two approaches to FF development, we propose a "recipe" of terms describing the physics of inter and intramolecular interactions whose inclusion in FFs would significantly improve our understanding of the energetics and dynamics of biomolecular systems resulting from molecular dynamics and other energy based simulations.

Biomolecular force fields: where have we been, where are we now, where do we need to go and how do we get there?

In this perspective, we review the theory and methodology of the derivation of force fields (FFs), and their validity, for molecular simulations, from their inception in the second half of the twentieth century to the improved representations at the end of the century. We examine the representations of the physics embodied in various force fields, their accuracy and deficiencies. The early days in the 1950s and 60s saw FFs first introduced to analyze vibrational spectra. The advent of computers was soon followed by the first molecular mechanics machine calculations. From the very first papers it was recognized that the accuracy with which the FFs represented the physics was critical if meaningful calculated structural and thermodynamic properties were to be achieved. We discuss the rigorous methodology formulated by Lifson, and later Allinger to derive molecular FFs, not only obtain optimal parameters but also uncover deficiencies in the representation of the physics and improve the functional form to account for this physics. In this context, the known coupling between valence coordinates and the importance of coupling terms to describe the physics of this coupling is evaluated. Early simplified, truncated FFs introduced to allow simulations of macromolecular systems are reviewed and their subsequent improvement assessed. We examine in some depth: the basis of the reformulation of the H-bond to its current description; the early introduction of QM in FF development methodology to calculate partial charges and rotational barriers; the powerful and abundant information provided by crystal structure and energetic observables to derive and test all aspects of a FF including both nonbond and intramolecular functional forms; the combined use of QM, along with crystallography and lattice energy calculations to derive rotational barriers about ɸ and ψ; the development and results of methodologies to derive "QM FFs" by sampling the QM energy surface, either by calculating energies at hundreds of configurations, or by describing the energy surface by energies, first and second derivatives sampled over the surface; and the use of the latter to probe the validity of the representations of the physics, reveal flaws and assess improved functional forms. Research demonstrating significant effects of the flaws in the use of the improper torsion angle to represent out of plane deformations, and the standard Lorentz-Berthelot combining rules for nonbonded interactions, and the more accurate descriptions presented are also reviewed. Finally, we discuss the thorough studies involved in deriving the 2nd generation all-atom versions of the CHARMm, AMBER and OPLS FFs, and how the extensive set of observables used in these studies allowed, in the spirit of Lifson, a characterization of both the abilities, but more importantly the deficiencies in the diagonal 12-6-1 FFs used. The significant contribution made by the extensive set of observables compiled in these papers as a basis to test improved forms is noted. In the following paper, we discuss the progress in improving the FFs and representations of the physics that have been investigated in the years following the research described above.

Mechanism of androgen receptor antagonism by bicalutamide in the treatment of prostate cancer.

The androgen receptor (AR) plays a key role in regulating gene expression in a variety of tissues, including the prostate. In that role, it is one of the primary targets in the development of new chemotherapeutics for treatment of prostate cancer and the target of the most widely prescribed current drug, bicalutamide (Bcu), for this disease. In view of its importance, and the absence of a crystal structure for any antagonist--AR complex, we have conducted a series of molecular dynamics-based simulations of the AR--Bcu complex and quantum mechanical (QM) calculations of Bcu, to elucidate the structural basis for antagonism of this key target. The structures that emerge show that bicalutamide antagonizes AR by accessing an additional binding pocket (B-site) adjacent to the hormone binding site (HBS), induced by displacing helix 12. This distorts the coactivator binding site and results in the inactivation of transcription. An alternative equienergetic conformational state of bicalutamide was found to bind in an expanded hormone pocket without materially perturbing either helix 12 or the coactivator binding site. Thus, both the structural basis of antagonism and the mechanism underlying agonist properties displayed by bicalutamide in different environments may be rationalized in terms of these structures. In addition, the antagonist structure and especially the induced second site (B-site) provide a structural framework for the design of novel antiandrogens.

Monte Carlo simulation of water behavior around the dipeptide N-acetylalanyl-N-methylamide.

Applications of Monte Carlo and molecular dynamics computer simulation techniques indicate that they are potentially powerful tools for understanding biological systems at the molecular level. The Monte Carlo technique can be used to study the solvent structure around a small peptide and the effect of the aqueous environment on the conformational equilibria of the peptide.

Dynamics and conformational energetics of a peptide hormone: vasopressin.

A theoretical methodology for use in conjunction with experiment was applied to the neurohypophyseal hormone lysine vasopressin for elucidation of its accessible molecular conformations and associated flexibility, conformational transitions, and dynamics. Molecular dynamics and energy minimization techniques make possible a description of the conformational properties of a peptide in terms of the precise positions of atoms, their fluctuations in time, and the interatomic forces acting on them. Analysis of the dynamic trajectory of lysine vasopressin shows the ability of a flexible peptide hormone to undergo spontaneous conformational transitions. The excursions of an individual phenylalanine residue exemplify the dynamic flexibility and multiple conformational states available to small peptide hormones and their component residues, even within constraints imposed by a cyclic hexapeptide ring.

Keeping it in the family: folding studies of related proteins.

Investigators have recently turned to studies of protein families to shed light on the mechanism of protein folding. In small proteins for which detailed analysis has been performed, recent studies show that transition-state structure is generally conserved. The number and structures of populated folding intermediates have been found to vary in homologous families of larger (greater than 100-residue) proteins, reflecting a balance of local and global interactions.

Quantum Derivative Fitting and Biomolecular Force Fields: Functional Form, Coupling Terms, Charge Flux, Nonbond Anharmonicity, and Individual Dihedral Potentials.

Computer simulations are increasingly prevalent, complementing experimental studies in all fields of biophysics, chemistry, and materials. Their utility, however, is critically dependent on the validity of the underlying force fields employed. In this Perspective we review the ability of quantum mechanics, and in particular analytical ab initio derivatives, to inform on the nature of intra- and intermolecular interactions. The power inherent in the exploitation of forces and second derivatives (Hessians) to derive force fields for a variety of compound types, including inorganic, organic, and biomolecules, is explored. We discuss the use of these quantities along with QM energies and geometries to determine force constants, including nonbond and electrostatic parameters, and to assess the functional form of the energy surface. The latter includes the optimal form of out-of-plane interactions and the necessity for anharmonicity, and terms to account for coupling between internals, to adequately represent the energy of intramolecular deformations. In addition, individual second derivatives of the energy with respect to selected interaction coordinates, such as interatomic distances or individual dihedral angles, have been shown to select out for the corresponding interactions, annihilating other interactions in the potential expression. Exploitation of these quantities allows one to probe the individual interaction and explore phenomena such as, for example, anisotropy of atom-atom nonbonded interactions, charge flux, or the functional form of isolated dihedral angles, e.g., a single dihedral X-C-C-Y about a tetrahedral C-C bond.

Conformational analysis of a highly potent dicyclic gonadotropin-releasing hormone antagonist by nuclear magnetic resonance and molecular dynamics.

Structural analysis of constrained (monocyclic) analogues of gonadotropin-releasing hormone (GnRH) has led to the development of a model for the receptor-bound conformation of GnRH and to the design of highly potent, dicyclic GnRH antagonists. This is one of the first cases where a dicyclic backbone has been introduced into analogues of a linear peptide hormone with retention of high biological activity. Here we present a conformational analysis of dicyclo(4-10,5-8)[Ac-D-2Nal1-D-pClPhe2-D-Trp3-Asp4+ ++-Glu5-D-Arg6-Leu7-Lys8- Pro9-Dpr10]-NH2 (I), using two-dimensional nuclear magnetic resonance (NMR) spectroscopy and molecular dynamics simulation. Compound I inhibits ovulation in the rat at a dose of 5-10 micrograms (Rivier et al. In Peptides: Chemistry, Structure ad Biology; Rivier, J. E., Marshall, G. R., Eds.; ESCOM: Leiden, The Netherlands, 1990; pp 33-37). The backbone conformation of the 4-10 cycle in this dicyclic compound is very similar to that found previously for a parent monocyclic (4-10) GnRH antagonist (Rizo et al. J. Am. Chem. Soc. 1992, 114, 2852-2859; ibid. 2860-2871), which gives strong support to the hypothesis that GnRH adopts a similar conformation upon binding to its receptor. In this conformation, residues 5-8 form a "beta-hairpin-like" structure that includes two transannular hydrogen bonds and a Type II' beta turn around residues D-Arg6-Leu7. The "tail" of the molecule formed by residues 1-3 is somewhat structured but does not populate a single major conformation. However, the orientation of the tail on the same side of the 4-10 cycle as the 5-8 bridge favors interactions between this bridge and the tail residues. These observations correlate with results obtained previously for the parent monocyclic (4-10) antagonist, and have led to the design of a series of new dicyclic GnRH antagonists with bridges between the tail residues and residues 5 or 8.

Computer simulation of the solvent structure around biological macromolecules.

The structure and energetics of the solvent, both ordered and disordered, in a small peptide crystal and in the triclinic lysozyme crystal have been simulated using the Monte Carlo method. The results are in good agreement with the experimental data and provide insight into the role of solvent structure around biological macromolecules in solution.

On the formation of protein tertiary structure on a computer.

In this paper we carry out computer simulation studies of some of the factors responsible for protein tertiary structure. We show that it is possible to obtain (fold) a compact globular conformation from a sequence of amino acids consisting of only glycines and alanines. Our results indicate that glycines play a central role in stabilizing globular structures by facilitating the formation of turns and by destabilizing helical structures. Using this simple two-amino-acid representation, which serves as a control experiment, we are able to obtain a conformation that resembles the native structure of pancreatic trypsin inhibitor, as closely as any obtained previously in folding studies. However, careful examination reveals that the true chain topology has not been reproduced here or in previous studies. We suggest that the discrepancies between calculated and observed structures are more significant than the similarities. The implications of these results for the validity of models for protein folding, the use of pancreatic trypsin inhibitor in folding studies, and the possible role of glycine in the evolution of protein structure are discussed.

Theoretical studies of the structure and molecular dynamics of a peptide crystal.

Energy minimizations and molecular dynamics simulations have been performed on the cyclic peptide cyclo-(Ala-Pro-D-Phe)2 in both the isolated and crystal states. The results of these calculations have been analyzed, both to investigate our ability to reproduce experimental data (structure and vibrational and NMR spectra) and to investigate the effects of environment on the energy, structure, and dynamics of peptides. Comparison of the minimized and time-averaged crystal systems with the experimental peptide structure shows that the calculations have closely reproduced the experimental structure. Molecular dynamics of the isolated molecule has led to a new conformation, which is approximately equal to 8.5 kcal/mol more stable than the conformation that exists in the crystal, the latter conformation being stabilized by intermolecular (packing) forces. This illustrates the considerable effect that environment can have on the conformation of peptides. The crystal environment has also been shown to significantly reduce the dynamic conformational fluctuations seen for the isolated molecule. The behavior of the peptide during the isolated simulation also supports the experimental NMR observation of a symmetric structure that differs from the asymmetric, instantaneous structures which characterize the molecule during the dynamics. Calculations of vibrational frequencies of the peptide in the crystal and isolated states show the expected shifts in bond-stretching frequencies due to intermolecular interactions. Finally, we have calculated NMR coupling constants from the dynamics simulation of the isolated peptide and have compared these with the experimental values. This has led to a possible reinterpretation of the experimental data.

Catalysis of a rotational transition in a peptide by crystal forces.

Detailed examination of the dynamics trajectories of the isolated cyclic peptide cyclo-(Ala-Pro-D-Phe)2 and of the molecule in its crystalline environment led to the unexpected observation that the methyl groups of the alanine residues rotated more frequently during a simulation in the crystal environment than in a simulation of the isolated peptide. In effect, the crystal environment is "catalyzing" the rotational isomerization of the methyl groups. In order to understand how the crystal forces increase the rate of this rotation, and to explore any possible analogy to the inducing of strained conformations of ligands by enzymes, the barriers to rotation in the two environments were studied by using the torsion angle forcing method. The crystal forces induce a different, higher energy, conformation of the peptide than is found for the isolated molecule, and the different rates of rotation have been explained in terms of the resulting specific intramolecular interactions that, it turns out, give rise to the lower rotational barrier. Molecular dynamics simulations of the peptide were also run at higher temperatures in order to calculate the barriers to rotation through the use of Arrhenius plots. The barriers obtained in this way agree well with the barriers obtained through an adiabatic reaction path derived by rotating the methyl through the barrier while minimizing all remaining degrees of freedom. The rates of rotation calculated from these adiabatic barriers also agree well with the rates observed during the 300 K simulations.

Predicted three-dimensional structure of the protease inhibitor domain of the Alzheimer's disease beta-amyloid precursor.

Alzheimer's disease is characterized by the deposition of amyloid beta-protein as plaques and tangles in the brains of its victims. The amyloid precursor can be expressed with or without the inclusion of a protease inhibitor domain, the potential role of which in amyloidogenesis has prompted the generation of a model of its three-dimensional structure based on the known structure of a related inhibitor. The model structure predicts that the mutated residues are almost entirely on the surface of the inhibitor domain, while conserved residues constitute the hydrophobic core. In addition, several pairs of structurally complementary, or concerted, mutations are seen. These structural features provide strong evidence for the validity of the modeled structure, and it is suggested that the presence of complementary mutations may be used as a criterion for evaluating protein structures built by homology, in addition to the (spatial) location of the mutations. The terminal residues delimiting the domain are among those furthest from the protease binding site and are in close proximity to one another, thus suggesting the ability of the domain to function as a structural cassette within the context of a larger protein. The electrostatic potentials of the inhibitor and of the related bovine pancreatic trypsin inhibitor reveal how two inhibitors with very different net charges can bind with approximately the same binding constant to trypsin and suggest a mutation of trypsin that might selectively enhance the binding of the amyloid inhibitor domain. The model provides a structural basis for understanding the functional roles of residues in the domain and for designing simpler molecules to test as pharmacologic agents for intervention in Alzheimer's disease.

Influence of environment on the antifolate drug trimethoprim: energy minimization studies.

Environmental effects on trimethoprim (TMP), an inhibitor of bacterial dihydrofolate reductase (DHFR), were investigated with energy minimizations in vacuo, in the crystal, and in aqueous solution. The conformations, harmonic dynamics, and energetics of the antibacterial drug calculated in these environments were compared with each other and with those of two enzyme-bound drugs. Valence and torsion angles and their energies and overall intra- and intermolecular energies compensated one another in the minimized TMP structures. The conformations of the isolated and aqueous molecules were similar to that of TMP bound to chicken liver DHFR, while the structures from the TMP crystal and from the Escherichia coli DHFR complex were unique. Since neither the small-molecule crystal nor a local minimum of the isolated molecule gave the conformation of TMP bound to the bacterial enzyme, a combination of several experimental and theoretical techniques may be necessary to probe accessible conformations of a molecule.

Derivation of force fields for molecular mechanics and dynamics from ab initio energy surfaces.

We present a technique for addressing the problem of deriving potential energy functions for the simulation of organic, polymeric, and biopolymeric systems, as well as for modeling vibrational spectroscopic properties. This method is designed to address three major objectives: deriving and comparing optimal functional forms for describing the energies of molecular deformations and interactions, developing a technique to rapidly and objectively determine reasonable force constants for intermolecular and intramolecular interactions, and determining the transferability of these potential forms and constants. The first two of these objectives are addressed in this paper, while the latter problem will be treated elsewhere. The technique uses ab initio molecular energy surfaces, which are described by the energy and its first and second derivatives with respect to coordinates. As an example, application to a small model compound (i.e., the formate anion) is given. A variety of analytical forms for the potential are tested against the data, to find which forms are best. The importance of anharmonicity and cross terms in accounting for structure and energy, as well as for dynamics, is demonstrated and a more accurate representation of the out-of-plane deformation for a trigonal center is derived from the energy surfaces.

Theoretical studies on the dihydrofolate reductase mechanism: electronic polarization of bound substrates.

We have applied local density functional theory, an ab initio quantum mechanical method, to study the shift in the spatial electron density of the substrate dihydrofolate that accompanies binding to the enzyme dihydrofolate reductase. The results shed light on fundamental electronic effects due to the enzyme that may contribute to catalysis. In particular, the enzyme induces a long-range polarization of the substrate that perturbs its electron density distribution in a specific and selective way in the vicinity of the bond that is reduced by the enzyme. Examination of the electron density changes that occur in folate reveals that a similar effect is seen but this time specifically at the bond that is reduced in this substrate. This suggests that the polarization effect may be implicated in the reaction mechanism and may play a role in determining the sequence whereby the 7,8-bond in folate is reduced first, followed by reduction of the 5,6-bond in the resulting dihydro compound.

Computer simulation of the free energy of peptides with the local states method: analogues of gonadotropin releasing hormone in the random coil and stable states.

The Helmholtz free energy F (rather than the energy) is the correct criterion for stability; therefore, calculation of F is important for peptides and proteins that can populate a large number of metastable states. The local states (LS) method proposed by H. Meirovitch [(1977) Chemical Physics Letters, Vol. 45, p. 389] enables one to obtain upper and lower bounds of the conformational free energy, FB (b,l) and FA (b,l), respectively, from molecular dynamics (MD) or Monte Carlo samples. The correlation parameter b is the number of consecutive dihedral or valence angles along the chain that are taken into account explicitly. The continuum angles are approximated by a discretization parameter l; the larger are b and l, the better the approximations; while FA can be estimated efficiently, it is more difficult to estimate FB. The method is further developed here by applying it to MD trajectories of a relatively large molecule (188 atoms), the potent "Asp4-Dpr10" antagonist [cyclo(4/10)-(Ac-delta 3Pro1-D-pFPhe2-D-Trp3-Asp4-Tyr5-D-Nal6-Leu7-Arg8 -Pro9- Dpr10-NH2)] of gonadotropin releasing hormone (GnRH). The molecule was simulated in vacuo at T = 300 K in two conformational states, previously investigated [J. Rizo et al. Journal of the American Chemical Society, (1992) Vol. 114, p. 2860], which differ by the orientation of the N-terminal tail, above (tail up, TU) and below (tail down, TD) the cyclic heptapeptide ring. As in previous applications of the LS method, we have found the following: (1) While FA is a crude approximation for the correct F, results for the difference, delta FA = FA (TD)-FA (TU) converge rapidly to 5.6 (1) kcal/mole as the approximation is improved (i.e., as b and l are increased), which suggests that this is the correct value for delta F; therefore TD is more stable than TU. (The corresponding difference in entropy, T delta SA = 1.3(2) kcal/mole, is equal to the value obtained by the harmonic approximation.) (2) The lowest approximation, which has the minimal number of local states, i.e., based on b = 0 (no correlations) and l = 1 (the angle values are distributed homogeneously), also leads to the correct value of delta F, within the error bars. This is important since the lowest approximation can be applied even to large proteins. (3) The method enables one to define the entropy of a part of the molecule and thus to measure the flexibility of this part.(ABSTRACT TRUNCATED AT 400 WORDS)