Implicit Solvent Models in Molecular Dynamics Simulations: A Brief Overview

Predicting solvation free energies with an implicit solvent machine learning potential

Machine learning (ML) potentials are a powerful tool in molecular modeling, enabling ab initio accuracy for comparably small computational costs. Nevertheless, all-atom simulations employing best-performing graph neural network architectures are still too expensive for applications requiring extensive sampling, such as free energy computations. Implicit solvent models could provide the necessary speed-up due to reduced degrees of freedom and faster dynamics. Here, we introduce a Solvation Free Energy Path Reweighting (ReSolv) framework to parameterize an implicit solvent ML potential for small organic molecules that accurately predicts the hydration free energy, an essential parameter in drug design and pollutant modeling. Learning on a combination of experimental hydration free energy data and ab initio data of molecules in vacuum, ReSolv bypasses the need for intractable ab initio data of molecules in an explicit bulk solvent and does not have to resort to less accurate data-generating models. On the FreeSolv dataset, ReSolv achieves a mean absolute error close to average experimental uncertainty, significantly outperforming standard explicit solvent force fields. Compared to the explicit solvent ML potential, ReSolv offers a computational speedup of four orders of magnitude and attains closer agreement with experiments. The presented framework paves the way for deep molecular models that are more accurate yet computationally more cost-effective than classical atomistic models.

Advanced computational approaches to understand protein aggregation

Protein aggregation is a widespread phenomenon implicated in debilitating diseases like Alzheimer's, Parkinson's, and cataracts, presenting complex hurdles for the field of molecular biology. In this review, we explore the evolving realm of computational methods and bioinformatics tools that have revolutionized our comprehension of protein aggregation. Beginning with a discussion of the multifaceted challenges associated with understanding this process and emphasizing the critical need for precise predictive tools, we highlight how computational techniques have become indispensable for understanding protein aggregation. We focus on molecular simulations, notably molecular dynamics (MD) simulations, spanning from atomistic to coarse-grained levels, which have emerged as pivotal tools in unraveling the complex dynamics governing protein aggregation in diseases such as cataracts, Alzheimer's, and Parkinson's. MD simulations provide microscopic insights into protein interactions and the subtleties of aggregation pathways, with advanced techniques like replica exchange molecular dynamics, Metadynamics (MetaD), and umbrella sampling enhancing our understanding by probing intricate energy landscapes and transition states. We delve into specific applications of MD simulations, elucidating the chaperone mechanism underlying cataract formation using Markov state modeling and the intricate pathways and interactions driving the toxic aggregate formation in Alzheimer's and Parkinson's disease. Transitioning we highlight how computational techniques, including bioinformatics, sequence analysis, structural data, machine learning algorithms, and artificial intelligence have become indispensable for predicting protein aggregation propensity and locating aggregation-prone regions within protein sequences. Throughout our exploration, we underscore the symbiotic relationship between computational approaches and empirical data, which has paved the way for potential therapeutic strategies against protein aggregation-related diseases. In conclusion, this review offers a comprehensive overview of advanced computational methodologies and bioinformatics tools that have catalyzed breakthroughs in unraveling the molecular basis of protein aggregation, with significant implications for clinical interventions, standing at the intersection of computational biology and experimental research.

DFT mechanistic study of the chemical fixation of CO2 by aziridine derivatives

Using density functional theory (DFT), we treat the reaction of coupling of CO2 with aziridine in gas phase, in the presence of water and of a green catalyst (NaBr). Computations show that, in gas phase, this ring-opening conversions to oxazolidinones initiates by coordinating a CO2 molecule to the nitrogen atom of the aziridine. Then, a nucleophilic interaction between one oxygen atom of the coordinated CO2 and the carbon atom of the aziridine occurs. For methyl substituted aziridine, two pathways are proposed leading either to 4-oxazolidinone or to 5-oxazolidinone. Besides, we show that the activation energy of this reaction reduces in aqueous solution, in the presence of a water molecule explicitly or NaBr catalyst. In addition, the corresponding reaction mechanisms and regioselectivity associated with this ring-opening conversions to oxazolidinones, in the presence of carbon dioxide are found to be influenced by solvent and catalyst. The present findings should allow better designing regioisomer oxazolidinones relevant for organic chemistry, medicinal and pharmacological applications.

Molecular Dynamics Simulations for Rationalizing Polymer Bioconjugation Strategies: Challenges, Recent Developments, and Future Opportunities

