Search

The anomeric ratio for D-xylose in aqueous solution was investigated by applying standard molecular dynamics (MD) and free energy simulations for the two anomers of D-xylopyranose in aqueous solution. All simulations employed the CHARMM molecular mechanics program. The simulations predicted a small free energy difference between the anomers, which results from the cancellation of an intramolecular term that favors the alpha anomer, indicating that the experimental value favors the beta form. The validity of the results were evaluated by performing a direct evaluation of the energy difference and its components from standard MD simulations of each anomer.

The effects of subnanosecond internal motions on proton-proton cross relaxation are investigated using the cyclic decapeptide antamanide (-Val1-Pro2-Pro3-Ala4-Phe5-Phe6-Pro7-Pro8-Phe9-Phe10-) in chloroform. An 800-ps molecular dynamics simulation of the antamanide shows that most of these effects may be illustrated by using a product approximation which distinguishes between radial and angular order parameters and an internal correlation time.

A study was conducted to develop novel protein parameters for the all-atom empirical energy function in the CHARMM program. The parameter evaluation focused on a self-consistent approach designed to realize a balance between the internal and interaction terms of the force field and interactions associated with solvents and solutes. The parameters were then obtained by fitting an extended set of experimental and ab initio results. In addition, internal parametrizations were selected to duplicate geometries from crystal structures, Raman spectroscopic data and ab initio computations.

The SPASIBA force field is parametrized for a large variety of chemical groups. It is shown that the refinement of the force field parameters using the spectroscopic data leads to a better prediction of various physical properties leading particularly to confidence when used in molecular dynamics simulations of molecules with biological interest including various chemical groups such as peptido-glycans or phospholipids.

Since its inception nearly a half century ago, CHARMM has been playing a central role in computational biochemistry and biophysics. Commensurate with the developments in experimental research and advances in computer hardware, the range of methods and applicability of CHARMM have also grown. This review summarizes major developments that occurred after 2009 when the last review of CHARMM was published. They include the following: new faster simulation engines, accessible user interfaces for convenient workflows, and a vast array of simulation and analysis methods that encompass quantum mechanical, atomistic, and coarse-grained levels, as well as extensive coverage of force fields. In addition to providing the current snapshot of the CHARMM development, this review may serve as a starting point for exploring relevant theories and computational methods for tackling contemporary and emerging problems in biomolecular systems. CHARMM is freely available for academic and nonprofit research at https://academiccharmm.org/program.

Infection by sarbecoviruses begins with the attachment of the homotrimeric viral "spike" protein to the angiotensin-converting enzyme 2 receptor on the surface of mammalian cells. This requires one or more receptor-binding domains (RBDs) to be in the open (up) position. Here, we present the results of long molecular dynamics simulations with umbrella sampling (US) to compute a one-dimensional free energy profile of RBD opening/closing and the associated transition times. After ≃3.58 μs of simulation time per US window (∼229 μs in total), which was required to approach trajectory decorrelation, the computed free energy profile was found to be without large barriers. This suggests that the RBD diffuses between the open and closed positions without significant energetic hindrance. This interpretation appears consistent with experiments but is at odds with some previous simulations. Modeling the RBD motion as diffusive dynamics along the computed free energy profile, we find that the overall time required for the transition is only about 2 μs, which is 5 orders of magnitude shorter than experimentally measured transition times. We speculate that the most likely reason for the transition time mismatch is our use of very short glycans, which was required to make the simulations performed here feasible. Despite the long simulation times, the final free energy profile is not fully converged with statistical errors of ≃1.16 kcal/mol, which were found to be consistent with the slow time decay in the autocorrelation of the conformational motions of the protein. The simulation lengths that would be required to obtain fully converged results remain unknown, but the present calculations would benefit from at least an order-of-magnitude extension.

The water dynamics, as characterized by the local hydrophobicity (LH), is investigated for tetrameric hemoglobin (Hb) and dimeric melittin. For the T 0 to R 0 transition in Hb, it is found that LH provides additional molecular-level insight into the Perutz mechanism, i.e., the breaking and formation of salt bridges at the α 1 /β 2 and α 2 /β 1 interface is accompanied by changes in LH. For Hb in cubic water boxes with 90 and 120 Å edge length it is observed that following a decrease in LH as a consequence of reduced water density or change of water orientation at the protein/water interface the α/β interfaces are destabilized; this is a hallmark of the Perutz stereochemical model for the T to R transition in Hb. The present work thus provides a dynamical view of the classical structural model relevant to the molecular foundations of Hb function. For dimeric melittin, earlier results by Cheng and Rossky [ Nature 1998, 392, 696-699] are confirmed and interpreted on the basis of LH from simulations in which the protein structure is frozen. For the flexible melittin dimer, the changes in the local hydration can be as much as 30% greater than for the rigid dimer, reflecting the fact that protein and water dynamics are coupled.

We describe the development and implementation of a quasi-equilibrium hydration shell model of biomolecular solvation with adaptive boundaries. Applying the model to microsecond-long molecular dynamics simulations of several protein systems of varying complexity, we find that the model simulation results are of comparable quality to those obtained from simulations of fully solvated systems, but at a reduced computational cost. We discuss the dominant sources of error in the model and outline directions for future improvements.

