AmberTools.

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.

Why Na+ has higher propensity than K+ to condense DNA in a crowded environment.

Experimentally, in the presence of the crowding agent polyethylene glycol (PEG), sodium ions compact double-stranded DNA more readily than potassium ions. Here, we have used molecular dynamics simulations and the "ion binding shells model" of DNA condensation to provide an explanation for the observed variations in condensation of short DNA duplexes in solutions containing different monovalent cations and PEG; several predictions are made. According to the model we use, externally bound ions contribute the most to the ion-induced aggregation of DNA duplexes. The simulations reveal that for two adjacent DNA duplexes, the number of externally bound Na+ ions is larger than the number of K+ ions over a wide range of chloride concentrations in the presence of PEG, providing a qualitative explanation for the higher propensity of sodium ions to compact DNA under crowded conditions. The qualitative picture is confirmed by an estimate of the corresponding free energy of DNA aggregation that is at least 0.2kBT per base pair more favorable in solution with NaCl than with KCl at the same ion concentration. The estimated attraction free energy of DNA duplexes in the presence of Na+ depends noticeably on the DNA sequence; we predict that AT-rich DNA duplexes are more readily condensed than GC-rich ones in the presence of Na+. Counter-intuitively, the addition of a small amount of a crowding agent with high affinity for the specific condensing ion may lead to the weakening of the ion-mediated DNA-DNA attraction, shifting the equilibrium away from the DNA condensed phase.

Explicit ions/implicit water generalized Born model for nucleic acids.

The ion atmosphere around highly charged nucleic acid molecules plays a significant role in their dynamics, structure, and interactions. Here we utilized the implicit solvent framework to develop a model for the explicit treatment of ions interacting with nucleic acid molecules. The proposed explicit ions/implicit water model is based on a significantly modified generalized Born (GB) model and utilizes a non-standard approach to define the solute/solvent dielectric boundary. Specifically, the model includes modifications to the GB interaction terms for the case of multiple interacting solutes-disconnected dielectric boundary around the solute-ion or ion-ion pairs. A fully analytical description of all energy components for charge-charge interactions is provided. The effectiveness of the approach is demonstrated by calculating the potential of mean force for Na+-Cl- ion pair and by carrying out a set of Monte Carlo (MC) simulations of mono- and trivalent ions interacting with DNA and RNA duplexes. The monovalent (Na+) and trivalent (CoHex3+) counterion distributions predicted by the model are in close quantitative agreement with all-atom explicit water molecular dynamics simulations used as reference. Expressed in the units of energy, the maximum deviations of local ion concentrations from the reference are within k B T. The proposed explicit ions/implicit water GB model is able to resolve subtle features and differences of CoHex distributions around DNA and RNA duplexes. These features include preferential CoHex binding inside the major groove of the RNA duplex, in contrast to CoHex biding at the "external" surface of the sugar-phosphate backbone of the DNA duplex; these differences in the counterion binding patters were earlier shown to be responsible for the observed drastic differences in condensation propensities between short DNA and RNA duplexes. MC simulations of CoHex ions interacting with the homopolymeric poly(dA·dT) DNA duplex with modified (de-methylated) and native thymine bases are used to explore the physics behind CoHex-thymine interactions. The simulations suggest that the ion desolvation penalty due to proximity to the low dielectric volume of the methyl group can contribute significantly to CoHex-thymine interactions. Compared to the steric repulsion between the ion and the methyl group, the desolvation penalty interaction has a longer range and may be important to consider in the context of methylation effects on DNA condensation.

Spermine Condenses DNA, but Not RNA Duplexes.

Interactions between the polyamine spermine and nucleic acids drive important cellular processes. Spermine condenses DNA and some RNAs, such as poly(rA):poly(rU). A large fraction of the spermine present in cells is bound to RNA but apparently does not condense it. Here, we study the effect of spermine binding to short duplex RNA and DNA, and compare our findings with predictions of molecular-dynamics simulations. When small numbers of spermine are introduced, RNA with a designed sequence containing a mixture of 14 GC pairs and 11 AU pairs resists condensation relative to DNA of an equivalent sequence or to 25 bp poly(rA):poly(rU) RNA. A comparison of wide-angle x-ray scattering profiles with simulation results suggests that spermine is sequestered deep within the major groove of mixed-sequence RNA. This prevents condensation by limiting opportunities to bridge to other molecules and stabilizes the RNA by locking it into a particular conformation. In contrast, for DNA, simulations suggest that spermine binds externally to the duplex, offering opportunities for intermolecular interaction. The goal of this study is to explain how RNA can remain soluble and available for interaction with other molecules in the cell despite the presence of spermine at concentrations high enough to precipitate DNA.