The covalent modification of proteins with polymers is a well-established method for improving the pharmacokinetic properties of therapeutically valuable biologics. The conjugated polymer chains of the resulting hybrid represent highly flexible macromolecular structures. As the dynamics of such systems remain rather elusive for established experimental techniques from the field of protein structure elucidation, molecular dynamics simulations have proven as a valuable tool for studying such conjugates at an atomistic level, thereby complementing experimental studies. With a focus on new developments, this review aims to provide researchers from the polymer bioconjugation field with a concise and up to date overview of such approaches. After introducing basic principles of molecular dynamics simulations, as well as methods for and potential pitfalls in modeling bioconjugates, the review illustrates how these computational techniques have contributed to the understanding of bioconjugates and bioconjugation strategies in the recent past and how they may lead to a more rational design of novel bioconjugates in the future.

Physics-Guided Deep Generative Model for New Ligand Discovery

Structure-based drug discovery aims to identify small molecules that can attach to a specific target protein and change its functionality. Recently, deep learning has shown great promise in generating drug-like molecules with specific biochemical features and conditioned with structural features. However, they usually fail to incorporate an essential factor: the underlying physics which guides molecular formation and binding in real-world scenarios. In this work, we describe a physics-guided deep generative model for new ligand discovery, conditioned not only on the binding site but also on physics-based features that describe the binding mechanism between a receptor and a ligand. The proposed hybrid model has been tested on large protein-ligand complexes and small host-guest systems. Using the top-N methodology, on average more than 75% of the generated structures by our hybrid model were stronger binders than the original reference ligand. All of them had higher ΔGbind (affinity) values than the ones generated by the previous state-of-the-art method by an average margin of 1.88 kcal/mol. The visualization of the top-5 ligands generated by the proposed physics-guided model and the reference deep learning model demonstrate more feasible conformations and orientations by the former. The future directions include training and testing the hybrid model on larger datasets, adding more relevant physics-based features, and interpreting the deep learning outcomes from biophysical perspectives.

Molecular Dynamics Simulation of Protein Cages

Molecular dynamics (MD) simulations enable the description of the physical movement of the system over time based on classical mechanics at various scales depending on the models. Protein cages are a particular group of different-size proteins with hollow, spherical structures and are widely found in nature, which have vast applications in numerous fields. The MD simulation of cage proteins is particularly important as a powerful tool to unveil their structures and dynamics for various properties, assembly behavior, and molecular transport mechanisms. Here, we describe how to conduct MD simulations for cage proteins, especially technical details, and analyze some of the properties of interest using GROMACS/NAMD packages.

Probing the effect of NEK7 and cofactor interactions on dynamics of NLRP3 monomer using molecular simulation

The NLRP3 inflammasome is a cytoplasmic complex that regulates the activation of inflammatory cytokines and, given its implication in a range of diseases, is an important therapeutic target. The cofactor ATP and the centrosomal kinase NEK7 are important for NLRP3 activation. Here we have constructed and simulated computational models of full‐length monomeric NLRP3 to shed light on the importance of NEK7 and cofactor interactions for its conformation and dynamics in aqueous solution. We find that molecular dynamics simulation reproduces well the features of the recently published cryo‐EM structure of the ADP‐bound NLRP3–NEK7 complex; on the removal of NEK7, the NLRP3 molecule adopts a more compact closed form during simulations. Replacement of ADP by ATP promotes a rearrangement of hydrogen‐bonding interactions, domain interfaces, and a degree of opening of the NLRP3 conformation. We also examine the dynamics of an acidic loop of the LRR domain of NLRP3, which samples in a region observed in the NEK7‐bound cryo‐EM structure but not in an oligomeric form of inactive NLRP3. During the molecular dynamics simulations of NLRP3, we find some plasticity in its topology that suggests access routes for ATP to the cofactor pocket not immediately evident from the existing NEK7‐bound cryo‐EM structure. These computed dynamical trajectories of NLRP3 provide insight into coordinates of deformation that may be key for cofactor binding and inflammasome activation.

Modeling electrokinetic flows with the discrete ion stochastic continuum overdamped solvent algorithm

In this article we develop an algorithm for the efficient simulation of electrolytes in the presence of physical boundaries. In previous work the discrete ion stochastic continuum overdamped solvent (DISCOS) algorithm was derived for triply periodic domains, and was validated through ion-ion pair correlation functions and Debye-Hückel-Onsager theory for conductivity, including the Wien effect for strong electric fields. In extending this approach to include an accurate treatment of physical boundaries we must address several important issues. First, the modifications to the spreading and interpolation operators necessary to incorporate interactions of the ions with the boundary are described. Next we discuss the modifications to the electrostatic solver to handle the influence of charges near either a fixed potential or dielectric boundary. An additional short-ranged potential is also introduced to represent interaction of the ions with a solid wall. Finally, the dry diffusion term is modified to account for the reduced mobility of ions near a boundary, which introduces an additional stochastic drift correction. Several validation tests are presented confirming the correct equilibrium distribution of ions in a channel. Additionally, the methodology is demonstrated using electro-osmosis and induced-charge electro-osmosis, with comparison made to theory and other numerical methods. Notably, the DISCOS approach achieves greater accuracy than a continuum electrostatic simulation method. We also examine the effect of under-resolving hydrodynamic effects using a "dry diffusion" approach, and find that considerable computational speedup can be achieved with a negligible impact on accuracy.

