Parallelization and improvements of the generalized born model with a simple sWitching function for modern graphics processors

A generalized Kirkwood implicit solvent for the polarizable AMOEBA protein model

Computational simulation of biomolecules can provide important insights into protein design, protein-ligand binding interactions, and ab initio biomolecular folding, among other applications. Accurate treatment of the solvent environment is essential in such applications, but the use of explicit solvents can add considerable cost. Implicit treatment of solvent effects using a dielectric continuum model is an attractive alternative to explicit solvation since it is able to describe solvation effects without the inclusion of solvent degrees of freedom. Previously, we described the development and parameterization of implicit solvent models for small molecules. Here, we extend the parameterization of the generalized Kirkwood (GK) implicit solvent model for use with biomolecules described by the AMOEBA force field via the addition of corrections to the calculation of effective radii that account for interstitial spaces that arise within biomolecules. These include element-specific pairwise descreening scale factors, a short-range neck contribution to describe the solvent-excluded space between pairs of nearby atoms, and finally tanh-based rescaling of the overall descreening integral. We then apply the AMOEBA/GK implicit solvent to a set of ten proteins and achieve an average coordinate root mean square deviation for the experimental structures of 2.0 Å across 500 ns simulations. Overall, the continued development of implicit solvent models will help facilitate the simulation of biomolecules on mechanistically relevant timescales.

CHARMM-GUI Implicit Solvent Modeler for Various Generalized Born Models in Different Simulation Programs

Implicit solvent models are widely used because they are advantageous to speed up simulations by drastically decreasing the number of solvent degrees of freedom, which allows one to achieve long simulation time scales for large system sizes. CHARMM-GUI, a web-based platform, has been developed to support the setup of complex multicomponent molecular systems and prepare input files. This study describes an Implicit Solvent Modeler (ISM) in CHARMM-GUI for various generalized Born (GB) implicit solvent simulations in different molecular dynamics programs such as AMBER, CHARMM, GENESIS, NAMD, OpenMM, and Tinker. The GB models available in ISM include GB-HCT, GB-OBC, GB-neck, GBMV, and GBSW with the CHARMM and Amber force fields for protein, DNA, RNA, glycan, and ligand systems. Using the system and input files generated by ISM, implicit solvent simulations of protein, DNA, and RNA systems produce similar results for different simulation packages with the same input information. Protein-ligand systems are also considered to further validate the systems and input files generated by ISM. Simple ligand root-mean-square deviation (RMSD) and molecular mechanics generalized Born surface area (MM/GBSA) calculations show that the performance of implicit simulations is better than docking and can be used for early stage ligand screening. These reasonable results indicate that ISM is a useful and reliable tool to provide various implicit solvent simulation applications.

A Guide to the Continuous Constant pH Molecular Dynamics Methods in Amber and CHARMM [Article v1.0]

Like temperature and pressure, solution pH is an important environmental variable in biomolecular simulations. Virtually all proteins depend on pH to maintain their structure and function. In conventional molecular dynamics (MD) simulations of proteins, pH is implicitly accounted for by assigning and fixing protonation states of titratable sidechains. This is a significant limitation, as the assigned protonation states may be wrong and they may change during dynamics. In this tutorial, we guide the reader in learning and using the various continuous constant pH MD methods in Amber and CHARMM packages, which have been applied to predict pK a values and elucidate proton-coupled conformational dynamics of a variety of proteins including enzymes and membrane transporters.

Discovery of Benzopyridone-Based Transient Receptor Potential Vanilloid 1 Agonists and Antagonists and the Structural Elucidation of Their Activity Shift

Among a series of benzopyridone-based scaffolds investigated as human transient receptor potential vanilloid 1 (TRPV1) ligands, two isomeric benzopyridone scaffolds demonstrated a consistent and distinctive functional profile in which 2-oxo-1,2-dihydroquinolin-5-yl analogues (e.g., 2) displayed high affinity and potent antagonism, whereas 1-oxo-1,2-dihydroisoquinolin-5-yl analogues (e.g., 3) showed full agonism with high potency. Our computational models provide insight into the agonist-antagonist boundary of the analogues suggesting that the Arg557 residue in the S4-S5 linker might be important for sensing the agonist binding and transmitting signals. These results provide structural insights into the TRPV1 and the protein-ligand interactions at a molecular level.