Partially Assembled Nucleosome Structures at Atomic Detail.

The evidence is now overwhelming that partially assembled nucleosome states (PANS) are as important as the canonical nucleosome structure for the understanding of how accessibility to genomic DNA is regulated in cells. We use a combination of molecular dynamics simulation and atomic force microscopy to deliver, in atomic detail, structural models of three key PANS: the hexasome (H2A·H2B)·(H3·H4)2, the tetrasome (H3·H4)2, and the disome (H3·H4). Despite fluctuations of the conformation of the free DNA in these structures, regions of protected DNA in close contact with the histone core remain stable, thus establishing the basis for the understanding of the role of PANS in DNA accessibility regulation. On average, the length of protected DNA in each structure is roughly 18 basepairs per histone protein. Atomistically detailed PANS are used to explain experimental observations; specifically, we discuss interpretation of atomic force microscopy, Förster resonance energy transfer, and small-angle x-ray scattering data obtained under conditions when PANS are expected to exist. Further, we suggest an alternative interpretation of a recent genome-wide study of DNA protection in active chromatin of fruit fly, leading to a conclusion that the three PANS are present in actively transcribing regions in a substantial amount. The presence of PANS may not only be a consequence, but also a prerequisite for fast transcription in vivo.

Grid-Based Surface Generalized Born Model for Calculation of Electrostatic Binding Free Energies.

Fast and accurate calculation of solvation free energies is central to many applications, such as rational drug design. In this study, we present a grid-based molecular surface implementation of "R6" flavor of the generalized Born (GB) implicit solvent model, named GBNSR6. The speed, accuracy relative to numerical Poisson-Boltzmann treatment, and sensitivity to grid surface parameters are tested on a set of 15 small protein-ligand complexes and a set of biomolecules in the range of 268 to 25099 atoms. Our results demonstrate that the proposed model provides a relatively successful compromise between the speed and accuracy of computing polar components of the solvation free energies (ΔGpol) and binding free energies (ΔΔGpol). The model tolerates a relatively coarse grid size h = 0.5 Å, where the grid artifact error in computing ΔΔGpol remains in the range of kBT ∼ 0.6 kcal/mol. The estimated ΔΔGpols are well correlated (r2 = 0.97) with the numerical Poisson-Boltzmann reference, while showing virtually no systematic bias and RMSE = 1.43 kcal/mol. The grid-based GBNSR6 model is available in Amber (AmberTools) package of molecular simulation programs.

The Role of Correlation and Solvation in Ion Interactions with B-DNA.

The ionic atmospheres around nucleic acids play important roles in biological function. Large-scale explicit solvent simulations coupled to experimental assays such as anomalous small-angle x-ray scattering can provide important insights into the structure and energetics of such atmospheres but are time- and resource intensive. In this article, we use classical density functional theory to explore the balance among ion-DNA, ion-water, and ion-ion interactions in ionic atmospheres of RbCl, SrCl2, and CoHexCl3 (cobalt hexamine chloride) around a B-form DNA molecule. The accuracy of the classical density functional theory calculations was assessed by comparison between simulated and experimental anomalous small-angle x-ray scattering curves, demonstrating that an accurate model should take into account ion-ion correlation and ion hydration forces, DNA topology, and the discrete distribution of charges on the DNA backbone. As expected, these calculations revealed significant differences among monovalent, divalent, and trivalent cation distributions around DNA. Approximately half of the DNA-bound Rb(+) ions penetrate into the minor groove of the DNA and half adsorb on the DNA backbone. The fraction of cations in the minor groove decreases for the larger Sr(2+) ions and becomes zero for CoHex(3+) ions, which all adsorb on the DNA backbone. The distribution of CoHex(3+) ions is mainly determined by Coulomb and steric interactions, while ion-correlation forces play a central role in the monovalent Rb(+) distribution and a combination of ion-correlation and hydration forces affect the Sr(2+) distribution around DNA. This does not imply that correlations in CoHex solutions are weaker or stronger than for other ions. Steric inaccessibility of the grooves to large CoHex ions leads to their binding at the DNA surface. In this binding mode, first-order electrostatic interactions (Coulomb) dominate the overall binding energy as evidenced by low sensitivity of ionic distribution to the presence or absence of second-order electrostatic correlation interactions.