A Physics-Guided Neural Network for Predicting Protein-Ligand Binding Free Energy: From Host-Guest Systems to the PDBbind Database

Calculation of protein-ligand binding affinity is a cornerstone of drug discovery. Classic implicit solvent models, which have been widely used to accomplish this task, lack accuracy compared to experimental references. Emerging data-driven models, on the other hand, are often accurate yet not fully interpretable and also likely to be overfitted. In this research, we explore the application of Theory-Guided Data Science in studying protein-ligand binding. A hybrid model is introduced by integrating Graph Convolutional Network (data-driven model) with the GBNSR6 implicit solvent (physics-based model). The proposed physics-data model is tested on a dataset of 368 complexes from the PDBbind refined set and 72 host-guest systems. Results demonstrate that the proposed Physics-Guided Neural Network can successfully improve the "accuracy" of the pure data-driven model. In addition, the "interpretability" and "transferability" of our model have boosted compared to the purely data-driven model. Further analyses include evaluating model robustness and understanding relationships between the physical features.

Simultaneous Assignment and Structure Determination of Proteins From Sparsely Labeled NMR Datasets

Sparsely labeled NMR samples provide opportunities to study larger biomolecular assemblies than is traditionally done by NMR. This requires new computational tools that can handle the sparsity and ambiguity in the NMR datasets. The MELD (modeling employing limited data) Bayesian approach was assessed to be the best performing in predicting structures from sparsely labeled NMR data in the 13th edition of the Critical Assessment of Structure Prediction (CASP) event—and limitations of the methodology were also noted. In this report, we evaluate the nature and difficulty in modeling unassigned sparsely labeled NMR datasets and report on an improved methodological pipeline leading to higher-accuracy predictions. We benchmark our methodology against the NMR datasets provided by CASP 13.

Conformational Transitions of Amyloid-β: A Langevin and Generalized Langevin Dynamics Simulation Study

The dynamics of conformational transitions of the disordered protein, amyloid-β, is studied via Langevin and generalized Langevin dynamics simulations. The transmission coefficient for the unfold-misfold transition of amyloid-β is calculated from multiple independent trajectories that originate at the transition state with different initial velocities and are directly correlated to Kramers and Grote-Hynes theories. For lower values of the frictional coefficient, a well-defined rate constant is obtained, whereas, for higher values, the transmission coefficient decays with time, indicating a breakdown of the Kramers and Grote-Hynes theories and the emergence of a dynamic disorder, which demonstrates the presence of multiple local minima in the misfolding potential energy surface. The calculated free energy profile describes a two-state transition of amyloid-β in the energy landscape. The transition path time distribution computed from these simulations is compared with the related experimental and theoretical results for the unfold-misfold transition of amyloid-β. The high free energy barrier for this transition confirms the misfolding of amyloid-β. These findings offer an insight into the dynamics of the unfold-misfold transition of this protein.

Actual Symmetry of Symmetric Molecular Adducts in the Gas Phase, Solution and in the Solid State

This review discusses molecular adducts, whose composition allows a symmetric structure. Such adducts are popular model systems, as they are useful for analyzing the effect of structure on the property selected for study since they allow one to reduce the number of parameters. The main objectives of this discussion are to evaluate the influence of the surroundings on the symmetry of these adducts, steric hindrances within the adducts, competition between different noncovalent interactions responsible for stabilizing the adducts, and experimental methods that can be used to study the symmetry at different time scales. This review considers the following central binding units: hydrogen (proton), halogen (anion), metal (cation), water (hydrogen peroxide).

Epitopes of Protein Binders Are Related to the Structural Flexibility of a Target Protein Surface