A method is developed to obtain simultaneously free energy profiles and diffusion constants from restrained molecular simulations in diffusive systems. The method is based on low-order expansions of the free energy and diffusivity as functions of the reaction coordinate. These expansions lead to simple analytical relationships between simulation statistics and model parameters. The method is tested on 1D and 2D model systems; its accuracy is found to be comparable to or better than that of the existing alternatives, which are briefly discussed. An important aspect of the method is that the free energy is constructed by integrating its derivatives, which can be computed without need for overlapping sampling windows. The implementation of the method in any molecular simulation program that supports external umbrella potentials (e.g., CHARMM) requires modification of only a few lines of code. As a demonstration of its applicability to realistic biomolecular systems, the method is applied to model the α-helix ↔ β-sheet transition in a 16-residue peptide in implicit solvent, with the reaction coordinate provided by the string method. Possible modifications of the method are briefly discussed; they include generalization to multidimensional reaction coordinates [in the spirit of the model of Ermak and McCammon (Ermak, D. L.; McCammon, J. A. J. Chem. Phys. 1978, 69, 1352-1360)], a higher-order expansion of the free energy surface, applicability in nonequilibrium systems, and a simple test for Markovianity. In view of the small overhead of the method relative to standard umbrella sampling, we suggest its routine application in the cases where umbrella potential simulations are appropriate.

Simulations of first-passage folding of the antiparallel β-sheet miniprotein beta3s, which has been intensively studied under equilibrium conditions by A. Caflisch and co-workers, show that the kinetics and dynamics are significantly different from those for equilibrium folding. Because the folding of a protein in a living system generally corresponds to the former (i.e., the folded protein is stable and unfolding is a rare event), the difference is of interest. In contrast to equilibrium folding, the Ch-curl conformations become very rare because they contain unfavorable parallel β-strand arrangements, which are difficult to form dynamically due to the distant N- and C-terminal strands. At the same time, the formation of helical conformations becomes much easier (particularly in the early stage of folding) due to short-range contacts. The hydrodynamic descriptions of the folding reaction have also revealed that while the equilibrium flow field presented a collection of local vortices with closed "streamlines", the first-passage folding is characterized by a pronounced overall flow from the unfolded states to the native state. The flows through the locally stable structures Cs-or and Ns-or, which are conformationally close to the native state, are negligible due to detailed balance established between these structures and the native state. Although there are significant differences in the general picture of the folding process from the equilibrium and first-passage folding simulations, some aspects of the two are in agreement. The rate of transitions between the clusters of characteristic protein conformations in both cases decreases approximately exponentially with the distance between the clusters in the hydrogen bond distance space of collective variables, and the folding time distribution in the first-passage segments of the equilibrium trajectory is in good agreement with that for the first-passage folding simulations.

The Bohr effect in hemoglobin, which refers to the dependence of the oxygen affinity on the pH, plays an important role in its cooperativity and physiological function. The dominant contribution to the Bohr effect arises from the difference in the pKa values of His residues of the unliganded (deoxy) and liganded (carbonmonoxy) structures. Using recent high resolution structures, the residue pKa values corresponding to the two structures are calculated. The method is based on determining the electrostatic interactions between residues in the protein, relative to those of the residue in solution, by use of the linearized finite difference Poisson-Boltzmann equation and Monte Carlo sampling of protonation states. Given that good agreement is obtained with the available experimental values for the contribution of His residues in HbA to the Bohr effect, the calculated results are used to determine the atomic origin of the pKa shift between deoxy and carbonmonoxy HbA. The contributions to the pKa shift calculated by means of the linear response approximation show that the salt bridge involving His146 plays an important role in the alkaline Bohr effect, as suggested by Perutz but that other interactions are significant as well. A corresponding analysis is made for the contribution of His143 to the acid Bohr effect for which there is no proposed explanation. The method used is summarized and the program by which it is implemented is described in the Appendix .

A new analysis of the 20 μs equilibrium folding/unfolding molecular dynamics simulations of the three-stranded antiparallel β-sheet miniprotein (beta3s) in implicit solvent is presented. The conformation space is reduced in dimensionality by introduction of linear combinations of hydrogen bond distances as the collective variables making use of a specially adapted principal component analysis (PCA); i.e., to make structured conformations more pronounced, only the formed bonds are included in determining the principal components. It is shown that a three-dimensional (3D) subspace gives a meaningful representation of the folding behavior. The first component, to which eight native hydrogen bonds make the major contribution (four in each beta hairpin), is found to play the role of the reaction coordinate for the overall folding process, while the second and third components distinguish the structured conformations. The representative points of the trajectory in the 3D space are grouped into conformational clusters that correspond to locally stable conformations of beta3s identified in earlier work. A simplified kinetic network based on the three components is constructed, and it is complemented by a hydrodynamic analysis. The latter, making use of "passive tracers" in 3D space, indicates that the folding flow is much more complex than suggested by the kinetic network. A 2D representation of streamlines shows there are vortices which correspond to repeated local rearrangement, not only around minima of the free energy surface but also in flat regions between minima. The vortices revealed by the hydrodynamic analysis are apparently not evident in folding pathways generated by transition-path sampling. Making use of the fact that the values of the collective hydrogen bond variables are linearly related to the Cartesian coordinate space, the RMSD between clusters is determined. Interestingly, the transition rates show an approximate exponential correlation with distance in the hydrogen bond subspace. Comparison with the many published studies shows good agreement with the present analysis for the parts that can be compared, supporting the robust character of our understanding of this "hydrogen atom" of protein folding.

A simple and robust formulation of the path-independent confinement method for the calculation of free energies is presented. The simplified confinement method (SCM) does not require matrix diagonalization or switching off the molecular force field, and has a simple convergence criterion. The method can be readily implemented in molecular dynamics programs with minimal or no code modifications. Because the confinement method is a special case of thermodynamic integration, it is trivially parallel over the integration variable. The accuracy of the method is demonstrated using a model diatomic molecule, for which exact results can be computed analytically. The method is then applied to the alanine dipeptide in vacuum, and to the α-helix ↔ β-sheet transition in a 16-residue peptide modeled in implicit solvent. The SCM requires less effort for the calculation of free energy differences than previous formulations because it does not require computing normal modes. The SCM has a diminished advantage for determining absolute free energy values, because it requires decreasing the MD integration step to obtain accurate results. An approximate confinement procedure is introduced, which can be used to estimate directly the configurational entropy difference between two macrostates, without the need for additional computation of the difference in the free energy or enthalpy. The approximation has convergence properties similar to those of the standard confinement method for the calculation of free energies. The use of the approximation requires about 5 times less wall-clock simulation time than that needed to compute enthalpy differences to similar precision from an MD trajectory. For the biomolecular systems considered in this study, the errors in the entropy approximation are under 10%. Practical applications of the methods to proteins are currently limited to implicit solvent simulations.