Opposing Effects of Multivalent Ions on the Flexibility of DNA and RNA.

Increasing the concentration of counterions (salt) is known to reduce the bending persistence length of DNA. Here we use atomistic molecular dynamics simulations to predict that multivalent counterions have the opposite effect on double-stranded RNA, increasing its bending rigidity by at least 30%. This counterintuitive effect is observed for various tri- and tetravalent ions alike, and is robust to methodological details and the RNA sequence. In contrast to DNA, multivalent counterions bind inside the RNA major groove, causing significant contraction of the molecule along its helical axis-as a result, its further deformation due to bending becomes energetically more expensive compared to bending without bound multivalent ions. Thus, the relationship between mechanical properties of a charged polymer and its ionic atmosphere may be richer than previously thought.

Accuracy limit of rigid 3-point water models.

Classical 3-point rigid water models are most widely used due to their computational efficiency. Recently, we introduced a new approach to constructing classical rigid water models [S. Izadi et al., J. Phys. Chem. Lett. 5, 3863 (2014)], which permits a virtually exhaustive search for globally optimal model parameters in the sub-space that is most relevant to the electrostatic properties of the water molecule in liquid phase. Here we apply the approach to develop a 3-point Optimal Point Charge (OPC3) water model. OPC3 is significantly more accurate than the commonly used water models of same class (TIP3P and SPCE) in reproducing a comprehensive set of liquid bulk properties, over a wide range of temperatures. Beyond bulk properties, we show that OPC3 predicts the intrinsic charge hydration asymmetry (CHA) of water - a characteristic dependence of hydration free energy on the sign of the solute charge - in very close agreement with experiment. Two other recent 3-point rigid water models, TIP3PFB and H2ODC, each developed by its own, completely different optimization method, approach the global accuracy optimum represented by OPC3 in both the parameter space and accuracy of bulk properties. Thus, we argue that an accuracy limit of practical 3-point rigid non-polarizable models has effectively been reached; remaining accuracy issues are discussed.

Multi-shell model of ion-induced nucleic acid condensation.

We present a semi-quantitative model of condensation of short nucleic acid (NA) duplexes induced by trivalent cobalt(iii) hexammine (CoHex) ions. The model is based on partitioning of bound counterion distribution around single NA duplex into "external" and "internal" ion binding shells distinguished by the proximity to duplex helical axis. In the aggregated phase the shells overlap, which leads to significantly increased attraction of CoHex ions in these overlaps with the neighboring duplexes. The duplex aggregationfree energy is decomposed into attractive and repulsive components in such a way that they can be represented by simple analytical expressions with parameters derived from molecular dynamic simulations and numerical solutions of Poisson equation. The attractive term depends on the fractions of bound ions in the overlapping shells and affinity of CoHex to the "external" shell of nearly neutralized duplex. The repulsive components of the free energy are duplex configurational entropy loss upon the aggregation and the electrostatic repulsion of the duplexes that remains after neutralization by bound CoHex ions. The estimates of the aggregationfree energy are consistent with the experimental range of NA duplex condensation propensities, including the unusually poor condensation of RNA structures and subtle sequence effects upon DNAcondensation. The model predicts that, in contrast to DNA, RNA duplexes may condense into tighter packed aggregates with a higher degree of duplex neutralization. An appreciable CoHex mediated RNA-RNA attraction requires closer inter-duplex separation to engage CoHex ions (bound mostly in the "internal" shell of RNA) into short-range attractive interactions. The model also predicts that longer NA fragments will condense more readily than shorter ones. The ability of this model to explain experimentally observed trends in NAcondensation lends support to proposed NAcondensation picture based on the multivalent "ion binding shells."