Protein binders including antibodies are known not to bind to random sites of target proteins, and their functional effectiveness mainly depends on the binding region, called the epitope. For the development of protein binders with desired functions, it is thus critical to understand which surface region protein binders prefer (or do not prefer) to bind. The current methods for epitope prediction focus on static indicators such as structural geometry or amino acid propensity, whereas protein binding events are in fact a consequence of dynamic interactions. Here, we demonstrate that the preference for a binding site by protein binders is strongly related to the structural flexibility of a target protein surface. Molecular dynamics simulations on unbound forms of antigen structures revealed that the antigen surface in direct contact with antibodies is less flexible than the rest of the surface. This tendency was shown to be similar in other non-antibody protein binders such as affibody, DARPin, monobody, and repebody. We also found that the relatedness of epitopes to the structural flexibility of a target protein surface is dependent on the secondary structure elements of paratopes. Monobody and repebody, whose binding sites are composed of β-strands, distinctively prefer to bind to a relatively more rigid region of a target protein. These observations enabled us to develop a simple epitope prediction method which shows a comparable performance to the commonly used ones.

DFT investigation of solvent, substituent, and catalysis effects on the intramolecular Diels-Alder reaction

In this study, we report on a DFT investigation of two intramolecular Diels-Alder furan reactions. Optimizations of the studied structures, TS and IRC calculations, were carried out at B3LYP/6-31G(d) level. We have studied the effect of substituent, solvent and Lewis acid catalyst on cyclization-retrocyclization equilibria, activation energies, and stability of the desired products. The analysis of orbital coefficients, IRC curves, and Wiberg indices have proved that both reactions are under orbital control. We have found that for the reaction I (2↔4 + 5), where R = H, the exo attack is favored by hydrogen bond interaction, while for R = t-Bu, the steric hindrance leads to the endo attack. For the reaction II (3 → 6 + 7), the t-Bu-substituted products are the most stable ones. At another level, we have found that it is recommended to use polar organic solvents as DMSO with Lewis acid catalyst BF₃. The latest leads to accelerate the reaction II with stabilization of the desired products.

Water in Nanopores and Biological Channels: A Molecular Simulation Perspective

This Review explores the dynamic behavior of water within nanopores and biological channels in lipid bilayer membranes. We focus on molecular simulation studies, alongside selected structural and other experimental investigations. Structures of biological nanopores and channels are reviewed, emphasizing those high-resolution crystal structures, which reveal water molecules within the transmembrane pores, which can be used to aid the interpretation of simulation studies. Different levels of molecular simulations of water within nanopores are described, with a focus on molecular dynamics (MD). In particular, models of water for MD simulations are discussed in detail to provide an evaluation of their use in simulations of water in nanopores. Simulation studies of the behavior of water in idealized models of nanopores have revealed aspects of the organization and dynamics of nanoconfined water, including wetting/dewetting in narrow hydrophobic nanopores. A survey of simulation studies in a range of nonbiological nanopores is presented, including carbon nanotubes, synthetic nanopores, model peptide nanopores, track-etched nanopores in polymer membranes, and hydroxylated and functionalized nanoporous silica. These reveal a complex relationship between pore size/geometry, the nature of the pore lining, and rates of water transport. Wider nanopores with hydrophobic linings favor water flow whereas narrower hydrophobic pores may show dewetting. Simulation studies over the past decade of the behavior of water in a range of biological nanopores are described, including porins and β-barrel protein nanopores, aquaporins and related polar solute pores, and a number of different classes of ion channels. Water is shown to play a key role in proton transport in biological channels and in hydrophobic gating of ion channels. An overall picture emerges, whereby the behavior of water in a nanopore may be predicted as a function of its hydrophobicity and radius. This informs our understanding of the functions of diverse channel structures and will aid the design of novel nanopores. Thus, our current level of understanding allows for the design of a nanopore which promotes wetting over dewetting or vice versa. However, to design a novel nanopore, which enables fast, selective, and gated flow of water de novo would remain challenging, suggesting a need for further detailed simulations alongside experimental evaluation of more complex nanopore systems.

Electrostatic control of calcineurin's intrinsically-disordered regulatory domain binding to calmodulin

Calcineurin (CaN) is a serine/threonine phosphatase that regulates a variety of physiological and pathophysiological processes in mammalian tissue. The calcineurin (CaN) regulatory domain (RD) is responsible for regulating the enzyme's phosphatase activity, and is believed to be highly-disordered when inhibiting CaN, but undergoes a disorder-to-order transition upon diffusion-limited binding with the regulatory protein calmodulin (CaM). The prevalence of polar and charged amino acids in the regulatory domain (RD) suggests electrostatic interactions are involved in mediating calmodulin (CaM) binding, yet the lack of atomistic-resolution data for the bound complex has stymied efforts to probe how the RD sequence controls its conformational ensemble and long-range attractions contribute to target protein binding. In the present study, we investigated via computational modeling the extent to which electrostatics and structural disorder facilitate CaM/CaN association kinetics. Specifically, we examined several RD constructs that contain the CaM binding region (CAMBR) to characterize the roles of electrostatics versus conformational diversity in controlling diffusion-limited association rates, via microsecond-scale molecular dynamics (MD) and Brownian dynamic (BD) simulations. Our results indicate that the RD amino acid composition and sequence length influence both the dynamic availability of conformations amenable to CaM binding, as well as long-range electrostatic interactions to steer association. These findings provide intriguing insight into the interplay between conformational diversity and electrostatically-driven protein-protein association involving CaN, which are likely to extend to wide-ranging diffusion-limited processes regulated by intrinsically-disordered proteins.

