Optimization of Gaussian surface calculations and extension to solvent-accessible surface areas

Insight of the molecular mechanism of inhibitors located at different allosteric sites regulating the activity of wild type and mutant KRAS (G12)

As the important hub of many cellular signaling networks, KRAS (Kirsten rat sarcoma viral oncogene homologue) has been identified as a tumor biomarker. It is the frequently mutated oncogene in human cancers, and KRAS protein activation caused by mutations, such as G12D, has been found in many human tumors tissues. Although, there are two specific allosteric sites (AS1 and AS2) on the KRAS protein that can be used as the targets for inhibitor development, the difference of regulatory mechanisms between two individual allosteric sites still not be reported. Here, using molecular dynamics simulations combined with molecular mechanics generalized born surface area (MM/GBSA) analysis, we found that both of the inhibitors, located at AS1 and AS2, were able to reduce the binding free energy between wild type, mutant KRAS (G12/D/V/S/C) and GTP remarkably, however the effect of inhibitors on the binding free energy between wild type, mutant KRAS and GDP was limited. In addition, the degree of decrease of binding free energy between KRAS and GTP caused by inhibitors at AS2 was significantly greater than that caused by inhibitors at AS1. Further analysis revealed that both inhibitors at AS1 and AS2 were able to regulate the fluctuation of Switch Ⅰ and Switch Ⅱ to expand the pocket of the orthosteric site (GTP binding site), thereby reducing the binding of KRAS to GTP. Noteworthy there was significant differences in the regulatory preferences on Switch Ⅰ and Switch Ⅱ between two type inhibitor. The inhibitor at AS2 mainly regulated Switch Ⅱ to affect the pocket of the orthosteric site, while the inhibitor at AS1 mainly expand the pocket of the orthosteric site by regulating the fluctuation of Switch Ⅰ. Our study compared the differences between two type inhibitors in regulating the KRAS protein activity and revealed the advantages of the AS2 as the small molecule drug target, aiming to provide theoretical guidance for the research of novel KRAS protein inhibitors.

An Unfitted Finite Element Poisson-Boltzmann Solver with Automatic Resolving of Curved Molecular Surface

So far, the existing Poisson-Boltzmann (PB) solvers that accurately take into account the interface jump conditions need a pregenerated body-fitted mesh (molecular surface mesh). However, qualified biomolecular surface meshing and its implementation into numerical methods remains a challenging and laborious issue, which practically hinders the progress of further developments and applications of a bunch of numerical methods in this field. In addition, even with a molecular surface mesh, it is only a low-order approximation of the original curved surface. In this article, an interface-penalty finite element method (IPFEM), which is a typical unfitted finite element method, is proposed to solve the Poisson-Boltzmann equation (PBE) without requiring the user to generate a molecular surface mesh. The Gaussian molecular surface is used to represent the molecular surface and can be automatically resolved with a high-order approximation within our method. Theoretical convergence rates of the IPFEM for the linear PB equation have been provided and are well validated on a benchmark problem with an analytical solution (we also noticed from numerical examples that the IPFEM has similar convergence rates for the nonlinear PBE). Numerical results on a set of different-sized biomolecules demonstrate that the IPFEM is numerically stable and accurate in the calculation of biomolecular electrostatic solvation energy.

Primary Nucleation of Polymorphic α-Synuclein Dimers Depends on Copper Concentrations and Definite Copper-Binding Site

The primary nucleation process of α-synuclein (AS) that forms toxic oligomeric species is the early stage of the pathological cause of Parkinson's disease. It is well-known that copper influences this primary nucleation process. While significant efforts have been made to solve the structures of polymorphic AS fibrils, the structures of AS oligomers and the copper-bound AS oligomers at the molecular level and the effect of copper concentrations on the primary nucleation are elusive. Here, we propose and demonstrate new molecular mechanism pathways of primary nucleation of AS that are tuned by distinct copper concentrations and by a specific copper-binding site. We present the polymorphic AS dimers bound to different copper-binding sites at the atomic resolution in high- and low-copper concentrations, using extensive molecular dynamics simulations. Our results show the complexity of the primary nucleation pathways that rely on the copper concentrations and the copper binding site. From a broader perspective, our study proposes a new strategy to control the primary nucleation of other toxic amyloid oligomers in other neurodegenerative diseases.

