SIMS: computation of a smooth invariant molecular surface

All-atom simulations of bent liquid crystal dimers: the twist-bend nematic phase and insights into conformational chirality

The liquid crystal dimer 1,7-bis-4-(4'-cyanobiphenyl)heptane (CB7CB) is known to exhibit a nematic-nematic phase transition, with the lower temperature phase identified as the twist-bend nematic (N_TB) phase. Despite the achiral nature of the mesogen, the N_TB phase demonstrates emergent chirality through the spontaneous formation of a helical structure. We present extensive molecular dynamics simulations of CB7CB using an all-atom force field. The N_TB phase is observed in this model and, upon heating, shows phase transitions into the nematic (N) and isotropic phases. The simulated N_TB phase returns a pitch of 8.35 nm and a conical tilt angle of 29°. Analysis of the bend angle between the mesogenic units reveals an average angle of 127°, which is invariant to the simulated phase. We have calculated distributions of the chirality order parameter, χ, for the ensemble of conformers in the N_TB and N phases. These distributions elucidate that CB7CB is statistically achiral but can adopt chiral conformers with no preference for a specific handedness. Furthermore, there is no change in the extent of conformational chirality between the N_TB and N phases. Using single-molecule stochastic dynamics simulations in the gas phase, we study the dimer series CBnCB (where n = 6, 7, 8 or 9) and CBX(CH₂)₅YCB (where X/Y = CH₂, O or S) in terms of the bend angle and conformational chirality. We confirm that the bent molecular shape determines the ability of a dimer to exhibit the N_TB phase rather than its potential to assume chiral conformers; as |χ|_max increases with the spacer length, but the even-membered dimers have a linear shape in contrast to the bent nature of dimers with spacers of odd parity. For CBX(CH₂)₅YCB, it is found that |χ|_max increases as the bend angle of the dimer decreases, while the flexibility of the dimers remains unchanged through the series.

An effective molecular blocker of ion channel of M2 protein as anti-influenza a drug

Design of a drug compound that can effectively bind to the M2 ion channel and block the diffusion of hydrogen ions (H⁺) through and inhibit influenza A virus replication is an important task. Known anti-influenza drugs amantadine and rimantadine have a weak effect on influenza A virus. A new class of positively charged, +2 e.u., molecules is proposed here to block diffusion of H⁺ ion through the M2 channel. Several drug candidates, derivatives of a lead compound (diazabicyclooctane), were proposed and investigated. Molecular dynamics of thermal fluctuations of M2 protein structure and ionization-conformation coupling of all the ionizable residues were simulated at physiological pH. The influence of the most probable mutations of key drug-binding amino acid residues in the M2 ion channel were investigated too. It is shown that the suggested new blocker has high binding affinity for the M2 channel. There are two in-channel binding sites of high affinity, the first one has H-bonds with two of four serine residues Ser-31A (B) or Ser-31C(D), and the second one has H-bonds with two of four histidine residues His-37A (B), or His-37C(D). The main advantage of the new drug molecule is the positive charge, +2 e.u., which creates a positive electrostatic potential barrier (in addition to a steric one) for a transfer of H⁺ ion through M2 channel and may serve as an effective anti-influenza A virus drug.Communicated by Ramaswamy H. Sarma.

Molecular Graphics: Bridging Structural Biologists and Computer Scientists

Visualization of molecular structures is one of the most common tasks carried out by structural biologists, typically using software, such as Chimera, COOT, PyMOL, or VMD. In this Perspective article, we outline how past developments in computer graphics and data visualization have expanded the understanding of biomolecular function, and we summarize recent advances that promise to further transform structural biology. We also highlight how progress in molecular graphics has been impeded by communication barriers between two communities: the computer scientists driving these advances, and the structural and computational biologists who stand to benefit. By pointing to canonical papers and explaining technical progress underlying new graphical developments in simple terms, we aim to improve communication between these communities; this, in turn, would help shape future developments in molecular graphics.

DG-GL: Differential geometry-based geometric learning of molecular datasets

MOTIVATION

Despite its great success in various physical modeling, differential geometry (DG) has rarely been devised as a versatile tool for analyzing large, diverse, and complex molecular and biomolecular datasets because of the limited understanding of its potential power in dimensionality reduction and its ability to encode essential chemical and biological information in differentiable manifolds.

RESULTS

We put forward a differential geometry-based geometric learning (DG-GL) hypothesis that the intrinsic physics of three-dimensional (3D) molecular structures lies on a family of low-dimensional manifolds embedded in a high-dimensional data space. We encode crucial chemical, physical, and biological information into 2D element interactive manifolds, extracted from a high-dimensional structural data space via a multiscale discrete-to-continuum mapping using differentiable density estimators. Differential geometry apparatuses are utilized to construct element interactive curvatures in analytical forms for certain analytically differentiable density estimators. These low-dimensional differential geometry representations are paired with a robust machine learning algorithm to showcase their descriptive and predictive powers for large, diverse, and complex molecular and biomolecular datasets. Extensive numerical experiments are carried out to demonstrate that the proposed DG-GL strategy outperforms other advanced methods in the predictions of drug discovery-related protein-ligand binding affinity, drug toxicity, and molecular solvation free energy.

AVAILABILITY AND IMPLEMENTATION

http://weilab.math.msu.edu/DG-GL/ Contact: wei@math.msu.edu.

Frontiers in biomolecular mesh generation and molecular visualization systems

With the development of biomolecular modeling and simulation, especially implicit solvent modeling, higher requirements are set for the stability, efficiency and mesh quality of molecular mesh generation software. In this review, we summarize the recent works in biomolecular mesh generation and molecular visualization. First, we introduce various definitions of molecular surface and corresponding meshing software. Second, as the mesh quality significantly influences biomolecular simulation, we investigate some remeshing methods in the fields of computer graphics and molecular modeling. Then, we show the application of biomolecular mesh in the boundary element method (BEM) and the finite element method (FEM). Finally, to conveniently visualize the numerical results based on the mesh, we present two types of molecular visualization systems.

A comprehensive analysis of the computed tautomer fractions of the imidazole ring of histidines in Loligo vulgaris

A recently introduced electrostatic-based method to determine the pKa values of ionizable residues and fractions of ionized and tautomeric forms of histidine (His) and acid residues in proteins, at a given fixed pH, is applied here to the analysis of a His-rich protein, namely Loligo vulgaris (pdb id 1E1A), a 314-residue all-β protein. The average tautomeric fractions for the imidazole ring of each of the six histidines in the sequence were computed using an approach that includes, but is not limited to, molecular dynamic simulations coupled with calculations of the ionization states for all 94 ionizable residues of protein 1E1A in water at pH 6.5 and 300 K. The electrostatic-calculated tautomeric fractions of the imidazole ring of His were compared with predictions obtained from an existent NMR-based methodology. Our results indicate that: (i) the averaged electrostatic-based tautomeric predictions for the imidazole ring of all histidines of Loligo vulgaris are dominated by the Nε2-H rather than the Nδ1-H form, although such preferences from the NMR-based methodology are not so well defined; (ii) the computed average absolute difference between the electrostatic- and the NMR-based tautomeric predictions among all six histidines vary among 0% to 17%; (iii) for the His showing the largest fraction of the neutral form (81%), the absolute difference between the NMR- and electrostatic-based computed tautomeric predictions is only 3%; and (iv) the tautomeric predictions for the imidazole ring of His computed with the NMR-based methodology are stable within a certain, well-defined, range of variations of a tautomer-related parameter.

Impact of the green tea ingredient epigallocatechin gallate and a short pentapeptide (Ile-Ile-Ala-Glu-Lys) on the structural organization of mixed micelles and the related uptake of cholesterol

BACKGROUND

High levels of blood cholesterol are conventionally linked to an increased risk of developing cardiovascular disease (Grundy, 1986). Here we examine the molecular mode of action of natural products with known cholesterol-lowering activity, such as for example the green tea ingredient epigallocatechin gallate and a short pentapeptide, Ile-Ile-Ala-Glu-Lys.

METHODS

Molecular Dynamics simulations are used to gain insight into the formation process of mixed micelles and, correspondingly, how active agents epigallocatechin gallate and Ile-Ile-Ala-Glu-Lys could possibly interfere with it.

RESULTS

Self-assembly of physiological micelles occurs on the order of 35-50 ns; most of the structural properties of mixed micelles are unaffected by epigallocatechin gallate or Ile-Ile-Ala-Glu-Lys which integrate into the micellar surface; the diffusive motion of constituting lipids palmitoyl-oleoyl-phosphatidylcholine and cholesterol is significantly down-regulated by both epigallocatechin gallate and Ile-Ile-Ala-Glu-Lys; CONCLUSIONS: The molecular mode of action of natural compounds epigallocatechin gallate and Ile-Ile-Ala-Glu-Lys is a significant down-regulation of the diffusive motion of micellar lipids.

GENERAL SIGNIFICANCE

Natural compounds like the green tea ingredient epigallocatechin gallate and a short pentapeptide, Ile-Ile-Ala-Glu-Lys, lead to a significant down-regulation of the diffusive motion of micellar lipids thereby modulating cholesterol absorption into physiological micelles.

Coupled molecular dynamics and continuum electrostatic method to compute the ionization pKa's of proteins as a function of pH. Test on a large set of proteins

A computational method, to predict the pKa values of the ionizable residues Asp, Glu, His, Tyr, and Lys of proteins, is presented here. Calculation of the electrostatic free-energy of the proteins is based on an efficient version of a continuum dielectric electrostatic model. The conformational flexibility of the protein is taken into account by carrying out molecular dynamics simulations of 10 ns in implicit water. The accuracy of the proposed method of calculation of pKa values is estimated from a test set of experimental pKa data for 297 ionizable residues from 34 proteins. The pKa-prediction test shows that, on average, 57, 86, and 95% of all predictions have an error lower than 0.5, 1.0, and 1.5 pKa units, respectively. This work contributes to our general understanding of the importance of protein flexibility for an accurate computation of pKa, providing critical insight about the significance of the multiple neutral states of acid and histidine residues for pKa-prediction, and may spur significant progress in our effort to develop a fast and accurate electrostatic-based method for pKa-predictions of proteins as a function of pH.

Effect of terminal chain length on the helical twisting power in achiral bent-core molecules doped in a cholesteric liquid crystal

We prepared a homologous series of achiral bent-core (BC) liquid crystals with different terminal alkoxy chain lengths, n (BC-n), and evaluated the helical twisting power (HTP) of the BC-n doped in a cholesteric liquid crystal. The BC-n molecules with longer terminal chains showed larger HTPs. To interpret this striking phenomenon, a stochastic dynamics simulation was performed to determine the distribution of the chirality order parameters (χ) for BC molecules with n = 8-16. The distribution of χ for each simulated conformation varied with n, and the variation tendency was different for molecules with n < 12 and n > 12 despite the linear relationship between HTP and n in the experiment.

Comparison Study of Polar and Nonpolar Contributions to Solvation Free Energy

In this study, we compared the contributions of polar and nonpolar interactions to the solvation free energy of a solute in solvent, which is decomposed into four different terms based on the nature of interactions: (i) electrostatic solvation free energy term counting for the work done to move solute charges from fixed points in some reference environment to their configuration positions in solvent; (ii) solute-solvent van der Waals dispersion interactions; (iii) change on solvent-solvent interactions and solvent entropy due to reorganization of solvent around solute cavity in solvent; and (iv) compensation of electrostatic forces acting on the dielectric surface boundary between solvent and solute. We compared these contributions to each other for a data set of 573 proteins, which were prepared using CHARMM22 and AMBER force fields. In addition, we compared the calculated with experimental hydration free energies for a data set of 642 small molecules, which were prepared using the general AMBER force field. Our results indicated the significance of each term to the total solvation free energy.

ESES: Software for Eulerian solvent excluded surface

Solvent excluded surface (SES) is one of the most popular surface definitions in biophysics and molecular biology. In addition to its usage in biomolecular visualization, it has been widely used in implicit solvent models, in which SES is usually immersed in a Cartesian mesh. Therefore, it is important to construct SESs in the Eulerian representation for biophysical modeling and computation. This work describes a software package called Eulerian solvent excluded surface (ESES) for the generation of accurate SESs in Cartesian grids. ESES offers the description of the solvent and solute domains by specifying all the intersection points between the SES and the Cartesian grid lines. Additionally, the interface normal at each intersection point is evaluated. Furthermore, for a given biomolecule, the ESES software not only provides the whole surface area, but also partitions the surface area according to atomic types. Homology theory is utilized to detect topological features, such as loops and cavities, on the complex formed by the SES. The sizes of loops and cavities are measured based on persistent homology with an evolutionary partial differential equation-based filtration. ESES is extensively validated by surface visualization, electrostatic solvation free energy computation, surface area and volume calculations, and loop and cavity detection and their size estimation. We used the Amber PBSA test set in our electrostatic solvation energy, area, and volume validations. Our results are either calibrated by analytical values or compared with those from the MSMS software. © 2017 Wiley Periodicals, Inc.

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.

Effects of mono- and divalent metal ions on DNA binding and catalysis of human apurinic/apyrimidinic endonuclease 1

Here, we used stopped-flow fluorescence techniques to conduct a comparative kinetic analysis of the conformational transitions in human apurinic/apyrimidinic endonuclease 1 (APE1) and in DNA containing an abasic site in the course of their interaction. Effects of monovalent (K(+)) and divalent (Mg(2+), Mn(2+), Ca(2+), Zn(2+), Cu(2+), and Ni(2+)) metal ions on DNA binding and catalytic stages were studied. It was shown that the first step of substrate binding (corresponding to formation of a primary enzyme-substrate complex) does not depend on the concentration (0.05-5.0 mM) or the nature of divalent metal ions. In contrast, the initial DNA binding efficiency significantly decreased at a high concentration (5-250 mM) of monovalent K(+) ions, indicating the involvement of electrostatic interactions in this stage. It was also shown that Cu(2+) ions abrogated the DNA binding ability of APE1, possibly, due to a strong interaction with DNA bases and the sugar-phosphate backbone. In the case of Ca(2+) ions, the catalytic activity of APE1 was lost completely with retention of binding potential. Thus, the enzymatic activity of APE1 is increased in the order Zn(2+) < Ni(2+) < Mn(2+) < Mg(2+). Circular dichroism spectra and calculation of the contact area between APE1 and DNA reveal that Mg(2+) ions stabilize the protein structure and the enzyme-substrate complex.

Parameterization for molecular Gaussian surface and a comparison study of surface mesh generation

The molecular Gaussian surface has been frequently used in the field of molecular modeling and simulation. Typically, the Gaussian surface is defined using two controlling parameters; the decay rate and isovalue. Currently, there is a lack of studies in which a systematic approach in the determination of optimal parameterization according to the geometric features has been done. In this paper, surface area, volume enclosed by the surface and Hausdorff distance are used as three criteria for the parameterization to make the Gaussian surface approximate the solvent excluded surface (SES) well. For each of these three criteria, a search of the parameter space is carried out in order to determine the optimal parameter values. The resulted parameters are close to each other and result in similar calculated molecular properties. Approximation of the VDW surface is also done by analyzing the explicit expressions of the Gaussian surface and VDW surface, which analysis and parameters can be similarly applied to the solvent accessible surface (SAS) due to its geometric similarity to the VDW surface. Once the optimal parameters are obtained, we compare the performance of our Gaussian surface generation software TMSmesh with other commonly used software programs, focusing primarily on mesh quality and fidelity. Additionally, the Poisson-Boltzmann solvation energies based on the surface meshes generated by TMSmesh and those generated by other software programs are calculated and compared for a set of molecules with different sizes. The results of these comparisons validate both the accuracy and the applicability of the parameterized Gaussian surface.

Thermodynamics of the DNA damage repair steps of human 8-oxoguanine DNA glycosylase

Human 8-oxoguanine DNA glycosylase (hOGG1) is a key enzyme responsible for initiating the base excision repair of 7,8-dihydro-8-oxoguanosine (oxoG). In this study a thermodynamic analysis of the interaction of hOGG1 with specific and non-specific DNA-substrates is performed based on stopped-flow kinetic data. The standard Gibbs energies, enthalpies and entropies of specific stages of the repair process were determined via kinetic measurements over a temperature range using the van’t Hoff approach. The three steps which are accompanied with changes in the DNA conformations were detected via 2-aminopurine fluorescence in the process of binding and recognition of damaged oxoG base by hOGG1. The thermodynamic analysis has demonstrated that the initial step of the DNA substrates binding is mainly governed by energy due to favorable interactions in the process of formation of the recognition contacts, which results in negative enthalpy change, as well as due to partial desolvation of the surface between the DNA and enzyme, which results in positive entropy change. Discrimination of non-specific G base versus specific oxoG base is occurring in the second step of the oxoG-substrate binding. This step requires energy consumption which is compensated by the positive entropy contribution. The third binding step is the final adjustment of the enzyme/substrate complex to achieve the catalytically competent state which is characterized by large endothermicity compensated by a significant increase of entropy originated from the dehydration of the DNA grooves.

Between algorithm and model: different Molecular Surface definitions for the Poisson-Boltzmann based electrostatic characterization of biomolecules in solution

The definition of a molecular surface which is physically sound and computationally efficient is a very interesting and long standing problem in the implicit solvent continuum modeling of biomolecular systems as well as in the molecular graphics field. In this work, two molecular surfaces are evaluated with respect to their suitability for electrostatic computation as alternatives to the widely used Connolly-Richards surface: the blobby surface, an implicit Gaussian atom centered surface, and the skin surface. As figures of merit, we considered surface differentiability and surface area continuity with respect to atom positions, and the agreement with explicit solvent simulations. Geometric analysis seems to privilege the skin to the blobby surface, and points to an unexpected relationship between the non connectedness of the surface, caused by interstices in the solute volume, and the surface area dependence on atomic centers. In order to assess the ability to reproduce explicit solvent results, specific software tools have been developed to enable the use of the skin surface in Poisson-Boltzmann calculations with the DelPhi solver. Results indicate that the skin and Connolly surfaces have a comparable performance from this last point of view.

Triangulated manifold meshing method preserving molecular surface topology

Generation of manifold mesh is an urgent issue in mathematical simulations of biomolecule using boundary element methods (BEM) or finite element method (FEM). Defects, such as not closed mesh, intersection of elements and missing of small structures, exist in surface meshes generated by most of the current meshing method. Usually the molecular surface meshes produced by existing methods need to be revised carefully by third party software to ensure the surface represents a continuous manifold before being used in a BEM and FEM calculations. Based on the trace technique proposed in our previous work, in this paper, we present an improved meshing method to avoid intersections and preserve the topology of the molecular Gaussian surface. The new method divides the whole Gaussian surface into single valued pieces along each of x, y, z directions by tracing the extreme points along the fold curves on the surface. Numerical test results show that the surface meshes produced by the new method are manifolds and preserve surface topologies. The result surface mesh can also be directly used in surface conforming volume mesh generation for FEM type simulation.

Reducing grid-dependence in finite-difference Poisson-Boltzmann calculations

Grid dependence in numerical reaction field energies and solvation forces is a well-known limitation in the finite-difference Poisson-Boltzmann methods. In this study we have investigated several numerical strategies to overcome the limitation. Specifically, we have included trimer arc dots during analytical molecular surface generation to improve the convergence of numerical reaction field energies and solvation forces. We have also utilized the level set function to trace the molecular surface implicitly to simplify the numerical mapping of the grid-independent solvent excluded surface. We have further explored to combine the weighted harmonic averaging of boundary dielectrics with a charge-based approach to improve the convergence and stability of numerical reaction field energies and solvation forces. Our test data show that the convergence and stability in both numerical energies and forces can be improved significantly when the combined strategy is applied to either the Poisson equation or the full Poisson-Boltzmann equation.

Thermodynamics of the multi-stage DNA lesion recognition and repair by formamidopyrimidine-DNA glycosylase using pyrrolocytosine fluorescence--stopped-flow pre-steady-state kinetics

Formamidopyrimidine-DNA glycosylase, Fpg protein from Escherichia coli, initiates base excision repair in DNA by removing a wide variety of oxidized lesions. In this study, we perform thermodynamic analysis of the multi-stage interaction of Fpg with specific DNA-substrates containing 7,8-dihydro-8-oxoguanosine (oxoG), or tetrahydrofuran (THF, an uncleavable abasic site analog) and non-specific (G) DNA-ligand based on stopped-flow kinetic data. Pyrrolocytosine, highly fluorescent analog of the natural nucleobase cytosine, is used to record multi-stage DNA lesion recognition and repair kinetics over a temperature range (10–30°C). The kinetic data were used to obtain the standard Gibbs energy, enthalpy and entropy of the specific stages using van’t Hoff approach. The data suggest that not only enthalpy-driven exothermic oxoG recognition, but also the desolvation-accompanied entropy-driven enzyme-substrate complex adjustment into the catalytically active state play equally important roles in the overall process.

Design of [(2-pyrimidinylthio)acetyl]benzenesulfonamides as inhibitors of human carbonic anhydrases

A series of [(2-pyrimidinylthio)acetyl]benzenesulfonamides were designed and synthesized. Their binding affinities as inhibitors of several recombinant human carbonic anhydrase (CA) isozymes were determined by isothermal titration calorimetry (ITC) and thermal shift assay (TSA). A group of compounds containing a chlorine atom in the benzenesulfonamide ring were found to exhibit higher selectivity but lower binding affinity toward tested CAs. The crystal structures of selected compounds in complex with CA II were determined to atomic resolution. Docking studies were performed to compare the binding modes of experimentally determined crystallographic structures with computational prediction of the pyrimidine derivative binding to CA II. Several compounds bound to select CAs with single-digit nanomolar affinities and could be used as leads for inhibitor development toward a select CA isozyme.

Structural features of aquaporin 4 supporting the formation of arrays and junctions in biomembranes

A limited class of aquaporins has been described to form regular arrays and junctions in membranes. The biological significance of these structures, however, remains uncertain. Here we analyze the underlying physical principles with the help of a computational procedure that takes into account protein-protein as well as protein-membrane interactions. Experimentally observed array/junction structures are systematically (dis)assembled and major driving forces identified. Aquaporin 4 was found to be markedly different from the non-junction forming aquaporin 1. The environmental stabilization resulting from embedding into the biomembrane was identified as the main driving force. This highlights the role of protein-membrane interactions in aquaporin 4. Analysis of the type presented here can help to decipher the biological role of membrane arrays and junctions formed by aquaporin.

Potential of mean force of water-proton bath and molecular dynamic simulation of proteins at constant pH

An advanced implicit solvent model of water-proton bath for protein simulations at constant pH is presented. The implicit water-proton bath model approximates the potential of mean force of a protein in water solvent in a presence of hydrogen ions. Accurate and fast computational implementation of the implicit water-proton bath model is developed using the continuum electrostatic Poisson equation model for calculation of ionization equilibrium and the corrected MSR6 generalized Born model for calculation of the electrostatic atom-atom interactions and forces. Molecular dynamics (MD) method for protein simulation in the potential of mean force of water-proton bath is developed and tested on three proteins. The model allows to run MD simulations of proteins at constant pH, to calculate pH-dependent properties and free energies of protein conformations. The obtained results indicate that the developed implicit model of water-proton bath provides an efficient way to study thermodynamics of biomolecular systems as a function of pH, pH-dependent ionization-conformation coupling, and proton transfer events.

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.

A proposal for the revision of molecular boundary typology

Molecular shape is essential in understanding molecular function, and understanding molecular shape requires definition of molecular boundaries. In this paper, we review the conceptual evolution of three molecular boundary types: the van der Waals surface, the Connolly surface, and the Lee-Richards (accessible) surface. Then, we point out the confusion among the names of these surfaces existing in the literature. Since it is desirable to have a well-defined terminology in a discipline, we propose the standard names of the three molecular boundary types and their corresponding volumes in order to maximize consistency among researchers, respect the first individual who defined or computed a surface type, and promote collaboration between biologists and geometers.

Blind docking method combining search of low-resolution binding sites with ligand pose refinement by molecular dynamics-based global optimization

This study describes the development of a new blind hierarchical docking method, bhDock, its implementation, and accuracy assessment. The bhDock method uses two-step algorithm. First, a comprehensive set of low-resolution binding sites is determined by analyzing entire protein surface and ranked by a simple score function. Second, ligand position is determined via a molecular dynamics-based method of global optimization starting from a small set of high ranked low-resolution binding sites. The refinement of the ligand binding pose starts from uniformly distributed multiple initial ligand orientations and uses simulated annealing molecular dynamics coupled with guided force-field deformation of protein-ligand interactions to find the global minimum. Assessment of the bhDock method on the set of 37 protein-ligand complexes has shown the success rate of predictions of 78%, which is better than the rate reported for the most cited docking methods, such as AutoDock, DOCK, GOLD, and FlexX, on the same set of complexes.

Current performance gains from utilizing the GPU or the ASIC MDGRAPE-3 within an enhanced Poisson Boltzmann approach

Scientific applications do frequently suffer from limited compute performance. In this article, we investigate the suitability of specialized computer chips to overcome this limitation. An enhanced Poisson Boltzmann program is ported to the graphics processing unit and the application specific integrated circuit MDGRAPE-3 and resulting execution times are compared to the conventional performance obtained on a modern central processing unit. Speed Up factors are measured and an analysis of numerical accuracy is provided. On both specialized architectures the improvement is increasing with problem size and reaches up to a Speed Up factor of 39 x for the largest problem studied. This type of alternative high performance computing can significantly improve the performance of demanding scientific applications.

Proton-linked dimerization of a retroviral capsid protein initiates capsid assembly

In mature retroviral particles, the capsid protein (CA) forms a shell encasing the viral replication complex. Human immunodeficiency virus (HIV) CA dimerizes in solution, through its C-terminal domain (CTD), and this interaction is important for capsid assembly. In contrast, other retroviral capsid proteins, including that of Rous sarcoma virus (RSV), do not dimerize with measurable affinity. Here we show, using X-ray crystallography and other biophysical methods, that acidification causes RSV CA to dimerize in a fashion analogous to HIV CA, and that this drives capsid assembly in vitro. A pair of aspartic acid residues, located within the CTD dimer interface, explains why dimerization is linked to proton binding. Our results show that despite overarching structural similarities, the intermolecular forces responsible for forming and stabilizing the retroviral capsid differ markedly across retroviral genera. Our data further suggest that proton binding may regulate RSV capsid assembly, or modulate stability of the assembled capsid.

Discretization of the induced-charge boundary integral equation

Boundary-element methods (BEMs) for solving integral equations numerically have been used in many fields to compute the induced charges at dielectric boundaries. In this paper, we consider a more accurate implementation of BEM in the context of ions in aqueous solution near proteins, but our results are applicable more generally. The ions that modulate protein function are often within a few angstroms of the protein, which leads to the significant accumulation of polarization charge at the protein-solvent interface. Computing the induced charge accurately and quickly poses a numerical challenge in solving a popular integral equation using BEM. In particular, the accuracy of simulations can depend strongly on seemingly minor details of how the entries of the BEM matrix are calculated. We demonstrate that when the dielectric interface is discretized into flat tiles, the qualocation method of Tausch [IEEE Trans Comput.-Comput.-Aided Des. 20, 1398 (2001)] to compute the BEM matrix elements is always more accurate than the traditional centroid-collocation method. Qualocation is not more expensive to implement than collocation and can save significant computational time by reducing the number of boundary elements needed to discretize the dielectric interfaces.

FAMBE-pH: a fast and accurate method to compute the total solvation free energies of proteins

A fast and accurate method to compute the total solvation free energies of proteins as a function of pH is presented. The method makes use of a combination of approaches, some of which have already appeared in the literature; (i) the Poisson equation is solved with an optimized fast adaptive multigrid boundary element (FAMBE) method; (ii) the electrostatic free energies of the ionizable sites are calculated for their neutral and charged states by using a detailed model of atomic charges; (iii) a set of optimal atomic radii is used to define a precise dielectric surface interface; (iv) a multilevel adaptive tessellation of this dielectric surface interface is achieved by using multisized boundary elements; and (v) 1:1 salt effects are included. The equilibrium proton binding/release is calculated with the Tanford-Schellman integral if the proteins contain more than approximately 20-25 ionizable groups; for a smaller number of ionizable groups, the ionization partition function is calculated directly. The FAMBE method is tested as a function of pH (FAMBE-pH) with three proteins, namely, bovine pancreatic trypsin inhibitor (BPTI), hen egg white lysozyme (HEWL), and bovine pancreatic ribonuclease A (RNaseA). The results are (a) the FAMBE-pH method reproduces the observed pK a's of the ionizable groups of these proteins within an average absolute value of 0.4 p K units and a maximum error of 1.2 p K units and (b) comparison of the calculated total pH-dependent solvation free energy for BPTI, between the exact calculation of the ionization partition function and the Tanford-Schellman integral method, shows agreement within 1.2 kcal/mol. These results indicate that calculation of total solvation free energies with the FAMBE-pH method can provide an accurate prediction of protein conformational stability at a given fixed pH and, if coupled with molecular mechanics or molecular dynamics methods, can also be used for more realistic studies of protein folding, unfolding, and dynamics, as a function of pH.

Dispersion terms and analysis of size- and charge dependence in an enhanced Poisson-Boltzmann approach

We implement a well-established concept to consider dispersion effects within a Poisson-Boltzmann approach of continuum solvation of proteins. The theoretical framework is particularly suited for boundary element methods. Free parameters are determined by comparison to experimental data as well as high-level quantum mechanical reference calculations. The method is general and can be easily extended in several directions. The model is tested on various chemical substances and found to yield good-quality estimates of the solvation free energy without obvious indication of any introduced bias. Once optimized, the model is applied to a series of proteins, and factors such as protein size or partial charge assignments are studied.

Systematic study of the boundary composition in Poisson Boltzmann calculations

We describe a three-stage procedure to analyze the dependence of Poisson Boltzmann calculations on the shape, size and geometry of the boundary between solute and solvent. Our study is carried out within the boundary element formalism, but our results are also of interest to finite difference techniques of Poisson Boltzmann calculations. At first, we identify the critical size of the geometrical elements for discretizing the boundary, and thus the necessary resolution required to establish numerical convergence. In the following two steps we perform reference calculations on a set of dipeptides in different conformations using the Polarizable Continuum Model and a high-level Density Functional as well as a high-quality basis set. Afterwards, we propose a mechanism for defining appropriate boundary geometries. Finally, we compare the classic Poisson Boltzmann description with the Quantum Chemical description, and aim at finding appropriate fitting parameters to get a close match to the reference data. Surprisingly, when using default AMBER partial charges and the rigorous geometric parameters derived in the initial two stages, no scaling of the partial charges is necessary and the best fit against the reference set is obtained automatically.

Molecular surface generation using a variable-radius solvent probe

Protein-ligand binding occurs through interactions at the molecular surface. Hence, a proper description of this surface is essential to our understanding of the process of molecular recognition. Recent studies have noted the inadequacy of using a fixed 1.4 A solvent probe radius to generate the molecular surface. This assumes that water molecules approach all surface atoms at an equal distance irrespective of polarity, which is not the case. To adequately model the protein-water boundary requires that the solvent probe radius change according to the polarity of its contacting atoms, smaller near polar atoms and larger near apolar atoms. To our knowledge, no method currently exists to generate the molecular surface of a protein in this manner. Using a modification of the marching tetrahedra algorithm, we present a method to generate molecular surfaces using a variable radius solvent probe. The resulting surface lacks many of the unrealistic small crevices in nonpolar regions that are found when utilizing an invariant 1.4 A solvent probe, while maintaining the fine detail of the surface at polar regions. On application of the method on a test set of 20 protein structures taken from the Protein Data Bank (PDB), we also find far fewer empty unsolvated cavities that are present when using only a 1.4 A solvent probe, while the majority of solvated and polar cavities is retained. This suggests that the majority of empty cavities previously observed in protein structures might simply be artifacts of the surfacing method. We also find that the variable probe surface can have significant effects on electrostatic calculations by generating a better tuned description of the protein-water boundary. We also examined the binding interfaces of a diverse set of 55 protein-protein complexes. We find that using a variable probe results in an increase in perceived shape complementarity at these sites compared to using a 1.4 A solvent probe. The molecular volume and surface area are geometric values that determine various important properties for macromolecules, and the altered description afforded by a variable solvent probe molecular surface can have significant implications in protein recognition, energetics, folding, and stability calculations.

A Single Disulfide Bond Differentiates Aggregation Pathways of ß2-Microglobulin

Deposition of wild-type beta2-microglobulin (beta2m) into amyloid fibrils is a complication in patients undergoing long-term hemodialysis. The native beta-sandwich fold of beta2m has a highly conserved disulfide bond linking Cys25 and Cys80. Oxidized beta2m forms needle-like amyloid fibrils at pH 2.5 in vitro, whereas reduced beta2m, at acid pH, in which the intra-chain disulfide bond is disrupted, cannot form typical fibrils. Instead, reduced beta2m forms thinner and more flexible filaments. To uncover the difference in molecular mechanisms underlying the aggregation of the oxidized and reduced beta2m, we performed molecular dynamics simulations of beta2m oligomerization under oxidized and reduced conditions. We show that, consistent with experimental observations, the oxidized beta2m forms domain-swapped dimer, in which the two proteins exchange their N-terminal segments complementing each other. In contrast, both dimers and trimers, formed by reduced beta2m, are comprised of parallel beta-sheets between monomers and stabilized by the hydrogen bond network along the backbone. The oligomerized monomers are in extended conformations, capable of further aggregation. We find that both reduced and oxidized dimers are thermodynamically less stable than their corresponding monomers, indicating that beta2m oligomerization is not accompanied by the formation of a thermodynamically stable dimer. Our studies suggest that the different aggregation pathways of oxidized and reduced beta2m are dictated by the formation of distinct precursor oligomeric species that are modulated by Cys25-Cys80 disulfide-bonds. We propose that the propagation of domain swapping is the aggregation mechanism for the oxidized beta2m, while "parallel stacking" of partially unfolded beta2m is the aggregation mechanism for the reduced beta2m.

Simple models for hydrophobic hydration

This tutorial review compares models that describe DeltaG(cavitation). Their qualitative agreement suggests the use of the simple, time-honored Pierotti equation. Its coefficients, fine-tuned with atomistic simulations, give a revised Pierotti approach, rPA. A discussion of the extension of the rPA model to non-spherical solutes is presented and the different roles of molecular volume and surface area of the solute are brought together. The tutorial review is aimed at experimentalists and theoreticians interested in the description of solvent effects.