Search

Ribonucleic acid (RNA) molecules can adopt a variety of secondary and tertiary structures in solution, with stem-loops being one of the more common motifs. Here, we present a systematic analysis of 15 RNA stem-loop sequences simulated with molecular dynamics simulations in an implicit solvent environment. Analysis of RNA cluster ensembles showed that the stem-loop structures can generally adopt the A-form RNA in the stem region. Loop structures are more sensitive, and experimental structures could only be reproduced with modification of CH···O interactions in the force field, combined with an implicit solvent nonpolar correction to better model base stacking interactions. Accurately modeling RNA with current atomistic physics-based models remains challenging, but the RNA systems studied herein may provide a useful benchmark set for testing other RNA modeling methods in the future.

One of the key challenges of k-means clustering is the seed selection or the initial centroid estimation since the clustering result depends heavily on this choice. Alternatives such as k-means++ have mitigated this limitation by estimating the centroids using an empirical probability distribution. However, with high-dimensional and complex data sets such as those obtained from molecular simulation, k-means++ fails to partition the data in an optimal manner. Furthermore, stochastic elements in all flavors of k-means++ will lead to a lack of reproducibility. K-means N-Ary Natural Initiation (NANI) is presented as an alternative to tackle this challenge by using efficient n-ary comparisons to both identify high-density regions in the data and select a diverse set of initial conformations. Centroids generated from NANI are not only representative of the data and different from one another, helping k-means to partition the data accurately, but also deterministic, providing consistent cluster populations across replicates. From peptide and protein folding molecular simulations, NANI was able to create compact and well-separated clusters as well as accurately find the metastable states that agree with the literature. NANI can cluster diverse data sets and be used as a standalone tool or as part of our MDANCE clustering package.

Virtual screening (VS) involves generation of poses for a library of ligands and ranking using simplified energy functions and limited flexibility. Top-scored poses are used to rank and prioritize ligands. Here, we adapt the reservoir replica exchange molecular dynamics (res-REMD) method to rerank poses generated through VS. REMD simulations are carried out but with occasional Monte Carlo jumps to alternate VS-generated poses using a Metropolis criterion. The simulations converge within 10 ns for all systems, generating populations of alternate poses in the context of fully flexible ligand and protein side chains. The protocol is applied to four model protein–ligand complexes, where DOCK resulted in two successes and two scoring failures. In all four systems, the most populated cluster from the final ensemble exhibits high similarity to the crystallographic pose with ligand RMSD values under 2.0 Å. Both DOCK failures were rescued. For one DOCK success, the protocol identified the correct pose but also sampled an alternate pose at equal probability. Opportunities for future improvements and extensions are discussed.

Generating precise ensembles is commonly a prerequisite to understand the energetics of biological processes using Molecular Dynamics (MD) simulations. Previously, we have shown how unweighted reservoirs built from high temperature MD simulations can accelerate convergence of Boltzmann-weighted ensembles by at least 10× with the Reservoir Replica Exchange MD (RREMD) method. Therefore, in this work, we explore whether an unweighted structure reservoir generated with one Hamiltonian (solute force field plus solvent model) can be reused to quickly generate accurately weighted ensembles for Hamiltonians other than the one that was used to generate the reservoir. We also extended this methodology to rapidly estimate the effects of mutations on peptide stability by using a reservoir of diverse structures obtained from wild-type simulations. These results suggest that structures generated via fast methods such as coarse-grained models or structures predicted by Rosetta or deep learning approaches could be integrated into a reservoir to accelerate generation of ensembles using more accurate representations.

RNA is a key participant in many biological processes, but studies of RNA using computer simulations lag behind those of proteins, largely due to less-developed force fields and the slow dynamics of RNA. Generating converged RNA ensembles for force field development and other studies remains a challenge. In this study, we explore the ability of replica exchange molecular dynamics to obtain well-converged conformational ensembles for two RNA hairpin systems in an implicit solvent. Even for these small model systems, standard REMD remains computationally costly, but coupling to a pre-generated structure library using the reservoir REMD approach provides a dramatic acceleration of ensemble convergence for both model systems. Such precise ensembles could facilitate RNA force field development and validation and applications of simulation to more complex RNA systems. The advantages and remaining challenges of applying R-REMD to RNA are investigated in detail.