Synergistic effect of conformational changes in phosphoglycerate kinase 1 product release

In the glycolysis pathway, phosphoglycerate kinase 1 (PGK1) transfers one phosphoryl-group from 1,3-diphosphoglycerate (1,3BPG) to ADP to product 3-phosphoglycerate (3PG) and ATP. The catalytic process is accompanied with the conversion between the open conformation and the closed conformation of PGK1. However, the dynamic collaboration mechanism between the PGK1 conformation transition and the products releasing process remains poorly understood. Here using molecular dynamics simulations combined with molecular mechanics generalized born surface area (MM/GBSA) analysis, we demonstrated that PGK1 in the closed conformation first releases the product ATP to reach a semi-open conformation, and releases the product 3PG to achieve the full open conformation, which could accept new substrates ADP and 1,3BPG for the next cycle. It is noteworthy that the phosphorylation of PGK1 at T243 causes the loop region (residues L248-E260) flip outside the protein, and the phosphorylation of Y324 leads PGK1 become looser. Both modifications cause the exposure of the ADP/ATP binding site, which was beneficial for the substrates/products binding/releasing of PGK1. In addition, the other post translational modifications (PTMs) were also able to regulate the ligands binding/releasing with different effects. Our results revealed the dynamic cooperative molecular mechanism of PGK1 conformational transition with products releasing, as well as the influence of PTMs, which would contribute to the understanding of PGK1 substrates/products conversion process and the development of small molecule drugs targeting PGK1.Communicated by Ramaswamy H. Sarma.

The molecular mechanism of Y473 phosphorylation of UGDH relieves the inhibition effect of UDP-glucose on HuR

Uridine diphosphate glucose (UDP-Glc) is able to accelerate the decay of snail family transcriptional repressor 1 (SNAI1) mRNA by inhibiting Hu antigen R (HuR, an RNA-binding protein), thereby preventing cancer invasiveness and drug resistance. Nevertheless, the phosphorylation of tyrosine 473 (Y473) of UDP-glucose dehydrogenase (UGDH is capable of converting UDP-Glc to uridine diphosphate glucuronic acid (UDP-GlcUA)) weakens the inhibition of UDP-Glc to HuR, thus initiating the epithelial-mesenchymal transformation of tumor cells and promoting tumor cell migration and metastasis. To address the mechanism, we performed molecular dynamics simulations combined with molecular mechanics generalized Born surface area (MM/GBSA) analysis on wild-type and Y473 phosphorylated UGDH and HuR, UDP-Glc, UDP-GlcUA complexes. We demonstrated that Y473 phosphorylation was able to enhance the binding between UGDH and the HuR/UDP-Glc complex. Compared with HuR, UGDH has a stronger binding ability with UDP-Glc; therefore, UDP-Glc was inclined to bind to UGDH and then was catalyzed to UDP-GlcUA by UGDH, which relieved the inhibition of UDP-Glc to HuR. In addition, the binding ability of HuR for UDP-GlcUA was lower than its affinity for UDP-Glc, significantly reducing the inhibition of HuR. Hence, HuR bound to SNAI1 mRNA more easily to increase the stability of mRNA. Our results revealed the micromolecular mechanism of Y473 phosphorylation of UGDH regulating the interaction between UGDH and HuR as well as relieving the inhibition of UDP-Glc on HuR, which contributed to understanding the role of UGDH and HuR in tumor metastasis and developing small molecule drugs targeting the interaction between UGDH and HuR.

Molecular recognition and binding of CcrA from Bacteroides fragilis with cefotaxime and ceftazidime by fluorescence spectra and molecular docking

In search of new super-bacterial inhibitor agents, the recognition and binding mechanism of the B1 subclass MβL CcrA from Bacteroides fragilis with cefotaxime (CTX) and ceftazidime (CAZ) were studied using spectroscopy analysis and molecular docking. The results showed that the fluorescence quenching of CcrA induced by CTX and CAZ were all due to the complex formation, which belonged to static quenching and was forced by hydrogen bonds and Van der Waals forces, despite the greater binding ability of CTX with CcrA than CAZ. Upon recognizing CTX or CAZ, the CcrA opened its binding pocket by the microenvironmental and conformational of three loops changing to promote an induced-fit of the freshly introduced antibiotics. In addition, the whole antibiotic molecule ultimately entered the active pocket of CcrA with its original carbonate replaced by the carboxyl oxygen of the hexatomic ring adjacent to the β-lactam ring in CTX or CAZ, forming a new tetrahedral coordination structure at the Zn2 site. Moreover, the difference in steric hindrance and electrostatic effects of the side chain affected the binding ability of the two antibiotics to the CcrA. This work showed the refined procedures of the antibiotics binding to CcrA and might provide useful information hint for the new strategy of developing the novel and innovative super-bacterial antibiotics.

