Efficient approximate all-atom solvent accessible surface area method parameterized for folded and denatured protein conformations

Continuing advances in computer hardware and software are permitting atomic-resolution molecular simulations for longer time scales and on larger systems. Despite these advances, routinely performing atomistic simulations with explicit water for even small proteins, which reach the folding time of such proteins, remains intractable for the foreseeable future. An implicit approximation of the solvent environment using a solvent accessible surface area (SASA) term in a molecular mechanics potential function allows exclusion of the explicit water molecules in protein simulations. This reduces the number of particles by approximately an order of magnitude. We present a fast and acceptably accurate approximate all-atom SASA method parameterized using a set of folded and heat-denatured conformations of globular proteins. The parameters are shown to be transferable to folded and heat-denatured conformations for another set of proteins. Calculation of the approximate SASA and the associated derivatives with respect to atomic positions for a 4644 atom protein requires only 1/11th the CPU time required for calculation of the nonbonded interactions for this system. On a per atom basis, this algorithm is three times faster than the fastest previously published approximate SASA method and achieves the same level of accuracy.     

A rapid method for calculating derivatives of solvent accessible surface areas of molecules

A novel method to calculate the derivatives of solvent accessible surface areas is presented. Unlike earlier analytic methods, which require the molecular topology and the use of global Gauss-Bonnet theorem, this method requires only the fractional accessibilities of surface arcs. We developed an efficient numerical algorithm to calculate the surface arcs by creating a uniform set of points on the circles of intersection between surface atoms. A hierarchical point density doubling scheme led to a logarithmic dependence of Central Processing Unit (CPU) time on the number of points used. This algorithm calculated area derivatives for a 1000-atom protein in 1.5 s on an SGI INDIGO² which were within 2% of the analytic area derivatives calculated with the program ANAREA. This algorithm scales linearly with the number of atoms for large molecules and is easily parallelizable. © 1995 by John Wiley & Sons, Inc.     

Maintaining solvent accessible surface area under rotamer substitution for protein design

Although quantities derived from solvent accessible surface areas (SASA) are useful in many applications in protein design and structural biology, the computational cost of accurate SASA calculation makes SASA-based scores difficult to integrate into commonly used protein design methodologies. We demonstrate a method for maintaining accurate SASA during a Monte Carlo search of sequence and rotamer space for a fixed protein backbone. We extend the fast Le Grand and Merz algorithm (Le Grand and Merz, J Comput Chem, 14, 349), which discretizes the solvent accessible surface for each atom by placing dots on a sphere and combines Boolean masks to determine which dots are exposed. By replacing semigroup operations with group operations (from Boolean logic to counting dot coverage) we support SASA updates. Our algorithm takes time proportional to the number of atoms affected by rotamer substitution, rather than the number of atoms in the protein. For design simulations with a one hundred residue protein our approach is approximately 145 times faster than performing a Le Grand and Merz SASA calculation from scratch following each rotamer substitution. To demonstrate practical effectiveness, we optimize a SASA-based measure of protein packing in the complete redesign of a large set of proteins and protein-protein interfaces.     

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.     

Fast,approximate algorithm for detection of solvent-inaccessible atoms

Up to about half of the atoms in biopolymers are inaccessible to solvents. If such atoms can be rapidly identified, time can be saved in the subsequent computation of atomic surface areas. A quick, approximate method, termed buried atom elimination (BAE), was developed for the detection of such atoms. Following the literature, the method makes use of a Gaussian function to calculate the neighbor density in four tetrahedral directions in 3-dimensional space, sometimes twice with different orientations. In macromolecules, our method detects between 63 and 81% of the buried atoms but also incorrectly classifies 2–8% of the exposed atoms as buried. These misidentified atoms all have small solvent-exposed (accessible) surface areas (SASAs): their surfaces sum to a maximum of 0.5% of the molecular SASA, and their maximum atomic SASA is 5.1 Ų. Using our recently reported LCPO method for computing atomic surfaces, which is one of the fastest available, the use of BAE increases the overall speed of computing the atomic SASAs by a factor of up to 1.6 for surfaces only and 1.9 when first and second derivatives are computed. BAE decreases the LCPO average absolute atomic error from about 2.3 Ų to about 1.7 Ų (average for larger compounds). BAE was introduced into the MacroModel molecular modeling package and tests show that it increases the efficiency of first- and second-derivative energy minimizations and molecular dynamics simulations without adversely affecting the stability or accuracy of the calculations. BAE parameters were developed for the most important atom types in biopolymers, based on a parameterization set of 18 compounds of different size (33–4346 atoms) and class (organics, proteins, DNA, and various complexes), consisting of a total of 23,186 atoms. ©1999 John Wiley & Sons, Inc. J Comput Chem 20, 586–596, 1999     