The popular targeted molecular dynamics (TMD) method for generating transition paths in complex biomolecular systems is revisited. In a typical TMD transition path, the large-scale changes occur early and the small-scale changes tend to occur later. As a result, the order of events in the computed paths depends on the direction in which the simulations are performed. To identify the origin of this bias, and to propose a method in which the bias is absent, variants of TMD in the restraint formulation are introduced and applied to the complex open ↔ closed transition in the protein calmodulin. Due to the global best-fit rotation that is typically part of the TMD method, the simulated system is guided implicitly along the lowest-frequency normal modes, until the large spatial scales associated with these modes are near the target conformation. The remaining portion of the transition is described progressively by higher-frequency modes, which correspond to smaller-scale rearrangements. A straightforward modification of TMD that avoids the global best-fit rotation is the locally restrained TMD (LRTMD) method, in which the biasing potential is constructed from a number of TMD potentials, each acting on a small connected portion of the protein sequence. With a uniform distribution of these elements, transition paths that lack the length-scale bias are obtained. Trajectories generated by steered MD in dihedral angle space (DSMD), a method that avoids best-fit rotations altogether, also lack the length-scale bias. To examine the importance of the paths generated by TMD, LRTMD, and DSMD in the actual transition, we use the finite-temperature string method to compute the free energy profile associated with a transition tube around a path generated by each algorithm. The free energy barriers associated with the paths are comparable, suggesting that transitions can occur along each route with similar probabilities. This result indicates that a broad ensemble of paths needs to be calculated to obtain a full description of conformational changes in biomolecules. The breadth of the contributing ensemble suggests that energetic barriers for conformational transitions in proteins are offset by entropic contributions that arise from a large number of possible paths.

The response of a protein to variation of a specific coordinate can provide insights into the role of the overall architecture in the structural change. Given that the calculated potential of mean force governing the fluctuation of an electron transfer donor-acceptor distance in the NAD(P)H:Flavin oxidoreductase (Fre)/FAD complex was shown to agree with experiment, an analysis of the structural response of the rest of the protein to that distance change was made. Significant displacements are found throughout much of the protein, and the coupling pathway resulting in the structural changes was determined. A covariance analysis based on the quasiharmonic modes of the unperturbed protein was used to provide information concerning how the residue motions are correlated. It is found that, of the three regions identified as moving together in an NMR study, two undergo significant structural changes when the electron donor-acceptor distance is varied, and the third does not.

To determine the factors involved in the specific recognition function of a bacterial 8-oxoguanine (oxoG) DNA glycosylase MutM, a series of potentials of mean force and thermodynamic integration simulations were performed with the wild type and a single-point E3Q mutant of MutM bound to oxoG and G-containing DNA, respectively. Interestingly, the mutation of the catalytically important Glu3 (E3) residue to Gln (Q) significantly changes the free-energy surface so that oxoG can bind stably in the active site of the enzyme. Free-energy simulations with the protonated and deprotonated E3 residue further showed that the protonation of the catalytically important E3 residue plays a key role in distinguishing oxoG versus G in the active site by lowering the free energy of oxoG preferentially in the active site. The results suggest that MutM utilizes the thermodynamic recognition mechanism for stable binding of the lesion base in the active site of the enzyme in addition to kinetic discrimination at the early stage of the base extrusion for facilitated extrusion of oxoG.

Discrete molecular dynamics has been used to study the folding of a SH3 domain with a Calpha-based Go-model at a temperature within the native state stability region. A standard analysis of the folding process, based on consideration of the mean-force (free energy) surfaces, contact maps and folding time distributions, is complemented by a "hydrodynamic" description of folding flows (Chekmarev et al., PRL, 2008, 018107) using two and three collective variables. Two types of folding trajectories (fast and slow) follow essentially different routes in the final stage of folding. The hydrodynamic description makes possible the calculation of folding flows corresponding to these routes. The results show that the probability flows do not correspond to the free energy surface and that vortex formation is involved in the slow trajectories. Comparison of the simulation results with the experimental data suggests that the two-state kinetics observed for Fyn and Src SH3 domain folding are associated with the slow trajectories, in which a partly formed N- and C-terminal beta sheet hinders the RT-loop from attaching to the protein core; the fast trajectories are not observed because they are in the dead time (1 ms) of the experiments.

We introduce the Gaussian-mixture umbrella sampling method (GAMUS) , a biased molecular dynamics technique based on adaptive umbrella sampling that efficiently escapes free energy minima in multidimensional problems. The prior simulation data are reweighted with a maximum likelihood formulation, and the new approximate probability density is fit to a Gaussian-mixture model, augmented by information about the unsampled areas. The method can be used to identify free energy minima in multidimensional reaction coordinates. To illustrate GAMUS , we apply it to the alanine dipeptide (2D reaction coordinate) and tripeptide (4D reaction coordinate).