Design of Thioether Cyclic Peptide Scaffolds with Passive Permeability and Oral Exposure

Advances in the design of permeable peptides and in the synthesis of large arrays of macrocyclic peptides with diverse amino acids have evolved on parallel but independent tracks. Less precedent combines their respective attributes, thereby limiting the potential to identify permeable peptide ligands for key targets. Herein, we present novel 6-, 7-, and 8-mer cyclic peptides (MW 774-1076 g·mol-1) with passive permeability and oral exposure that feature the amino acids and thioether ring-closing common to large array formats, including DNA- and RNA-templated synthesis. Each oral peptide herein, selected from virtual libraries of partially N-methylated peptides using in silico methods, reflects the subset consistent with low energy conformations, low desolvation penalties, and passive permeability. We envision that, by retaining the backbone N-methylation pattern and consequent bias toward permeability, one can generate large peptide arrays with sufficient side chain diversity to identify permeability-biased ligands to a variety of protein targets.

Molecular Sparse Representation by a 3D Ellipsoid Radial Basis Function Neural Network via L1 Regularization

The three-dimensional structures and shapes of biomolecules provide essential information about their interactions and functions. Unfortunately, the computational cost of biomolecular shape representation is an active challenge which increases rapidly as the number of atoms increase. Recent developments in sparse representation and deep learning have shown significant improvements in terms of time and space. A sparse representation of molecular shape is also useful in various other applications, such as molecular structure alignment, docking, and coarse-grained molecular modeling. We have developed an ellipsoid radial basis function neural network (ERBFNN) and an algorithm for sparsely representing molecular shape. To evaluate a sparse representation model of molecular shape, the Gaussian density map of the molecule is approximated using ERBFNN with a relatively small number of neurons. The deep learning models were trained by optimizing a nonlinear loss function with L1 regularization. Experimental results reveal that our algorithm can represent the original molecular shape with a relatively higher accuracy and fewer scale of ERBFNN. Our network in principle is applicable to the multiresolution sparse representation of molecular shape and coarse-grained molecular modeling. Executable files are available at https://github.com/SGUI-LSEC/SparseGaussianMolecule. The program was implemented in PyTorch and was run on Linux.

A grid-based algorithm in conjunction with a gaussian-based model of atoms for describing molecular geometry

A novel grid-based method is presented, which in conjunction with a smooth Gaussian-based model of atoms, is used to compute molecular volume (MV) and surface area (MSA). The MV and MSA are essential for computing nonpolar component of free energies. The objective of our grid-based approach is to identify solute atom pairs that share overlapping volumes in space. Once completed, this information is used to construct a rooted tree using depth-first method to yield the final volume and SA by using the formulations of the Gaussian model described by Grant and Pickup (J. Phys Chem, 1995, 99, 3503). The method is designed to function uninterruptedly with the grid-based finite-difference method implemented in Delphi, a popular and open-source package used for solving the Poisson-Boltzmann equation (PBE). We demonstrate the time efficacy of the method while also validating its performance in terms of the effect of grid-resolution, positioning of the solute within the grid-map and accuracy in identification of overlapping atom pairs. We also explore and discuss different aspects of the Gaussian model with key emphasis on its physical meaningfulness. This development and its future release with the Delphi package are intended to provide a physically meaningful, fast, robust and comprehensive tool for MM/PBSA based free energy calculations. © 2019 Wiley Periodicals, Inc.

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.

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.

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.

Electrostatics of proteins in dielectric solvent continua. II. Hamiltonian reaction field dynamics