Inspired by the role of cell-surface glycoproteins as coreceptors for pathogens, we report the development of GlycoGrip: a glycopolymer-based lateral flow assay for detecting SARS-CoV-2 and its variants. GlycoGriputilizes glycopolymers for primary capture and antispike antibodies labeled with gold nanoparticles for signal-generating detection. A lock-step integration between experiment and computation has enabled efficient optimization of GlycoGriptest strips which can selectively, sensitively, and rapidly detect SARS-CoV-2 and its variants in biofluids. Employing the power of the glycocalyx in a diagnostic assay has distinct advantages over conventional immunoassays as glycopolymers can bind to antigens in a multivalent capacity and are highly adaptable for mutated strains. As new variants of SARS-CoV-2 are identified, GlycoGripwill serve as a highly reconfigurable biosensor for their detection. Additionally, via extensive ensemble-based docking simulations which incorporate protein and glycan motion, we have elucidated important clues as to how heparan sulfate and other glycocalyx components may bind the spike glycoprotein during SARS-CoV-2 host-cell infection. GlycoGripis a promising and generalizable alternative to costly, labor-intensive RT-PCR, and we envision it will be broadly useful, including for rural or low-income populations that are historically undertested and under-reported in infection statistics.

The SARS-CoV-2 coronavirus is an enveloped, positive-sense single-stranded RNA virus that is responsible for the COVID-19 pandemic. The spike is a class I viral fusion glycoprotein that extends from the viral surface and is responsible for viral entry into the host cell and is the primary target of neutralizing antibodies. The receptor binding domain (RBD) of the spike samples multiple conformations in a compromise between evading immune recognition and searching for the host-cell surface receptor. Using atomistic simulations of the glycosylated wild-type spike in the closed and 1-up RBD conformations, we map the free energy landscape for RBD opening and identify interactions in an allosteric pocket that influence RBD dynamics. The results provide an explanation for experimental observation of increased antibody binding for a clinical variant with a substitution in this pocket. Our results also suggest the possibility of allosteric targeting of the RBD equilibrium to favor open states via binding of small molecules to the hinge pocket. In addition to potential value as experimental probes to quantify RBD conformational heterogeneity, small molecules that modulate the RBD equilibrium could help explore the relationship between RBD opening and S1 shedding.

In silicoscreening of drug–target interactions is a key part of the drug discovery process. Changes in the drug scaffold via contraction or expansion of rings, the breaking of rings, and the introduction of cyclic structures from acyclic structures are commonly applied by medicinal chemists to improve binding affinity and enhance favorable properties of candidate compounds. These processes, commonly referred to as scaffold hopping, are challenging to model computationally. Although relative binding free energy (RBFE) calculations have shown success in predicting binding affinity changes caused by perturbing R-groups attached to a common scaffold, applications of RBFE calculations to modeling scaffold hopping are relatively limited. Scaffold hopping inevitably involves breaking and forming bond interactions of quadratic functional forms, which is highly challenging. A novel method for handling ring opening/closure/contraction/expansion and linker contraction/expansion is presented here. To the best of our knowledge, RBFE calculations on linker contraction/expansion have not been previously reported. The method uses auxiliary restraints to hold the atoms at the ends of a bond in place during the breaking and forming of the bonds. The broad applicability of the method was demonstrated by examining perturbations involving small-molecule macrocycles and mutations of proline in proteins. High accuracy was obtained using the method for most of the perturbations studied. The rigor of the method was isolated from the force field by validating the method using relative and absolute hydration free energy calculations compared to standard simulation results. Unlike other methods that rely on λ-dependent functional forms for bond interactions, the method presented here can be employed using modern molecular dynamics software without modification of codes or force field functions.