Molecular Dynamics Simulations of Protein-Drug Complexes: A Computational Protocol for Investigating the Interactions of Small-Molecule Therapeutics with Biological Targets and Biosensors

MD simulations provide a powerful tool for the investigation of protein-drug complexes. The following chapter uses the aryl acylamidase-acetaminophen system as an example to describe a general protocol for preparing and running simulations of protein-drug complexes, complete with a step-by-step tutorial. The described approach is broadly applicable toward the study of drug interactions in the context of both biological targets and biosensing enzymes.

Looking "Under the Hood" of Cellular Mechanotransduction with Computational Tools: A Systems Biomechanics Approach across Multiple Scales

Signal modulation has been developed in living cells throughout evolution to promote utilizing the same machinery for multiple cellular functions. Chemical and mechanical modules of signal transmission and transduction are interconnected and necessary for organ development and growth. However, due to the high complexity of the intercommunication of physical intracellular connections with biochemical pathways, there are many missing details in our overall understanding of mechanotransduction processes, i.e., the process by which mechanical signals are converted to biochemical cascades. Cell-matrix adhesions are mechanically coupled to the nucleus through the cytoskeleton. This modulated and tightly integrated network mediates the transmission of mechanochemical signals from the extracellular matrix to the nucleus. Various experimental and computational techniques have been utilized to understand the basic mechanisms of mechanotransduction, yet many aspects have remained elusive. Recently, in silico experiments have made important contributions to the field of mechanobiology. Herein, computational modeling efforts devoted to understanding integrin-mediated mechanotransduction pathways are reviewed, and an outlook is presented for future directions toward using suitable computational approaches and developing novel techniques for addressing important questions in the field of mechanotransduction.

Parameterization of the Hamiltonian Dielectric Solvent (HADES) Reaction-Field Method for the Solvation Free Energies of Amino Acid Side-Chain Analogs

Optimization of the Hamiltonian dielectric solvent (HADES) method for biomolecular simulations in a dielectric continuum is presented with the goal of calculating accurate absolute solvation free energies while retaining the model's accuracy in predicting conformational free-energy differences. The solvation free energies of neutral and polar amino acid side-chain analogs calculated by using HADES, which may optionally include nonpolar contributions, were optimized against experimental data to reach a chemical accuracy of about 0.5 kcal mol(-1). The new parameters were evaluated for charged side-chain analogs. The HADES results were compared with explicit-solvent, generalized Born, Poisson-Boltzmann, and QM-based methods. The potentials of mean force (PMFs) between pairs of side-chain analogs obtained by using HADES and explicit-solvent simulations were used to evaluate the effects of the improved parameters optimized for solvation free energies on intermolecular potentials.

Structure Based Modeling of Small Molecules Binding to the TLR7 by Atomistic Level Simulations

Toll-Like Receptors (TLR) are a large family of proteins involved in the immune system response. Both the activation and the inhibition of these receptors can have positive effects on several diseases, including viral pathologies and cancer, therefore prompting the development of new compounds. In order to provide new indications for the design of Toll-Like Receptor 7 (TLR7)-targeting drugs, the mechanism of interaction between the TLR7 and two important classes of agonists (imidazoquinoline and adenine derivatives) was investigated through docking and Molecular Dynamics simulations. To perform the computational analysis, a new model for the dimeric form of the receptors was necessary and therefore created. Qualitative and quantitative differences between agonists and inactive compounds were determined. The in silico results were compared with previous experimental observations and employed to define the ligand binding mechanism of TLR7.

Molecular modeling of nucleic Acid structure: electrostatics and solvation

