The role of molecular structure of sugar-phosphate backbone and nucleic acid bases in the formation of single-stranded and double-stranded DNA structures.

Our previous DFT computations of deoxydinucleoside monophosphate complexes with Na(+)-ions (dDMPs) have demonstrated that the main characteristics of Watson-Crick (WC) right-handed duplex families are predefined in the local energy minima of dDMPs. In this work, we study the mechanisms of contribution of chemically monotonous sugar-phosphate backbone and the bases into the double helix irregularity. Geometry optimization of sugar-phosphate backbone produces energy minima matching the WC DNA conformations. Studying the conformational variability of dDMPs in response to sequence permutation, we found that simple replacement of bases in the previously fully optimized dDMPs, e.g. by constructing Pyr-Pur from Pur-Pyr, and Pur-Pyr from Pyr-Pur sequences, while retaining the backbone geometry, automatically produces the mutual base position characteristic of the target sequence. Based on that, we infer that the directionality and the preferable regions of the sugar-phosphate torsions, combined with the difference of purines from pyrimidines in ring shape, determines the sequence dependence of the structure of WC DNA. No such sequence dependence exists in dDMPs corresponding to other DNA conformations (e.g., Z-family and Hoogsteen duplexes). Unlike other duplexes, WC helix is unique by its ability to match the local energy minima of the free single strand to the preferable conformations of the duplex.

Quantum mechanical binding free energy calculation for phosphopeptide inhibitors of the Lck SH2 domain.

The accurate and efficient calculation of binding free energies is essential in computational biophysics. We present a linear-scaling quantum mechanical (QM)-based end-point method termed MM/QM-COSMO to calculate binding free energies in biomolecular systems, with an improved description of entropic changes. Molecular dynamics trajectories are re-evaluated using a semiempirical Hamiltonian and a continuum solvent model; translational and rotational entropies are calculated using configurational integrals, and internal entropy is calculated using the harmonic oscillator approximation. As an application, we studied the binding of a series of phosphotyrosine tetrapeptides to the human Lck SH2 domain, a key component in intracellular signal transduction, modulation of which can have therapeutic relevance in the treatment of cancer, osteoporosis, and autoimmune diseases. Calculations with molecular mechanics Poisson-Boltzmann, and generalized Born surface area methods showed large discrepancies with experimental data stemming from the enthalpic component, in agreement with an earlier report. The empirical force field-based solvent interaction energy scoring function yielded improved results, with average unsigned error of 3.6 kcal/mol, and a better ligand ranking. The MM/QM-COSMO method exhibited the best agreement both for absolute (average unsigned error = 0.7 kcal/mol) and relative binding free energy calculations. These results show the feasibility and promise of a full QM-based end-point method with an adequate balance of accuracy and computational efficiency.

Computational and experimental studies of the interaction between phospho-peptides and the C-terminal domain of BRCA1.

The C-terminal domain of BRCA1(BRCT) is involved in the DNA repair pathway by recognizing the pSXXF motif in interacting proteins. It has been reported that short peptides containing this motif bind to BRCA1(BRCT) in the micromolar range with high specificity. In this work, the binding of pSXXF peptides has been studied computationally and experimentally in order to characterize their interaction with BRCA1(BRCT). Elucidation of the contacts that drive the protein-ligand interaction is critical for the development of high affinity small-molecule BRCA1 inhibitors. Molecular dynamics (MD) simulations revealed the key role of threonine at the peptide P+2 position in providing structural rigidity to the ligand in the bound state. The mutation at P+1 had minor effects. Peptide extension at the N-terminal position with the naphthyl amino acid exhibited a modest increase in binding affinity, what could be explained by the dispersion interaction of the naphthyl side-chain with a hydrophobic patch. Three in silico end-point methods were considered for the calculation of binding free energy. The Molecular Mechanics Poisson-Boltzmann Surface Area and the Solvated Interaction Energy methods gave reasonable agreement with experimental data, exhibiting a Pearlman predictive index of 0.71 and 0.78, respectively. The MM-quantum mechanics-surface area method yielded improved results, which was characterized by a Pearlman index of 0.78. The correlation coefficients were 0.59, 0.61 and 0.69, respectively. The ability to apply a QM level of theory within an end-point binding free energy protocol may provide a way for a consistent improvement of accuracy in computer-aided drug design.