In Paper I of this work [S. Bauer, G. Mathias, and P. Tavan, J. Chem. Phys. 140, 104102 (2014)] we have presented a reaction field (RF) method, which accurately solves the Poisson equation for proteins embedded in dielectric solvent continua at a computational effort comparable to that of polarizable molecular mechanics (MM) force fields. Building upon these results, here we suggest a method for linearly scaling Hamiltonian RF/MM molecular dynamics (MD) simulations, which we call "Hamiltonian dielectric solvent" (HADES). First, we derive analytical expressions for the RF forces acting on the solute atoms. These forces properly account for all those conditions, which have to be self-consistently fulfilled by RF quantities introduced in Paper I. Next we provide details on the implementation, i.e., we show how our RF approach is combined with a fast multipole method and how the self-consistency iterations are accelerated by the use of the so-called direct inversion in the iterative subspace. Finally we demonstrate that the method and its implementation enable Hamiltonian, i.e., energy and momentum conserving HADES-MD, and compare in a sample application on Ac-Ala-NHMe the HADES-MD free energy landscape at 300 K with that obtained in Paper I by scanning of configurations and with one obtained from an explicit solvent simulation.

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.

MOLECULAR MODELING OF THE INTERACTIONS BETWEEN HISTONE DEACETYLASE 8 AND INHIBITORS

Inhibitors of histone deacetylases (HDACs) have become an attractive class of anticancer agent. To understand the interaction between HDAC8 and inhibitors, including "pan-" inhibitors that inhibit many HDACs isoforms and selective inhibitors with no linker domain, docking and molecular dynamics simulation were conducted. Docking results showed the presence of π-π interactions between "linkerless" inhibitors and the aromatic amino acid residues of HDAC8 in the active site. Binding between HDAC8 and inhibitors was also stabilized by hydrogen bond and hydrophobic interaction. In molecular dynamics simulations, the zinc ion was shown to coordinate one more atom of HDAC8-"linkerless" inhibitor complexes than HDAC8-"pan-" inhibitor complexes. Persistent hydrogen bonds also existed between Tyr306 of HDAC8 and some inhibitors. When inhibitors with large cap groups bound to the active pocket of HDAC8, Phe152 and Met274 shifted from their initial positions and the entrance of the active pocket became more open, resulting in the formation of sub-pocket. Hydrophobic interactions contributed most favorably to the binding free energy between HDAC8 and inhibitors. Lys33, Asp178, Asp267, Tyr306 and Leu308 of HDAC8 were favorable for binding with all inhibitors.

A novel series of l-2-benzyloxycarbonylamino-8-(2-pyridyl)-disulfidyloctanoic acid derivatives as histone deacetylase inhibitors: design, synthesis and molecular modeling study

Histone deacetylases inhibitors (HDACIs) have become an attractive class of anticancer agents. In order to find some novel potent HDACIs, we designed and synthesized a series of l-2-benzyloxycarbonylamino-8-(2-pyridyl)-disulfidyloctanoic acid derivatives. All compounds exhibited potent HDAC-inhibitory activity, and two of them had similar potency to TSA. The introduction of 2-amino-4-phenylthiazole or 9-methyleneoxy-fluorenyl group at the surface recognize domain of these HDACIs could greatly increase their HDAC-inhibitory activity. Molecular modeling studies indicated that coordination of the zinc ion by these inhibitors, and formation of hydrogen bond and hydrophobic interaction between inhibitors and HDACs were essential for the HDAC-inhibitory activities of these inhibitors. Asp181, Asp269, Leu276 and Tyr308 in the active site of HDAC2 gave favorable contributions for binding with all compounds.

MOLECULAR FACE — A NOVEL MOLECULAR MODEL SHOWING BOTH MOLECULAR SPATIAL CONTOUR AND FRONTIER ELECTRON DENSITY

The spatial knowledge is the first one of all information about an object. Molecular shape and size, molecular van der Waals surface and/or solvent-accessible surface etc. have been widely studied and applied. This paper is to show that a molecular face (MF) for a molecule may be defined uniquely and intrinsically via the molecular intrinsic characteristic contours (MICC) with coding the molecular electron density (ED) as the fourth dimension. The significant feature of an MF provides both molecular spatial appearance and its frontier electron density, being an intuitive picture as a molecular fingerprint or face. With simple examples, the physical significance of an MF is then demonstrated.

Geometric measures of large biomolecules: surface, volume, and pockets