The conformational space of a 20-residue three-stranded antiparallel beta-sheet peptide (double hairpin) was sampled by equilibrium folding/unfolding molecular dynamics simulations for a total of 20 micros. The resulting one-dimensional free-energy profiles (FEPs) provide a detailed description of the free-energy basins and barriers for the folding reaction. The similarity of the FEPs obtained using the probability of folding before unfolding (pfold) or the mean first passage time supports the robustness of the procedure. The folded state and the most populated free-energy basins in the denatured state are described by the one-dimensional FEPs, which avoid the overlap of states present in the usual one- or two-dimensional projections. Within the denatured state, a basin with fluctuating helical conformations and a heterogeneous entropic state are populated near the melting temperature at about 11% and 33%, respectively. Folding pathways from the helical basin or enthalpic traps (with only one of the two hairpins formed) reach the native state through the entropic state, which is on-pathway and is separated by a low barrier from the folded state. A simplified equilibrium kinetic network based on the FEPs shows the complexity of the folding reaction and indicates, as augmented by additional analyses, that the basins in the denatured state are connected primarily by the native state. The overall folding kinetics shows single-exponential behavior because barriers between the non-native basins and the folded state have similar heights.

We derive a nonequilibrium thermodynamics identity (the "differential fluctuation theorem") that connects forward and reverse joint probabilities of nonequilibrium work and of arbitrary generalized coordinates corresponding to states of interest. This identity allows us to estimate the free energy difference between domains of these states. Our results follow from a general symmetry relation between averages over nonequilibrium forward and backward path functions derived by Crooks [Crooks, G. E. Phys. Rev. E 2000, 61, 2361-2366]. We show how several existing nonequilibrium thermodynamic identities can be obtained directly from the differential fluctuation theorem. We devise an approach for measuring conformational free energy differences, and we demonstrate its applicability to the analysis of molecular dynamics simulations by estimating the free energy difference between two conformers of the alanine dipeptide model system. We anticipate that these developments can be applied to the analysis of laboratory experiments.

We present a kinetic model for the walking of myosin V on actin under conditions of zero external force. The model includes three pathways and the termination of the processivity. Experimentally measured kinetic parameters are used in the model to obtain quantitative results. Using the model and associated parameters, we compute the proportion of the pathway containing an intermediate state, as well as the walking velocities and run lengths at various concentrations of ATP and ADP. The resulting trends agree with experimental data. The model explains the surprising experimental finding that myosin walks at a faster speed but for a shorter distance as the ATP concentration increases in the absence of ADP. It also suggests that under physiological condition ([ADP] approximately 12-50 microM), myosin walks with a higher speed and for longer distances when ATP is more abundant.

We have designed a model lattice protein that has two stable folded states, the lower free energy native state and a latent state of somewhat higher energy. The two states have a sizable part of their structures in common (two "alpha-helices") and differ in the content of "alpha-helices" and "beta-strands" in the rest of their structures; i.e. for the native state, this part is alpha-helical, and for the latent state it is composed of beta-strands. Thus, the lattice protein free energy surface mimics that of amyloidogenic proteins that form well organized fibrils under appropriate conditions. A Go-like potential was used and the folding process was simulated with a Monte Carlo method. To gain insight into the equilibrium free energy surface and the folding kinetics, we have combined standard approaches (reduced free energy surfaces, contact maps, time-dependent populations of the characteristic states, and folding time distributions) with a new approach. The latter is based on a principal coordinate analysis of the entire set of contacts, which makes possible the introduction of unbiased reaction coordinates and the construction of a kinetic network for the folding process. The system is found to have four characteristic basins, namely a semicompact globule, an on-pathway intermediate (the bifurcation basin), and the native and latent states. The bifurcation basin is shallow and consists of the structure common to the native and latent states, with the rest disorganized. On the basis of the simulation results, a simple kinetic model describing the transitions between the characteristic states was developed, and the rate constants for the essential transitions were estimated. During the folding process the system dwells in the bifurcation basin for a relatively short time before it proceeds to the native or latent state. We suggest that such a bifurcation may occur generally for proteins in which native and latent states have a sizable part of their structures in common. Moreover, there is the possibility of introducing changes in the system (e.g., mutations), which guide the system toward the native or misfolded state.

Measurements of time-resolved Stokes shifts on picosecond to nanosecond time scales have been used to probe the polar solvation dynamics of biological systems. Since it is difficult to decompose the measurements into protein and solvent contributions, computer simulations are useful to aid in understanding the details of the molecular behavior. Here we report the analysis of simulations of the electrostatic interactions of the rest of the protein and the solvent with 11 residues of the immunoglobulin binding domain B1 of protein G. It is shown that the polar solvation dynamics are position-dependent and highly heterogeneous. The contributions due to interactions with the protein and with the solvent are determined. The solvent contributions are found to vary from negligible after a few picoseconds to dominant on a scale of hundreds of picoseconds. The origin for the latter is found to involve coupled hydration and protein conformational dynamics. The resulting microscopic picture demonstrates that a wide range of possibilities have to be considered in the interpretation of time-resolved Stokes shift measurements.

A theoretical study showed that the base excision repair enzyme uracil-DNA glycosylase (UDG) exploits electrostatic interactions with backbone phosphate groups in the substrate for catalysis. Although experiments performed to test the calculated results confirmed the predicted importance of the -2, -1, and +1 phosphate groups, there was an apparent disagreement with regard to the +2 phosphate group. The calculations indicated that it made an important contribution, while experimentally, the effect of its deletion or neutralization was small. The +2 phosphate group interacts directly with an active site histidine (H148 in humans) in the crystal structure of UDG in complex with double-stranded (ds) DNA. We previously calculated that H148 has a strong anticatalytic effect due to its protonation, and here we use alchemical free energy simulations to estimate its site-specific pKa. We find that it is positively charged over the entire experimental pH range (4-10), so its deprotonation cannot compensate for deletion or neutralization of the +2 phosphate group. The free energy simulations are facilitated by an efficient charge-scaling procedure that allows quantitative correction for the implicit treatment of solvent far from the active site; improvements are made to that method to account carefully for differences in the truncation of electrostatic interactions in the contributing molecular-mechanical and continuum-electrostatic approaches. Additional simulations are used to demonstrate that the +2 phosphate group is fully solvent exposed in complexes with single-stranded DNA substrates like those used in the experiments. In contrast, it is well-structured and buried in the dsDNA complex used in the original simulations. Differences in solvent shielding thus account for the apparent lack of an effect observed experimentally upon neutralization or deletion of this group.

