CHARMM additive and polarizable force fields for biophysics and computer-aided drug design.

BACKGROUND: Molecular Mechanics (MM) is the method of choice for computational studies of biomolecular systems owing to its modest computational cost, which makes it possible to routinely perform molecular dynamics (MD) simulations on chemical systems of biophysical and biomedical relevance. SCOPE OF REVIEW: As one of the main factors limiting the accuracy of MD results is the empirical force field used, the present paper offers a review of recent developments in the CHARMM additive force field, one of the most popular biomolecular force fields. Additionally, we present a detailed discussion of the CHARMM Drude polarizable force field, anticipating a growth in the importance and utilization of polarizable force fields in the near future. Throughout the discussion emphasis is placed on the force fields' parametrization philosophy and methodology. MAJOR CONCLUSIONS: Recent improvements in the CHARMM additive force field are mostly related to newly found weaknesses in the previous generation of additive force fields. Beyond the additive approximation is the newly available CHARMM Drude polarizable force field, which allows for MD simulations of up to 1µs on proteins, DNA, lipids and carbohydrates. GENERAL SIGNIFICANCE: Addressing the limitations ensures the reliability of the new CHARMM36 additive force field for the types of calculations that are presently coming into routine computational reach while the availability of the Drude polarizable force fields offers an inherently more accurate model of the underlying physical forces driving macromolecular structures and dynamics. This article is part of a Special Issue entitled "Recent developments of molecular dynamics".

Automation of the CHARMM General Force Field (CGenFF) II: assignment of bonded parameters and partial atomic charges.

Molecular mechanics force fields are widely used in computer-aided drug design for the study of drug candidates interacting with biological systems. In these simulations, the biological part is typically represented by a specialized biomolecular force field, while the drug is represented by a matching general (organic) force field. In order to apply these general force fields to an arbitrary drug-like molecule, functionality for assignment of atom types, parameters, and partial atomic charges is required. In the present article, algorithms for the assignment of parameters and charges for the CHARMM General Force Field (CGenFF) are presented. These algorithms rely on the existing parameters and charges that were determined as part of the parametrization of the force field. Bonded parameters are assigned based on the similarity between the atom types that define said parameters, while charges are determined using an extended bond-charge increment scheme. Charge increments were optimized to reproduce the charges on model compounds that were part of the parametrization of the force field. A "penalty score" is returned for every bonded parameter and charge, allowing the user to quickly and conveniently assess the quality of the force field representation of different parts of the compound of interest. Case studies are presented to clarify the functioning of the algorithms and the significance of their output data.

Automation of the CHARMM General Force Field (CGenFF) I: bond perception and atom typing.

Molecular mechanics force fields are widely used in computer-aided drug design for the study of drug-like molecules alone or interacting with biological systems. In simulations involving biological macromolecules, the biological part is typically represented by a specialized biomolecular force field, while the drug is represented by a matching general (organic) force field. In order to apply these general force fields to an arbitrary drug-like molecule, functionality for assignment of atom types, parameters, and charges is required. In the present article, which is part I of a series of two, we present the algorithms for bond perception and atom typing for the CHARMM General Force Field (CGenFF). The CGenFF atom typer first associates attributes to the atoms and bonds in a molecule, such as valence, bond order, and ring membership among others. Of note are a number of features that are specifically required for CGenFF. This information is then used by the atom typing routine to assign CGenFF atom types based on a programmable decision tree. This allows for straightforward implementation of CGenFF's complicated atom typing rules and for equally straightforward updating of the atom typing scheme as the force field grows. The presented atom typer was validated by assigning correct atom types on 477 model compounds including in the training set as well as 126 test-set molecules that were constructed to specifically verify its different components. The program may be utilized via an online implementation at https://www.paramchem.org/ .

Development of the CHARMM Force Field for Lipids.

