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.

Improving the Accuracy of Physics-Based Hydration-Free Energy Predictions by Machine Learning the Remaining Error Relative to the Experiment.

The accuracy of computational models of water is key to atomistic simulations of biomolecules. We propose a computationally efficient way to improve the accuracy of the prediction of hydration-free energies (HFEs) of small molecules: the remaining errors of the physics-based models relative to the experiment are predicted and mitigated by machine learning (ML) as a postprocessing step. Specifically, the trained graph convolutional neural network attempts to identify the "blind spots" in the physics-based model predictions, where the complex physics of aqueous solvation is poorly accounted for, and partially corrects for them. The strategy is explored for five classical solvent models representing various accuracy/speed trade-offs, from the fast analytical generalized Born (GB) to the popular TIP3P explicit solvent model; experimental HFEs of small neutral molecules from the FreeSolv set are used for the training and testing. For all of the models, the ML correction reduces the resulting root-mean-square error relative to the experiment for HFEs of small molecules, without significant overfitting and with negligible computational overhead. For example, on the test set, the relative accuracy improvement is 47% for the fast analytical GB, making it, after the ML correction, almost as accurate as uncorrected TIP3P. For the TIP3P model, the accuracy improvement is about 39%, bringing the ML-corrected model's accuracy below the 1 kcal/mol threshold. In general, the relative benefit of the ML corrections is smaller for more accurate physics-based models, reaching the lower limit of about 20% relative accuracy gain compared with that of the physics-based treatment alone. The proposed strategy of using ML to learn the remaining error of physics-based models offers a distinct advantage over training ML alone directly on reference HFEs: it preserves the correct overall trend, even well outside of the training set.

Strong interactions between highly dynamic lamina-associated domains and the nuclear envelope stabilize the 3D architecture of Drosophila interphase chromatin.

BACKGROUND: Interactions among topologically associating domains (TADs), and between the nuclear envelope (NE) and lamina-associated domains (LADs) are expected to shape various aspects of three-dimensional (3D) chromatin structure and dynamics; however, relevant genome-wide experiments that may provide statistically significant conclusions remain difficult. RESULTS: We have developed a coarse-grained dynamical model of D. melanogaster nuclei at TAD resolution that explicitly accounts for four distinct epigenetic classes of TADs and LAD-NE interactions. The model is parameterized to reproduce the experimental Hi-C map of the wild type (WT) nuclei; it describes time evolution of the chromatin over the G1 phase of the interphase. The simulations include an ensemble of nuclei, corresponding to the experimentally observed set of several possible mutual arrangements of chromosomal arms. The model is validated against multiple structural features of chromatin from several different experiments not used in model development. Predicted positioning of all LADs at the NE is highly dynamic-the same LAD can attach, detach and move far away from the NE multiple times during interphase. The probabilities of LADs to be in contact with the NE vary by an order of magnitude, despite all having the same affinity to the NE in the model. These probabilities are mostly determined by a highly variable local linear density of LADs along the genome, which also has the same strong effect on the predicted positioning of individual TADs -- higher probability of a TAD to be near NE is largely determined by a higher linear density of LADs surrounding this TAD. The distribution of LADs along the chromosome chains plays a notable role in maintaining a non-random average global structure of chromatin. Relatively high affinity of LADs to the NE in the WT nuclei substantially reduces sensitivity of the global radial chromatin distribution to variations in the strength of TAD-TAD interactions compared to the lamin depleted nuclei, where a small (0.5 kT) increase of cross-type TAD-TAD interactions doubles the chromatin density in the central nucleus region. CONCLUSIONS: A dynamical model of the entire fruit fly genome makes multiple genome-wide predictions of biological interest. The distribution of LADs along the chromatin chains affects their probabilities to be in contact with the NE and radial positioning of highly mobile TADs, playing a notable role in creating a non-random average global structure of the chromatin. We conjecture that an important role of attractive LAD-NE interactions is to stabilize global chromatin structure against inevitable cell-to-cell variations in TAD-TAD interactions.

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.

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.

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.

Accurate evaluation of charge asymmetry in aqueous solvation.

Charge hydration asymmetry (CHA)-a characteristic dependence of hydration free energy on the sign of the solute charge-quantifies the asymmetric response of water to electric field at microscopic level. Accurate estimates of CHA are critical for understanding hydration effects ubiquitous in chemistry and biology. However, measuring hydration energies of charged species is fraught with significant difficulties, which lead to unacceptably large (up to 300%) variation in the available estimates of the CHA effect. We circumvent these difficulties by developing a framework which allows us to extract and accurately estimate the intrinsic propensity of water to exhibit CHA from accurate experimental hydration free energies of neutral polar molecules. Specifically, from a set of 504 small molecules we identify two pairs that are analogous, with respect to CHA, to the K(+) /F(-) pair-a classical probe for the effect. We use these "CHA-conjugate" molecule pairs to quantify the intrinsic charge-asymmetric response of water to the microscopic charge perturbations: the asymmetry of the response is strong, ∼50% of the average hydration free energy of these molecules. The ability of widely used classical water models to predict hydration energies of small molecules correlates with their ability to predict CHA.

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.