We show that the balanced minimum-cut procedure introduced in PNAS 2004, 101, 14766 can be reinterpreted as a method for solving the constrained optimization problem of finding the minimum cut among the cuts with a particular value of an additive function of the nodes on either side of the cut. Such an additive function (e.g., the partition function of the reactant region) can be used as a progress coordinate to determine a one-dimensional profile (FEP) of the free-energy surface of the protein-folding reaction as well as other complex reactions. The algorithm is based on the network (obtained from an equilibrium molecular dynamics simulation) that represents the calculated reaction behavior. The resulting FEP gives the exact values of the free energy as a function of the progress coordinate; i.e., at each value of the progress coordinate, the profile is obtained from the surface with the minimal partition function among the surfaces that divide the full free-energy surface between two chosen end points. In many cases, the balanced minimum-cut procedure gives results for only a limited set of points. An approximate method based on p(fold) is shown to provide the profile for a more complete set of values of the progress coordinate. Applications of the approach to model problems and to realistic systems (beta-hairpin of protein G, LJ38 cluster) are presented.

A density functional study of the hydrolysis reaction of phosphodiesters with a series of attacking nucleophiles in the gas phase and in solution is presented. The nucleophiles HOH, HO-, CH3OH, and CH3O- were studied in reactions with ethylene phosphate, 2'3'-ribose cyclic phosphate and in their neutral (protonated) and monoanionic forms. Stationary-point geometries for the reactions were determined at the density functional B3LYP/6-31++G(d,p) level followed by energy refinement at the B3LYP/6-311++G(3df,2p) level. Solvation effects were estimated by using a dielectric approximation with the polarizable continuum model (PCM) at the gas-phase optimized geometries. This series of reactions characterizes factors that influence the intrinsic reactivity of the model phosphate compounds, including the effect of nucleophile, protonation state, cyclic structure, and solvent. The present study of the in-line mechanism for phosphodiester hydrolysis, a reaction of considerable biological importance, has implications for enzymatic mechanisms. The analysis generally supports the associative mechanism for phosphate ester hydrolysis. The results highlight the importance for the reaction barrier of charge neutralization resulting from the protonation of the nonbridging phosphoryl oxygens and the role of internal hydrogen transfer in the gas-phase mechanism. It also shows that solvent stabilization has a profound influence on the relative barrier heights for the dianionic, monoanionic, and neutral reactions. The calculations provide a comprehensive data set for the in-line hydrolysis mechanisms that can be used for the development of improved semiempirical quantum models for phosphate hydrolysis reactions.

Dynamic disorder in proteins, as demonstrated by variations in single-molecule electron transfer rates, is investigated by molecular dynamics simulations. The potential of mean force for the fluctuating donor-acceptor distance is calculated for the NAD(P)H:flavin oxidoreductase (Fre) complex with flavin adenine dinucleotide (FAD) and is found to be in agreement with that estimated from electron transfer experiments. The calculated autocorrelation function of the distance fluctuations has a simple exponential behavior at low temperatures and stretched exponential behavior at higher temperatures on femtosecond to nanosecond time scales. This indicates that the calculated dynamic disorder arises from a wide range of trapping times in potential wells on the protein energy landscape and suggests a corresponding origin for the stretched exponential behavior observed experimentally on longer time scales.

The ubiquitin mutant UbG folding experiments of Sabelko et al., in which "strange kinetics" were observed, are interpreted in terms of a simple kinetic model. A minimal set of states consisting of a semicompact globule, two off-pathway traps, and the native state are included; the fully unfolded state is not considered because folding to the semicompact globule is fast. Both the low- and the high-temperature experiments of Sabelko et al. are fitted by a system of kinetic equations determining the transitions between these states. It is possible that cold- and heat-denaturated states of UbG are the basis of the off-pathway traps. The fits of the kinetic model to the experimental results provides an estimate of the rate constants for the various reaction channels and show how their contributions vary with temperature. Introduction of an on-pathway intermediate instead of one of the off-pathway traps does not lead to agreement with the experiments.

The hammerhead ribozyme is one of the best studied ribozymes, but it still presents challenges for our understanding of RNA catalysis. It catalyzes a transesterification reaction that converts a 5',3' diester to a 2',3' cyclic phosphate diester via an S(N)2 mechanism. Thus, the overall reaction corresponds to that catalyzed by bovine pancreatic ribonuclease. However, an essential distinguishing aspect is that metal ions are not involved in RNase catalysis but appear to be important in ribozymes. Although various techniques have been used to assign specific functions to metals in the hammerhead ribozyme, their number and roles in catalysis is not clear. Two recent theoretical studies on RNA catalysis examined the reaction mechanism of a single-metal-ion model. A two-metal-ion model, which is supported by experiment and based on ab initio and density functional theory calculations, is described here. The proposed mechanism of the reaction has four chemical steps with three intermediates and four transition states along the reaction pathway. Reaction profiles are calculated in the gas phase and in solution. The early steps of the reaction are found to be fast (with low activation barriers), and the last step, corresponding to the departure of the leaving group, is rate limiting. This two-metal-ion model differs from the models proposed previously in that the two metal ions function not only as Lewis acids but also as general acids/bases. Comparison with experiment shows good agreement with thermodynamic and kinetic data. A detailed analysis based on natural bond orbitals (NBOs) and natural energy decomposition (NEDA) provides insights into the role of metal ions and other factors important for catalysis.