The development of the CHARMM additive all-atom lipid force field (FF) is traced from the early 1990's to the most recent version (C36) published in 2010. Though simulations with early versions yielded useful results, they failed to reproduce two important quantities: a zero surface tension at the experimental bilayer surface area, and the signature splitting of the deuterium order parameters in the glycerol and upper chain carbons. Systematic optimization of parameters based on high level quantum mechanical data and free energy simulations have resolved these issues, and bilayers with a wide range of lipids can be simulated in tensionless ensembles using C36. Issues associated with other all-atom lipid FFs, success and limitations in the C36 FF and ongoing developments are also discussed.

Targeting of an interrupted polypurine:polypyrimidine sequence in mammalian cells by a triplex-forming oligonucleotide containing a novel base analogue.

The DNA triple helix consists of a third strand of nucleic acid lying in the major groove of an intact DNA duplex. The most stable triplexes form on polypurine:polypyrimidine sequences, and pyrimidine interruptions in the purine strand are destabilizing. Sequence stringency is imparted by specific Hoogsteen hydrogen bonds between third strand bases and the purine bases in the duplex. Appropriate base and sugar modifications of triple helix-forming oligonucleotides (TFOs) confer chromosome targeting activity in living cells. However, broad utilization of TFOs as gene targeting reagents in mammalian cells has been limited by the requirement for homopurine target sequences. Although there have been a number of base analogues described that appear to be promising as candidates for triplex target expansion, none has been examined in a biological system. We have employed a postsynthetic strategy to prepare a collection of TFOs with base analogues at a defined position. Following assessment of affinity for a triplex target with a single C:G inversion, TFOs with a second generation of analogues were synthesized. One of these, TFO-5a, with 2'-OMe-guanidinylethyl-5-methylcytosine at the position corresponding to the C:G interruption in the target sequence, was further modified to confer bioactivity. The activity of this TFO, linked to psoralen, was measured in a mammalian cell line that was engineered by directed sequence conversion to carry a triplex target with a single C:G interruption. TFO-5a was active against this target and inactive against the corresponding target with an uninterrupted polypurine:polypyrimidine sequence.

CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields.

The widely used CHARMM additive all-atom force field includes parameters for proteins, nucleic acids, lipids, and carbohydrates. In the present article, an extension of the CHARMM force field to drug-like molecules is presented. The resulting CHARMM General Force Field (CGenFF) covers a wide range of chemical groups present in biomolecules and drug-like molecules, including a large number of heterocyclic scaffolds. The parametrization philosophy behind the force field focuses on quality at the expense of transferability, with the implementation concentrating on an extensible force field. Statistics related to the quality of the parametrization with a focus on experimental validation are presented. Additionally, the parametrization procedure, described fully in the present article in the context of the model systems, pyrrolidine, and 3-phenoxymethylpyrrolidine will allow users to readily extend the force field to chemical groups that are not explicitly covered in the force field as well as add functional groups to and link together molecules already available in the force field. CGenFF thus makes it possible to perform "all-CHARMM" simulations on drug-target interactions thereby extending the utility of CHARMM force fields to medicinally relevant systems.

CHARMM: the biomolecular simulation program.

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

Inhibition of (cytosine C5)-methyltransferase by oligonucleotides containing flexible (cyclopentane) and conformationally constrained (bicyclo[3.1.0]hexane) abasic sites.

Pseudorotationally locked sugar analogues based on bicyclo[3.1.0]-hexane templates were placed in DNA duplexes as abasic target sites in the M. HhaI recognition sequence. The binding affinity of the enzyme increases when the abasic site is constrained to the South conformation and decreases when it is constrained to the North conformation. A structural understanding of these differences is provided.

Reevaluation of stereoelectronic contributions to the conformational properties of the phosphodiester and N3'-phosphoramidate moieties of nucleic acids.

