Computational methods for biomolecular electrostatics

Underlying features for the enhanced electrostatic strength of the extremophilic malate dehydrogenase interface salt-bridge compared to the mesophilic one

Malate dehydrogenase (MDH) exists in multimeric form in normal and extreme solvent conditions where residues of the interface are involved in specific interactions. The interface salt-bridge (ISB) and its microenvironment (ME) residues may have a crucial role in the stability and specificity of the interface. To gain insight into this, we have analyzed 218 ISBs from 42 interfaces of 15 crystal structures along with their sequences. Comparative analyses demonstrate that the ISB strength is ∼30 times greater in extremophilic cases than that of the normal one. To this end, the interface residue propensity, ISB design and pair selection, and ME-residue's types, i.e., type-I and type-II, are seen to be intrinsically involved. Although Type-I is a common type, Type-II appears to be extremophile-specific, where the net ME-residue count is much lower with an excessive net ME-energy contribution, which seems to be a novel interface compaction strategy. Furthermore, the interface strength can be enhanced by selecting the desired mutant from the net-energy profile of all possible mutations of an unfavorable ME-residue. The study that applies to other similar systems finds applications in protein-protein interaction and protein engineering.

Insight on the interaction between the scorpion toxin blocker Discrepin on potassium voltage-gated channel Kv4.3 by molecular dynamics simulations

Discrepin is a 38-residue α-toxin extracted from the venom of the Venezuelan scorpion Tityus discrepans, which inhibits ionic transit in the voltage-dependent potassium channels (Kv) of A-type current. The effect of specific residues on the IC₅₀ between Discrepine and Kv4.3, the main component of A-type currents, is known; however, the molecular details of the toxin-channel interaction are not known. In this work, we present interaction models between Discrepin (wt) and two peptide variants (V6K/D20K and K13A) on the pore-forming domain of the Kv4.3 channel obtained from homology, docking, and molecular dynamics modeling techniques. The free energy calculations in these models correspond to the order of the experimentally determined IC₅₀ values. Our studies shed light on the role of the K13 residue as responsible for occluding the Kv4.3 selectivity filter and the importance of the V6K mutation in the approach and stabilization of toxin-channel complex interactions.Communicated by Ramaswamy H. Sarma.

Efficient evaluation of electrostatic potential with computerized optimized code

The evaluation of molecular electrostatic potential (ESP) is a performance bottleneck for many computational chemical tasks like restrained ESP charge fitting or quantum mechanics/molecular mechanics simulations. In this paper, an efficient algorithm for the evaluation of ESP is proposed. It regroups the expression in terms of primitive Gaussian type orbitals (GTOs) with identical angular momentum types and nuclei centers. Each term is calculated using a computerized optimized code. This algorithm was integrated into the wavefunction analysis program Multiwfn and was tested on several large systems. In the cases of dopamine and remdesivir, the performance of this algorithm was comparable to or better than some popular state-of-the-art codes. For meta1-organic framework-5, where the number of GTOs and ESP points is 4840 and 259 262, respectively, our code could finish the evaluation in 1874 seconds on ordinary hardware. It also exhibits good parallelization scaling. The source code of this algorithm is freely available and can become a useful tool for computational chemists.

Testing the Limitations of MD-Based Local Electric Fields Using the Vibrational Stark Effect in Solution: Penicillin G as a Test Case

Noncovalent interactions underlie nearly all molecular processes in the condensed phase from solvation to catalysis. Their quantification within a physically consistent framework remains challenging. Experimental vibrational Stark effect (VSE)-based solvatochromism can be combined with molecular dynamics (MD) simulations to quantify the electrostatic forces in solute-solvent interactions for small rigid molecules and, by extension, when these solutes bind in enzyme active sites. While generalizing this approach toward more complex (bio)molecules, such as the conformationally flexible and charged penicillin G (PenG), we were surprised to observe inconsistencies in MD-based electric fields. Combining synthesis, VSE spectroscopy, and computational methods, we provide an intimate view on the origins of these discrepancies. We observe that the electric fields are correlated to conformation-dependent effects of the flexible PenG side chain, including both the local solvation structure and solute conformational sampling in MD. Additionally, we identified that MD-based electric fields are consistently overestimated in three-point water models in the vicinity of charged groups; this cannot be entirely ameliorated using polarizable force fields (AMOEBA) or advanced water models. This work demonstrates the value of the VSE as a direct method for experiment-guided refinements of MD force fields and establishes a general reductionist approach to calibrating vibrational probes for complex (bio)molecules.

g_elpot: A Tool for Quantifying Biomolecular Electrostatics from Molecular Dynamics Trajectories

Electrostatic forces drive a wide variety of biomolecular processes by defining the energetics of the interaction between biomolecules and charged substances. Molecular dynamics (MD) simulations provide trajectories that contain ensembles of structural configurations sampled by biomolecules and their environment. Although this information can be used for high-resolution characterization of biomolecular electrostatics, it has not yet been possible to calculate electrostatic potentials from MD trajectories in a way allowing for quantitative connection to energetics. Here, we present g_elpot, a GROMACS-based tool that utilizes the smooth particle mesh Ewald method to quantify the electrostatics of biomolecules by calculating potential within water molecules that are explicitly present in biomolecular MD simulations. g_elpot can extract the global distribution of the electrostatic potential from MD trajectories and measure its time course in functionally important regions of a biomolecule. To demonstrate that g_elpot can be used to gain biophysical insights into various biomolecular processes, we applied the tool to MD trajectories of the P2X3 receptor, TMEM16 lipid scramblases, the secondary-active transporter Glt_Ph, and DNA complexed with cationic polymers. Our results indicate that g_elpot is well suited for quantifying electrostatics in biomolecular systems to provide a deeper understanding of its role in biomolecular processes.

Predicting Monovalent Ion Correlation Effects in Nucleic Acids

Ion correlation and fluctuation can play a potentially significant role in metal ion-nucleic acid interactions. Previous studies have focused on the effects for multivalent cations. However, the correlation and fluctuation effects can be important also for monovalent cations around the nucleic acid surface. Here, we report a model, gMCTBI, that can explicitly treat discrete distributions of both monovalent and multivalent cations and can account for the correlation and fluctuation effects for the cations in the solution. The gMCTBI model enables investigation of the global ion binding properties as well as the detailed discrete distributions of the bound ions. Accounting for the ion correlation effect for monovalent ions can lead to more accurate predictions, especially in a mixed monovalent and multivalent salt solution, for the number and location of the bound ions. Furthermore, although the monovalent ion-mediated correlation does not show a significant effect on the number of bound ions, the correlation may enhance the accumulation of monovalent ions near the nucleic acid surface and hence affect the ion distribution. The study further reveals novel ion correlation-induced effects in the competition between the different cations around nucleic acids.