Ligands such as CO, O(2), or NO are involved in the biological function of myoglobin. Here we investigate the energetics and dynamics of NO interacting with the Fe(II) heme group in native myoglobin using ab initio and molecular dynamics simulations. At the global minimum of the ab initio potential energy surface (PES), the binding energy of 23.4 kcal/mol and the Fe-NO structure compare well with the experimental results. Interestingly, the PES is found to exhibit two minima: There exists a metastable, linear Fe-O-N minimum in addition to the known, bent Fe-N-O global minimum conformation. Moreover, the T-shaped configuration is found to be a saddle point, in contrast to the corresponding minimum for NO interacting with Fe(III). To use the ab initio results for finite temperature molecular dynamics simulations, an analytical function was fitted to represent the Fe-NO interaction. The simulations show that the secondary minimum is dynamically stable up to 250 K and has a lifetime of several hundred picoseconds at 300 K. The difference in the topology of the heme-NO PES from that assumed previously (one deep, single Fe-NO minimum) suggests that it is important to use the full PES for a quantitative understanding of this system. Why the metastable state has not been observed in the many spectroscopic studies of myoglobin interacting with NO is discussed, and possible approaches to finding it are outlined.

The molecular chaperone HSP90 is an attractive target for chemotherapy because its activity is required for the functional maturation of a number of oncogenes. Among the known inhibitors, radicicol, a 14-member macrolide, stands out as the most potent. A molecular dynamics/minimization of radicicol showed that there were three low energy conformers of the macrocycle. The lowest of these is the bioactive conformation observed in the cocrystal structure of radicicol with HSP90. Corresponding conformational analyses of several known analogues gave a good correlation between the bioactivity and the energy of the bioactive conformer, relative to other conformers. Based on this observation, a number of proposed analogues were analyzed for their propensity to adopt the bioactive conformation prior to synthesis. This led to the identification of pochonin D, a recently isolated secondary metabolite of Pochonia chlamydosporia, as a potential inhibitor of HSP90. Pochonin D was synthesized using polymer-bound reagents and shown to be nearly as potent an HSP90 inhibitor as radicicol.

A 27-residue lattice heteropolymer subject to Monte Carlo dynamics on a simple cubic lattice is studied over a range of temperatures. Folding time distributions are used to obtain information concerning the details of folding kinetics. The results are compared with those from methods based on mean force surfaces expressed in terms of a reduced set of variables and on a disconnectivity graph for the same system. A detailed analysis of the folding trajectories is given, and the importance of dead-end traps in determining the folding time is demonstrated. We show that the calculated folding kinetics can be modeled by a system of kinetic equations, with the essential rate constants determined from the Monte Carlo simulations and the resulting folding time distributions. The kinetic equations make possible an analysis of the variation of the importance of different channels with temperature. In particular, we show that the presence of intermediates may be masked in the folding time distributions, with the mean folding time being independent of the height of the barrier between the intermediates and collapsed globule state of the system. This and other results demonstrate that care has to be used in interpreting experimental folding data in terms of the underlying kinetics. Correspondingly, simulations are shown to have to satisfy certain requirements to obtain proper sampling of the dead-end traps.

Previous work by Wand et al. (Nature 2001, 411, 501-504) showed that the NMR order parameters characterizing the amplitude of motion of protein side chains seemed to form a multimodal distribution. At the time, no detailed explanation of this at the molecular level was offered, yet three "classes" of motion were inferred. We have analyzed a larger published data set and found that, although the distribution is multimodal, the evidence for three classes is weak. More significantly, we have been able to provide a simple physical explanation for the distributions based on the results of molecular dynamics simulations. This result will aid in the interpretation of data from NMR dynamics experiments.

We present a method for determining the ensembles of native protein structures that result from the large fluctuations of low probability revealed by hydrogen-exchange experiments. The measured protection factors are used to bias Monte Carlo simulations to sample the structures of the exchange competent species. The approach is illustrated by its application to the case of alpha-lactalbumin.

A conformational analysis of a stereochemically complete set of peptide analogues based on a cis-enediol unit is presented. The cis-enediol unit, which can replace a two or a three amino acid segment of a peptide, contains two "side chains", four asymmetrical carbon atoms, and six free dihedral angles. To determine the accessible conformational space, the molecules are divided into three fragments, each containing two free dihedral angles. The energy surfaces are computed for all dihedral angle values, and the possible conformations of the cis-enediol unit analogues are built using all combinations of the surface minima. Such a "build-up" procedure, which is very fast, is able to reproduce 75% of the minima obtained from a full dihedral angle exploration of the conformational space. The cis-enediol unit minima are compared with the corresponding di- and tripeptide minima; all peptide minima can be closely matched by a cis-enediol unit minimum of low energy (less than 2.2 kcal/mol above the lowest energy conformer). However, there are low energy minima of the cis-enediol unit that have no corresponding minima in peptides. The results are shown to depend strongly on the chirality of the analogues. The ability of each of the stereoisomers to mimic natural peptides, evaluated by the present approach, is correlated with its experimental activity in a renin inhibition assay.

There has recently been considerable interest in using NMR spectroscopy to identify ligand binding sites of macromolecules. In particular, a modular approach has been put forward by Fesik et al. (Shuker, S. B.; Hajduk, P. J.; Meadows, R. P.; Fesik, S. W. Science 1996, 274, 1531-1534) in which small ligands that bind to a particular target are identified in a first round of screening and subsequently linked together to form ligands of higher affinity. Similar strategies have also been proposed for in silico drug design, where the binding sites of small chemical groups are identified, and complete ligands are subsequently assembled from different groups that have favorable interactions with the macromolecular target. In this paper, we compare experimental and computational results on a selected target (FKBP12). The binding sites of three small ligands ((2S)1-acetylprolinemethylester, 1-formylpiperidine, 1-piperidinecarboxamide) in FKBP12 were identified independently by NMR and by computational methods. The subsequent comparison of the experimental and computational data showed that the computational method identified and ranked favorably ligand positions that satisfy the experimental NOE constraints.