QM/QM docking method based on the variational finite localized molecular orbital approximation.

We present a derivation of the semiempirical variational finite localized molecular orbital (VFL) approximation, which was introduced by Anikin et al. (J Chem Phys 2004, 121, 1266). On the basis of VFL approximation, we developed the novel semiempirical (quantum mechanical) QM/QM method in which a part of the system, including the ligand and protein active site, are treated self-consistently, while the protein bulk is considered as carrying a frozen electronic density matrix. The developed method is applied toward the QM docking study for the p56 LCK SH2 domain. The virtual search has predicted 10 most potent inhibitors by searching through the database of 200,000 empirically docked poses of 20,000 drug-like molecules. Energy score calculation of each complex roughly consisting of 1700 atoms took 14.54 s of single-CPU time at the NDDO AM1 level. The entire computation performed on a 32-CPU cluster would be accomplished in 1 day. Flexible ligand QM docking studies, performed on a subset of 10,000 poses, required 153.03 s of single-CPU time per complex. The entire calculation performed on the 32-CPU cluster would be finished in half-day.

Understanding the dielectric properties of liquid amides from a polarizable force field.

The role played by electronic polarization in the dielectric properties of liquid N-methyl acetamide (NMA) is examined using molecular dynamics simulations with a polarizable force field based on classical Drude oscillators. The model presented is the first force field shown to reproduce the anomalously large dielectric constant of liquid NMA. Details of the molecular polarizability are found to be important. For instance, all elements of the polarizability tensor, rather then just the trace, impact on the condensed phase properties. Two factors related to electronic polarizability are found to contribute to this large dielectric constant. First is the significant enhancement of the mean amide molecular dipole magnitude, which is 50% larger in the liquid than in the gas phase. Second is the consequent strong hydrogen bonding between molecular neighbors that enhances the orientational alignment of the molecular dipoles. Polarizable models of amide compounds that have two (acetamide) and zero (N,N-dimethyl acetamide) polar hydrogen-bond donor atoms are also investigated. Experimentally, the neat liquid dielectric constants at 373 K are 100 for NMA, 66 for acetamide and 26 for N,N-dimethyl acetamide. The polarizable models replicate this trend, predicting a dielectric constant of 92+/-5 for NMA, 66+/-3 for acetamide and 23+/-1 for N,N-dimethyl acetamide.

DFT study of B-like conformations of deoxydinucleoside monophosphates containing Gua and/or Cyt and their complexes with Na+ cation.

B-like minimum energy conformations of deoxydinucleoside monophosphate anions (dDMPs) containing Gua and/or Cyt and their Na+ complexes have been studied by the DFT PW91PW91/DZVP method. The optimized geometry of the dDMPs is in close agreement with experimental observations and the obtained minimum energy conformations are consistent with purine-purine, purine-pyrimidine, and pyrimidine-purine arrangements in crystals of B-DNA duplexes. All the studied systems are characterized by pyramidalization of the amino groups, which participate in the formation of unusual hydrogen bond between the carbonyl oxygen of the second base in the dGpdC, dCpdG dDMPs, and their Na+ complexes. In all the obtained structures the bases assume a nearly parallel disposition to each other and this effect is independent on the degree of their spatial superposition. From this it is concluded that the parallel disposition of the bases in the B-like single-stranded conformations is dictated by the sugar-phosphate backbone. Correspondingly, the base-base interactions attain a secondary role in the formation of these spatial structures. The formation of a weak C6-H6...O5' hydrogen bond between cytosine and the phosphate oxygen is reported, in agreement with experimental observations.

Cooperative binding of DNA and CBFbeta to the Runt domain of the CBFalpha studied via MD simulations.