Geometry plays a major role in our attempts to understand the activity of large molecules. For example, surface area and volume are used to quantify the interactions between these molecules and the water surrounding them in implicit solvent models. In addition, the detection of pockets serves as a starting point for predictive studies of biomolecule-ligand interactions. The alpha shape theory provides an exact and robust method for computing these geometric measures. Several implementations of this theory are currently available. We show however that these implementations fail on very large macromolecular systems. We show that these difficulties are not theoretical; rather, they are related to the architecture of current computers that rely on the use of cache memory to speed up calculation. By rewriting the algorithms that implement the different steps of the alpha shape theory such that we enforce locality, we show that we can remediate these cache problems; the corresponding code, UnionBall has an apparent O(n) behavior over a large range of values of n (up to tens of millions), where n is the number of atoms. As an example, it takes 136 sec with UnionBall to compute the contribution of each atom to the surface area and volume of a viral capsid with more than five million atoms on a commodity PC. UnionBall includes functions for computing analytically the surface area and volume of the intersection of two, three and four spheres that are fully detailed in an appendix. UnionBall is available as an OpenSource software.

Investigation of the molecular surface area and volume: Defined and calculated by the molecular face theory

Based on the molecular face (MF) theory, the molecular face surface area (MFSA) and molecular face volume (MFV) are defined. For a variety of organic molecules and several inorganic molecules, the MFSA and MFV have been studied and calculated in terms of an algorithm of our own via the Matlab package. The MFV shows a very good linear relationship with the experimentally measured critical molar volume. It is also found that the MFSA and MFV have significant linear correlations with those of the commonly used hard-sphere model and the electron density isosurface.

Poisson-Nernst-Planck Equations for Simulating Biomolecular Diffusion-Reaction Processes I: Finite Element Solutions

In this paper we developed accurate finite element methods for solving 3-D Poisson-Nernst-Planck (PNP) equations with singular permanent charges for electrodiffusion in solvated biomolecular systems. The electrostatic Poisson equation was defined in the biomolecules and in the solvent, while the Nernst-Planck equation was defined only in the solvent. We applied a stable regularization scheme to remove the singular component of the electrostatic potential induced by the permanent charges inside biomolecules, and formulated regular, well-posed PNP equations. An inexact-Newton method was used to solve the coupled nonlinear elliptic equations for the steady problems; while an Adams-Bashforth-Crank-Nicolson method was devised for time integration for the unsteady electrodiffusion. We numerically investigated the conditioning of the stiffness matrices for the finite element approximations of the two formulations of the Nernst-Planck equation, and theoretically proved that the transformed formulation is always associated with an ill-conditioned stiffness matrix. We also studied the electroneutrality of the solution and its relation with the boundary conditions on the molecular surface, and concluded that a large net charge concentration is always present near the molecular surface due to the presence of multiple species of charged particles in the solution. The numerical methods are shown to be accurate and stable by various test problems, and are applicable to real large-scale biophysical electrodiffusion problems.

A Fast Variational Method for the Construction of Resolution Adaptive C-Smooth Molecular Surfaces

We present a variational approach to smooth molecular (proteins, nucleic acids) surface constructions, starting from atomic coordinates, as available from the protein and nucleic-acid data banks. Molecular dynamics (MD) simulations traditionally used in understanding protein and nucleic-acid folding processes, are based on molecular force fields, and require smooth models of these molecular surfaces. To accelerate MD simulations, a popular methodology is to employ coarse grained molecular models, which represent clusters of atoms with similar physical properties by psuedo- atoms, resulting in coarser resolution molecular surfaces. We consider generation of these mixed-resolution or adaptive molecular surfaces. Our approach starts from deriving a general form second order geometric partial differential equation in the level-set formulation, by minimizing a first order energy functional which additionally includes a regularization term to minimize the occurrence of chemically infeasible molecular surface pockets or tunnel-like artifacts. To achieve even higher computational efficiency, a fast cubic B-spline C(2) interpolation algorithm is also utilized. A narrow band, tri-cubic B-spline level-set method is then used to provide C(2) smooth and resolution adaptive molecular surfaces.

Shape complementarity of protein-protein complexes at multiple resolutions