In recent years, molecular dynamics simulations of biomolecular free energy differences have benefited from significant methodological advances and increased computer power. Applications to molecular recognition provide an understanding of the interactions involved that goes beyond, and is an important complement to, experimental studies. Poisson-Boltzmann electrostatic models provide a faster and simpler free energy method in cases where electrostatic interactions are important. We illustrate both molecular dynamics and Poisson-Boltzmann methods with a detailed study of amino acid recognition by aspartyl-tRNA synthetase, whose specificity is important for maintaining the integrity of the genetic code.

Knowledge of the pK(a) of phosphoranes is important for the interpretation of phosphate ester hydrolysis. Calculated pK(a)'s of the model phosphorane, ethylene phosphorane, are reported. The method of calculation is based on the use of dimethyl phosphate as a reference state for evaluating relative pK(a) values, and on the optimization of the oxygen and acidic hydrogen van der Waals radii to give reasonable pK(1)(a), pK(2)(a), and pK(3)(a) for phosphoric acid in solution. Density functional theory is employed to calculate the gas-phase protonation energies, and continuum dielectric methods are used to determine the solvation corrections. The calculated pK(1)(a) and p(2)(a) for the model phosphorane are 7.9 and 14.3, respectively. These values are within the range of proposed experimental values, 6.5-11.0 for pK(1)(a), and 11.3-15.0 for pK(2)(a). The mechanistic implications of the calculated pK(a)'s are discussed.

The role of tunneling for two proton-transfer steps in the reactions catalyzed by triosephosphate isomerase (TIM) has been studied. One step is the rate-limiting proton transfer from Calpha in the substrate to Glu 165, and the other is an intrasubstrate proton transfer proposed for the isomerization of the enediolate intermediate. The latter, which is not important in the wild-type enzyme but is a useful model system because of its simplicity, has also been examined in the gas phase and in solution. Variational transition-state theory with semiclassical ground-state tunneling was used for the calculation with potential energy surface determined by an AM1 method specifically parametrized for the TIM system. The effect of tunneling on the reaction rate was found to be less than a factor of 10 at room temperature; the tunneling becomes more important at lower temperature, as expected. The imaginary frequency (barrier) mode and modes that have large contributions to the reaction path curvature are localized on the atoms in the active site, within 4 A of the substrate. This suggests that only a small number of atoms that are close to the substrate and their motions (e.g., donor-acceptor vibration) directly determine the magnitude of tunneling. Atoms that are farther away influence the effect of tunneling indirectly by modulating the energetics of the proton transfer. For the intramolecular proton transfer, tunneling was found to be most important in the gas phase, to be similar in the enzyme, and to be the smallest in water. The major reason for this trend is that the barrier frequency is substantially lower in solution than in the gas phase and enzyme; the broader solution barrier is caused by the strong electrostatic interaction between the highly charged solute and the polar solvent molecules. Analysis of isotope effects showed that the conventional Arrenhius parameters are more useful as experimental criteria for determining the magnitude of tunneling than the widely used Swain-Schaad exponent (SSE). For the primary SSE, although values larger than the transition-state theory limit (3.3) occur when tunneling is included, there is no clear relationship between the calculated magnitudes of tunneling and the SSE. Also, the temperature dependence of the primary SSE is rather complex; the value of SSE tends to decrease as the temperature is lowered (i.e., when tunneling becomes more significant). For the secondary SSE, the results suggest that it is more relevant for evaluating the "coupled motion" between the secondary hydrogen and the reaction coordinate than the magnitude of tunneling. Although tunneling makes a significant contribution to the rate of proton transfer, it appears not to be a major aspect of the catalysis by TIM at room temperature; i.e., the tunneling factor of 10 is "small" relative to the overall rate acceleration by 10(9). For the intramolecular proton transfer, the tunneling in the enzyme is larger by a factor of 5 than in solution.

Computational methods were used to design structure-based combinatorial libraries of antipicornaviral capsid-binding ligands. The multiple copy simultaneous search (MCSS) program was employed to calculate functionality maps for many diverse functional groups for both the poliovirus and rhinovirus capsid structures in the region of the known drug binding pocket. Based on the results of the MCSS calculations, small combinatorial libraries consisting of 10s or 100s of three-monomer compounds were designed and synthesized. Ligand binding was demonstrated by a noncell-based mass spectrometric assay, a functional immuno-precipitation assay, and crystallographic analysis of the complexes of the virus with two of the candidate ligands. The P1/Mahoney poliovirus strain was used in the experimental studies. A comparison showed that the MCSS calculations had correctly identified the observed binding site for all three monomer units in one ligand and for two out of three in the other ligand. The correct central monomer position in the second ligand was reproduced in calculations in which the several key residues lining the pocket were allowed to move. This study validates the computational methodology. It also illustrates that subtle changes in protein structure can lead to differences in docking results and points to the importance of including target flexibility, as well as ligand flexibility, in the design process.

The mechanism of the alkaline hydrolysis of phosphate and sulfate esters is of great interest. Ab initio quantum mechanical calculations and dielectric continuum methods are used to investigate the effect of the solvent on the associative/dissociative and the in-line/sideways character of the hydrolysis reaction of ethylene sulfate (ES) and ethylene phosphate (EP(-)), and their acyclic counterparts, dimethyl sulfate (DMS) and dimethyl phosphate (DMP(-)). The gas-phase reaction coordinates are determined by Hartree-Fock and density functional theory. For ES, the reaction coordinate in solution is determined; for the other three reactions only the transition state in solution is obtained. The alterations in the reaction induced by solvent are interpreted by use of the Hammond and anti-Hammond postulates.

Three mechanisms proposed for the triosephosphate isomerase (TIM) catalyzed reactions were studied with the QM/MM approach using B3LYP/6-31+G(d,p) as the QM method. The two pathways that involve an enediol species were found to give similar values for the barriers and the calculated rates are in satisfactory agreement with experiment. By contrast, the mechanism that involves intramolecular proton transfer in the enediolate was found to be energetically unfavorable due to electrostatic interactions with His 95, a conserved residue in TIM from different organisms. A perturbation analysis was used to determine the residues that make the major contribution to catalysis.