Understanding nucleic acid structural changes by comparing wide-angle x-ray scattering (WAXS) experiments to molecular dynamics simulations.

Wide-angle x-ray scattering (WAXS) is emerging as a powerful tool for increasing the resolution of solution structure measurements of biomolecules. Compared to its better known complement, small angle x-ray scattering (SAXS), WAXS targets higher scattering angles and can enhance structural studies of molecules by accessing finer details of solution structures. Although the extension from SAXS to WAXS is easy to implement experimentally, the computational tools required to fully harness the power of WAXS are still under development. Currently, WAXS is employed to study structural changes and ligand binding in proteins; however, the methods are not as fully developed for nucleic acids. Here, we show how WAXS can qualitatively characterize nucleic acid structures as well as the small but significant structural changes driven by the addition of multivalent ions. We show the potential of WAXS to test all-atom molecular dynamics (MD) simulations and to provide insight into understanding how the trivalent ion cobalt(III) hexammine (CoHex) affects the structure of RNA and DNA helices. We find that MD simulations capture the RNA structural change that occurs due to addition of CoHex.

Why double-stranded RNA resists condensation.

The addition of small amounts of multivalent cations to solutions containing double-stranded DNA leads to inter-DNA attraction and eventual condensation. Surprisingly, the condensation is suppressed in double-stranded RNA, which carries the same negative charge as DNA, but assumes a different double helical form. Here, we combine experiment and atomistic simulations to propose a mechanism that explains the variations in condensation of short (25 base-pairs) nucleic acid (NA) duplexes, from B-like form of homopolymeric DNA, to mixed sequence DNA, to DNA:RNA hybrid, to A-like RNA. Circular dichroism measurements suggest that duplex helical geometry is not the fundamental property that ultimately determines the observed differences in condensation. Instead, these differences are governed by the spatial variation of cobalt hexammine (CoHex) binding to NA. There are two major NA-CoHex binding modes--internal and external--distinguished by the proximity of bound CoHex to the helical axis. We find a significant difference, up to 5-fold, in the fraction of ions bound to the external surfaces of the different NA constructs studied. NA condensation propensity is determined by the fraction of CoHex ions in the external binding mode.

Two-phase stretching of molecular chains.

Although stretching of most polymer chains leads to rather featureless force-extension diagrams, some, notably DNA, exhibit nontrivial behavior with a distinct plateau region. Here, we propose a unified theory that connects force-extension characteristics of the polymer chain with the convexity properties of the extension energy profile of its individual monomer subunits. Namely, if the effective monomer deformation energy as a function of its extension has a nonconvex (concave up) region, the stretched polymer chain separates into two phases: the weakly and strongly stretched monomers. Simplified planar and 3D polymer models are used to illustrate the basic principles of the proposed model. Specifically, we show rigorously that, when the secondary structure of a polymer is mostly caused by weak noncovalent interactions, the stretching is two phase, and the force-stretching diagram has the characteristic plateau. We then use realistic coarse-grained models to confirm the main findings and make direct connection to the microscopic structure of the monomers. We show in detail how the two-phase scenario is realized in the α-helix and DNA double helix. The predicted plateau parameters are consistent with single-molecules experiments. Detailed analysis of DNA stretching shows that breaking of Watson-Crick bonds is not necessary for the existence of the plateau, although some of the bonds do break as the double helix extends at room temperature. The main strengths of the proposed theory are its generality and direct microscopic connection.

Speed of conformational change: comparing explicit and implicit solvent molecular dynamics simulations.