"Big Data" Fast Chemoinformatics Model to Predict Generalized Born Radius and Solvent Accessibility as a Function of Geometry

The Generalized Born (GB) solvent model is offering the best accuracy/computing effort ratio yet requires drastic simplifications to estimate of the Effective Born Radii (EBR) in bypassing a too expensive volume integration step. EBRs are a measure of the degree of burial of an atom and not very sensitive to small changes of geometry: in molecular dynamics, the costly EBR update procedure is not mandatory at every step. This work however aims at implementing a GB model into the Sampler for Multiple Protein-Ligand Entities (S4MPLE) evolutionary algorithm with mandatory EBR updates at each step triggering arbitrarily large geometric changes. Therefore, a quantitative structure-property relationship has been developed in order to express the EBRs as a linear function of both the topological neighborhood and geometric occupancy of the space around atoms. A training set of 810 molecular systems, starting from fragment-like to drug-like compounds, proteins, host-guest systems, and ligand-protein complexes, has been compiled. For each species, S4MPLE generated several hundreds of random conformers. For each atom in each geometry of each species, its "standard" EBR was calculated by numeric integration and associated to topological and geometric descriptors of the atom neighborhood. This training set (EBR, atom descriptors) involving >5 M entries was subjected to a boot-strapping multilinear regression process with descriptor selection. In parallel, the strategy was repurposed to also learn atomic solvent-accessible areas (SA) based on the same descriptors. Resulting linear equations were challenged to predict EBR and SA values for a similarly compiled external set of >2000 new molecular systems. Solvation energies calculated with estimated EBR and SA match "standard" energies within the typical error of a force-field-based approach (a few kilocalories per mole). Given the extreme diversity of molecular systems covered by the model, this simple EBR/SA estimator covers a vast applicability domain.

Computer simulations of protein-membrane systems

The interactions between proteins and membranes play critical roles in signal transduction, cell motility, and transport, and they are involved in many types of diseases. Molecular dynamics (MD) simulations have greatly contributed to our understanding of protein-membrane interactions, promoted by a dramatic development of MD-related software, increasingly accurate force fields, and available computer power. In this chapter, we present available methods for studying protein-membrane systems with MD simulations, including an overview about the various all-atom and coarse-grained force fields for lipids, and useful software for membrane simulation setup and analysis. A large set of case studies is discussed.

Accelerating the Generalized Born with Molecular Volume and Solvent Accessible Surface Area Implicit Solvent Model Using Graphics Processing Units

The generalized Born with molecular volume and solvent accessible surface area (GBMV2/SA) implicit solvent model provides an accurate description of molecular volume and has the potential to accurately describe the conformational equilibria of structured and disordered proteins. However, its broader application has been limited by the computational cost and poor scaling in parallel computing. Here, we report an efficient implementation of both the electrostatic and nonpolar components of GBMV2/SA on graphics processing unit (GPU) within the CHARMM/OpenMM module. The GPU-GBMV2/SA is numerically equivalent to the original CPU-GBMV2/SA. The GPU acceleration offers ~60- to 70-fold speedup on a single NVIDIA TITAN X (Pascal) graphics card for molecular dynamic simulations of both folded and unstructured proteins of various sizes. The current implementation can be further optimized to achieve even greater acceleration with minimal reduction on the numerical accuracy. The successful development of GPU-GBMV2/SA greatly facilitates its application to biomolecular simulations and paves the way for further development of the implicit solvent methodology. © 2019 Wiley Periodicals, Inc.

CHARMM-GUI Glycan Modeler for modeling and simulation of carbohydrates and glycoconjugates