This unit presents an overview of computer simulation techniques as applied to nucleic acid systems, ranging from simple in vacuo molecular modeling techniques to more complete all-atom molecular dynamics treatments that include an explicit representation of the environment. The third in a series of four units, this unit focuses on critical issues in solvation and the treatment of electrostatics. UNITS 7.5 & 7.8 introduced the modeling of nucleic acid structure at the molecular level. This included a discussion of how to generate an initial model, how to evaluate the utility or reliability of a given model, and ultimately how to manipulate this model to better understand its structure, dynamics, and interactions. Subject to an appropriate representation of the energy, such as a specifically parameterized empirical force field, the techniques of minimization and Monte Carlo simulation, as well as molecular dynamics (MD) methods, were introduced as a way of sampling conformational space for a better understanding of the relevance of a given model. This discussion highlighted the major limitations with modeling in general. When sampling conformational space effectively, difficult issues are encountered, such as multiple minima or conformational sampling problems, and accurately representing the underlying energy of interaction. In order to provide a realistic model of the underlying energetics for nucleic acids in their native environments, it is crucial to include some representation of solvation (by water) and also to properly treat the electrostatic interactions. These subjects are discussed in detail in this unit.

Many-body effect in ion binding to RNA

Ion-mediated electrostatic interactions play an important role in RNA folding stability. For a RNA in a solution with higher Mg(2+) ion concentration, more counterions in the solution can bind to the RNA, causing a strong many-body coupling between the bound ions. The many-body effect can change the effective potential of mean force between the tightly bound ions. This effect tends to dampen ion binding and lower RNA folding stability. Neglecting the many-body effect leads to a systematic error (over-estimation) of RNA folding stability at high Mg(2+) ion concentrations. Using the tightly bound ion model combined with a conformational ensemble model, we investigate the influence of the many-body effect on the ion-dependent RNA folding stability. Comparisons with the experimental data indicate that including the many-body effect led to much improved predictions for RNA folding stability at high Mg(2+) ion concentrations. The results suggest that the many-body effect can be important for RNA folding in high concentrations of multivalent ions. Further investigation showed that the many-body effect can influence the spatial distribution of the tightly bound ions and the effect is more pronounced for compact RNA structures and structures prone to the formation of local clustering of ions.

Modified Poisson–Boltzmann equations for characterizing biomolecular solvation

Methods for computing electrostatic interactions often account implicitly for the solvent, due to the much smaller number of degrees of freedom involved. In the Poisson–Boltzmann (PB) approach the electrostatic potential is obtained by solving the Poisson–Boltzmann equation (PBE), where the solvent region is modeled as a homogeneous medium with a high dielectric constant. PB however is not exempt of problems. It does not take into account for example the sizes of the ions in the atmosphere surrounding the solute, nor does it take into account the inhomogeneous dielectric response of water due to the presence of a highly charged surface. In this paper we review two major modifications of PB that circumvent these problems, namely the size-modified PB (SMPB) equation and the Dipolar Poisson–Boltzmann Langevin (DPBL) model. In SMPB, steric effects between ions are accounted for with a lattice gas model. In DPBL, the solvent region is no longer modeled as a homogeneous dielectric media but rather as an assembly of self-orienting interacting dipoles of variable density. This model results in a dielectric profile that transits smoothly from the solute to the solvent region as well as in a variable solvent density that depends on the charges of the solute. We show successful applications of the DPBL formalism to computing the solvation free energies of isolated ions in water. Further developments of more accurately modified PB models are discussed.

Design and application of implicit solvent models in biomolecular simulations

•

Implicit solvent replaces explicit water by a potential of mean force.

•

Popular models are SASA, VOL and Generalized Born.

•

Implicit solvent is used in MD, protein modelling, folding, design, prediction and drug screening.

•

Large-scale simulations allow for parametrisation via force matching.

•

Application to nucleic acids and membranes is challenging.

We review implicit solvent models and their parametrisation by introducing the concepts and recent devlopments of the most popular models with a focus on parametrisation via force matching. An overview of recent applications of the solvation energy term in protein dynamics, modelling, design and prediction is given to illustrate the usability and versatility of implicit solvation in reproducing the physical behaviour of biomolecular systems. Limitations of implicit modes are discussed through the example of more challenging systems like nucleic acids and membranes.

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.

Improving Generalized Born Models by Exploiting Connections to Polarizable Continuum Models. II. Corrections for Salt Effects

A previous analytical investigation of the generalized Born (GB) implicit solvation model is extended to solvents of nonzero ionic strength. The GB model with salt effects (GB-SE) is shown to resemble the Debye-Hückel-like screening model (DESMO), a polarizable continuum model (PCM) that we have recently developed for salty solutions. DESMO may be regarded either as a generalization of the conductor-like PCM (C-PCM) that extends C-PCM to electrolyte solutions or alternatively as a generalization of Debye-Hückel theory to arbitrary cavity shapes. The connection between GB-SE and DESMO suggests how the former can be modified to account for the exclusion of mobile ions from the cavity interior, an effect that is typically absent in GB-SE models. We propose two simple GB-SE models that are exact for a point charge in a spherical cavity and that introduce the ability to account, albeit approximately, for the finite size of the mobile ions. The accuracy of these new models is demonstrated by applications to both model systems and real proteins. These tests also demonstrate the accuracy of the DESMO approach, as compared to more sophisticated PCMs developed for electrolyte solutions.