Proteins fold on relatively smooth free energy landscapes which are biased toward the native state, but even simple topologies which fold rapidly can experience roughness on their free energy landscape. The details of these interactions are difficult to elucidate experimentally. Closely related to the problem of deciphering the details of the free energy landscape is the problem of defining the interactions in the denatured state ensemble (DSE) which is populated under native conditions, that is, under conditions where the native state is stable. The DSE of many proteins deviates from random coil models, but quantifying and defining the energetics of the transiently populated interactions in this ensemble is extremely challenging. Characterization of the DSE of proteins which fold to compact structures is also relevant to studies of intrinsically disordered proteins (IDPs) since interactions in the dynamic ensemble populated by IDPs can modulate their behavior. Here we show how experimental thermodynamic and pKameasurements can be combined with computational thermodynamic integration to quantify interactions in the DSE. We show that non-native side chain interactions can stabilize native backbone structure in the DSE and demonstrate that that even rapidly folding proteins can form energetically significant non-native interactions in their DSE. As an example, we characterize a non-native salt bridge that stabilizes local native backbone structure in the DSE of a widely studied fast-folding protein, the villin headpiece helical domain. The combined computational experimental approach is applicable to other protein unfolded states and provides insight that is impossible to obtain with either method alone.

The coronavirus disease 2019 (COVID-19) pandemic has swept over the world in the past months, causing significant loss of life and consequences to human health. Although numerous drug and vaccine development efforts are underway, there are many outstanding questions on the mechanism of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) viral association to angiotensin-converting enzyme 2 (ACE2), its main host receptor, and host cell entry. Structural and biophysical studies indicate some degree of flexibility in the viral extracellular spike glycoprotein and at the receptor-binding domain (RBD)-receptor interface, suggesting a role in infection. Here, we perform explicitly solvated, all-atom, molecular dynamics simulations of the glycosylated, full-length, membrane-bound ACE2 receptor in both an apo and spike RBD-bound state to probe the intrinsic dynamics of the ACE2 receptor in the context of the cell surface. A large degree of fluctuation in the full-length structure is observed, indicating hinge bending motions at the linker region connecting the head to the transmembrane helix while still not disrupting the ACE2 homodimer or ACE2-RBD interfaces. This flexibility translates into an ensemble of ACE2 homodimer conformations that could sterically accommodate binding of the spike trimer to more than one ACE2 homodimer and suggests a mechanical contribution of the host receptor toward the large spike conformational changes required for cell fusion. This work presents further structural and functional insights into the role of ACE2 in viral infection that can potentially be exploited for the rational design of effective SARS-CoV-2 therapeutics.

Temperature replica exchange molecular dynamics (REMD) is a widely used enhanced sampling method for accelerating biomolecular simulations. During the past 2 decades, several variants of REMD have been developed to further improve the rate of conformational sampling of REMD. One such variant, reservoir REMD (RREMD), was shown to improve the rate of conformational sampling by around 5–20×. Despite the significant increase in the sampling speed, RREMD methods have not been widely used because of the difficulties in building the reservoir and also because of the code not being available on the graphics processing units (GPUs). In this work, we ported the Amber RREMD code onto GPUs making it 20× faster than the central processing unit code. Then, we explored protocols for building Boltzmann-weighted reservoirs as well as non-Boltzmann reservoirs and tested how each choice affects the accuracy of the resulting RREMD simulations. We show that, using the recommended protocols outlined here, RREMD simulations can accurately reproduce Boltzmann-weighted ensembles obtained by much more expensive conventional temperature-based REMD simulations, with at least 15× faster convergence rates even for larger proteins (>50 amino acids) compared to conventional REMD.