Fast summation boundary element method for calculating solvation free energies of macromolecules

We present a boundary element method (BEM) for calculating the reaction field energy of a macromolecule embedded in a high-dielectric medium such as water. In a BEM calculation, the key computational task is the calculation of the induced surface charge distribution at the dielectric boundary. This is obtained by solving a system of linear equations whose dimension can run into the tens of thousands for a macromolecule. In this work, we use a fast summation hierarchical multipole method to solve for the induced surface charge densities. By careful analysis of the levels of approximation required for the various terms in the calculation, we avoid the unnecessary computation of terms that contribute negligibly to the final outcome and, consequently, achieve high computational efficiency. For a protein such as BPTI with 890 atoms, the calculation of the induced surface charge density distribution and the reaction field energy was completed in 7.9 s on an SGI workstation with an R10000 CPU. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1494–1504, 1998     

Evaluation of the protein solvent-accessible surface using reduced representations in terms of critical points of the electron density

The aim of our study is the development of a method for calculating the interface of dimerization of protein-protein complexes based on simplified medium-resolution structures. In particular, we wished to evaluate if the existing concepts for the computation of the Solvent-Accessible Surface Area (SASA) of macromolecules could be applied to medium-resolution models. Therefore, we selected a set of 140 protein chains and computed their reduced representations by topological analysis of their electron density maps at 2.85 A crystallographic resolution. This procedure leads to a limited number of critical points (CPs) that can be identified and associated to backbone and side-chain parts. To evaluate the SASA and interfaces of dimerization of the reduced representations, we chose and modified two existing programs that calculate the SASA of atomic representations, and tested (1) several radii tables of amino acids, (2) the influence of the backbone and side-chain points, and (3) the radius of the solvent molecule, which rolls over the surface. The results are shown in terms of relative error compared to the values calculated on the corresponding atomic representations of the proteins.     

Combining the polarizable Drude force field with a continuum electrostatic Poisson–Boltzmann implicit solvation model

In this work, we have combined the polarizable force field based on the classical Drude oscillator with a continuum Poisson–Boltzmann/solvent‐accessible surface area (PB/SASA) model. In practice, the positions of the Drude particles experiencing the solvent reaction field arising from the fixed charges and induced polarization of the solute must be optimized in a self‐consistent manner. Here, we parameterized the model to reproduce experimental solvation free energies of a set of small molecules. The model reproduces well‐experimental solvation free energies of 70 molecules, yielding a root mean square difference of 0.8 kcal/mol versus 2.5 kcal/mol for the CHARMM36 additive force field. The polarization work associated with the solute transfer from the gas‐phase to the polar solvent, a term neglected in the framework of additive force fields, was found to make a large contribution to the total solvation free energy, comparable to the polar solute–solvent solvation contribution. The Drude PB/SASA also reproduces well the electronic polarization from the explicit solvent simulations of a small protein, BPTI. Model validation was based on comparisons with the experimental relative binding free energies of 371 single alanine mutations. With the Drude PB/SASA model the root mean square deviation between the predicted and experimental relative binding free energies is 3.35 kcal/mol, lower than 5.11 kcal/mol computed with the CHARMM36 additive force field. Overall, the results indicate that the main limitation of the Drude PB/SASA model is the inability of the SASA term to accurately capture non‐polar solvation effects. © 2018 Wiley Periodicals, Inc.     