Biological complexes typically exhibit intermolecular interfaces of high shape complementarity. Many computational docking approaches use this surface complementarity as a guide in the search for predicting the structures of protein-protein complexes. Proteins often undergo conformational changes to create a highly complementary interface when associating. These conformational changes are a major cause of failure for automated docking procedures when predicting binding modes between proteins using their unbound conformations. Low resolution surfaces in which high frequency geometric details are omitted have been used to address this problem. These smoothed, or blurred, surfaces are expected to minimize the differences between free and bound structures, especially those that are due to side chain conformations or small backbone deviations. Despite the fact that this approach has been used in many docking protocols, there has yet to be a systematic study of the effects of such surface smoothing on the shape complementarity of the resulting interfaces. Here we investigate this question by computing shape complementarity of a set of 66 protein-protein complexes represented by multiresolution blurred surfaces. Complexed and unbound structures are available for these protein-protein complexes. They are a subset of complexes from a nonredundant docking benchmark selected for rigidity (i.e. the proteins undergo limited conformational changes between their bound and unbound states). In this work, we construct the surfaces by isocontouring a density map obtained by accumulating the densities of Gaussian functions placed at all atom centers of the molecule. The smoothness or resolution is specified by a Gaussian fall-off coefficient, termed "blobbyness." Shape complementarity is quantified using a histogram of the shortest distances between two proteins' surface mesh vertices for both the crystallographic complexes and the complexes built using the protein structures in their unbound conformation. The histograms calculated for the bound complex structures demonstrate that medium resolution smoothing (blobbyness = -0.9) can reproduce about 88% of the shape complementarity of atomic resolution surfaces. Complexes formed from the free component structures show a partial loss of shape complementarity (more overlaps and gaps) with the atomic resolution surfaces. For surfaces smoothed to low resolution (blobbyness = -0.3), we find more consistency of shape complementarity between the complexed and free cases. To further reduce bad contacts without significantly impacting the good contacts we introduce another blurred surface, in which the Gaussian densities of flexible atoms are reduced. From these results we discuss the use of shape complementarity in protein-protein docking.

Joint neighbors approximation of macromolecular solvent accessible surface area

A new method for approximate analytical calculations of solvent accessible surface area (SASA) for arbitrary molecules and their gradients with respect to their atomic coordinates was developed. This method is based on the recursive procedure of pairwise joining of neighboring atoms. Unlike other available methods of approximate SASA calculations, the method has no empirical parameters, and therefore can be used with comparable accuracy in calculations of SASA in folded and unfolded conformations of macromolecules of any chemical nature. As shown by tests with globular proteins in folded conformations, average errors in absolute atomic surface area is around 1 A2, while for unfolded protein conformations it varies from 1.65 to 1.87 A2. Computational times of the method are comparable with those by GETAREA, one of the fastest exact analytical methods available today.

AGBNP: an analytic implicit solvent model suitable for molecular dynamics simulations and high-resolution modeling

We have developed an implicit solvent effective potential (AGBNP) that is suitable for molecular dynamics simulations and high-resolution modeling. It is based on a novel implementation of the pairwise descreening Generalized Born model for the electrostatic component and a new nonpolar hydration free energy estimator. The nonpolar term consists of an estimator for the solute-solvent van der Waals dispersion energy designed to mimic the continuum solvent solute-solvent van der Waals interaction energy, in addition to a surface area term corresponding to the work of cavity formation. AGBNP makes use of a new parameter-free algorithm to calculate the scaling coefficients used in the pairwise descreening scheme to take into account atomic overlaps. The same algorithm is also used to calculate atomic surface areas. We show that excellent agreement is achieved for the GB self-energies and surface areas in comparison to accurate, but much more expensive, numerical evaluations. The parameter-free approach used in AGBNP and the sensitivity of the AGBNP model with respect to large and small conformational changes makes the model suitable for high-resolution modeling of protein loops and receptor sites as well as high-resolution prediction of the structure and thermodynamics of protein-ligand complexes. We present illustrative results for these kinds of benchmarks. The model is fully analytical with first derivatives and is computationally efficient. It has been incorporated into the IMPACT molecular simulation program.

Method and algorithm of obtaining the molecular intrinsic characteristic contours (MICCs) of organic molecules