Adequate sampling of conformation space remains challenging in atomistic simulations, especially if the solvent is treated explicitly. Implicit-solvent simulations can speed up conformational sampling significantly. We compare the speed of conformational sampling between two commonly used methods of each class: the explicit-solvent particle mesh Ewald (PME) with TIP3P water model and a popular generalized Born (GB) implicit-solvent model, as implemented in the AMBER package. We systematically investigate small (dihedral angle flips in a protein), large (nucleosome tail collapse and DNA unwrapping), and mixed (folding of a miniprotein) conformational changes, with nominal simulation times ranging from nanoseconds to microseconds depending on system size. The speedups in conformational sampling for GB relative to PME simulations, are highly system- and problem-dependent. Where the simulation temperatures for PME and GB are the same, the corresponding speedups are approximately onefold (small conformational changes), between ∼1- and ∼100-fold (large changes), and approximately sevenfold (mixed case). The effects of temperature on speedup and free-energy landscapes, which may differ substantially between the solvent models, are discussed in detail for the case of miniprotein folding. In addition to speeding up conformational sampling, due to algorithmic differences, the implicit solvent model can be computationally faster for small systems or slower for large systems, depending on the number of solute and solvent atoms. For the conformational changes considered here, the combined speedups are approximately twofold, ∼1- to 60-fold, and ∼50-fold, respectively, in the low solvent viscosity regime afforded by the implicit solvent. For all the systems studied, 1) conformational sampling speedup increases as Langevin collision frequency (effective viscosity) decreases; and 2) conformational sampling speedup is mainly due to reduction in solvent viscosity rather than possible differences in free-energy landscapes between the solvent models.

Protonation and pK changes in protein-ligand binding.

Formation of protein-ligand complexes causes various changes in both the receptor and the ligand. This review focuses on changes in pK and protonation states of ionizable groups that accompany protein-ligand binding. Physical origins of these effects are outlined, followed by a brief overview of the computational methods to predict them and the associated corrections to receptor-ligand binding affinities. Statistical prevalence, magnitude and spatial distribution of the pK and protonation state changes in protein-ligand binding are discussed in detail, based on both experimental and theoretical studies. While there is no doubt that these changes occur, they do not occur all the time; the estimated prevalence varies, both between individual complexes and by method. The changes occur not only in the immediate vicinity of the interface but also sometimes far away. When receptor-ligand binding is associated with protonation state change at particular pH, the binding becomes pH dependent: we review the interplay between sub-cellular characteristic pH and optimum pH of receptor-ligand binding. It is pointed out that there is a tendency for protonation state changes upon binding to be minimal at physiologically relevant pH for each complex (no net proton uptake/release), suggesting that native receptor-ligand interactions have evolved to reduce the energy cost associated with ionization changes. As a result, previously reported statistical prevalence of these changes - typically computed at the same pH for all complexes - may be higher than what may be expected at optimum pH specific to each complex. We also discuss whether proper account of protonation state changes appears to improve practical docking and scoring outcomes relevant to structure-based drug design. An overview of some of the existing challenges in the field is provided in conclusion.

H++ 3.0: automating pK prediction and the preparation of biomolecular structures for atomistic molecular modeling and simulations.