The anomeric effect in the phosphodiester backbone of nucleic acids is a stereoelectronic effect that has conventionally been linked to interactions between lone pairs on the O(ester) atoms and P-O(ester) antibonding orbitals. The present study demonstrates that the anomeric effect in the phosphodiester backbone is significantly more complex than portrayed by this description. The presence of multiple lone pairs and antibonding orbitals around the phosphorus atom leads to additional contributions to the anomeric effect, especially involving the anionic oxygen lone pairs. On the basis of the structural changes and Natural Bond Orbital analysis it is shown that a complex balance between stereoelectronic effects involving both the ester and anionic oxygen lone pairs governs the conformational properties of the phosphodiester backbone. The N3'-phosphoramidate DNA backbone differs from the phosphodiester backbone due to the N3'-H moiety having only a single lone pair instead of the two lone pairs present on the O3' atom substituted. The present study uses N3'-phosphoramidate as a control to understand the changes in stereoelectronic effects as a result of changes in the structure and conformation. Two previously uncharacterized properties of the N3'-phosphoramidate backbone are also observed and explained through the complex balance of the postulated electronic delocalizations. The first observation is that the N3'-H moiety in N3'-phosphoramidate is a flexible moiety that can change the orientation of its hydrogen through inversion without a significant energetic penalty in both the gas phase and the aqueous phase. The second observation is that the stabilization of the C3'-endo conformation in N3'-phosphoramidate is primarily due to aqueous solvation rather than intrinsic gas-phase effects involving the reduced electronegativity of the 3'-substituent.

Re-examination of the intrinsic, dynamic and hydration properties of phosphoramidate DNA.

Intrinsic energetic and solvation factors contributing to the unusual structural and biochemical properties of N3'-phosphoramidate DNA analogs have been re-examined using a combination of quantum mechanical and molecular dynamics methods. Evaluation of the impact of the N3'-H substitution was performed via comparison of N3'-phosphoramidate DNA starting from both A- and B-form structures, B-form DNA and A-form RNA. The N3'-H group is shown to be flexible, undergoing reversible inversion transitions associated with motion of the hydrogen atom attached to the N3' atom. The inversion process is correlated with both sugar pucker characteristics as well as other local backbone torsional dynamics, yielding increased dihedral flexibility over DNA. Solvation of N3'-phosphoramidate DNA is shown to be similar to RNA, consistent with thermodynamic data on the two species. A previously unobserved intrinsic conformational perturbation caused by the N5'-phosphoramidate substitution is identified and suggested to be linked to the differences in the properties of N3'- and N5'-phosphoramidate oligonucleotide analogs.

Ab initio conformational analysis of nucleic acid components: intrinsic energetic contributions to nucleic acid structure and dynamics.

In recent years, the use of high-level ab initio calculations has allowed for the intrinsic conformational properties of nucleic acid building blocks to be revisited. This has provided new insights into the intrinsic conformational energetics of these compounds and its relationship to nucleic acids structure and dynamics. In this article we review recent developments and present new results. New data include comparison of various levels of theory on conformational properties of nucleic acid building blocks, calculations on the abasic sugar, known to occur in vivo in DNA, on the TA conformation of DNA observed in the complex with the TATA box binding protein, and on inosine. Tests of the Hartree-Fock (HF), second-order Møller-Plesset (MP2), and Density Functional Theory/Becke3, Lee, Yang and Par (DFT/B3LYP) levels of theory show the overall shape of backbone torsional energy profiles (for gamma, epsilon, and chi) to be similar for the different levels, though some systematic differences are identified between the MP2 and DFT/B3LYP profiles. The east pseudorotation energy barrier in deoxyribonucleosides is also sensitive to the level of theory, with the HF and DFT/B3LYP east barriers being significantly lower (approximately 2.5 kcal/mol) than the MP2 counterpart (approximately 4.0 kcal/mol). Additional calculations at various levels of theory suggest that the east barrier in deoxyribonucleosides is between 3.0 and 4.0 kcal/mol. In the abasic sugar, the west pseudorotation energy barrier is found to be slightly lower than the east barrier and the south pucker is favored more than in standard nucleosides. Results on the TA conformation suggest that, at the nucleoside level, this conformation is significantly destabilized relative to the global energy minimum, or relative to the A- and B-DNA conformations. Deoxyribocytosine would destabilize the TA conformation more than other bases relative to the A-DNA conformation, but not relative to the B-DNA conformation.