The Runt domain (RD) is the DNA-binding region of the Runx genes. A related protein, known as core binding factor beta (CBFbeta) also binds to the RD to enhance RD-DNA interaction by 6- to 10-fold. Here, we report results from molecular dynamics (MD) simulations of RD alone, as a dimer in complexes with DNA and CBFbeta and in a ternary complex with DNA and CBFbeta. Consistent with the experimental findings, in the presence of CBFbeta the estimated free energy of binding of RD to the DNA is more favorable, which is shown to be due to more favorable intermolecular interactions and desolvation contributions. Also contributing to the enhanced binding are favorable intramolecular interactions between the 'wing' residues (RD residues 139-145) and the 'wing1' residues (RD residues 104-116). The simulation studies also indicate that the RD-CBFbeta binding is more favorable in the presence of DNA due to a more favorable RD-CBFbeta interaction energy. In addition, it is predicted that long-range interactions involving ionic residues contribute to binding cooperativity. Results from the MD calculations are used to interpret a variety of experimental mutagenesis data. A novel role for RD Glu116 to the RD-CBFbeta interaction is predicted.

Polarizable empirical force field for alkanes based on the classical Drude oscillator model.

Recent extensions of potential energy functions used in empirical force field calculations have involved the inclusion of electronic polarizability. To properly include this extension into a potential energy function it is necessary to systematically and rigorously optimize the associated parameters based on model compounds for which extensive experimental data are available. In the present work, optimization of parameters for alkanes in a polarizable empirical force field based on a classical Drude oscillator is presented. Emphasis is placed on the development of parameters for CH3, CH2, and CH moieties that are directly transferable to long chain alkanes, as required for lipids and other biomolecules. It is shown that a variety of quantum mechanical and experimental target data are reproduced by the polarizable model. Notable is the proper treatment of the dielectric constant of pure alkanes by the polarizable force field, a property essential for the accurate treatment of, for example, hydrophobic solvation in lipid bilayers. The present alkane force field will act as the basis for the aliphatic moieties in an extensive empirical force field for biomolecules that includes the explicit treatment of electronic polarizability.

Post Hartree-Fock studies of the canonical Watson-Crick DNA base pairs: molecular structure and the nature of stability.

Gas-phase gradient optimization was carried out on the canonical Watson-Crick DNA base pairs using the second-order Møller-Plesset perturbation method at the 6-31G(d) and 6-31G(d,p) basis sets. It is detected that full geometry optimization at the MP2 level leads to an intrinsically nonplanar propeller-twisted and buckled geometry of G-C and A-T base pairs; while HF and DFT methods predict perfect planar or almost planar geometry of the base pairs. Supposedly the nonplanarity of the pairs is caused by pyramidalization of the amino nitrogen atoms, which is underestimated by the HF and DFT methods. This justifies the importance of geometry optimization at the MP2 level for obtaining reliable prediction of the charge distribution, molecular dipole moments and geometrical structure of the base pairs. The Morokuma-Kitaura and the Reduced Variational Space methods of the decomposition for molecular HF interaction energies were used for investigation of the hydrogen bonding in the Watson-Crick base pairs. It is shown that the HF stability of the hydrogen-bonded DNA base pairs originates mainly from electrostatic interactions. At the same time, the calculated magnitude of the second order intramolecular correlation correction to the Coulomb energy showed that electron correlation reduces the contribution of the electrostatic term to the attractive interaction for the A-T and G-C base pairs. Polarization, charge transfer and dispersion interactions also make considerable contribution to the attraction energy of bases.

Optimization of the Coupled Cluster Implementation in NWChem on Petascale Parallel Architectures.