Fast Analytical Methods for Macroscopic Electrostatic Models in Biomolecular Simulations

We review recent developments of fast analytical methods for macroscopic electrostatic calculations in biological applications, including the Poisson-Boltzmann (PB) and the generalized Born models for electrostatic solvation energy. The focus is on analytical approaches for hybrid solvation models, especially the image charge method for a spherical cavity, and also the generalized Born theory as an approximation to the PB model. This review places much emphasis on the mathematical details behind these methods.

A simple polarizable continuum solvation model for electrolyte solutions

We propose a Debye-Hückel-like screening model (DESMO) that generalizes the familiar conductor-like screening model (COSMO) to solvents with non-zero ionic strength and furthermore provides a numerical generalization of the Debye-Hückel model that is applicable to non-spherical solute cavities. The numerical implementation of DESMO is based upon the switching/Gaussian (SWIG) method for smooth cavity discretization, which we have recently introduced in the context of polarizable continuum models (PCMs). This approach guarantees that the potential energy is a smooth function of the solute geometry and analytic gradients for DESMO are reported here. The SWIG formalism also facilitates analytic implementation of two other PCMs that are based on a screened Coulomb potential: the "integral equation formalism" (IEF-PCM) and the "surface and simulation of volume polarization for electrostatics" [SS(V)PE] method. Fully analytic implementations of these screened PCMs are reported here for the first time. Numerical results, for model systems where an exact solution of the linearized Poisson-Boltzmann equation is available, demonstrate that these screened PCMs are highly accurate. In realistic test cases, they are as accurate as the best available three-dimensional finite-difference methods. In polar solvents, DESMO is nearly as accurate as more sophisticated screened PCMs, but is significantly simpler and more efficient.

A strategy for reducing gross errors in the generalized Born models of implicit solvation

The "canonical" generalized Born (GB) formula [C. Still, A. Tempczyk, R. C. Hawley, and T. Hendrickson, J. Am. Chem. Soc. 112, 6127 (1990)] is known to provide accurate estimates for total electrostatic solvation energies ΔG(el) of biomolecules if the corresponding effective Born radii are accurate. Here we show that even if the effective Born radii are perfectly accurate, the canonical formula still exhibits significant number of gross errors (errors larger than 2k(B)T relative to numerical Poisson equation reference) in pairwise interactions between individual atomic charges. Analysis of exact analytical solutions of the Poisson equation (PE) for several idealized nonspherical geometries reveals two distinct spatial modes of the PE solution; these modes are also found in realistic biomolecular shapes. The canonical GB Green function misses one of two modes seen in the exact PE solution, which explains the observed gross errors. To address the problem and reduce gross errors of the GB formalism, we have used exact PE solutions for idealized nonspherical geometries to suggest an alternative analytical Green function to replace the canonical GB formula. The proposed functional form is mathematically nearly as simple as the original, but depends not only on the effective Born radii but also on their gradients, which allows for better representation of details of nonspherical molecular shapes. In particular, the proposed functional form captures both modes of the PE solution seen in nonspherical geometries. Tests on realistic biomolecular structures ranging from small peptides to medium size proteins show that the proposed functional form reduces gross pairwise errors in all cases, with the amount of reduction varying from more than an order of magnitude for small structures to a factor of 2 for the largest ones.

Treecode-based generalized Born method

We have developed a treecode-based O(N log N) algorithm for the generalized Born (GB) implicit solvation model. Our treecode-based GB (tGB) is based on the GBr6 [J. Phys. Chem. B 111, 3055 (2007)], an analytical GB method with a pairwise descreening approximation for the R6 volume integral expression. The algorithm is composed of a cutoff scheme for the effective Born radii calculation, and a treecode implementation of the GB charge-charge pair interactions. Test results demonstrate that the tGB algorithm can reproduce the vdW surface based Poisson solvation energy with an average relative error less than 0.6% while providing an almost linear-scaling calculation for a representative set of 25 proteins with different sizes (from 2815 atoms to 65456 atoms). For a typical system of 10k atoms, the tGB calculation is three times faster than the direct summation as implemented in the original GBr6 model. Thus, our tGB method provides an efficient way for performing implicit solvent GB simulations of larger biomolecular systems at longer time scales.

Advances in implicit models of water solvent to compute conformational free energy and molecular dynamics of proteins at constant pH