Identification of HIV-1 integrase inhibitors via three-dimensional database searching using ASV and HIV-1 integrases as targets.

Integration of viral DNA into the host cell genome is a critical step in the life cycle of HIV. This essential reaction is catalyzed by integrase (IN) through two steps, 3'-processing and DNA strand transfer. Integrase is an attractive target for drug design because there is no known cellular analogue and integration is essential for successful replication of HIV. A computational three-dimensional (3-D) database search was used to identify novel HIV-1 integrase inhibitors. Starting from the previously identified Y3 (4-acetylamino-5-hydroxynaphthalene-2,7-disulfonic acid) binding site on the avian sarcoma virus integrase (ASV IN), a preliminary search of all compounds in the nonproprietary, open part of the National Cancer Institute 3-D database yielded a collection of 3100 compounds. A more rigorous scoring method was used to rescreen the 3100 compounds against both ASV IN and HIV-1 IN. Twenty-two of those compounds were selected for inhibition assays against HIV-1 IN. Thirteen of the 22 showed inhibitory activity against HIV-1 IN at concentrations less than 200 microM and three of them showed antiviral activities in HIV-1 infected CEM cells with effective concentrations (EC50) ranging from 0.8 to 200 microM. Analysis of the computer-generated binding modes of the active compounds to HIV-1 IN showed that simultaneous interaction with the Y3 site and the catalytic site is possible. In addition, interactions between the active compounds and the flexible loop involved in the binding of DNA by IN are indicated to occur. The structural details and the unique binding motif between the HIV-1 IN and its inhibitors identified in the present work may contribute to the future development of IN inhibitors.

New insights into the structure of abasic DNA from molecular dynamics simulations.

Abasic (AP) sites constitute a common form of DNA damage, arising from the spontaneous or enzymatic breakage of the N-glycosyl bond and the loss of a nucleotide base. To examine the effects of such damage on DNA structure, especially in the vicinity of the abasic sugar, four 1.5 ns molecular dynamics simulations of double-helical DNA dodecamers with and without a single abasic (tetrahydrofuran, X) lesion in a 5'-d(CXT) context have been performed and analyzed. The results indicate that the abasic site does not maintain a hole or gap in the DNA, but instead perturbs the canonical structure and induces additional flexibility close to the abasic site. In the apurinic simulations (i.e., when a pyrimidine is opposite the AP site), the abasic sugar flipped in and out of the minor groove, and the gap was water filled, except during the occurrence of a novel non-Watson-Crick C-T base pair across the abasic site. The apyrimidinic gap was not penetrated by water until the abasic sugar flipped out and remained extrahelical. Both AP helices showed kinks of 20-30 degrees at the abasic site. The Watson-Crick hydrogen bonds are more transient throughout the DNA double helices containing an abasic site. The abasic sugar displayed an unusually broad range of sugar puckers centered around the northern pucker. The increased motion of the bases and backbone near the abasic site appear to correlate with sequence-dependent helical stability. The data indicate that abasic DNA contorts more easily and in specific ways relative to unmodified DNA, an aspect likely to be important in abasic site recognition and hydrolysis.

Development and current status of the CHARMM force field for nucleic acids.

The CHARMM27 all-atom force field for nucleic acids represents a highly optimized model for investigations of nucleic acids via empirical force field calculations. The force field satisfactorily treats the A, B, and Z forms of DNA as well as RNA, and it also useful for nucleosides and nucleotides. In addition, it is compatible with the CHARMM force fields for proteins and lipids, allowing for simulation studies of heterogeneous systems.

Structure, force, and energy of a double-stranded DNA oligonucleotide under tensile loads.