Characterizing glycans and glycoconjugates in the context of three-dimensional structures is important in understanding their biological roles and developing efficient therapeutic agents. Computational modeling and molecular simulation have become an essential tool complementary to experimental methods. Here, we present a computational tool, Glycan Modeler for in silico N-/O-glycosylation of the target protein and generation of carbohydrate-only systems. In our previous study, we developed Glycan Reader, a web-based tool for detecting carbohydrate molecules from a PDB structure and generation of simulation system and input files. As integrated into Glycan Reader in CHARMM-GUI, Glycan Modeler (Glycan Reader & Modeler) enables to generate the structures of glycans and glycoconjugates for given glycan sequences and glycosylation sites using PDB glycan template structures from Glycan Fragment Database (http://glycanstructure.org/fragment-db). Our benchmark tests demonstrate the universal applicability of Glycan Reader & Modeler to various glycan sequences and target proteins. We also investigated the structural properties of modeled glycan structures by running 2-μs molecular dynamics simulations of HIV envelope protein. The simulations show that the modeled glycan structures built by Glycan Reader & Modeler have the similar structural features compared to the ones solved by X-ray crystallography. We also describe the representative examples of glycoconjugate modeling with video demos to illustrate the practical applications of Glycan Reader & Modeler. Glycan Reader & Modeler is freely available at http://charmm-gui.org/input/glycan.

Molecular Docking and Molecular Dynamics (MD) Simulation of Human Anti-Complement Factor H (CFH) Antibody Ab42 and CFH Polypeptide

An understanding of the interaction between the antibody and its targeted antigen and knowing of the epitopes are critical for the development of monoclonal antibody drugs. Complement factor H (CFH) is implied to play a role in tumor growth and metastasis. An autoantibody to CHF is associated with anti-tumor cell activity. The interaction of a human monoclonal antibody Ab42 that was isolated from a cancer patient with CFH polypeptide (pCFH) antigen was analyzed by molecular docking, molecular dynamics (MD) simulation, free energy calculation, and computational alanine scanning (CAS). Experimental alanine scanning (EAS) was then carried out to verify the results of the theoretical calculation. Our results demonstrated that the Ab42 antibody interacts with pCFH by hydrogen bonds through the Tyr315, Ser100, Gly33, and Tyr53 residues on the complementarity-determining regions (CDRs), respectively, with the amino acid residues of Pro441, Ile442, Asp443, Asn444, Ile447, and Thr448 on the pCFH antigen. In conclusion, this study has explored the mechanism of interaction between Ab42 antibody and its targeted antigen by both theoretical and experimental analysis. Our results have important theoretical significance for the design and development of relevant antibody drugs.

Generalized Born Implicit Solvent Models for Biomolecules

It would often be useful in computer simulations to use an implicit description of solvation effects, instead of explicitly representing the individual solvent molecules. Continuum dielectric models often work well in describing the thermodynamic aspects of aqueous solvation and can be very efficient compared to the explicit treatment of the solvent. Here, we review a particular class of so-called fast implicit solvent models, generalized Born (GB) models, which are widely used for molecular dynamics (MD) simulations of proteins and nucleic acids. These approaches model hydration effects and provide solvent-dependent forces with efficiencies comparable to molecular-mechanics calculations on the solute alone; as such, they can be incorporated into MD or other conformational searching strategies in a straightforward manner. The foundations of the GB model are reviewed, followed by examples of newer, emerging models and examples of important applications. We discuss their strengths and weaknesses, both for fidelity to the underlying continuum model and for the ability to replace explicit consideration of solvent molecules in macromolecular simulations.

CHARMM-GUI Membrane Builder for Complex Biological Membrane Simulations with Glycolipids and Lipoglycans

Glycolipids (such as glycoglycerolipids, glycosphingolipids, and glycosylphosphatidylinositol) and lipoglycans (such as lipopolysaccharides (LPS), lipooligosaccharides (LOS), mycobacterial lipoarabinomannan, and mycoplasma lipoglycans) are typically found on the surface of cell membranes and play crucial roles in various cellular functions. Characterizing their structure and dynamics at the molecular level is essential to understand their biological roles, but systematic generation of glycolipid and lipoglycan structures is challenging because of great variations in lipid structures and glycan sequences (i.e., carbohydrate types and their linkages). To facilitate the generation of all-atom glycolipid/LPS/LOS structures, we have developed Glycolipid Modeler and LPS Modeler in CHARMM-GUI ( http://www.charmm-gui.org ), a web-based interface that simplifies building of complex biological simulation systems. In addition, we have incorporated these modules into Membrane Builder so that users can readily build a complex symmetric or asymmetric biological membrane system with various glycolipids and LPS/LOS. These tools are expected to be useful in innovative and novel glycolipid/LPS/LOS modeling and simulation research by easing tedious and intricate steps in modeling complex biological systems and shall provide insight into structures, dynamics, and underlying mechanisms of complex glycolipid-/LPS-/LOS-containing biological membrane systems.

Peptide and Protein Structure Prediction with a Simplified Continuum Solvent Model

A continuum solvent model based on screened Coulomb potentials has been simplified and parametrized to sample native-like structures in replica-exchange simulations of each of six different peptides and miniproteins. Low-energy, native, and non-native structures were used to iteratively refine 11 parameter values. The centroid of the largest cluster of structures sampled in simulations initiated from an extended conformation represents the predicted structure. The main-chain rms deviation of this prediction from the experimental structure was 0.47 Å for the 12-residue Trp-zip2, 0.86 Å for the 14-residue MBH12, 2.53 Å for the 17-residue U(1-17)T9D, 2.03 Å for the 20-residue BS1, 1.08 Å for the 20-residue Trp-cage, and 3.64 Å for the 35-residue villin headpiece subdomain HP35. The centroid of the sixth largest cluster sampled for HP35 deviated by 0.91 Å. The CHARMM22/CMAP force field was used, with an additional ψ torsion term for residues other than glycine and proline. Six parameters govern the dielectric response of the continuum solvent, and four values of surface tension approximate nonpolar effects. An atom's self-energy and interaction energies are screened independently, each depending on whether the atom is part of a charged group, a neutral hydrogen-bonding main-chain group, or any other neutral group. The parameters inferred result in strong main-chain hydrogen bonds, consistent with the view that protein folding is dominated by the formation of these bonds. (1,2) Conformations of MBH12 and BS1 were excluded from the energy-function refinement, suggesting the parameters, referred to as SCP18, are transferable. An efficient estimate of solvent-accessible surface area is also described.

In silico site-directed mutagenesis of neutralizing mAb 4C4 and analysis of its interaction with G-H loop of VP1 to explore its therapeutic applications against FMD

Investigating the behaviour of bio-molecules through computational mutagenesis is gaining interest to facilitate the development of new therapeutic solutions for infectious diseases. The antigenetically variant genotypes of foot and mouth disease virus (FMDV) and their subsequent infections are challenging to tackle with traditional vaccination. In such scenario, neutralizing antibodies might provide an alternate solution to manage the FMDV infection. Thus, we have analysed the interaction of the mAb 4C4 with a synthetic G-H loop of FMDV-VP1 through in silico mutagenesis and molecular modelling. Initially, a set of 25,434 mutants were designed and the mutants having better energetic stability than 4C4 were clustered based on sequence identity. The best mutant representing each cluster was selected and evaluated for its binding affinity with the antigen in terms of docking scores, interaction energy and binding energy. Six mutants have confirmed better binding affinities towards the antigen than 4C4. Further, interaction of these mutants with the natural G-H loop that is bound to mAb SD6 was also evaluated. One 4C4 variant having mutations at the positions 2034(N→L), 2096(N→C), 2098(D→Y), 2532(T→K) and 2599(A→G) has revealed better binding affinities towards both the synthetic and natural G-H loops than 4C4 and SD6, respectively. A molecular dynamic simulation for 50 ns was conducted for mutant and wild-type antibody structures which supported the pre-simulation results. Therefore, these mutations on mAb 4C4 are believed to provide a better antibody-based therapeutic option for FMD. Communicated by Ramaswamy H. Sarma.

Glycan Reader is improved to recognize most sugar types and chemical modifications in the Protein Data Bank

Motivation

Glycans play a central role in many essential biological processes. Glycan Reader was originally developed to simplify the reading of Protein Data Bank (PDB) files containing glycans through the automatic detection and annotation of sugars and glycosidic linkages between sugar units and to proteins, all based on atomic coordinates and connectivity information. Carbohydrates can have various chemical modifications at different positions, making their chemical space much diverse. Unfortunately, current PDB files do not provide exact annotations for most carbohydrate derivatives and more than 50% of PDB glycan chains have at least one carbohydrate derivative that could not be correctly recognized by the original Glycan Reader.

Results

Glycan Reader has been improved and now identifies most sugar types and chemical modifications (including various glycolipids) in the PDB, and both PDB and PDBx/mmCIF formats are supported. CHARMM-GUI Glycan Reader is updated to generate the simulation system and input of various glycoconjugates with most sugar types and chemical modifications. It also offers a new functionality to edit the glycan structures through addition/deletion/modification of glycosylation types, sugar types, chemical modifications, glycosidic linkages, and anomeric states. The simulation system and input files can be used for CHARMM, NAMD, GROMACS, AMBER, GENESIS, LAMMPS, Desmond, OpenMM, and CHARMM/OpenMM. Glycan Fragment Database in GlycanStructure.Org is also updated to provide an intuitive glycan sequence search tool for complex glycan structures with various chemical modifications in the PDB.

Availability and implementation

http://www.charmm-gui.org/input/glycan and http://www.glycanstructure.org.

Contact

wonpil@lehigh.edu.

Supplementary information

Supplementary data are available at Bioinformatics online.

Hybrid All-Atom/Coarse-Grained Simulations of Proteins by Direct Coupling of CHARMM and PRIMO Force Fields

Hybrid all-atom/coarse-grained (AA/CG) simulations of proteins offer a computationally efficient compromise where atomistic details are only applied to biologically relevant regions while benefiting from the speedup of treating the remaining parts of a given system at the CG level. The recently developed CG model, PRIMO, allows a direct coupling with an atomistic force field with no additional modifications or coupling terms and the ability to carry out dynamic simulations without any restraints on secondary or tertiary structures. A hybrid AA/CG scheme based on combining all-atom CHARMM and coarse-grained PRIMO representations was validated via molecular dynamics and replica exchange simulations of soluble and membrane proteins. The AA/CG scheme was also tested in the calculation of the free energy profile for the transition from the closed to the open state of adenylate kinase via umbrella sampling molecular dynamics method. The overall finding is that the AA/CG scheme generates dynamics and energetics that are qualitatively and quantitatively comparable to AA simulations while offering the computational advantages of coarse-graining. This model opens the door to challenging applications where high accuracy is required only in parts of large biomolecular complexes.

Growth-incompetent monomers of human calcitonin lead to a noncanonical direct relationship between peptide concentration and aggregation lag time

The role of the peptide hormone calcitonin in skeletal protection has led to its use as a therapeutic for osteoporosis. However, calcitonin aggregation into amyloid fibrils limits its therapeutic efficacy, necessitating a modification of calcitonin's aggregation kinetics. Here, we report a direct relationship between human calcitonin (hCT) concentration and aggregation lag time. This kinetic trend was contrary to the conventional understanding of amyloid aggregation and persisted over a range of aggregation conditions, as confirmed by thioflavin-T kinetics assays, CD spectroscopy, and transmission EM. Dynamic light scattering, 1H NMR experiments, and seeded thioflavin-T assay results indicated that differences in initial peptide species contribute to this trend more than variations in the primary nucleus formation rate. On the basis of kinetics modeling results, we propose a mechanism whereby a structural conversion of hCT monomers is needed before incorporation into the fibril. Our kinetic mechanism recapitulates the experimentally observed relationship between peptide concentration and lag time and represents a novel mechanism in amyloid aggregation. Interestingly, hCT at low pH and salmon calcitonin (sCT) exhibited the canonical inverse relationship between concentration and lag time. Comparative studies of hCT and sCT with molecular dynamics simulations and CD indicated an increased α-helical structure in sCT and low-pH hCT monomers compared with neutral-pH hCT, suggesting that α-helical monomers represent a growth-competent species, whereas unstructured random coil monomers represent a growth-incompetent species. Our finding that initial monomer concentration is positively correlated with lag time in hCT aggregation could help inform future efforts for improving therapeutic applications of CT.

CHARMM-GUI 10 years for biomolecular modeling and simulation

CHARMM-GUI, http://www.charmm-gui.org, is a web-based graphical user interface that prepares complex biomolecular systems for molecular simulations. CHARMM-GUI creates input files for a number of programs including CHARMM, NAMD, GROMACS, AMBER, GENESIS, LAMMPS, Desmond, OpenMM, and CHARMM/OpenMM. Since its original development in 2006, CHARMM-GUI has been widely adopted for various purposes and now contains a number of different modules designed to set up a broad range of simulations: (1) PDB Reader & Manipulator, Glycan Reader, and Ligand Reader & Modeler for reading and modifying molecules; (2) Quick MD Simulator, Membrane Builder, Nanodisc Builder, HMMM Builder, Monolayer Builder, Micelle Builder, and Hex Phase Builder for building all-atom simulation systems in various environments; (3) PACE CG Builder and Martini Maker for building coarse-grained simulation systems; (4) DEER Facilitator and MDFF/xMDFF Utilizer for experimentally guided simulations; (5) Implicit Solvent Modeler, PBEQ-Solver, and GCMC/BD Ion Simulator for implicit solvent related calculations; (6) Ligand Binder for ligand solvation and binding free energy simulations; and (7) Drude Prepper for preparation of simulations with the CHARMM Drude polarizable force field. Recently, new modules have been integrated into CHARMM-GUI, such as Glycolipid Modeler for generation of various glycolipid structures, and LPS Modeler for generation of lipopolysaccharide structures from various Gram-negative bacteria. These new features together with existing modules are expected to facilitate advanced molecular modeling and simulation thereby leading to an improved understanding of the structure and dynamics of complex biomolecular systems. Here, we briefly review these capabilities and discuss potential future directions in the CHARMM-GUI development project. © 2016 Wiley Periodicals, Inc.

Efficient gaussian density formulation of volume and surface areas of macromolecules on graphical processing units

We present an algorithm to efficiently compute accurate volumes and surface areas of macromolecules on graphical processing unit (GPU) devices using an analytic model which represents atomic volumes by continuous Gaussian densities. The volume of the molecule is expressed by means of the inclusion-exclusion formula, which is based on the summation of overlap integrals among multiple atomic densities. The surface area of the molecule is obtained by differentiation of the molecular volume with respect to atomic radii. The many-body nature of the model makes a port to GPU devices challenging. To our knowledge, this is the first reported full implementation of this model on GPU hardware. To accomplish this, we have used recursive strategies to construct the tree of overlaps and to accumulate volumes and their gradients on the tree data structures so as to minimize memory contention. The algorithm is used in the formulation of a surface area-based non-polar implicit solvent model implemented as an open source plug-in (named GaussVol) for the popular OpenMM library for molecular mechanics modeling. GaussVol is 50 to 100 times faster than our best optimized implementation for the CPUs, achieving speeds in excess of 100 ns/day with 1 fs time-step for protein-sized systems on commodity GPUs. © 2017 Wiley Periodicals, Inc.

Membrane insertion of fusion peptides from Ebola and Marburg viruses studied by replica-exchange molecular dynamics simulations

This work presents replica-exchange molecular dynamics simulations of inserting a 16-residue Ebola virus fusion peptide into a membrane bilayer. A computational approach is applied for modeling the peptide at the explicit all-atom level and the membrane-aqueous bilayer by a generalized Born continuum model with a smoothed switching function (GBSW). We provide an assessment of the model calculations in terms of three metrics: (1) the ability to reproduce the NMR structure of the peptide determined in the presence of SDS micelles and comparable structural data on other fusion peptides; (2) determination of the effects of the mutation Trp-8 to Ala and sequence discrimination of the homologous Marburg virus; and (3) calculation of potentials of mean force for estimating the partitioning free energy and their comparison to predictions from the Wimley-White interfacial hydrophobicity scale. We found the GBSW implicit membrane model to produce results of limited accuracy in conformational properties of the peptide when compared to the NMR structure, yet the model resolution is sufficient to determine the effect of sequence differentiation on peptide-membrane integration. © 2016 Wiley Periodicals, Inc.

Coarse-grained models using local-density potentials optimized with the relative entropy: Application to implicit solvation

Bottom-up multiscale techniques are frequently used to develop coarse-grained (CG) models for simulations at extended length and time scales but are often limited by a compromise between computational efficiency and accuracy. The conventional approach to CG nonbonded interactions uses pair potentials which, while computationally efficient, can neglect the inherently multibody contributions of the local environment of a site to its energy, due to degrees of freedom that were coarse-grained out. This effect often causes the CG potential to depend strongly on the overall system density, composition, or other properties, which limits its transferability to states other than the one at which it was parameterized. Here, we propose to incorporate multibody effects into CG potentials through additional nonbonded terms, beyond pair interactions, that depend in a mean-field manner on local densities of different atomic species. This approach is analogous to embedded atom and bond-order models that seek to capture multibody electronic effects in metallic systems. We show that the relative entropy coarse-graining framework offers a systematic route to parameterizing such local density potentials. We then characterize this approach in the development of implicit solvation strategies for interactions between model hydrophobes in an aqueous environment.

Accuracy comparison of several common implicit solvent models and their implementations in the context of protein-ligand binding

In this study several commonly used implicit solvent models are compared with respect to their accuracy of estimating solvation energies of small molecules and proteins, as well as desolvation penalty in protein-ligand binding. The test set consists of 19 small proteins, 104 small molecules, and 15 protein-ligand complexes. We compared predicted hydration energies of small molecules with their experimental values; the results of the solvation and desolvation energy calculations for small molecules, proteins and protein-ligand complexes in water were also compared with Thermodynamic Integration calculations based on TIP3P water model and Amber12 force field. The following implicit solvent (water) models considered here are: PCM (Polarized Continuum Model implemented in DISOLV and MCBHSOLV programs), GB (Generalized Born method implemented in DISOLV program, S-GB, and GBNSR6 stand-alone version), COSMO (COnductor-like Screening Model implemented in the DISOLV program and the MOPAC package) and the Poisson-Boltzmann model (implemented in the APBS program). Different parameterizations of the molecules were examined: we compared MMFF94 force field, Amber12 force field and the quantum-chemical semi-empirical PM7 method implemented in the MOPAC package. For small molecules, all of the implicit solvent models tested here yield high correlation coefficients (0.87-0.93) between the calculated solvation energies and the experimental values of hydration energies. For small molecules high correlation (0.82-0.97) with the explicit solvent energies is seen as well. On the other hand, estimated protein solvation energies and protein-ligand binding desolvation energies show substantial discrepancy (up to 10kcal/mol) with the explicit solvent reference. The correlation of polar protein solvation energies and protein-ligand desolvation energies with the corresponding explicit solvent results is 0.65-0.99 and 0.76-0.96 respectively, though this difference in correlations is caused more by different parameterization and less by methods and indicates the need for further improvement of implicit solvent models parameterization. Within the same parameterization, various implicit methods give practically the same correlation with results obtained in explicit solvent model for ligands and proteins: e.g. correlation values of polar ligand solvation energies and the corresponding energies in the frame of explicit solvent were 0.953-0.966 for the APBS program, the GBNSR6 program and all models used in the DISOLV program. The DISOLV program proved to be on a par with the other used programs in the case of proteins and ligands solvation energy calculation. However, the solution of the Poisson-Boltzmann equation (APBS program) and Generalized Born method (implemented in the GBNSR6 program) proved to be the most accurate in calculating the desolvation energies of complexes.

Efficient implementation of constant pH molecular dynamics on modern graphics processors

The treatment of pH sensitive ionization states for titratable residues in proteins is often omitted from molecular dynamics (MD) simulations. While static charge models can answer many questions regarding protein conformational equilibrium and protein-ligand interactions, pH-sensitive phenomena such as acid-activated chaperones and amyloidogenic protein aggregation are inaccessible to such models. Constant pH molecular dynamics (CPHMD) coupled with the Generalized Born with a Simple sWitching function (GBSW) implicit solvent model provide an accurate framework for simulating pH sensitive processes in biological systems. Although this combination has demonstrated success in predicting pKa values of protein structures, and in exploring dynamics of ionizable side-chains, its speed has been an impediment to routine application. The recent availability of low-cost graphics processing unit (GPU) chipsets with thousands of processing cores, together with the implementation of the accurate GBSW implicit solvent model on those chipsets (Arthur and Brooks, J. Comput. Chem. 2016, 37, 927), provide an opportunity to improve the speed of CPHMD and ionization modeling greatly. Here, we present a first implementation of GPU-enabled CPHMD within the CHARMM-OpenMM simulation package interface. Depending on the system size and nonbonded force cutoff parameters, we find speed increases of between one and three orders of magnitude. Additionally, the algorithm scales better with system size than the CPU-based algorithm, thus allowing for larger systems to be modeled in a cost effective manner. We anticipate that the improved performance of this methodology will open the door for broad-spread application of CPHMD in its modeling pH-mediated biological processes. © 2016 Wiley Periodicals, Inc.