Predicting RNA-Metal Ion Binding with Ion Dehydration Effects

Metal ions play essential roles in nucleic acids folding and stability. The interaction between metal ions and nucleic acids can be highly complicated because of the interplay between various effects such as ion correlation, fluctuation, and dehydration. These effects may be particularly important for multivalent ions such as Mg²⁺ ions. Previous efforts to model ion correlation and fluctuation effects led to the development of the Monte Carlo tightly bound ion model. Here, by incorporating ion hydration/dehydration effects into the Monte Carlo tightly bound ion model, we develop a, to our knowledge, new approach to predict ion binding. The new model enables predictions for not only the number of bound ions but also the three-dimensional spatial distribution of the bound ions. Furthermore, the new model reveals several intriguing features for the bound ions such as the mutual enhancement/inhibition in ion binding between the fully hydrated (diffuse) ions, the outer-shell dehydrated ions, and the inner-shell dehydrated ions and novel features for the monovalent-divalent ion interplay due to the hydration effect.

Theory Meets Experiment: Metal Ion Effects in HCV Genomic RNA Kissing Complex Formation

The long-range base pairing between the 5BSL3. 2 and 3′X domains in hepatitis C virus (HCV) genomic RNA is essential for viral replication. Experimental evidence points to the critical role of metal ions, especially Mg²⁺ ions, in the formation of the 5BSL3.2:3′X kissing complex. Furthermore, NMR studies suggested an important ion-dependent conformational switch in the kissing process. However, for a long time, mechanistic understanding of the ion effects for the process has been unclear. Recently, computational modeling based on the Vfold RNA folding model and the partial charge-based tightly bound ion (PCTBI) model, in combination with the NMR data, revealed novel physical insights into the role of metal ions in the 5BSL3.2-3′X system. The use of the PCTBI model, which accounts for the ion correlation and fluctuation, gives reliable predictions for the ion-dependent electrostatic free energy landscape and ion-induced population shift of the 5BSL3.2:3′X kissing complex. Furthermore, the predicted ion binding sites offer insights about how ion-RNA interactions shift the conformational equilibrium. The integrated theory-experiment study shows that Mg²⁺ ions may be essential for HCV viral replication. Moreover, the observed Mg²⁺-dependent conformational equilibrium may be an adaptive property of the HCV genomic RNA such that the equilibrium is optimized to the intracellular Mg²⁺ concentration in liver cells for efficient viral replication.

Predicting Ion Effects in an RNA Conformational Equilibrium

We develop a partial charge-based tightly bound ion (PCTBI) model for the ion effects in RNA folding. On the basis of the Monte Carlo tightly bound ion (MCTBI) approach, the model can account for ion fluctuation and correlation effects, and can predict the ion distribution around the RNA. Furthermore, unlike the previous coarse-grained RNA charge models, where negative charges are placed on the phosphates only, the current new model considers the detailed all-atom partial charge distribution on the RNA. Thus, the model not only keeps the advantage of the MCTBI model, but also has the potential to provide important detailed information unattainable by the previous MCTBI models. For example, the model predicts the reduction in ion binding upon protein binding and ion-induced conformational switches. For hepatitis C virus genomic RNA, the model predicts a Mg²⁺-induced stabilization of a kissing motif for a cis-acting regulatory element in the genomic RNA. Extensive theory-experiment comparisons support the reliability of the theoretical predictions. Therefore, the model may serve as a robust starting point for further development of an accurate method for ion effects in an RNA conformational equilibrium and RNA-cofactor interactions.

MCTBI: a web server for predicting metal ion effects in RNA structures