The end-to-end stretching of a duplex DNA oligonucleotide has been studied using potential of mean force (PMF) calculations based on molecular dynamics (MD) simulations and atomic force microscopy (AFM) experiments. Near quantitative agreement between the calculations and experiments was obtained for both the extension length and forces associated with strand separation. The PMF calculations show that the oligonucleotide extends without a significant energetic barrier from a length shorter than A-DNA to a length 2.4 times the contour length of B-DNA at which the barrier to strand separation is encountered. Calculated forces associated with the barrier are 0.09 +/- 0.03 nN, based on assumptions concerning tip and thermal-activated barrier crossing contributions to the forces. Direct AFM measurements show the oligonucleotide strands separating at 2.6 +/- 0.8 contour lengths with a force of 0.13 +/- 0.05 nN. Analysis of the energies from the MD simulations during extension reveals compensation between increases in the DNA-self energy and decreases in the DNA-solvent interaction energy, allowing for the barrierless extension of DNA beyond the canonical B form. The barrier to strand separation occurs when unfavorable DNA interstrand repulsion cannot be compensated for by favorable DNA-solvent interactions. The present combination of single molecule theoretical and experimental approaches produces a comprehensive picture of the free energy surface of biological macromolecular structural transitions.

Intrinsic conformational properties of deoxyribonucleosides: implicated role for cytosine in the equilibrium among the A, B, and Z forms of DNA.

Structural properties of biomolecules are dictated by their intrinsic conformational energetics in combination with environmental contributions. Calculations using high-level ab initio methods on the deoxyribonucleosides have been performed to investigate the influence of base on the intrinsic conformational energetics of nucleosides. Energy minima in the north and south ranges of the deoxyribose pseudorotation surfaces have been located, allowing characterization of the influence of base on the structures and energy differences between those minima. With all bases, chi values associated with the south energy minimum are lower than in canonical B-DNA, while chi values associated with the north energy minimum are close to those in canonical A-DNA. In deoxycytidine, chi adopts an A-DNA conformation in both the north and south energy minima. Energy differences between the A and B conformations of the nucleosides are <0.5 kcal/mol in the present calculations, except with deoxycytidine, where the A form is favored by 2.3 kcal/mol, leading the intrinsic conformational energetics of GC basepairs to favor the A form of DNA by 1.5 kcal/mol as compared with AT pairs. This indicates that the intrinsic conformational properties of cytosine at the nucleoside level contribute to the A form of DNA containing predominately GC-rich sequences. In the context of a B versus Z DNA equilibrium, deoxycytidine favors the Z form over the B form by 1.6 kcal/mol as compared with deoxythymidine, suggesting that the intrinsic conformational properties of cytosine also contribute to GC-rich sequences occurring in Z DNA with a higher frequency than AT-rich sequences. Results show that the east pseudorotation energy barrier involves a decrease in the furanose amplitude and is systematically lower than the inversion barrier, with the energy differences influenced by the base. Energy barriers going from the south (B form) sugar pucker to the east pseudorotation barrier are lower in pyrimidines as compared with purines, indicating that the intrinsic conformational properties associated with base may also influence the sugar pseudorotational population distribution seen in DNA crystal structures and the kinetics of B to A transitions. The present work provides evidence that base composition, in addition to base sequence, can influence DNA conformation.

TIT for TAT: the properties of inosine and adenosine in TATA box DNA.

The sequence dependent conformation, flexibility and hydration properties of DNA molecules constitute selectivity determinants in the formation of protein-DNA complexes. TATA boxes in which AT basepairs (bp) have been substituted by IC bp (TITI box) allow for probing these selectivity determinants for the complexation with the TATA box-binding protein (TBP) with different sequences but identical chemical surfaces. The reference promoter Adenovirus 2 Major Late Promoter (mlp) is formed by the apposition of two sequences with very different dynamic properties: an alternating TATA sequence and an A-tract. For a comparative study, we carried out molecular dynamics simulations of two DNA oligomers, one containing the mlp sequence (2 ns), and the other an analog where AT basepairs were substituted by IC basepairs (1 ns). The simulations, carried out with explicit solvent and counterinons, yield straight purine tracts, the A-tract being stiffer than the I-tract, an alternating structure for the YRYR tracts, and hydration patterns that differ between the purine tracts and the alternating sequence tracts. A detailed analysis of the proposed interactions responsible for the stiffness of the purine tracts indicates that the stacking between the bases bears the strongest correlation to stiffness. The hydration properties of the minor groove in the two oligomers are distinctly different. Such differences are likely to be responsible for the stronger binding of TBP to mlp over the inosine-substituted variant. The calculations were made possible by the development, described here, of a new set of forcefield parameters for inosine that complement the published CHARMM all-hydrogen nucleic acid parametrization.