Light-activated protein domains provide a convenient, modular, and genetically encodable sensor for optogenetics and optobiology. Although these domains have now been deployed in numerous systems, the precise mechanism of photoactivation and the accompanying structural dynamics that modulate output domain activity remain to be fully elucidated. In the C-terminal light–oxygen–voltage (LOV) domain of plant phototropins (LOV2), blue light activation leads to formation of an adduct between a conserved Cys residue and the embedded FMN chromophore, rotation of a conserved Gln (Q513), and unfolding of a helix (Jα-helix) which is coupled to the output domain. In the present work, we focus on the allosteric pathways leading to Jα helix unfolding in Avena sativaLOV2 (AsLOV2) using an interdisciplinary approach involving molecular dynamics simulations extending to 7 μs, time-resolved infrared spectroscopy, solution NMR spectroscopy, and in-cell optogenetic experiments. In the dark state, the side chain of N414 is hydrogen bonded to the backbone N–H of Q513. The simulations predict a lever-like motion of Q513 after Cys adduct formation resulting in a loss of the interaction between the side chain of N414 and the backbone C═O of Q513, and formation of a transient hydrogen bond between the Q513 and N414 side chains. The central role of N414 in signal transduction was evaluated by site-directed mutagenesis supporting a direct link between Jα helix unfolding dynamics and the cellular function of the Zdk2-AsLOV2 optogenetic construct. Through this multifaceted approach, we show that Q513 and N414 are critical mediators of protein structural dynamics, linking the ultrafast (sub-ps) excitation of the FMN chromophore to the microsecond conformational changes that result in photoreceptor activation and biological function.

Molecular dynamics (MD) simulations have become increasingly popular in studying the motions and functions of biomolecules. The accuracy of the simulation, however, is highly determined by the molecular mechanics (MM) force field (FF), a set of functions with adjustable parameters to compute the potential energies from atomic positions. However, the overall quality of the FF, such as our previously published ff99SB and ff14SB, can be limited by assumptions that were made years ago. In the updated model presented here (ff19SB), we have significantly improved the backbone profiles for all 20 amino acids. We fit coupled φ/ψ parameters using 2D φ/ψ conformational scans for multiple amino acids, using as reference data the entire 2D quantum mechanics (QM) energy surface. We address the polarization inconsistency during dihedral parameter fitting by using both QM and MM in aqueous solution. Finally, we examine possible dependency of the backbone fitting on side chain rotamer. To extensively validate ff19SB parameters, and to compare to results using other Amber models, we have performed a total of ∼5 ms MD simulations in explicit solvent. Our results show that after amino-acid-specific training against QM data with solvent polarization, ff19SB not only reproduces the differences in amino-acid-specific Protein Data Bank (PDB) Ramachandran maps better but also shows significantly improved capability to differentiate amino-acid-dependent properties such as helical propensities. We also conclude that an inherent underestimation of helicity is present in ff14SB, which is (inexactly) compensated for by an increase in helical content driven by the TIP3P bias toward overly compact structures. In summary, ff19SB, when combined with a more accurate water model such as OPC, should have better predictive power for modeling sequence-specific behavior, protein mutations, and also rational protein design. Of the explicit water models tested here, we recommend use of OPC with ff19SB.

Intrinsically disordered proteins (IDPs) play important roles in biology, but little is known about the energetics of their inter-residue interactions. Methods that have been successfully applied to analyze the energetics of globular proteins are not applicable to the fluctuating partially ordered ensembles populated by IDPs. A combined computational experimental strategy is introduced for analyzing the energetic role of individual residues in the free state of IDPs. The approach combines experimental measurements of the binding of wild-type and mutant IDPs to their partners with alchemical free energy calculations of the structured complexes. These data allow quantitative information to be deduced about the free state via a thermodynamic cycle. The approach is validated by the analysis of the effects of mutations upon the binding free energy of the ovomucoid inhibitor third binding domain to its partners and is applied to the C-terminal domain of the measles virus nucleoprotein, a 125-residue IDP involved in the RNA transcription and replication of measles virus. The analysis reveals significant inter-residue interactions in the unbound IDP and suggests a biological role for them. The work demonstrates that advances in force fields and computational hardware have now led to the point where it is possible to develop methods, which integrate experimental and computational techniques to reveal insights that cannot be studied using either technique alone.