Metal ions play critical roles in RNA structure and function. However, web servers and software packages for predicting ion effects in RNA structures are notably scarce. Furthermore, the existing web servers and software packages mainly neglect ion correlation and fluctuation effects, which are potentially important for RNAs. We here report a new web server, the MCTBI server (http://rna.physics.missouri.edu/MCTBI), for the prediction of ion effects for RNA structures. This server is based on the recently developed MCTBI, a model that can account for ion correlation and fluctuation effects for nucleic acid structures and can provide improved predictions for the effects of metal ions, especially for multivalent ions such as Mg²⁺ effects, as shown by extensive theory-experiment test results. The MCTBI web server predicts metal ion binding fractions, the most probable bound ion distribution, the electrostatic free energy of the system, and the free energy components. The results provide mechanistic insights into the role of metal ions in RNA structure formation and folding stability, which is important for understanding RNA functions and the rational design of RNA structures.

A New Method to Predict Ion Effects in RNA Folding

The strong interaction between metal ions in solution and highly charged RNA molecules is critical for RNA structure formation and stabilization. Metal ions binding to RNA can induce RNA structural changes that are important for RNA cellular functions. Therefore, quantitative modeling of the ion effects is essential for RNA structure prediction and RNA-based molecular design. Recently, inspired by the increasing experimental evidence that supports the importance of ion correlation and fluctuation in ion-RNA interactions, we developed a new computational model, Monte Carlo Tightly Bound Ion (MCTBI) model. The validity of the model is shown by the improved accuracy in the predictions for ion binding properties and ion-dependent free energies for RNA structures. In this chapter, using homodimeric tetraloop-receptor docking as an illustrative example, we showcase the MCTBI method for the computational prediction of the ion effects in RNA folding.

Accurate estimation of electrostatic binding energy with Poisson-Boltzmann equation solver DelPhi program

Poisson–Boltzmann (PB) model is a widely used implicit solvent approximation in biophysical modeling because of its ability to provide accurate and reliable PB electrostatic salvation free energies ([Formula: see text] as well as electrostatic binding free energy ([Formula: see text] estimations. However, a recent study has warned that the 0.5[Formula: see text]Å grid spacing which is normally adopted can produce unacceptable errors in [Formula: see text] estimation with the solvent excluded surface (SES) (Harris RC, Boschitsch AH and Fenley MO, Influence of grid spacing in Poisson–Boltzmann equation binding energy estimation, J Chem Theory Comput 19: 3677–3685, 2013). In this work, we investigate the grid dependence of the widely used PB solver DelPhi v6.2 with molecular surface (MS) for estimating both electrostatic solvation free energies and electrostatic binding free energies. Our results indicate that, for the molecular complex and components the absolute errors of [Formula: see text] are smaller than that of [Formula: see text], and grid spacing of 0.8[Formula: see text]Å with DelPhi program ensures the accuracy and reliability of [Formula: see text]; however, the accuracy of [Formula: see text] largely relies on the order of magnitude of [Formula: see text] itself rather than that of [Formula: see text] or [Formula: see text]. Our findings suggest that grid spacing of 0.5[Formula: see text]Å is enough to produce accurate [Formula: see text] for molecules whose [Formula: see text] are large, but finer grids are needed when [Formula: see text] is very small.

Electrostatic component of binding energy: Interpreting predictions from poisson-boltzmann equation and modeling protocols

Macromolecular interactions are essential for understanding numerous biological processes and are typically characterized by the binding free energy. Important component of the binding free energy is the electrostatics, which is frequently modeled via the solutions of the Poisson-Boltzmann Equations (PBE). However, numerous works have shown that the electrostatic component (ΔΔGelec ) of binding free energy is very sensitive to the parameters used and modeling protocol. This prompted some researchers to question the robustness of PBE in predicting ΔΔGelec . We argue that the sensitivity of the absolute ΔΔGelec calculated with PBE using different input parameters and definitions does not indicate PBE deficiency, rather this is what should be expected. We show how the apparent sensitivity should be interpreted in terms of the underlying changes in several numerous and physical parameters. We demonstrate that PBE approach is robust within each considered force field (CHARMM-27, AMBER-94, and OPLS-AA) once the corresponding structures are energy minimized. This observation holds despite of using two different molecular surface definitions, pointing again that PBE delivers consistent results within particular force field. The fact that PBE delivered ΔΔGelec values may differ if calculated with different modeling protocols is not a deficiency of PBE, but natural results of the differences of the force field parameters and potential functions for energy minimization. In addition, while the absolute ΔΔGelec values calculated with different force field differ, their ordering remains practically the same allowing for consistent ranking despite of the force field used. © 2016 Wiley Periodicals, Inc.

Monte Carlo Tightly Bound Ion Model: Predicting Ion-Binding Properties of RNA with Ion Correlations and Fluctuations

Experiments have suggested that ion correlation and fluctuation effects can be potentially important for multivalent ions in RNA folding. However, most existing computational methods for the ion electrostatics in RNA folding tend to ignore these effects. The previously reported tightly bound ion (TBI) model can treat ion correlation and fluctuation but its applicability to biologically important RNAs is severely limited by the low computational efficiency. Here, on the basis of Monte Carlo sampling for the many-body ion distribution, we develop a new computational model, the Monte Carlo tightly bound ion (MCTBI) model, for ion-binding properties around an RNA. Because of an enhanced sampling algorithm for ion distribution, the model leads to a significant improvement in computational efficiency. For example, for a 160-nt RNA, the model causes a more than 10-fold increase in the computational efficiency, and the improvement in computational efficiency is more pronounced for larger systems. Furthermore, unlike the earlier model that describes ion distribution using the number of bound ions around each nucleotide, the current MCTBI model is based on the three-dimensional coordinates of the ions. The higher efficiency of the model allows us to treat the ion effects for medium to large RNA molecules, RNA-ligand complexes, and RNA-protein complexes. This new model together with proper RNA conformational sampling and the energetics model may serve as a starting point for further development for the ion effects in RNA folding and conformational changes and for large nucleic acid systems.

MPBEC, a Matlab Program for Biomolecular Electrostatic Calculations

One of the most used and efficient approaches to compute electrostatic properties of biological systems is to numerically solve the Poisson-Boltzmann (PB) equation. There are several software packages available that solve the PB equation for molecules in aqueous electrolyte solutions. Most of these software packages are useful for scientists with specialized training and expertise in computational biophysics. However, the user is usually required to manually take several important choices, depending on the complexity of the biological system, to successfully obtain the numerical solution of the PB equation. This may become an obstacle for researchers, experimentalists, even students with no special training in computational methodologies. Aiming to overcome this limitation, in this article we present MPBEC, a free, cross-platform, open-source software that provides non-experts in the field an easy and efficient way to perform biomolecular electrostatic calculations on single processor computers. MPBEC is a Matlab script based on the Adaptative Poisson Boltzmann Solver, one of the most popular approaches used to solve the PB equation. MPBEC does not require any user programming, text editing or extensive statistical skills, and comes with detailed user-guide documentation. As a unique feature, MPBEC includes a useful graphical user interface (GUI) application which helps and guides users to configure and setup the optimal parameters and approximations to successfully perform the required biomolecular electrostatic calculations. The GUI also incorporates visualization tools to facilitate users pre- and post- analysis of structural and electrical properties of biomolecules.

Cation-Anion Interactions within the Nucleic Acid Ion Atmosphere Revealed by Ion Counting

The ion atmosphere is a critical structural, dynamic, and energetic component of nucleic acids that profoundly affects their interactions with proteins and ligands. Experimental methods that "count" the number of ions thermodynamically associated with the ion atmosphere allow dissection of energetic properties of the ion atmosphere, and thus provide direct comparison to theoretical results. Previous experiments have focused primarily on the cations that are attracted to nucleic acid polyanions, but have also showed that anions are excluded from the ion atmosphere. Herein, we have systematically explored the properties of anion exclusion, testing the zeroth-order model that anions of different identity are equally excluded due to electrostatic repulsion. Using a series of monovalent salts, we find, surprisingly, that the extent of anion exclusion and cation inclusion significantly depends on salt identity. The differences are prominent at higher concentrations and mirror trends in mean activity coefficients of the electrolyte solutions. Salts with lower activity coefficients exhibit greater accumulation of both cations and anions within the ion atmosphere, strongly suggesting that cation-anion correlation effects are present in the ion atmosphere and need to be accounted for to understand electrostatic interactions of nucleic acids. To test whether the effects of cation-anion correlations extend to nucleic acid kinetics and thermodynamics, we followed the folding of P4-P6, a domain of the Tetrahymena group I ribozyme, via single-molecule fluorescence resonance energy transfer in solutions with different salts. Solutions of identical concentration but lower activity gave slower and less favorable folding. Our results reveal hitherto unknown properties of the ion atmosphere and suggest possible roles of oriented ion pairs or anion-bridged cations in the ion atmosphere for electrolyte solutions of salts with reduced activity. Consideration of these new results leads to a reevaluation of the strengths and limitations of Poisson-Boltzmann theory and highlights the need for next-generation atomic-level models of the ion atmosphere.

Ultrafast protein structure-based virtual screening with Panther

Molecular docking is by far the most common method used in protein structure-based virtual screening. This paper presents Panther, a novel ultrafast multipurpose docking tool. In Panther, a simple shape-electrostatic model of the ligand-binding area of the protein is created by utilizing the protein crystal structure. The features of the possible ligands are then compared to the model by using a similarity search algorithm. On average, one ligand can be processed in a few minutes by using classical docking methods, whereas using Panther processing takes <1 s. The presented Panther protocol can be used in several applications, such as speeding up the early phases of drug discovery projects, reducing the number of failures in the clinical phase of the drug development process, and estimating the environmental toxicity of chemicals. Panther-code is available in our web pages (http://www.jyu.fi/panther) free of charge after registration.

Multidimensional persistence in biomolecular data

Persistent homology has emerged as a popular technique for the topological simplification of big data, including biomolecular data. Multidimensional persistence bears considerable promise to bridge the gap between geometry and topology. However, its practical and robust construction has been a challenge. We introduce two families of multidimensional persistence, namely pseudomultidimensional persistence and multiscale multidimensional persistence. The former is generated via the repeated applications of persistent homology filtration to high-dimensional data, such as results from molecular dynamics or partial differential equations. The latter is constructed via isotropic and anisotropic scales that create new simiplicial complexes and associated topological spaces. The utility, robustness, and efficiency of the proposed topological methods are demonstrated via protein folding, protein flexibility analysis, the topological denoising of cryoelectron microscopy data, and the scale dependence of nanoparticles. Topological transition between partial folded and unfolded proteins has been observed in multidimensional persistence. The separation between noise topological signatures and molecular topological fingerprints is achieved by the Laplace-Beltrami flow. The multiscale multidimensional persistent homology reveals relative local features in Betti-0 invariants and the relatively global characteristics of Betti-1 and Betti-2 invariants.

DNA under Force: Mechanics, Electrostatics, and Hydration

Quantifying the basic intra- and inter-molecular forces of DNA has helped us to better understand and further predict the behavior of DNA. Single molecule technique elucidates the mechanics of DNA under applied external forces, sometimes under extreme forces. On the other hand, ensemble studies of DNA molecular force allow us to extend our understanding of DNA molecules under other forces such as electrostatic and hydration forces. Using a variety of techniques, we can have a comprehensive understanding of DNA molecular forces, which is crucial in unraveling the complex DNA functions in living cells as well as in designing a system that utilizes the unique properties of DNA in nanotechnology.

Structural and electrostatic analysis of HLA B-cell epitopes: inference on immunogenicity and prediction of humoral alloresponses

PURPOSE OF REVIEW

The immunogenic capacity of donor human leukocyte antigen (HLA) to induce humoral immune responses is not an intrinsic property of the mismatched alloantigen but depends on the HLA phenotype of the recipient. In recent years, advances in molecular sequence technology and information from X-ray crystallography have enabled structural comparison of donor and recipient HLA type providing an opportunity for a more rational approach for determining HLA compatibility. In this article, we review studies investigating the molecular basis of antibody-antigen interactions and present computational approaches to determine the complex physiochemical and structural properties of B-cell epitopes.

RECENT FINDINGS

The relative immunogenicity of individual HLA mismatches may be predicted from analysis of polymorphic amino acids at continuous and discontinuous HLA sequence positions. The use of alloantigen sequence information alone, however, provides limited insight into key determinants of B-cell epitope immunogenicity, such as the orientation, accessibility and physiochemical properties of amino acid side chains. Advances in computational molecular modelling techniques now enable assessment of HLA-alloantibody interactions at the atomic level. Recent evidence supports a strong link between HLA B-cell epitope surface electrostatic potential and their immunogenicity.

SUMMARY

Assessment of the surface electrostatic properties of HLA alloantigens and computational analyses of HLA-alloantibody interactions represent a promising area for future research into the molecular basis of HLA immunogenicity and antigenicity.

Poisson-Boltzmann versus Size-Modified Poisson-Boltzmann Electrostatics Applied to Lipid Bilayers

Mean-field methods, such as the Poisson–Boltzmann equation (PBE), are often used to calculate the electrostatic properties of molecular systems. In the past two decades, an enhancement of the PBE, the size-modified Poisson–Boltzmann equation (SMPBE), has been reported. Here, the PBE and the SMPBE are reevaluated for realistic molecular systems, namely, lipid bilayers, under eight different sets of input parameters. The SMPBE appears to reproduce the molecular dynamics simulation results better than the PBE only under specific parameter sets, but in general, it performs no better than the Stern layer correction of the PBE. These results emphasize the need for careful discussions of the accuracy of mean-field calculations on realistic systems with respect to the choice of parameters and call for reconsideration of the cost-efficiency and the significance of the current SMPBE formulation.

Protein-protein docking with dynamic residue protonation states

Protein-protein interactions depend on a host of environmental factors. Local pH conditions influence the interactions through the protonation states of the ionizable residues that can change upon binding. In this work, we present a pH-sensitive docking approach, pHDock, that can sample side-chain protonation states of five ionizable residues (Asp, Glu, His, Tyr, Lys) on-the-fly during the docking simulation. pHDock produces successful local docking funnels in approximately half (79/161) the protein complexes, including 19 cases where standard RosettaDock fails. pHDock also performs better than the two control cases comprising docking at pH 7.0 or using fixed, predetermined protonation states. On average, the top-ranked pHDock structures have lower interface RMSDs and recover more native interface residue-residue contacts and hydrogen bonds compared to RosettaDock. Addition of backbone flexibility using a computationally-generated conformational ensemble further improves native contact and hydrogen bond recovery in the top-ranked structures. Although pHDock is designed to improve docking, it also successfully predicts a large pH-dependent binding affinity change in the Fc–FcRn complex, suggesting that it can be exploited to improve affinity predictions. The approaches in the study contribute to the goal of structural simulations of whole-cell protein-protein interactions including all the environmental factors, and they can be further expanded for pH-sensitive protein design.

Protein-protein interactions are fundamental for biological function and are strongly influenced by their local environment. Cellular pH is tightly controlled and is one of the critical environmental factors that regulates protein-protein interactions. Three-dimensional structures of the protein complexes can help us understand the mechanism of the interactions. Since experimental determination of the structures of protein-protein complexes is expensive and time-consuming, computational docking algorithms are helpful to predict the structures. However, none of the current protein-protein docking algorithms account for the critical environmental pH effects. So we developed a pH-sensitive docking algorithm that can dynamically pick the favorable protonation states of the ionizable amino-acid residues. Compared to our previous standard docking algorithm, the new algorithm improves docking accuracy and generates higher-quality predictions over a large dataset of protein-protein complexes. We also use a case study to demonstrate efficacy of the algorithm in predicting a large pH-dependent binding affinity change that cannot be captured by the other methods that neglect pH effects. In principle, the approaches in the study can be used for rational design of pH-dependent protein inhibitors or industrial enzymes that are active over a wide range of pH values.

Performance comparison of Poisson–Boltzmann equation solvers DelPhi and PBSA in calculation of electrostatic solvation energies

Many Poisson–Boltzmann equation (PBE) solvers have been developed to calculate electrostatic energy of biomolecules, and each PBE solver has its advantages. In this study two PBE solvers — DelPhi and PBSA module in AMBER — are compared in terms of calculation results, convergence and program-running time by calculating the electrostatic solvation energy for a test set composed of 4 Kirkwood models and another test set of 25 protein structures. The protein structures were pretreated by AMBER and AMBER99SB force field parameters were used in all calculations. It is found that (i) At fine grids, both PBE solvers can produce accurate results on test set 1 and consistent results with each other on test set 2, with differences between each other varying from several to ~ 10 kcal/mol. (ii) Under convergence criterion "absolute value of relative error is less than or equal to 2% (|RE| ≤ 2%)", both PBE solvers need very fine grids to produce convergent results on small and complex Kirkwood models, while grid spacings of ≤ 0.5 Å–0.6 Å are well enough for them to achieve good convergent results on various molecular structures. We recommend users to adopt such grid spacing in using of these PBE solvers so as to get good enough convergent results. (iii) In terms of time consumption, DelPhi appears to be more time-saving than PBSA. In summary, according to our comparison, DelPhi and PBSA are paralleled good PBE solvers. The convergence of PBSA is a little better than DelPhi, while DelPhi exceeds PBSA in running speed.

Ion competition in condensed DNA arrays in the attractive regime

Physical origin of DNA condensation by multivalent cations remains unsettled. Here, we report quantitative studies of how one DNA-condensing ion (Cobalt(3+) Hexammine, or Co(3+)Hex) and one nonDNA-condensing ion (Mg(2+)) compete within the interstitial space in spontaneously condensed DNA arrays. As the ion concentrations in the bath solution are systematically varied, the ion contents and DNA-DNA spacings of the DNA arrays are determined by atomic emission spectroscopy and x-ray diffraction, respectively. To gain quantitative insights, we first compare the experimentally determined ion contents with predictions from exact numerical calculations based on nonlinear Poisson-Boltzmann equations. Such calculations are shown to significantly underestimate the number of Co(3+)Hex ions, consistent with the deficiencies of nonlinear Poisson-Boltzmann approaches in describing multivalent cations. Upon increasing the concentration of Mg(2+), the Co(3+)Hex-condensed DNA array expands and eventually redissolves as a result of ion competition weakening DNA-DNA attraction. Although the DNA-DNA spacing depends on both Mg(2+) and Co(3+)Hex concentrations in the bath solution, it is observed that the spacing is largely determined by a single parameter of the DNA array, the fraction of DNA charges neutralized by Co(3+)Hex. It is also observed that only ∼20% DNA charge neutralization by Co(3+)Hex is necessary for spontaneous DNA condensation. We then show that the bath ion conditions can be reduced to one variable with a simplistic ion binding model, which is able to describe the variations of both ion contents and DNA-DNA spacings reasonably well. Finally, we discuss the implications on the nature of interstitial ions and cation-mediated DNA-DNA interactions.

A fast alternating direction implicit algorithm for geometric flow equations in biomolecular surface generation

In this paper, a new alternating direction implicit (ADI) method is introduced to solve potential driven geometric flow PDEs for biomolecular surface generation. For such PDEs, an extra factor is usually added to stabilize the explicit time integration. However, two existing implicit ADI schemes are also based on the scaled form, which involves nonlinear cross derivative terms that have to be evaluated explicitly. This affects the stability and accuracy of these ADI schemes. To overcome these difficulties, we propose a new ADI algorithm based on the unscaled form so that cross derivatives are not involved. Central finite differences are employed to discretize the nonhomogenous diffusion process of the geometric flow. The proposed ADI algorithm is validated through benchmark examples with analytical solutions, reference solutions, or literature results. Moreover, quantitative indicators of a biomolecular surface, including surface area, surface-enclosed volume, and solvation free energy, are analyzed for various proteins. The proposed ADI method is found to be unconditionally stable and more accurate than the existing ADI schemes in all tests. This enables the use of a large time increment in the steady state simulation so that the proposed ADI algorithm is very efficient for biomolecular surface generation.

Understanding nucleic acid-ion interactions

Ions surround nucleic acids in what is referred to as an ion atmosphere. As a result, the folding and dynamics of RNA and DNA and their complexes with proteins and with each other cannot be understood without a reasonably sophisticated appreciation of these ions' electrostatic interactions. However, the underlying behavior of the ion atmosphere follows physical rules that are distinct from the rules of site binding that biochemists are most familiar and comfortable with. The main goal of this review is to familiarize nucleic acid experimentalists with the physical concepts that underlie nucleic acid-ion interactions. Throughout, we provide practical strategies for interpreting and analyzing nucleic acid experiments that avoid pitfalls from oversimplified or incorrect models. We briefly review the status of theories that predict or simulate nucleic acid-ion interactions and experiments that test these theories. Finally, we describe opportunities for going beyond phenomenological fits to a next-generation, truly predictive understanding of nucleic acid-ion interactions.

Tertiary structure-based analysis of microRNA-target interactions

Current computational analysis of microRNA interactions is based largely on primary and secondary structure analysis. Computationally efficient tertiary structure-based methods are needed to enable more realistic modeling of the molecular interactions underlying miRNA-mediated translational repression. We incorporate algorithms for predicting duplex RNA structures, ionic strength effects, duplex entropy and free energy, and docking of duplex-Argonaute protein complexes into a pipeline to model and predict miRNA-target duplex binding energies. To ensure modeling accuracy and computational efficiency, we use an all-atom description of RNA and a continuum description of ionic interactions using the Poisson-Boltzmann equation. Our method predicts the conformations of two constructs of Caenorhabditis elegans let-7 miRNA-target duplexes to an accuracy of ∼3.8 Å root mean square distance of their NMR structures. We also show that the computed duplex formation enthalpies, entropies, and free energies for eight miRNA-target duplexes agree with titration calorimetry data. Analysis of duplex-Argonaute docking shows that structural distortions arising from single-base-pair mismatches in the seed region influence the activity of the complex by destabilizing both duplex hybridization and its association with Argonaute. Collectively, these results demonstrate that tertiary structure-based modeling of miRNA interactions can reveal structural mechanisms not accessible with current secondary structure-based methods.

Ion-mediated RNA structural collapse: effect of spatial confinement

RNAs are negatively charged molecules that reside in cellular environments with macromolecular crowding. Macromolecular confinement can influence the ion effects in RNA folding. In this work, using the recently developed tightly bound ion model for ion fluctuation and correlation, we investigate the effect of confinement on ion-mediated RNA structural collapse for a simple model system. We find that for both Na(+) and Mg(2+), the ion efficiencies in mediating structural collapse/folding are significantly enhanced by the structural confinement. This enhancement of ion efficiency is attributed to the decreased electrostatic free-energy difference between the compact conformation ensemble and the (restricted) extended conformation ensemble due to the spatial restriction.

Rapid calculation of protein pKa values using Rosetta

We developed a Rosetta-based Monte Carlo method to calculate the pK(a) values of protein residues that commonly exhibit variable protonation states (Asp, Glu, Lys, His, and Tyr). We tested the technique by calculating pK(a) values for 264 residues from 34 proteins. The standard Rosetta score function, which is independent of any environmental conditions, failed to capture pK(a) shifts. After incorporating a Coulomb electrostatic potential and optimizing the solvation reference energies for pK(a) calculations, we employed a method that allowed side-chain flexibility and achieved a root mean-square deviation (RMSD) of 0.83 from experimental values (0.68 after discounting 11 predictions with an error over 2 pH units). Additional degrees of side-chain conformational freedom for the proximal residues facilitated the capture of charge-charge interactions in a few cases, resulting in an overall RMSD of 0.85 pH units. The addition of backbone flexibility increased the overall RMSD to 0.93 pH units but improved relative pK(a) predictions for proximal catalytic residues. The method also captures large pK(a) shifts of lysine and some glutamate point mutations in staphylococcal nuclease. Thus, a simple and fast method based on the Rosetta score function and limited conformational sampling produces pK(a) values that will be useful when rapid estimation is essential, such as in docking, design, and folding.

Poisson-Boltzmann Calculations: van der Waals or Molecular Surface?

The Poisson-Boltzmann equation is widely used for modeling the electrostatics of biomolecules, but the calculation results are sensitive to the choice of the boundary between the low solute dielectric and the high solvent dielectric. The default choice for the dielectric boundary has been the molecular surface, but the use of the van der Waals surface has also been advocated. Here we review recent studies in which the two choices are tested against experimental results and explicit-solvent calculations. The assignment of the solvent high dielectric constant to interstitial voids in the solute is often used as a criticism against the van der Waals surface. However, this assignment may not be as unrealistic as previously thought, since hydrogen exchange and other NMR experiments have firmly established that all interior parts of proteins are transiently accessible to the solvent.

iAPBS: a programming interface to Adaptive Poisson-Boltzmann Solver (APBS)

The Adaptive Poisson-Boltzmann Solver (APBS) is a state-of-the-art suite for performing Poisson-Boltzmann electrostatic calculations on biomolecules. The iAPBS package provides a modular programmatic interface to the APBS library of electrostatic calculation routines. The iAPBS interface library can be linked with a FORTRAN or C/C++ program thus making all of the APBS functionality available from within the application. Several application modules for popular molecular dynamics simulation packages - Amber, NAMD and CHARMM are distributed with iAPBS allowing users of these packages to perform implicit solvent electrostatic calculations with APBS.

Dissipative electro-elastic network model of protein electrostatics

We propose a dissipative electro-elastic network model to describe the dynamics and statistics of electrostatic fluctuations at active sites of proteins. The model combines the harmonic network of residue beads with overdamped dynamics of the normal modes of the network characterized by two friction coefficients. The electrostatic component is introduced to the model through atomic charges of the protein force field. The overall effect of the electrostatic fluctuations of the network is recorded through the frequency-dependent response functions of the electrostatic potential and electric field at the protein active site. We also consider the dynamics of displacements of individual residues in the network and the dynamics of distances between pairs of residues. The model is tested against loss spectra of residue displacements and the electrostatic potential and electric field at the heme's iron from all-atom molecular dynamics simulations of three hydrated globular proteins.

Predicting ion-nucleic acid interactions by energy landscape-guided sampling

The recently developed Tightly Bound Ion (TBI) model offers improved predictions for ion effect in nucleic acid systems by accounting for ion correlation and fluctuation effects. However, further application of the model to larger systems is limited by the low computational efficiency of the model. Here, we develop a new computational efficient TBI model using free energy landscape-guided sampling method. The method leads to drastic reduction in the computer time by a factor of 50 for RNAs of 50-100 nucleotides long. The improvement in the computational efficiency would be more significant for larger structures. To test the new method, we apply the model to predict the free energies and the number of bound ions for a series of RNA folding systems. The validity of this new model is supported by the nearly exact agreement with the results from the original TBI model and the agreement with the experimental data. The method may pave the way for further applications of the TBI model to treat a broad range of biologically significant systems such as tetraloop-receptor and riboswitches.

Differential geometry based solvation model II: Lagrangian formulation

Solvation is an elementary process in nature and is of paramount importance to more sophisticated chemical, biological and biomolecular processes. The understanding of solvation is an essential prerequisite for the quantitative description and analysis of biomolecular systems. This work presents a Lagrangian formulation of our differential geometry based solvation models. The Lagrangian representation of biomolecular surfaces has a few utilities/advantages. First, it provides an essential basis for biomolecular visualization, surface electrostatic potential map and visual perception of biomolecules. Additionally, it is consistent with the conventional setting of implicit solvent theories and thus, many existing theoretical algorithms and computational software packages can be directly employed. Finally, the Lagrangian representation does not need to resort to artificially enlarged van der Waals radii as often required by the Eulerian representation in solvation analysis. The main goal of the present work is to analyze the connection, similarity and difference between the Eulerian and Lagrangian formalisms of the solvation model. Such analysis is important to the understanding of the differential geometry based solvation model. The present model extends the scaled particle theory of nonpolar solvation model with a solvent-solute interaction potential. The nonpolar solvation model is completed with a Poisson-Boltzmann (PB) theory based polar solvation model. The differential geometry theory of surfaces is employed to provide a natural description of solvent-solute interfaces. The optimization of the total free energy functional, which encompasses the polar and nonpolar contributions, leads to coupled potential driven geometric flow and PB equations. Due to the development of singularities and nonsmooth manifolds in the Lagrangian representation, the resulting potential-driven geometric flow equation is embedded into the Eulerian representation for the purpose of computation, thanks to the equivalence of the Laplace-Beltrami operator in the two representations. The coupled partial differential equations (PDEs) are solved with an iterative procedure to reach a steady state, which delivers desired solvent-solute interface and electrostatic potential for problems of interest. These quantities are utilized to evaluate the solvation free energies and protein-protein binding affinities. A number of computational methods and algorithms are described for the interconversion of Lagrangian and Eulerian representations, and for the solution of the coupled PDE system. The proposed approaches have been extensively validated. We also verify that the mean curvature flow indeed gives rise to the minimal molecular surface and the proposed variational procedure indeed offers minimal total free energy. Solvation analysis and applications are considered for a set of 17 small compounds and a set of 23 proteins. The salt effect on protein-protein binding affinity is investigated with two protein complexes by using the present model. Numerical results are compared to the experimental measurements and to those obtained by using other theoretical methods in the literature.

Differential geometry based solvation model II: Lagrangian formulation

Solvation is an elementary process in nature and is of paramount importance to more sophisticated chemical, biological and biomolecular processes. The understanding of solvation is an essential prerequisite for the quantitative description and analysis of biomolecular systems. This work presents a Lagrangian formulation of our differential geometry based solvation models. The Lagrangian representation of biomolecular surfaces has a few utilities/advantages. First, it provides an essential basis for biomolecular visualization, surface electrostatic potential map and visual perception of biomolecules. Additionally, it is consistent with the conventional setting of implicit solvent theories and thus, many existing theoretical algorithms and computational software packages can be directly employed. Finally, the Lagrangian representation does not need to resort to artificially enlarged van der Waals radii as often required by the Eulerian representation in solvation analysis. The main goal of the present work is to analyze the connection, similarity and difference between the Eulerian and Lagrangian formalisms of the solvation model. Such analysis is important to the understanding of the differential geometry based solvation model. The present model extends the scaled particle theory of nonpolar solvation model with a solvent-solute interaction potential. The nonpolar solvation model is completed with a Poisson-Boltzmann (PB) theory based polar solvation model. The differential geometry theory of surfaces is employed to provide a natural description of solvent-solute interfaces. The optimization of the total free energy functional, which encompasses the polar and nonpolar contributions, leads to coupled potential driven geometric flow and PB equations. Due to the development of singularities and nonsmooth manifolds in the Lagrangian representation, the resulting potential-driven geometric flow equation is embedded into the Eulerian representation for the purpose of computation, thanks to the equivalence of the Laplace-Beltrami operator in the two representations. The coupled partial differential equations (PDEs) are solved with an iterative procedure to reach a steady state, which delivers desired solvent-solute interface and electrostatic potential for problems of interest. These quantities are utilized to evaluate the solvation free energies and protein-protein binding affinities. A number of computational methods and algorithms are described for the interconversion of Lagrangian and Eulerian representations, and for the solution of the coupled PDE system. The proposed approaches have been extensively validated. We also verify that the mean curvature flow indeed gives rise to the minimal molecular surface and the proposed variational procedure indeed offers minimal total free energy. Solvation analysis and applications are considered for a set of 17 small compounds and a set of 23 proteins. The salt effect on protein-protein binding affinity is investigated with two protein complexes by using the present model. Numerical results are compared to the experimental measurements and to those obtained by using other theoretical methods in the literature.

Structure determination in "shiftless" solid state NMR of oriented protein samples

An efficient formalism for calculating protein structures from oriented-sample NMR data in the torsion-angle space is presented. Angular anisotropies of the NMR observables are treated by utilizing an irreducible spherical basis of rotations. An intermediate rotational transformation is introduced that greatly speeds up structural fitting by rendering the dependence on the torsion angles Φ and Ψ in a purely diagonal form. Back-calculation of the simulated solid-state NMR spectra of protein G involving 15N chemical shift anisotropy (CSA), and 1H-15N and 1Hα-13Cα dipolar couplings was performed by taking into account non-planarity of the peptide linkages and experimental uncertainty. Even a relatively small (to within 1 ppm) random variation in the CSA values arising from uncertainties in the tensor parameters yields the RMSD's of the back-calculated structures of more than 10 Å. Therefore, the 15N CSA has been substituted with heteronuclear dipolar couplings which are derived from the highly conserved bond lengths and bond angles associated with the amino-acid covalent geometry. Using the additional 13Cα-15N and 13C'-15N dipolar couplings makes it possible to calculate protein structures entirely from "shiftless" solid-state NMR data. With the simulated "experimental" uncertainty of 15 Hz for protein G and 120 Hz for a helical hairpin derived from bacteriorhodopsin, back-calculation of the synthetic dipolar NMR spectra yielded a converged set of solutions. The use of distance restraints dramatically improves structural convergence even if larger experimental uncertainties are assumed.

Web servers and services for electrostatics calculations with APBS and PDB2PQR

APBS and PDB2PQR are widely utilized free software packages for biomolecular electrostatics calculations. Using the Opal toolkit, we have developed a Web services framework for these software packages that enables the use of APBS and PDB2PQR by users who do not have local access to the necessary amount of computational capabilities. This not only increases accessibility of the software to a wider range of scientists, educators, and students but also increases the availability of electrostatics calculations on portable computing platforms. Users can access this new functionality in two ways. First, an Opal-enabled version of APBS is provided in current distributions, available freely on the web. Second, we have extended the PDB2PQR web server to provide an interface for the setup, execution, and visualization of electrostatic potentials as calculated by APBS. This web interface also uses the Opal framework which ensures the scalability needed to support the large APBS user community. Both of these resources are available from the APBS/PDB2PQR website: http://www.poissonboltzmann.org/.

Electrostatic clustering and free energy calculations provide a foundation for protein design and optimization

Electrostatic interactions are ubiquitous in proteins and dictate stability and function. In this review, we discuss several methods for the analysis of electrostatics in protein–protein interactions. We discuss alanine-scanning mutagenesis, Poisson–Boltzmann electrostatics, free energy calculations, electrostatic similarity distances, and hierarchical clustering of electrostatic potentials. Our recently developed computational framework, known as Analysis of Electrostatic Similarities Of Proteins (AESOP), incorporates these tools to efficiently elucidate the role of electrostatic potentials in protein interactions. We present the application of AESOP to several proteins and protein complexes, for which charge is purported to facilitate protein association. Specifically, we illustrate how recent work has shaped the formulation of electrostatic calculations, the correlation of electrostatic free energies and electrostatic potential clustering results with experimental binding and activity data, the pH dependence of protein stability and association, the design of mutant proteins with enhanced immunological activity, and how AESOP can expose deficiencies in structural models and experimental data. This integrative approach can be utilized to develop mechanistic models and to guide experimental studies by predicting mutations with desired physicochemical properties and function. Alteration of the electrostatic properties of proteins offers a basis for the design of proteins with optimized binding and activity.

Differential geometry based solvation model I: Eulerian formulation

This paper presents a differential geometry based model for the analysis and computation of the equilibrium property of solvation. Differential geometry theory of surfaces is utilized to define and construct smooth interfaces with good stability and differentiability for use in characterizing the solvent-solute boundaries and in generating continuous dielectric functions across the computational domain. A total free energy functional is constructed to couple polar and nonpolar contributions to the salvation process. Geometric measure theory is employed to rigorously convert a Lagrangian formulation of the surface energy into an Eulerian formulation so as to bring all energy terms into an equal footing. By minimizing the total free energy functional, we derive coupled generalized Poisson-Boltzmann equation (GPBE) and generalized geometric flow equation (GGFE) for the electrostatic potential and the construction of realistic solvent-solute boundaries, respectively. By solving the coupled GPBE and GGFE, we obtain the electrostatic potential, the solvent-solute boundary profile, and the smooth dielectric function, and thereby improve the accuracy and stability of implicit solvation calculations. We also design efficient second order numerical schemes for the solution of the GPBE and GGFE. Matrix resulted from the discretization of the GPBE is accelerated with appropriate preconditioners. An alternative direct implicit (ADI) scheme is designed to improve the stability of solving the GGFE. Two iterative approaches are designed to solve the coupled system of nonlinear partial differential equations. Extensive numerical experiments are designed to validate the present theoretical model, test computational methods, and optimize numerical algorithms. Example solvation analysis of both small compounds and proteins are carried out to further demonstrate the accuracy, stability, efficiency and robustness of the present new model and numerical approaches. Comparison is given to both experimental and theoretical results in the literature.

Binding thermodynamics of phosphorylated inhibitors to triosephosphate isomerase and the contribution of electrostatic interactions

Electrostatic interactions have a central role in some biological processes, such as recognition of charged ligands by proteins. We characterized the binding energetics of yeast triosephosphate isomerase (TIM) with phosphorylated inhibitors 2-phosphoglycollate (2PG) and phosphoglycolohydroxamate (PGH). We determined the thermodynamic parameters of the binding process (K(b), ΔG(b), ΔH(b), ΔS(b) and ΔC(p)) with different concentrations of NaCl, using fluorimetric and calorimetric titrations in the conventional mode of ITC and a novel method, multithermal titration calorimetry (MTC), which enabled us to measure ΔC(p) in a single experiment. We ruled out specific interactions of Na(+) and Cl(-) with the native enzyme and did not detect significant linked protonation effects upon the binding of inhibitors. Increasing ionic strength (I) caused K(b), ΔG(b) and ΔH(b) to become less favorable, while ΔS(b) became less unfavorable. From the variation of K(b) with I, we determined the electrostatic contribution of TIM-2PG and TIM-PGH to ΔG(b) at I=0.06 M and 25 °C to be 36% and 26%, respectively. The greater affinity of PGH for TIM is due to a more favorable ΔH(b) compared to 2PG (by 19-24 kJ mol(-1) at 25 °C). This difference is compatible with PGH establishing up to five more hydrogen bonds with TIM. Both binding ΔC(p)s were negative, and less negative with increasing ionic strength. ΔC(p)s at I=0.06 M were much more negative than predicted by surface area models. Water molecules trapped in the interface when ligands bind to protein could explain the highly negative ΔCps. Thermodynamic binding functions for TIM-2PG changed more with ionic strength than those for TIM-PGH. This greater dependence is consistent with linked, but compensated, protonation equilibriums yielding the dianionic species of 2PG that binds to TIM, process that is not required for PGH.

Differential geometry based multiscale models

Large chemical and biological systems such as fuel cells, ion channels, molecular motors, and viruses are of great importance to the scientific community and public health. Typically, these complex systems in conjunction with their aquatic environment pose a fabulous challenge to theoretical description, simulation, and prediction. In this work, we propose a differential geometry based multiscale paradigm to model complex macromolecular systems, and to put macroscopic and microscopic descriptions on an equal footing. In our approach, the differential geometry theory of surfaces and geometric measure theory are employed as a natural means to couple the macroscopic continuum mechanical description of the aquatic environment with the microscopic discrete atomistic description of the macromolecule. Multiscale free energy functionals, or multiscale action functionals are constructed as a unified framework to derive the governing equations for the dynamics of different scales and different descriptions. Two types of aqueous macromolecular complexes, ones that are near equilibrium and others that are far from equilibrium, are considered in our formulations. We show that generalized Navier-Stokes equations for the fluid dynamics, generalized Poisson equations or generalized Poisson-Boltzmann equations for electrostatic interactions, and Newton's equation for the molecular dynamics can be derived by the least action principle. These equations are coupled through the continuum-discrete interface whose dynamics is governed by potential driven geometric flows. Comparison is given to classical descriptions of the fluid and electrostatic interactions without geometric flow based micro-macro interfaces. The detailed balance of forces is emphasized in the present work. We further extend the proposed multiscale paradigm to micro-macro analysis of electrohydrodynamics, electrophoresis, fuel cells, and ion channels. We derive generalized Poisson-Nernst-Planck equations that are coupled to generalized Navier-Stokes equations for fluid dynamics, Newton's equation for molecular dynamics, and potential and surface driving geometric flows for the micro-macro interface. For excessively large aqueous macromolecular complexes in chemistry and biology, we further develop differential geometry based multiscale fluid-electro-elastic models to replace the expensive molecular dynamics description with an alternative elasticity formulation.

Coarse-grained ions without charges: reproducing the solvation structure of NaCl in water using short-ranged potentials

A coarse-grained model of NaCl in water is presented where the ions are modeled without charge to avoid computationally challenging electrostatics. A monatomic model of water [V. Molinero and E. B. Moore, J. Phys. Chem. B 113, 4008 (2009)] is used as the basis for this coarse-grain approach. The ability of Na(+) to disrupt the native tetrahedral arrangement of water molecules, and of Cl(-) to integrate within this organization, is preserved in this mW-ion model through parametrization focused on water's solvation of these ions. This model successfully reproduces the structural effect of ions on water, referenced to observations from experiments and atomistic molecular dynamics simulations, while using extremely short-ranged potentials. Without Coulomb interactions the model replicates details of the ion-water structure such as distinguishing contact and solvent-separated ion pairs and the free energy barriers between them. The approach of mimicking ionic effects with short-ranged interactions results in performance gains of two orders of magnitude compared to Ewald methods. Explored over a broad range of salt concentration, the model reproduces the solvation structure and trends of diffusion relative to atomistic simulations and experimental results. The functional form of the mW-ion model can be parametrized to represent other electrolytes. With increased computational efficiency and reliable structural fidelity, this model promises to be an asset for accessing significantly longer simulation time scales with an explicit solvent in a coarse-grained system involving, for example, polyelectrolytes such as proteins, nucleic acids, and fuel-cell membranes.