The coupled cluster singles and doubles (CCSD) algorithm in the NWChem software package has been optimized to alleviate the communication bottleneck. This optimization provided a 2-fold to 5-fold speedup in the CCSD iteration time depending on the problem size and available memory, and improved the CCSD scaling to 20 000 nodes of the NCSA Blue Waters supercomputer. On 20 000 XE6 nodes of Blue Waters, a complete conventional CCSD(T) calculation of a system encountering 1042 basis functions and 103 occupied correlated orbitals obtained a performance of 0.32 petaflop/s and took 5 h and 24 min to complete. The reported time and performance included all stages of the calculation from initialization to termination for iterative single and double excitations as well as perturbative triples correction. In perturbative triples alone, the computation sustained a rate of 1.18 petaflop/s. The CCSD and (T) phases took approximately (3)/4 and (1)/4 of the total time to solution, respectively, showing that CCSD is the most time-consuming part at the large scale. The MP2, CCSD, and CCSD(T) computations in 6-311++G** basis set performed on guanine-cytosine deoxydinucleotide monophosphate probed the conformational energy difference between the A- and B-conformations of single stranded DNA. Good agreement between MP2 and coupled cluster methods has been obtained, suggesting the utility of MP2 for conformational analysis in these systems. The study revealed a significant discrepancy between the quantum mechanical and classical force field predictions, suggesting a need to improve the dihedral parameters.

Charge transfer effects in the GroEL-GroES chaperonin tetramer in solution.

In this work, we present the results of a large-scale, semiempirical LocalSCF quantum mechanical study of GroEL-GroES chaperonin in solution containing 2,481,723 atoms. We find that large biological systems exhibit strong quantum mechanical character, the extent of which was not previously known. Our data show that protein transfers -743 electron units of charge to solvent, which is not described by classical force fields. Contrary to the commonly held belief, which is based on classical mechanics, our computational data suggest that the quantum mechanical effects of charge transfer increase with the size of biological systems. We show that the neglect of charge transfer in classical force fields leads to significant error in the electrostatic potential of the macromolecule. These findings illustrate that a quantum mechanical framework is necessary for a realistic description of electrostatic interactions in large biological systems.

Hydration free energies using semiempirical quantum mechanical Hamiltonians and a continuum solvent model with multiple atomic-type parameters.

To build the foundation for accurate quantum mechanical (QM) simulation of biomacromolecules in an aqueous environment, we undertook the optimization of the COnductor-like Screening MOdel (COSMO) atomic radii and atomic surface tension coefficients for different semiempirical Hamiltonians adhering to the same computational conditions recently followed in the simulation of biomolecular systems. This optimization was achieved by reproducing experimental hydration free energies of a set consisting of 507 neutral and 99 ionic molecules. The calculated hydration free energies were significantly improved by introducing a multiple atomic-type scheme that reflects different chemical environments. The nonpolar contribution was treated according to the scaled particle Claverie-Pierotti formalism. Separate radii and surface tension coefficient sets have been developed for AM1, PM3, PM5, and RM1 semiempirical Hamiltonians, with an average unsigned error for neutral molecules of 0.64, 0.66, 0.73, and 0.71 kcal/mol, respectively. Free energy calculation of each molecule took on average 0.5 s on a single processor. The new sets of parameters will enhance the quality of semiempirical QM calculations using COSMO in biomolecular systems. Overall, these results further extend the utility of QM methods to chemical and biological systems in the condensed phase.

Development of CHARMM polarizable force field for nucleic acid bases based on the classical Drude oscillator model.

A polarizable force field for nucleic acid bases based on the classical Drude oscillator model is presented. Parameter optimization was performed to reproduce crystallographic geometries, crystal unit cell parameters, heats of sublimation, vibrational frequencies and assignments, dipole moments, molecular polarizabilities and quantum mechanical base-base and base-water interaction energies. The training and validation data included crystals of unsubstituted and alkyl-substituted adenine, guanine, cytosine, uracil, and thymine bases, hydrated crystals, and hydrogen bonded base pairs. Across all compounds, the RMSD in the calculated heats of sublimation is 4.1%. This equates to an improvement of more than 2.5 kcal/mol in accuracy compared to the nonpolarizable CHARMM27 force field. However, the level of agreement with experimental molecular volume decreased from 1.7% to 2.1% upon moving from the nonpolarizable to the polarizable model. The representation of dipole moments is significantly improved with the Drude polarizable force field. Unlike in additive force fields, there is no requirement for the gas-phase dipole moments to be overestimated, illustrating the ability of the Drude polarizable force field to treat accurately differently dielectric environments and indicating the improvements in the electrostatic model. Validation of the model was performed on the basis of the calculation of the gas-phase binding enthalpies of base pairs obtained via potential of mean force calculations; the additive and polarizable models both performed satisfactorily with average differences of 0.2 and 0.9 kcal/mol, respectively, and rms differences of 1.3 and 1.7 kcal/mol, respectively. Overall, considering the number of significant improvements versus the additive CHARMM force field, the incorporation of explicit polarizability into the force field for nucleic acid bases represents an additional step toward accurate computational modeling of biological systems.