Modern implicit solvent models for macromolecular simulations in water-proton bath are considered. The fundamental quantity that implicit models approximate is the solute potential of mean force, which is obtained by averaging over solvent degrees of freedom. The implicit solvent models suggest practical ways to calculate free energies of macromolecular conformations taking into account equilibrium interactions with water solvent and proton bath, while the explicit solvent approach is unable to do that due to the need to account for a large number of solvent degrees of freedom. The most advanced realizations of the implicit continuum models by different research groups are discussed, their accuracy are examined, and some applications of the implicit solvent models to macromolecular modeling, such as free energy calculations, protein folding, and constant pH molecular dynamics are highlighted.

Multipole electrostatics in hydration free energy calculations

Hydration free energy (HFE) is generally used for evaluating molecular solubility, which is an important property for pharmaceutical and chemical engineering processes. Accurately predicting HFE is also recognized as one fundamental capability of molecular mechanics force field. Here, we present a systematic investigation on HFE calculations with AMOEBA polarizable force field at various parameterization and simulation conditions. The HFEs of seven small organic molecules have been obtained alchemically using the Bennett Acceptance Ratio method. We have compared two approaches to derive the atomic multipoles from quantum mechanical calculations: one directly from the new distributed multipole analysis and the other involving fitting to the electrostatic potential around the molecules. Wave functions solved at the MP2 level with four basis sets (6-311G*, 6-311++G(2d,2p), cc-pVTZ, and aug-cc-pVTZ) are used to derive the atomic multipoles. HFEs from all four basis sets show a reasonable agreement with experimental data (root mean square error 0.63 kcal/mol for aug-cc-pVTZ). We conclude that aug-cc-pVTZ gives the best performance when used with AMOEBA, and 6-311++G(2d,2p) is comparable but more efficient for larger systems. The results suggest that the inclusion of diffuse basis functions is important for capturing intermolecular interactions. The effect of long-range correction to van der Waals interaction on the hydration free energies is about 0.1 kcal/mol when the cutoff is 12Å, and increases linearly with the number of atoms in the solute/ligand. In addition, we also discussed the results from a hybrid approach that combines polarizable solute with fixed-charge water in the HFE calculation.

Accelerating electrostatic surface potential calculation with multi-scale approximation on graphics processing units

Tools that compute and visualize biomolecular electrostatic surface potential have been used extensively for studying biomolecular function. However, determining the surface potential for large biomolecules on a typical desktop computer can take days or longer using currently available tools and methods. Two commonly used techniques to speed-up these types of electrostatic computations are approximations based on multi-scale coarse-graining and parallelization across multiple processors. This paper demonstrates that for the computation of electrostatic surface potential, these two techniques can be combined to deliver significantly greater speed-up than either one separately, something that is in general not always possible. Specifically, the electrostatic potential computation, using an analytical linearized Poisson-Boltzmann (ALPB) method, is approximated using the hierarchical charge partitioning (HCP) multi-scale method, and parallelized on an ATI Radeon 4870 graphical processing unit (GPU). The implementation delivers a combined 934-fold speed-up for a 476,040 atom viral capsid, compared to an equivalent non-parallel implementation on an Intel E6550 CPU without the approximation. This speed-up is significantly greater than the 42-fold speed-up for the HCP approximation alone or the 182-fold speed-up for the GPU alone.

The AGBNP2 Implicit Solvation Model

The AGBNP2 implicit solvent model, an evolution of the Analytical Generalized Born plus Non-Polar (AGBNP) model we have previously reported, is presented with the aim of modeling hydration effects beyond those described by conventional continuum dielectric representations. A new empirical hydration free energy component based on a procedure to locate and score hydration sites on the solute surface is introduced to model first solvation shell effects, such as hydrogen bonding, which are poorly described by continuum dielectric models. This new component is added to the Generalized Born and non-polar AGBNP terms. Also newly introduced is an analytical Solvent Excluded Volume (SEV) model which improves the solute volume description by reducing the effect of spurious high-dielectric interstitial spaces present in conventional van der Waals representations. The AGBNP2 model is parametrized and tested with respect to experimental hydration free energies of small molecules and the results of explicit solvent simulations. Modeling the granularity of water is one of the main design principles employed for the the first shell solvation function and the SEV model, by requiring that water locations have a minimum available volume based on the size of a water molecule. It is shown that the new volumetric model produces Born radii and surface areas in good agreement with accurate numerical evaluations of these quantities. The results of molecular dynamics simulations of a series of mini-proteins show that the new model produces conformational ensembles in substantially better agreement with reference explicit solvent ensembles than the original AGBNP model with respect to both structural and energetics measures.