The selection of appropriate descriptors is an important step in the successful formulation of quantitative structure-activity relationships (QSARs). This paper compares a number of feature selection routines and mapping methods that are in current use. They include forward stepping regression (FSR), genetic function approximation (GFA), generalized simulated annealing (GSA), and genetic neural network (GNN). On the basis of a data set of steroids of known in vitro binding affinity to the progsterone receptor, a number of QSAR models are constructed. A comparison of the predictive qualities for both training and test compounds demonstrates that the GNN protocol achieves the best results among the 2D QSAR that are considered. Analysis of the choice of descriptors by the GNN method shows that the results are consistent with established SARs on this series of compounds.

The utility of genetic neural network (GNN) to obtain quantitative structure-activity relationships (QSAR) from molecular similarity matrices is described. In this application, the corticosteroid-binding globulin (CBG) binding affinity of the well-known steroid data set is examined. Excellent predictivity can be obtained through the use of either electrostatic or shape properties alone. Statistical validation using a standard randomization test indicates that the results are not due to chance correlations. Application of GNN on the combined electrostatic and shape matrix produces a six-descriptor model with a cross-validated r2 value of 0.94. The model is superior to those obtained from partial least-squares and genetic regressions, and it also compares favorably with the results for the same data set from other established 3D QSAR methods. The theoretical basis for the use of molecular similarity in QSAR is discussed.

Validation of a method that uses a genetic neural network with electrostatic and steric similarity matrices (SM/GNN) to obtain quantitative structure-activity relationships (QSARs) is performed with eight data sets. Biological and physicochemical properties from a broad range of chemical classes are correlated and predicted using this technique. Quantitatively the results compare favorably with the benchmarks obtained by a number of well-established QSAR methods; qualitatively the models are consistent with the published descriptions on the relative contribution of steric and electrostatic factors. The results demonstrate the general utility of this method in deriving QSARs. The implication of the importance of molecular alignment and possible methodological improvements are discussed.

A novel tool, called a genetic neural network (GNN), has been developed for obtaining quantitative structure-activity relationships (QSAR) for high-dimensional data sets (J. Med. Chem. 1996, 39, 1521-1530). The GNN method uses a neural network to correlate activity with descriptors that are preselected by a genetic algorithm. To provide an extended test of the GNN method, the data on 57 benzodiazepines given by Maddalena and Johnston (MJ; J. Med. Chem. 1995, 38, 715-724) have been examined with an enhanced version of GNN, and the results are compared with the excellent QSAR of MJ. The problematic steepest descent training has been replaced by the scaled conjugate gradient algorithm. This leads to a substantial gain in performance in both robustness of prediction and speed of computation. The cross-validation GNN simulation and the subsequent run based on an unbiased and more efficient protocol led to the discovery of other 10-descriptor QSARs that are superior to the best model of MJ based on backward elimination selection and neural network training. Results from a series of GNNs with a different number of inputs showed that a neural network with fewer inputs can produce QSARs as good as or even better than those with higher dimensions. The top-ranking models from a GNN simulation using only six input descriptors are presented, and the chemical significance of the chosen descriptors is discussed. The statistical significance of these GNN QSARs is validated. The best QSARs are used to provide a graphical tool that aids the design of new drug analogues. By replacing functional groups at the 7- and 2'-positions with ones that have optimal substituent parameters, a number of new benzodiazepines with high potency are predicted.

A new hybrid method (GNN) combining a genetic algorithm and an artificial neural network has been developed for quantitative structure-activity relationship (QSAR) studies. A suitable set of molecular descriptors are selected by a genetic algorithm. This set serves as input to a neural network, in which model-free mapping of multivariate data is performed. Multiple predictors are generated that are superior to results obtained from previous studies of the Selwood data set, which is used to test the method. The neural network technique provides a graphical description of the functional form of the descriptors that play an important role in determining drug activity. This can serve as an aid in future design of drug analogues. The effectiveness of GNN is tested by comparing its results with a benchmark obtained by exhaustive enumeration. Different fitness strategies that tune the evolution of genetic models are examined, and QSARs with higher predictiveness are found. From these results, a composite model is constructed by averaging predictions from several high-ranking models. The predictions of the resulting QSAR should be more reliable than those derived from a single predictor because it makes greater use of information and also permits error estimation. An analysis of the sets of descriptors selected by GNN shows that it is essential to have one each for the steric, electrostatic, and hydrophobic attributes of a drug candidate to obtain a satisfactory QSAR for this data set. This type of result is expected to be of general utility in designing and understanding QSAR.

A detailed mechanism for the catalysis of prolyl isomerization by the rotamase enzyme FKBP is proposed on the basis of a model constructed from the known structure of the FK506/FKBP complex. The model substrate is bound as a type VIa proline turn with the ends exposed to permit longer polypeptide chains (e.g., protein loops) to act as substrates. An ab initio potential for the isomerized imide bond is combined with a molecular mechanics representation of the rest of the system to calculate the reaction path. The resulting activation energy for the enzymatic cis-->trans isomerization is equal to about 6 kcal/mol, in good agreement with experiment. The lowering of the barrier relative to the solution value of 19 kcal/mol is found to arise from a combination of desolvation of the imide carbonyl, ground-state destabilization, substrate autocatalysis, and preferential transition-state binding. Minimal rearrangements are required in the enzyme and the substrate along the reaction path. The enzyme residues that participate in catalysis agree with the available mutation data. The type VIa turn model corresponds to a sequence-specific structural motif commonly found on the surface of proteins. It is likely to have a role in the formation of protein complexes with FKBP-like domains that function as foldases or chaperones.