Binding mode prediction and MD/MMPBSA-based free energy ranking for agonists of REV-ERBα/NCoR

The knowledge of the free energy of binding of small molecules to a macromolecular target is crucial in drug design as is the ability to predict the functional consequences of binding. We highlight how a molecular dynamics (MD)-based approach can be used to predict the free energy of small molecules, and to provide priorities for the synthesis and the validation via in vitro tests. Here, we study the dynamics and energetics of the nuclear receptor REV-ERBα with its co-repressor NCoR and 35 novel agonists. Our in silico approach combines molecular docking, molecular dynamics (MD), solvent-accessible surface area (SASA) and molecular mechanics poisson boltzmann surface area (MMPBSA) calculations. While docking yielded initial hints on the binding modes, their stability was assessed by MD. The SASA calculations revealed that the presence of the ligand led to a higher exposure of hydrophobic REV-ERB residues for NCoR recruitment. MMPBSA was very successful in ranking ligands by potency in a retrospective and prospective manner. Particularly, the prospective MMPBSA ranking-based validations for four compounds, three predicted to be active and one weakly active, were confirmed experimentally.     

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.     

vmdICE: a plug-in for rapid evaluation of molecular dynamics simulations using VMD

Molecular dynamics (MD) is a powerful in silico method to investigate the interactions between biomolecules. It solves Newton's equations of motion for atoms over a specified period of time and yields a trajectory file, containing the different spatial arrangements of atoms during the simulation. The movements and energies of each single atom are recorded. For evaluating of these simulation trajectories with regard to biomedical implications, several methods are available. Three well-known ones are the root mean square deviation (RMSD), the root mean square fluctuation (RMSF) and solvent accessible surface area (SASA). Herein, we present a novel plug-in for the software "visual molecular dynamics" (VMD) that allows an interactive 3D representation of RMSD, RMSF, and SASA, directly on the molecule. On the one hand, our plug-in is easy to handle for inexperienced users, and on the other hand, it provides a fast and flexible graphical impression of the spatial dynamics of a system for experts in the field.     

Using path sampling to build better Markovian state models: predicting the folding rate and mechanism of a tryptophan zipper beta hairpin

We propose an efficient method for the prediction of protein folding rate constants and mechanisms. We use molecular dynamics simulation data to build Markovian state models (MSMs), discrete representations of the pathways sampled. Using these MSMs, we can quickly calculate the folding probability (P(fold)) and mean first passage time of all the sampled points. In addition, we provide techniques for evaluating these values under perturbed conditions without expensive recomputations. To demonstrate this method on a challenging system, we apply these techniques to a two-dimensional model energy landscape and the folding of a tryptophan zipper beta hairpin.     

化合物分子表面积计算方法的研究

化合物的分子表面积是重要的物理化学性质参数．根据不同的应用领域及数学方法,提出了多种分子表面积算法[1-10],由此产生了“分子表面积法”,各方法的有效性在各具体应用领域都已被验证．本文将统计方法中的随机变量引入计算化合物分子表面积体系,由分子模型化技术得到化合物分子的原子坐标,不考虑化合物分子中原子的相互作用及分子间近似,直接计算分子表面积．使用该方法可以计算“净”分子表面积、“溶剂可及表面积”、甚至分子结构片段,分子结构中有交叉重叠片段及存在“空洞”的各种分子表面积,该算法及程序较简捷,适应范围广,计算结果较为满意．     

Study of protein fluctuation with an effective inter-Cα atomic potential derived from average distances between amino acids in proteins