Introducing Charge Hydration Asymmetry into the Generalized Born Model.

The effect of charge hydration asymmetry (CHA)-non-invariance of solvation free energy upon solute charge inversion-is missing from the standard linear response continuum electrostatics. The proposed charge hydration asymmetric-generalized Born (CHA-GB) approximation introduces this effect into the popular generalized Born (GB) model. The CHA is added to the GB equation via an analytical correction that quantifies the specific propensity of CHA of a given water model; the latter is determined by the charge distribution within the water model. Significant variations in CHA seen in explicit water (TIP3P, TIP4P-Ew, and TIP5P-E) free energy calculations on charge-inverted "molecular bracelets" are closely reproduced by CHA-GB, with the accuracy similar to models such as SEA and 3D-RISM that go beyond the linear response. Compared against reference explicit (TIP3P) electrostatic solvation free energies, CHA-GB shows about a 40% improvement in accuracy over the canonical GB, tested on a diverse set of 248 rigid small neutral molecules (root mean square error, rmse = 0.88 kcal/mol for CHA-GB vs 1.24 kcal/mol for GB) and 48 conformations of amino acid analogs (rmse = 0.81 kcal/mol vs 1.26 kcal/mol). CHA-GB employs a novel definition of the dielectric boundary that does not subsume the CHA effects into the intrinsic atomic radii. The strategy leads to finding a new set of intrinsic atomic radii optimized for CHA-GB; these radii show physically meaningful variation with the atom type, in contrast to the radii set optimized for GB. Compared to several popular radii sets used with the original GB model, the new radii set shows better transferability between different classes of molecules.

Charge hydration asymmetry: the basic principle and how to use it to test and improve water models.

Charge hydration asymmetry (CHA) manifests itself in the experimentally observed strong dependence of free energy of ion hydration on the sign of the ion charge. This asymmetry is not consistently accounted for by popular models of solvation; its magnitude varies greatly between the models. While it is clear that CHA is somehow related to charge distribution within a water molecule, the exact nature of this relationship is unknown. We propose a simple, yet general and rigorous criterion that relates rotational and charge inversion properties of a water molecule's charge distribution with its ability to cause CHA. We show which electric multipole components of a water molecule are key to explain its ability for asymmetric charge hydration. We then test several popular water models and explain why specific models show none, little, or strong CHA in simulations. We use the gained insight to derive an analogue of the Born equation that includes the missing physics necessary to account for CHA and does not rely on redefining the continuum dielectric boundary. The proposed formula is as simple as the original, does not contain any fitting parameters, and predicts hydration free energies and entropies of spherical cations and anions within experimental uncertainty. Our findings suggest that the gap between the practical continuum electrostatics framework and the more fundamental explicit solvent treatment may be reduced considerably by explicitly introducing CHA into the existing continuum framework.

Binding free energy and counterion release for adsorption of the antimicrobial peptide lactoferricin B on a POPG membrane.

Molecular dynamics (MD) simulations are used to study the interaction of an anionic palmitoyl-oleoyl-phosphatidylglycerol (POPG) bilayer with the cationic antimicrobial peptide bovine lactoferricin (LFCinB) in a 100 mM NaCl solution at 310 K. The interaction of LFCinB with a POPG bilayer is employed as a model system for studying the details of membrane adsorption selectivity of cationic antimicrobial peptides. Seventy eight 4 ns MD production run trajectories of the equilibrated system, with six restrained orientations of LFCinB at 13 different separations from the POPG membrane, are generated to determine the free energy profile for the peptide as a function of the distance between LFCinB and the membrane surface. To calculate the profile for this relatively large system, a variant of constrained MD and thermodynamic integration is used. A simplified method for relating the free energy profile to the LFCinB-POPG membrane binding constant is employed to predict a free energy of adsorption of -5.4+/-1.3 kcal/mol and a corresponding maximum adsorption binding force of about 58 pN. We analyze the results using Poisson-Boltzmann theory. We find the peptide-membrane attraction to be dominated by the entropy increase due to the release of counterions and polarized water from the region between the charged membrane and peptide, as the two approach each other. We contrast these results with those found earlier for adsorption of LFCinB on the mammalianlike palmitoyl-oleoyl-phosphatidylcholine membrane.

Prediction of binding free energy for adsorption of antimicrobial peptide lactoferricin B on a POPC membrane.