(-)-3 beta-Substituted ecgonine methyl esters as inhibitors for cocaine binding and dopamine uptake.

Ten 3 beta-ecgonine analogues were synthesized and characterized by 1H and 13C NMR, MS, and elemental analysis. The compounds were synthesized as (-)-stereoisomers from (-)-cocaine. These compounds were assessed for their ability to inhibit [3H]cocaine binding to rat striatal tissue and to inhibit [3H]DA uptake into rat striatal synaptosomes. In this series of compounds, the length of the spacer between the aryl group and the tropane skeleton ranged from 1 to 4 bond distances, and conformational flexibility of the linkage and orientation of the aryl ring system were controlled by various types of linkages. The most potent of the analogues was methyl-(1R-2-exo-3-exo)-8-methyl-3-(beta-styrenyl)-8-azabicyclo[3. 2.1] octane-2-carboxylate. One of the less potent compounds was found to inhibit [3H]cocaine binding and [3H]DA uptake with significantly different IC50 values, in contrast to 14 other 3 beta-substituted analogues. Molecular modeling and CoMFA analysis were used to obtain a rigorous structure-function relationship for the studied compounds. The results showed that the potencies of these 3 beta-substituted ecgonine methyl esters were dominated by steric effects and were acutely sensitive to the distance between the aryl ring and the tropane skeleton and to the orientation of the aryl ring system relative to the tropane skeleton. The current study provides a clearer picture of the shape and size of the putative hydrophobic binding pocket for the 3 beta substituent at the cocaine receptor as well as emphasizing the importance of a drug's free energy of solvation in obtaining structure-activity relationships.

All-atom empirical potential for molecular modeling and dynamics studies of proteins.

New protein parameters are reported for the all-atom empirical energy function in the CHARMM program. The parameter evaluation was based on a self-consistent approach designed to achieve a balance between the internal (bonding) and interaction (nonbonding) terms of the force field and among the solvent-solvent, solvent-solute, and solute-solute interactions. Optimization of the internal parameters used experimental gas-phase geometries, vibrational spectra, and torsional energy surfaces supplemented with ab initio results. The peptide backbone bonding parameters were optimized with respect to data for N-methylacetamide and the alanine dipeptide. The interaction parameters, particularly the atomic charges, were determined by fitting ab initio interaction energies and geometries of complexes between water and model compounds that represented the backbone and the various side chains. In addition, dipole moments, experimental heats and free energies of vaporization, solvation and sublimation, molecular volumes, and crystal pressures and structures were used in the optimization. The resulting protein parameters were tested by applying them to noncyclic tripeptide crystals, cyclic peptide crystals, and the proteins crambin, bovine pancreatic trypsin inhibitor, and carbonmonoxy myoglobin in vacuo and in crystals. A detailed analysis of the relationship between the alanine dipeptide potential energy surface and calculated protein φ, χ angles was made and used in optimizing the peptide group torsional parameters. The results demonstrate that use of ab initio structural and energetic data by themselves are not sufficient to obtain an adequate backbone representation for peptides and proteins in solution and in crystals. Extensive comparisons between molecular dynamics simulations and experimental data for polypeptides and proteins were performed for both structural and dynamic properties. Energy minimization and dynamics simulations for crystals demonstrate that the latter are needed to obtain meaningful comparisons with experimental crystal structures. The presented parameters, in combination with the previously published CHARMM all-atom parameters for nucleic acids and lipids, provide a consistent set for condensed-phase simulations of a wide variety of molecules of biological interest.