The accuracy of atomistic biomolecular modeling and simulation studies depend on the accuracy of the input structures. Preparing these structures for an atomistic modeling task, such as molecular dynamics (MD) simulation, can involve the use of a variety of different tools for: correcting errors, adding missing atoms, filling valences with hydrogens, predicting pK values for titratable amino acids, assigning predefined partial charges and radii to all atoms, and generating force field parameter/topology files for MD. Identifying, installing and effectively using the appropriate tools for each of these tasks can be difficult for novice and time-consuming for experienced users. H++ (http://biophysics.cs.vt.edu/) is a free open-source web server that automates the above key steps in the preparation of biomolecular structures for molecular modeling and simulations. H++ also performs extensive error and consistency checking, providing error/warning messages together with the suggested corrections. In addition to numerous minor improvements, the latest version of H++ includes several new capabilities and options: fix erroneous (flipped) side chain conformations for HIS, GLN and ASN, include a ligand in the input structure, process nucleic acid structures and generate a solvent box with specified number of common ions for explicit solvent MD.

Accuracy of continuum electrostatic calculations based on three common dielectric boundary definitions.

We investigate the influence of three common definitions of the solute/solvent dielectric boundary (DB) on the accuracy of the electrostatic solvation energy ΔGel computed within the Poisson Boltzmann and the generalized Born models of implicit solvation. The test structures include small molecules, peptides and small proteins; explicit solvent ΔGel are used as accuracy reference. For common atomic radii sets BONDI, PARSE (and ZAP9 for small molecules) the use of van der Waals (vdW) DB results, on average, in considerably larger errors in ΔGel than the molecular surface (MS) DB. The optimal probe radius ρw for which the MS DB yields the most accurate ΔGel varies considerably between structure types. The solvent accessible surface (SAS) DB becomes optimal at ρw ~ 0.2 Å (exact value is sensitive to the structure and atomic radii), at which point the average accuracy of ΔGel is comparable to that of the MS-based boundary. The geometric equivalence of SAS to vdW surface based on the same atomic radii uniformly increased by ρw gives the corresponding optimal vdW DB. For small molecules, the optimal vdW DB based on BONDI + 0.2 Å radii can yield ΔGel estimates at least as accurate as those based on the optimal MS DB. Also, in small molecules, pairwise charge-charge interactions computed with the optimal vdW DB are virtually equal to those computed with the MS DB, suggesting that in this case the two boundaries are practically equivalent by the electrostatic energy criteria. In structures other than small molecules, the optimal vdW and MS dielectric boundaries are not equivalent: the respective pairwise electrostatic interactions in the presence of solvent can differ by up to 5 kcal/mol for individual atomic pairs in small proteins, even when the total ΔGel are equal. For small proteins, the average decrease in pairwise electrostatic interactions resulting from the switch from optimal MS to optimal vdW DB definition can be mimicked within the MS DB definition by doubling of the solute dielectric constant. However, the use of the higher interior dielectric does not eliminate the large individual deviations between pairwise interactions computed within the two DB definitions. It is argued that while the MS based definition of the dielectric boundary is more physically correct in some types of practical calculations, the choice is not so clear in some other common scenarios.

Inclusion of Water Multipoles into the Implicit Solvation Framework Leads to Accuracy Gains.

The current practical "workhorses" of the atomistic implicit solvation─the Poisson-Boltzmann (PB) and generalized Born (GB) models─face fundamental accuracy limitations. Here, we propose a computationally efficient implicit solvation framework, the Implicit Water Multipole GB (IWM-GB) model, that systematically incorporates the effects of multipole moments of water molecules in the first hydration shell of a solute, beyond the dipole water polarization already present at the PB/GB level. The framework explicitly accounts for coupling between polar and nonpolar contributions to the total solvation energy, which is missing from many implicit solvation models. An implementation of the framework, utilizing the GAFF force field and AM1-BCC atomic partial charges model, is parametrized and tested against the experimental hydration free energies of small molecules from the FreeSolv database. The resulting accuracy on the test set (RMSE ∼ 0.9 kcal/mol) is 12% better than that of the explicit solvation (TIP3P) treatment, which is orders of magnitude slower. We also find that the coupling between polar and nonpolar parts of the solvation free energy is essential to ensuring that several features of the IWM-GB model are physically meaningful, including the sign of the nonpolar contributions.

The probability of chromatin to be at the nuclear lamina has no systematic effect on its transcription level in fruit flies.

BACKGROUND: Multiple studies have demonstrated a negative correlation between gene expression and positioning of genes at the nuclear envelope (NE) lined by nuclear lamina, but the exact relationship remains unclear, especially in light of the highly stochastic, transient nature of the gene association with the NE. RESULTS: In this paper, we ask whether there is a causal, systematic, genome-wide relationship between the expression levels of the groups of genes in topologically associating domains (TADs) of Drosophila nuclei and the probabilities of TADs to be found at the NE. To investigate the nature of this possible relationship, we combine a coarse-grained dynamic model of the entire Drosophila nucleus with genome-wide gene expression data; we analyze the TAD averaged transcription levels of genes against the probabilities of individual TADs to be in contact with the NE in the control and lamins-depleted nuclei. Our findings demonstrate that, within the statistical error margin, the stochastic positioning of Drosophila melanogaster TADs at the NE does not, by itself, systematically affect the mean level of gene expression in these TADs, while the expected negative correlation is confirmed. The correlation is weak and disappears completely for TADs not containing lamina-associated domains (LADs) or TADs containing LADs, considered separately. Verifiable hypotheses regarding the underlying mechanism for the presence of the correlation without causality are discussed. These include the possibility that the epigenetic marks and affinity to the NE of a TAD are determined by various non-mutually exclusive mechanisms and remain relatively stable during interphase. CONCLUSIONS: At the level of TADs, the probability of chromatin being in contact with the nuclear envelope has no systematic, causal effect on the transcription level in Drosophila. The conclusion is reached by combining model-derived time-evolution of TAD locations within the nucleus with their experimental gene expression levels.