We present a fast implementation of the nudged elastic band (NEB) method into the particle mesh Ewald molecular dynamics module of the Amber software package for both central processing units (CPU) and graphics processing units (GPU). The accuracy of the new implementation has been validated for three cases: a conformational change of alanine dipeptide, the α-helix to β-sheet transition in polyalanine, and a large conformational transition in the human 8-oxoguanine–DNA glycosylase with DNA complex (OGG1–DNA). Timing benchmark tests were performed on the explicitly solvated OGG1–DNA system containing ∼50 000 atoms. The GPU-optimized implementation of NEB achieves a more than two orders of magnitude speedup compared with the previous CPU implementation performed with a two-core CPU processor. The speed and scalable features of this implementation will enable NEB applications on larger and more complex systems.

Phosphorylated amino acids are involved in many cell regulatory networks; proteins containing these post-translational modifications are widely studied both experimentally and computationally. Simulations are used to investigate a wide range of structural and dynamic properties of biomolecules, such as ligand binding, enzyme-reaction mechanisms, and protein folding. However, the development of force field parameters for the simulation of proteins containing phosphorylated amino acids using the Amber program has not kept pace with the development of parameters for standard amino acids, and it is challenging to model these modified amino acids with accuracy comparable to proteins containing only standard amino acids. In particular, the popular ff14SB and ff19SB models do not contain parameters for phosphorylated amino acids. Here, the dihedral parameters for the side chains of the most common phosphorylated amino acids are trained against reference data from QM calculations adopting the ff14SB approach, followed by validation against experimental data. Library files and corresponding parameter files are provided, with versions that are compatible with both ff14SB and ff19SB.

The precise prediction of Major Histocompatibility Complex (MHC)-peptide complex structures is pivotal for understanding cellular immune responses and advancing vaccine design. In this study, we enhanced AlphaFold’s capabilities by fine-tuning it with a specialized dataset consisting of exclusively high-resolution Class I MHC-peptide crystal structures. This tailored approach aimed to address the generalist nature of AlphaFold’s original training, which, while broad-ranging, lacked the granularity necessary for the high-precision demands of Class I MHC-peptide interaction prediction. A comparative analysis was conducted against the homology-modeling-based method Pandora, as well as the AlphaFold multimer model. Our results demonstrate that our fine-tuned model outperforms others in terms of RMSD (median value for C-alpha atoms for peptides is 0.66 Å) and also provides enhanced predicted lDDT scores, offering a more reliable assessment of the predicted structures. These advances have substantial implications for computational immunology, potentially accelerating the development of novel therapeutics and vaccines by providing a more precise computational lens through which to view MHC-peptide interactions.

The role of ribonucleic acid (RNA) in biology continues to grow, but insight into important aspects of RNA behavior is lacking, such as dynamic structural ensembles in different environments, how flexibility is coupled to function, and how function might be modulated by small molecule binding. In the case of proteins, much progress in these areas has been made by complementing experiments with atomistic simulations, but RNA simulation methods and force fields are less mature. It remains challenging to generate stable RNA simulations, even for small systems where well-defined, thermostable structures have been established by experiments. Many different aspects of RNA energetics have been adjusted in force fields, seeking improvements that are transferable across a variety of RNA structural motifs. In this work, the role of weak CH···O interactions is explored, which are ubiquitous in RNA structure but have received less attention in RNA force field development. By comparing data extracted from high-resolution RNA crystal structures to energy profiles from quantum mechanics and force field calculations, it is shown that CH···O interactions are overly repulsive in the widely used Amber RNA force fields. A simple, targeted adjustment of CH···O repulsion that leaves the remainder of the force field unchanged was developed. Then, the standard and modified force fields were tested using molecular dynamics (MD) simulations with explicit water and salt, amassing over 300 μs of data for multiple RNA systems containing important features such as the presence of loops, base stacking interactions as well as canonical and noncanonical base pairing. In this work and others, standard force fields lead to reproducible unfolding of the NMR-based structures. Including a targeted CH···O adjustment in an otherwise identical protocol dramatically improves the outcome, leading to stable simulations for all RNA systems tested.