The conformational dynamics around the native structure of bovine pancreatic trypsin inhibitor (BPTI) in both the oxidized and reduced forms was analyzed by a Monte Carlo method using an approximate residue–residue potential derived from the statistics of average distances between C^α atoms of residues as proposed by the present author (T. Kikuchi, J Comput Chem 1996, 17, 226–237). The results from using this effective potential are similar to those from molecular dynamics simulations, taking all atoms into account, and are consistent with temperature factors from an X-ray analysis and disulfide formation from a kinetic experiment. This agreement suggests that the essential nature of the potential energy surface formed by the potential around the native structure closely mimics the actual energy landscape within the resolution of C^α atomic fluctuation. Furthermore, it is expected that the potential we found can describe the basic properties of folding kinetics. Examination of the fluctuation property of the native structure of BPTI threaded by a sequence from cytochrome b562 reveals differences specific to the sequence and this result also shows that the dynamical properties obtained in our calculations are not only ascribed to the geometrical constraints of the initial conformation but also the force field specifically produced by a sequence. ©1999 John Wiley & Sons, Inc. J Comput Chem 20: 713–719, 1999     

Topological analysis of complex molecular surfaces

An algorithm for the calculation of local and global curvatures of molecular surfaces is presented. The analysis is based on a surface representation as a set of points in 3-D space (“dotted surface” representation). The surface data are used to subdivide the surfac into domains with different curvatures. All domains are characterized by a reference point with a corresponding curvature profile specifiying the topological properties in its neighborhood. The curvature profiles provide a method for a systematic comparison of the shapes of different molecules. Such a strategy is important for the treatment of molecular recognition problems. The enzyme-inhibitor complex trypsin/BPTI was chosen to demonstrate the scopes of the method.     

Incorporation of Molecular Reorientation into Modeling Surface Pressure-Area Isotherms of Langmuir Monolayers

Langmuir monolayers can be assembled from molecules that change from a low-energy orientation occupying a large cross-sectional area to a high-energy orientation of small cross-sectional area as the lateral pressure grows. Examples include cyclosporin A, amphotericin B, nystatin, certain alpha-helical peptides, cholesterol oxydation products, dumbbell-shaped amphiphiles, organic–inorganic nanoparticles and hybrid molecular films. The transition between the two orientations leads to a shoulder in the surface pressure-area isotherm. We propose a theoretical model that describes the shoulder and can be used to extract the energy cost per molecule for the reorientation. Our two-state model is based on a lattice–sublattice approximation that hosts the two orientations and a corresponding free energy expression which we minimize with respect to the orientational distribution. Inter-molecular interactions other than steric repulsion are ignored. We provide an analysis of the model, including an analytic solution for one specific lateral pressure near a point of inflection in the surface pressure-area isotherm, and an approximate solution for the entire range of the lateral pressures. We also use our model to estimate energy costs associated with orientational transitions from previously reported experimental surface pressure-area isotherms.     

Analytical calculation of the solvent-accessible surface area and its nuclear gradients by stereographic projection: A general approach for molecules,polymers, nanotubes,helices, and surfaces

In this article, we explore an alternative to the analytical Gauss–Bonnet approach for computing the solvent-accessible surface area (SASA) and its nuclear gradients. These two key quantities are required to evaluate the nonelectrostatic contribution to the solvation energy and its nuclear gradients in implicit solvation models. We extend a previously proposed analytical approach for finite systems based on the stereographic projection technique to infinite periodic systems such as polymers, nanotubes, helices, or surfaces and detail its implementation in the Crystal code. We provide the full derivation of the SASA nuclear gradients, and introduce an iterative perturbation scheme of the atomic coordinates to stabilize the gradients calculation for certain difficult symmetric systems. An excellent agreement of computed SASA with reference analytical values is found for finite systems, while the SASA size-extensivity is verified for infinite periodic systems. In addition, correctness of the analytical gradients is confirmed by the excellent agreement obtained with numerical gradients and by the translational invariance achieved, both for finite and infinite periodic systems. Overall therefore, the stereographic projection approach appears as a general, simple, and efficient technique to compute the key quantities required for the calculation of the nonelectrostatic contribution to the solvation energy and its nuclear gradients in implicit solvation models applicable to both finite and infinite periodic systems.