Simulating Monovalent and Divalent Ions in Aqueous Solution Using a Drude Polarizable Force Field.

An accurate representation of ion solvation in aqueous solution is critical for meaningful computer simulations of a broad range of physical and biological processes. Polarizable models based on classical Drude oscillators are introduced and parametrized for a large set of monoatomic ions including cations of the alkali metals (Li(+), Na(+), K(+), Rb(+) and Cs(+)) and alkaline earth elements (Mg(2+), Ca(2+), Sr(2+) and Ba(2+)) along with Zn(2+) and halide anions (F(-), Cl(-), Br(-) and I(-)). The models are parameterized, in conjunction with the polarizable SWM4-NDP water model [Lamoureux et al., Chem. Phys. Lett. 418, 245 (2006)], to be consistent with a wide assortment of experimentally measured aqueous bulk thermodynamic properties and the energetics of small ion-water clusters. Structural and dynamic properties of the resulting ion models in aqueous solutions at infinite dilution are presented.

Re-evaluation of the reported experimental values of the heat of vaporization of N-methylacetamide.

The accuracy of empirical force fields is inherently related to the quality of the target data used for optimization of the model. With the heat of vaporization (ΔH(vap)) of N-methylacetamide (NMA), a range of values have been reported as target data for optimization of the nonbond parameters associated with the peptide bond in proteins. In the present work, the original experimental data and Antoine constants used for the determination of the ΔH(vap) of NMA are reanalyzed. Based on this analysis, the wide range of ΔH(vap) values reported in the literature are shown to be due to incorrect reporting of the temperatures at which the original values were extracted and limitations in the quality of experimental vapor pressure-temperature data over a wide range of temperatures. Taking these problems into account, a consistent ΔH(vap) value is extracted from three studies for which experimental data are available. This analysis suggests that the most reliable value for ΔH(vap) is 13.0±0.1 at 410 K for use in force field optimization studies. The present results also indicate that similar analyses, including analysis of Antoine constants alone, may be of utility when reported ΔH(vap) values are not consistent for a given neat liquid.

Polarizable empirical force field for the primary and secondary alcohol series based on the classical Drude model.

A polarizable empirical force field based on the classical Drude oscillator has been developed for the aliphatic alcohol series. The model is optimized with emphasis on condensed-phase properties and is validated against a variety of experimental data. Transferability of the developed parameters is emphasized by the use of a single electrostatic model for the hydroxyl group throughout the alcohol series. Aliphatic moiety parameters were transferred from the polarizable alkane parameter set, with only the Lennard-Jones parameters on the carbon in methanol optimized. The developed model yields good agreement with pure solvent properties with the exception of the heats of vaporization of 1-propanol and 1-butanol, which are underestimated by approximately 6%; special LJ parameters for the oxygen in these two molecules that correct for this limitation are presented. Accurate treatment of the free energies of aqueous solvation required the use of atom-type specific O(alcohol)-O(water) LJ interaction terms, with specific terms used for the primary and secondary alcohols. With respect to gas phase properties the polarizable model overestimates experimental dipole moments and quantum mechanical interaction energies with water by approximately 10 and 8 %, respectively, a significant improvement over 44 and 46 % overestimations of the corresponding properties in the CHARMM22 fixed-charge additive model. Comparison of structural properties of the polarizable and additive models for the pure solvents and in aqueous solution shows significant differences indicating atomic details of intermolecular interactions to be sensitive to the applied force field. The polarizable model predicts pure solvent and aqueous phase dipole moment distributions for ethanol centered at 2.4 and 2.7 D, respectively, a significant increase over the gas phase value of 1.8 D, whereas in a solvent of lower polarity, benzene, a value of 1.9 is obtained. The ability of the polarizable model to yield changes in dipole moment as well as the reproduction of a range of condensed phase properties indicates its utility in the study of the properties of alcohols in a variety of condensed phase environments as well as representing an important step in the development of a comprehensive force field for biological molecules.