In this paper, we present differentiable solvent-accessible surface area (dSASA), an exact geometric method to calculate SASA analytically along with atomic derivatives on GPUs. The atoms in a molecule are first assigned to tetrahedra in groups of four atoms by Delaunay tetrahedralization adapted for efficient GPU implementation, and the SASA values for atoms and molecules are calculated based on the tetrahedralization information and inclusion–exclusion method. The SASA values from the numerical icosahedral-based method can be reproduced with >98% accuracy for both proteins and RNAs. Having been implemented on GPUs and incorporated into AMBER, we can apply dSASA to implicit solvent molecular dynamics simulations with the inclusion of this nonpolar term. The current GPU version of GB/SA simulations has been accelerated up to nearly 20-fold compared to the CPU version, outperforming LCPO, a commonly used, fast algorithm for calculating SASA, as the system size increases. While we focus on the accuracy of the SASA calculations for proteins and nucleic acids, we also demonstrate stable GB/SA MD mini-protein simulations.

We propose a pairwise and readily parallelizable SASA-based nonpolar solvation approach for protein simulations, inspired by our previous pairwise GB polar solvation model development. In this work, we developed a novel function to estimate the atomic and molecular SASAs of proteins, which results in comparable accuracy as the LCPO algorithm in reproducing numerical icosahedral-based SASA values. Implemented in Amber software and tested on consumer GPUs, our pwSASA method reasonably reproduces LCPO simulation results, but accelerates MD simulations up to 30 times compared to the LCPO implementation, which is greatly desirable for protein simulations facing sampling challenges. The value of incorporating the nonpolar term in implicit solvent simulations is explored on a peptide fragment containing the hydrophobic core of HP36 and evaluating thermal stability profiles of four small proteins.

The enoyl-ACP reductase (ENR) catalyzes the last reaction in the elongation cycle of the bacterial type II fatty acid biosynthesis (FAS-II) pathway. While the FabI ENR is a well-validated drug target in organisms such as Mycobacterium tuberculosisand Staphylococcus aureus, alternate ENR isoforms have been discovered in other pathogens, including the FabV enzyme that is the sole ENR in Yersinia pestis(ypFabV). Previously, we showed that the prototypical ENR inhibitor triclosan was a poor inhibitor of ypFabV and that inhibitors based on the 2-pyridone scaffold were more potent [Hirschbeck, M. (2012) Structure20(1), 89–100]. These studies were performed with the T276S FabV variant. In the work presented here, we describe a detailed examination of the mechanism and inhibition of wild-type ypFabV and the T276S variant. The T276S mutation significantly reduces the affinity of diphenyl ether inhibitors for ypFabV (20-fold → 100-fold). In addition, while T276S ypFabV generally displays an affinity for 2-pyridone inhibitors higher than that of the wild-type enzyme, the 4-pyridone scaffold yields compounds with similar affinity for both wild-type and T276S ypFabV. T276 is located at the N-terminus of the helical substrate-binding loop, and structural studies coupled with site-directed mutagenesis reveal that alterations in this residue modulate the size of the active site portal. Subsequently, we were able to probe the mechanism of time-dependent inhibition in this enzyme family by extending the inhibition studies to include P142W ypFabV, a mutation that results in a gain of slow-onset inhibition for the 4-pyridone PT156.

AmberTools is a free and open-source collection of programs used to set up, run, and analyze molecular simulations. The newer features contained within AmberTools23 are briefly described in this Application note.