Molecular dynamics (MD) simulations are used to study the interaction of a zwitterionic palmitoyl-oleoyl-phosphatidylcholine (POPC) bilayer with the cationic antimicrobial peptide bovine lactoferricin (LFCinB) in a 100 mM NaCl solution at 310 K. The interaction of LFCinB with POPC is used as a model system for studying the details of membrane-peptide interactions, with the peptide selected because of its antimicrobial nature. Seventy-two 3 ns MD simulations, with six orientations of LFCinB at 12 different distances from a POPC membrane, are carried out to determine the potential of mean force (PMF) or free energy profile for the peptide as a function of the distance between LFCinB and the membrane surface. To calculate the PMF for this relatively large system a new variant of constrained MD and thermodynamic integration is developed. A simplified method for relating the PMF to the LFCinB-membrane binding free energy is described and used to predict a free energy of adsorption (or binding) of -1.05+/-0.39 kcal/mol , and corresponding maximum binding force of about 20 pN, for LFCinB-POPC. The contributions of the ions-LFCinB and the water-LFCinB interactions to the PMF are discussed. The method developed will be a useful starting point for future work simulating peptides interacting with charged membranes and interactions involved in the penetration of membranes, features necessary to understand in order to rationally design peptides as potential alternatives to traditional antibiotics.

Nonequilibrium molecular dynamics calculation of the conductance of the KcsA potassium ion channel.

We have developed an efficient method for the atomistic determination of the conductance of a biological ion channel model by applying an external field to the conducting ions only. The underlying theory is discussed and demonstrated in a simple test system consisting of two ions in a box of water. Finally, the theory is applied to the experimentally determined structure of the KcsA potassium channel from which a conductance in reasonable agreement with the experimental result is predicted.

Unified modeling of conductance kinetics for low- and high-conductance potassium ion channels.

A kinetics model is proposed for the description of ion conductance of low- and high-conductance potassium ion channels. The model describes ion permeation through the selectivity filter, which is assumed to be the only conductance determining part of the open channel. The filter occupancy can vary from zero to three ions, affecting the ion entry and exit rates. Ion motion between the binding sites inside the filter is assumed fast compared to the latter rates allowing averaging the equilibrium entry and exit rate constants over the possible ion configurations in the filter with a particular occupancy. Averaged rate constants related to a pair of adjacent occupancy states characterize a particular ion permeation mechanism. An expression for the channel conductance as a function of the symmetrical external ion concentration is derived. It comprises a sum of concentration independent conductance amplitudes for different ion permeation mechanisms weighted by the equilibrium filter occupancy probabilities. It is shown that each amplitude (i.e., maximum contribution to the channel conductance from each conductance mechanism) is proportional to an averaged exit rate constant and to quantities characterizing the effect of the applied electric field on the rate constants and the equilibrium ion distribution in the filter. The conductance expression derived provides a good description of the experimentally observed conductance-concentration curves for low-conductance (e.g., Kir2.1) and high-conductance (e.g., KcsA) potassium channels. It enables one to obtain equilibrium ion binding constants at different filter occupancies and to calculate the average number of ions in the selectivity filter for a given external ion concentration. For KcsA this number (2.0 at 200 mM) is in a good agreement with the available experimental value (2.1 at 200 mM). For the high-conductance potassium channels the net negative electrical charge around the selectivity filter increases the ion binding constants, thereby causing the larger occupancy probabilities to occur at smaller external ion concentrations compared to the low-conductance channels. This substantially increases the contributions of the two- and three-ion permeation mechanisms, with the larger conductance amplitudes leading to increased channel conductance compared to the low-conductance channels.

Non-Michaelis-Menten kinetics model for conductance of low-conductance potassium ion channels.

A reduced kinetics model is proposed for ion permeation in low-conductance potassium ion channels with zero net electrical charge in the selectivity filter region. The selectivity filter is assumed to be the only conductance-determining part of the channel. Ion entry and exit rate constants depend on the occupancy of the filter due to ion-ion interactions. The corresponding rates are assumed slow relative to the rates of ion motion between binding sites inside the filter, allowing a reduction of the kinetics model of the filter by averaging the entry and exit rate constants over the states with a particular occupancy number. The reduced kinetics model for low-conductance channels is described by only three states and two sets of effective rate constants characterizing transitions between these states. An explicit expression for the channel conductance as a function of symmetrical external ion concentration is derived under the assumption that the average electrical mobility of ions in the selectivity filter region in a limited range of ion concentrations does not depend on these concentrations. The simplified conductance model is shown to provide a good description of the experimentally observed conductance-concentration curve for the low-conductance potassium channel Kir2.1, and also predicts the mean occupancy of the selectivity filter of this channel. We find that at physiological external ion concentrations this occupancy is much lower than the value of two ions observed for one of the high-conductance potassium channels, KcsA.