Additive and Classical Drude Polarizable Force Fields for Linear and Cyclic Ethers.

Empirical force field parameters consistent with the CHARMM additive and classical Drude based polarizable force fields are presented for linear and cyclic ethers. Initiation of the optimization process involved validation of the aliphatic parameters based on linear alkanes and cyclic alkanes. Results showed the transfer to cyclohexane to yield satisfactory agreement with target data; however, in the case of cyclopentane direct transfer of the Lennard-Jones parameters was not sufficient due to ring strain, requiring additional optimization of these parameters for this molecule. Parameters for the ethers were then developed starting with the available aliphatic parameters, with the nonbond parameters for the oxygens optimized to reproduce both gas- and condensed-phase properties. Nonbond parameters for the polarizable model include the use of an anisotropic electrostatic model on the oxygens. Parameter optimization emphasized the development of transferable parameters between the ethers of a given class. The ether models are shown to be in satisfactory agreement with both pure solvent and aqueous solvation properties, and the resulting parameters are transferable to test molecules. The presented force field will allow for simulation studies of ethers in condensed phase and provides a basis for ongoing developments in both additive and polarizable force fields for biological molecules.

Validation of Linear Scaling Semiempirical LocalSCF Method.

The numerical accuracy of linear scaling semiempirical methods LocalSCF and MOZYME is analyzed in comparison to conventional matrix diagonalization with respect to a variety of molecular properties including conformational energy, dipole moment, atomic charges, and bond orders. Major semiempirical MNDO, AM1, PM3, and PM5 Hamiltonians were considered in the study. As the numerical tests demonstrate, both LocalSCF and MOZYME reasonably reproduce matrix diagonalization results with the deviations being below the accuracy of semiempirical methods. However, the economical LocalSCF memory consumption and faster calculations are more beneficial for the quantum-mechanical modeling of large biological systems. The computational performance of the LocalSCF method is tested on the conformational energy calculation of a series of molecular dynamics snapshots of insulin in a large box of water.

Atomic Level Anisotropy in the Electrostatic Modeling of Lone Pairs for a Polarizable Force Field Based on the Classical Drude Oscillator.

Electron pairs in the valence shell of an atom that do not participate in the bonding of a molecule ("lone pairs") give rise to a concentrated electron density away from the atom center. To account for the asymmetry in the electron charge density that arises from lone pairs, an electrostatic model is developed that is parametrically anisotropic at the atomic level. The model uses virtual interaction sites with partial charges that are associated but not coincident with the nuclei. In addition, the model incorporates anisotropic atomic polarizabilities. The protocol previously outlined in Anisimov et al. [J. Chem. Theory Comput. 2005, 1, 153] for parametrizing the electrostatic potential energy of a polarizable force field using classical Drude oscillators is extended to incorporate additional lone pair parameters. To probe the electrostatic environment around the lone pairs, the static (molecule alone) and perturbed (molecule in the presence of a test charge) electrostatic potential (ESP) are evaluated and compared to high level quantum mechanical (QM) electronic structure calculations. The parametrization of the virtual sites relies on data from the QM static ESP. The contribution to the perturbed ESP from the electronic polarization of the molecule is used to resolve the components of the atomic polarizability tensor. The model is tested in the case of four molecules: methanol, acetone, methylamine, and pyridine. Interaction energies with water and sodium are used to assess the accuracy of the model. The results are compared with simpler models placing all the charge on the nuclei as well as using only isotropic atomic polarizabilities. Analysis shows that the addition of virtual sites reduces the average error relative to the QM calculations. In contrast to models with atom centered charges, the virtual site models correctly predict the minimum energy conformation for acetone and methanol, with water, to be closely coordinated with the lone pair direction. Furthermore, addition of anisotropic atomic polarizabilities to the virtual site model allows for precise fitting to the local perturbed QM ESP.