In this application note, we describe a tool which we developed to help structural biologists who study the SARS-CoV-2 spike glycoprotein. There are more than 500 structures of this protein available in the Protein Data Bank. These structures are available in different flavors: wild type spike, different variants, 2P substitutions, structures with bound antibodies, structures with Receptor Binding Domains (RBD) in closed or open conformation, etc. Understanding differences between these structures could provide insight into how the spike structure changes in different variants or upon interaction with different molecules such as receptors or antibodies. However, inconsistencies among deposited structures, such as different chain or sequence numbering, hamper a straightforward comparison of all structures. The tool described in this note fixes those chain inconsistencies and calculates the distribution of the requested distance between any two atoms across all SARS-CoV-2 spike structures available in the Protein Data Bank (excluding PDB files with only spike fragments such as the RBD), with the option to filter by various selections. The tool provides a histogram and cumulative frequency of the calculated distribution, as the ability to download the results and corresponding PDB IDs.

Force fields, such as Amber’s ff12SB, can be fairly accurate models of the physical forces in proteins and other biomolecules. When coupled with accurate solvation models, force fields are able to bring insight into the conformational preferences, transitions, pathways, and free energies for these biomolecules. When computational speed/cost matters, implicit solvent is often used but at the cost of accuracy. We present an empirical grid-like correction term, in the spirit of cMAPs, to the combination of the ff12SB protein force field and the GBneck2 implicit-solvent model. Ff12SB-cMAP is parametrized on experimental helicity data. We provide validation on a set of peptides and proteins. Ff12SB-cMAP successfully improves the secondary structure biases observed in ff12SB + Gbneck2. Ff12SB-cMAP can be downloaded (https://github.com/laufercenter/Amap.git) and used within the Amber package. It can improve the agreement of force fields + implicit solvent with experiments.

The Generalized Born (GB) implicit solvent model has undergone significant improvements in accuracy for modeling of proteins and small molecules. However, GB still remains a less widely explored option for nucleic acid simulations, in part because fast GB models are often unable to maintain stable nucleic acid structures or they introduce structural bias in proteins, leading to difficulty in application of GB models in simulations of protein–nucleic acid complexes. Recently, GB-neck2 was developed to improve the behavior of protein simulations. In an effort to create a more accurate model for nucleic acids, a similar procedure to the development of GB-neck2 is described here for nucleic acids. The resulting parameter set significantly reduces absolute and relative energy error relative to Poisson–Boltzmann for both nucleic acids and nucleic acid–protein complexes, when compared to its predecessor GB-neck model. This improvement in solvation energy calculation translates to increased structural stability for simulations of DNA and RNA duplexes, quadruplexes, and protein–nucleic acid complexes. The GB-neck2 model also enables successful folding of small DNA and RNA hairpins to near native structures as determined from comparison with experiment. The functional form and all required parameters are provided here and also implemented in the AMBER software.

Molecular mechanics is powerful for its speed in atomistic simulations, but an accurate force field is required. The Amber ff99SB force field improved protein secondary structure balance and dynamics from earlier force fields like ff99, but weaknesses in side chain rotamer and backbone secondary structure preferences have been identified. Here, we performed a complete refit of all amino acid side chain dihedral parameters, which had been carried over from ff94. The training set of conformations included multidimensional dihedral scans designed to improve transferability of the parameters. Improvement in all amino acids was obtained as compared to ff99SB. Parameters were also generated for alternate protonation states of ionizable side chains. Average errors in relative energies of pairs of conformations were under 1.0 kcal/mol as compared to QM, reduced 35% from ff99SB. We also took the opportunity to make empirical adjustments to the protein backbone dihedral parameters as compared to ff99SB. Multiple small adjustments of φ and ψ parameters were tested against NMR scalar coupling data and secondary structure content for short peptides. The best results were obtained from a physically motivated adjustment to the φ rotational profile that compensates for lack of ff99SB QM training data in the β-ppII transition region. Together, these backbone and side chain modifications (hereafter called ff14SB) not only better reproduced their benchmarks, but also improved secondary structure content in small peptides and reproduction of NMR χ1 scalar coupling measurements for proteins in solution. We also discuss the Amber ff12SB parameter set, a preliminary version of ff14SB that includes most of its improvements.