The molecular intrinsic characteristic contour (MICC) is defined as the set of all the classical turning points of electron movement in a molecule. Studies on the MICCs of some medium organic molecules, such as dimethylether, acetone, and some homologues of alkanes, alkenes, and alkynes, as well as the electron density distributions on the MICCs, are shown for the first time. Results show that the MICC is an intrinsic approach to shape and size of a molecule. Unlike the van der Waals hard-sphere model, the MICC is a smooth contour, and it has a clear physical meaning. Detailed investigations on the cross-sections of MICCs have provided a kind of important information about atomic size changing in the process of forming molecules. Studies on electron density distribution on the MICC not only provide a new insight into molecular shape, but also show that the electron density distribution on the boundary surface relates closely with molecular properties and reactivities. For the homologues of alkanes, Rout(H), Dmin, and Dmax (the minimum and maximum of electron density on the MICC), all have very good linear relationships with minus of the molecular ionization potential. This work may serve as a basis for exploring a new reactivity indicator of chemical reactions and for studying molecular shape properties of large organic and biological molecules.

Application of the frozen atom approximation to the GB/SA continuum model for solvation free energy

The generalized Born/surface area (GB/SA) continuum model for solvation free energy is a fast and accurate alternative to using discrete water molecules in molecular simulations of solvated systems. However, computational studies of large solvated molecular systems such as enzyme-ligand complexes can still be computationally expensive even with continuum solvation methods simply because of the large number of atoms in the solute molecules. Because in such systems often only a relatively small portion of the system such as the ligand binding site is under study, it becomes less attractive to calculate energies and derivatives for all atoms in the system. To curtail computation while still maintaining high energetic accuracy, atoms distant from the site of interest are often frozen; that is, their coordinates are made invariant. Such frozen atoms do not require energetic and derivative updates during the course of a simulation. Herein we describe methodology and results for applying the frozen atom approach to both the generalized Born (GB) and the solvent accessible surface area (SASA) parts of the GB/SA continuum model for solvation free energy. For strictly pairwise energetic terms, such as the Coulombic and van-der-Waals energies, contributions from pairs of frozen atoms can be ignored. This leaves energetic differences unaffected for conformations that vary only in the positions of nonfrozen atoms. Due to the nonlocal nature of the GB analytical form, however, excluding such pairs from a GB calculation leads to unacceptable inaccuracies. To apply a frozen-atom scheme to GB calculations, a buffer region within the frozen-atom zone is generated based on a user-definable cutoff distance from the nonfrozen atoms. Certain pairwise interactions between frozen atoms in the buffer region are retained in the GB computation. This allows high accuracy in conformational GB comparisons to be maintained while achieving significant savings in computational time compared to the full (nonfrozen) calculation. A similar approach for using a buffer region of frozen atoms is taken for the SASA calculation. The SASA calculation is local in nature, and thus exact SASA energies are maintained. With a buffer region of 8 A for the frozen-atom cases, excellent agreement in differences in energies for three different conformations of cytochrome P450 with a bound camphor ligand are obtained with respect to the nonfrozen cases. For various minimization protocols, simulations run 2 to 10.5 times faster and memory usage is reduced by a factor of 1.5 to 5. Application of the frozen atom method for GB/SA calculations thus can render computationally tractable biologically and medically important simulations such as those used to study ligand-receptor binding conformations and energies in a solvated environment.

Approximate solvent-accessible surface areas from tetrahedrally directed neighbor densities

A fast analytical formula (TDND) has been derived for the calculation of approximate atomic and molecular solvent-accessible surface areas (SASA), as well as the first and second derivatives of these quantities with respect to atomic coordinates. Extending the work of Stouten et al. (Molecular Simulation, 1993, Vol. 10, pp. 97-120), as well as our own (Journal of Computational Chemistry, 1999, Vol. 20, pp. 586-596), the method makes use of a Gaussian function to calculate the neighbor density in four tetrahedral directions in three-dimensional space, sometimes twice with different orientations. SASA and first derivatives of the 2366 heavy atoms of penicillopepsin are computed in 0.13 s on an SGI R10000/194 MHz processor. When second derivatives are computed as well, the total time is 0.23 s. This is considerably faster than timings reported previously for other algorithms. Based on a parameterization set of nineteen compounds of different size (11-4346 atoms) and class (organics, proteins, DNA, and various complexes) consisting of a total 23,197 atoms, the method exhibits relative errors in the range 0.2-12.6% for total molecular surface areas and average absolute atomic surface area deviations in the range 1.7 to 3.6 A(2).