Free energy of solvation from molecular dynamics simulation applying Voronoi-Delaunay triangulation to the cavity creation

The free energy of solvation for a large number of representative solutes in various solvents has been calculated from the polarizable continuum model coupled to molecular dynamics computer simulation. A new algorithm based on the Voronoi-Delaunay triangulation of atom-atom contact points between the solute and the solvent molecules is presented for the estimation of the solvent-accessible surface surrounding the solute. The volume of the inscribed cavity is used to rescale the cavitational contribution to the solvation free energy for each atom of the solute atom within scaled particle theory. The computation of the electrostatic free energy of solvation is performed using the Voronoi-Delaunay surface around the solute as the boundary for the polarizable continuum model. Additional short-range contributions to the solvation free energy are included directly from the solute-solvent force field for the van der Waals-type interactions. Calculated solvation free energies for neutral molecules dissolved in benzene, water, CCl4, and octanol are compared with experimental data. We found an excellent correlation between the experimental and computed free energies of solvation for all the solvents. In addition, the employed algorithm for the cavity creation by Voronoi-Delaunay triangulation is compared with the GEPOL algorithm and is shown to predict more accurate free energies of solvation, especially in solvents composed by molecules with nonspherical molecular shapes.     

The effect of box shape on the dynamic properties of proteins simulated under periodic boundary conditions

The effect of the box shape on the dynamic behavior of proteins simulated under periodic boundary conditions is evaluated. In particular, the influence of simulation boxes defined by the near-densest lattice packing (NDLP) in conjunction with rotational constraints is compared to that of standard box types without these constraints. Three different proteins of varying size, shape, and secondary structure content were examined in the study. The statistical significance of differences in RMSD, radius of gyration, solvent-accessible surface, number of hydrogen bonds, and secondary structure content between proteins, box types, and the application or not of rotational constraints has been assessed. Furthermore, the differences in the collective modes for each protein between different boxes and the application or not of rotational constraints have been examined. In total 105 simulations were performed, and the results compared using a three-way multivariate analysis of variance (MANOVA) for properties derived from the trajectories and a three-way univariate analysis of variance (ANOVA) for collective modes. It is shown that application of roto-translational constraints does not have a statistically significant effect on the results obtained from the different simulations. However, the choice of simulation box was found to have a small (5-10%), but statistically significant effect on the behavior of two of the three proteins included in the study.     

A molecular dynamics study of the correlations between solvent-accessible surface, molecular volume, and folding state

We analyzed the correlations between molecular volume, solvent-accessible surface, and folding state (secondary structure content) for unfolded conformers of alpha (holo- and apomyoglobin) and beta (retinal-binding protein) proteins and a small water-soluble alanine-rich alpha-helical peptide. Conformers with different degrees of folding were obtained using molecular dynamics at constant temperature and pressure with implicit solvent (dielectric constant adjustment) for all four systems and with explicit solvent for the single helix peptide. Our results support the view that unfolded conformations are not necessary extended, that volume variation is not a good indication of folding state and that the simple model of water penetrating the interior of the protein does not explain the increase in volume upon unfolding.     

A Fast and Robust Poisson-Boltzmann Solver Based on Adaptive Cartesian Grids

An adaptive Cartesian grid (ACG) concept is presented for the fast and robust numerical solution of the 3D Poisson-Boltzmann Equation (PBE) governing the electrostatic interactions of large-scale biomolecules and highly charged multi-biomolecular assemblies such as ribosomes and viruses. The ACG offers numerous advantages over competing grid topologies such as regular 3D lattices and unstructured grids. For very large biological molecules and multi-biomolecule assemblies, the total number of grid-points is several orders of magnitude less than that required in a conventional lattice grid used in the current PBE solvers thus allowing the end user to obtain accurate and stable nonlinear PBE solutions on a desktop computer. Compared to tetrahedral-based unstructured grids, ACG offers a simpler hierarchical grid structure, which is naturally suited to multigrid, relieves indirect addressing requirements and uses fewer neighboring nodes in the finite difference stencils. Construction of the ACG and determination of the dielectric/ionic maps are straightforward, fast and require minimal user intervention. Charge singularities are eliminated by reformulating the problem to produce the reaction field potential in the molecular interior and the total electrostatic potential in the exterior ionic solvent region. This approach minimizes grid-dependency and alleviates the need for fine grid spacing near atomic charge sites. The technical portion of this paper contains three parts. First, the ACG and its construction for general biomolecular geometries are described. Next, a discrete approximation to the PBE upon this mesh is derived. Finally, the overall solution procedure and multigrid implementation are summarized. Results obtained with the ACG-based PBE solver are presented for: (i) a low dielectric spherical cavity, containing interior point charges, embedded in a high dielectric ionic solvent - analytical solutions are available for this case, thus allowing rigorous assessment of the solution accuracy; (ii) a pair of low dielectric charged spheres embedded in a ionic solvent to compute electrostatic interaction free energies as a function of the distance between sphere centers; (iii) surface potentials of proteins, nucleic acids and their larger-scale assemblies such as ribosomes; and (iv) electrostatic solvation free energies and their salt sensitivities - obtained with both linear and nonlinear Poisson-Boltzmann equation - for a large set of proteins. These latter results along with timings can serve as benchmarks for comparing the performance of different PBE solvers.     

Temperature dependence of the configurational free energy of a branched polymer

We consider a lattice model of branched polymers in dilute solution in which the polymer is modelled as an animal, weakly embeddable in the (simple cubic) lattice. In order to model the effect on the thermodynamic properties of changing the temperature or the quality of the solvent, we include an energy associated with the number of nearneighbour contacts between pairs of vertices of the animals. We show that the configurational free energy of the animal is a continuous function of the temperature and derive rigorous upper and lower bounds on the temperature dependence of the free energy. Finally, we comment on similarities between these results and corresponding ones for a model in which the energy is associated with the cyclomatic index of the animal.     

Atomic radii: Incorporation of solvation effects

Atomic radii used to define the solute cavity in continuum-based methods are determined by reproducing the solvent-accessible surface defined as the loci of minima in a potential (solvent interaction potential) between the solute and a probe. This potential includes electrostatic interaction (ion–dipole, ion–quadrupole, and ion-induced dipole) terms as well as a Lennard–Jones energy term. The method alleviates the need to distinguish solute atoms in different chemical environments. These radii, when used in the calculation of solvation free energies, are shown to be superior to fixed atom-specific radii or to radii obtained from the electron isodensity surface from quantum-mechanical calculations. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1482–1493, 1998     

二十面体壳层结构：纪念C60发现10周年

圆顶建筑、病毒衣壳和Ｃ６０分子笼的生成机制和尺寸差别很大，但都属二十面体壳层，遵循同一几何规律。圆顶建筑和病毒衣壳由近似一等的三角面构成，有１２个五键加３＾５顶，１０（Ｔ－１）个六键连３＾６顶，其中三角面数Ｔ＝ｈ＾２＋ｈｋ＋ｋ＾２是二十面体的一个三角面中的小三角面数，ｈ，ｋ是六角坐标系中３＾５顶的坐标。另一方面，全碳分子笼Ｃｎ（包括其原型Ｃ６０），由１２个五角面与１０（Ｔ－１）个六角面围成，全是三     

Generation and enumeration of compact conformations on the two-dimensional triangular and three-dimensional fcc lattices

We enumerated all compact conformations within simple geometries on the two-dimensional (2D) triangular and three-dimensional (3D) face centered cubic (fcc) lattice. These compact conformations correspond mathematically to Hamiltonian paths and Hamiltonian circuits and are frequently used as simple models of proteins. The shapes that were studied for the 2D triangular lattice included mxn parallelograms, regular equilateral triangles, and various hexagons. On the 3D fcc lattice we generated conformations for a limited class of skewed parallelepipeds. Symmetries of the shape were exploited to reduce the number of conformations. We compared surface to volume ratios against protein length for compact conformations on the 3D cubic lattice and for a selected set of real proteins. We also show preliminary work in extending the transfer matrix method, previously developed by us for the 2D square and the 3D cubic lattices, to the 2D triangular lattice. The transfer matrix method offers a superior way of generating all conformations within a given geometry on a lattice by completely avoiding attrition and reducing this highly complicated geometrical problem to a simple algebraic problem of matrix multiplication.     

From the Cube to the Dyck and Klein Tessellations: Implications for the Structures of Zeolite-like Carbon and Boron Nitride Allotropes

Geometrical transformations can be described which have the effect of multiplying the numbers of vertices in a trivalent polyhedron by three, four, or seven. Tripling the cube by the so-called leapfrog transformation gives the truncated octahedron. Quadrupling the cube followed by identifying the square faces to give a genus 3 surface gives the Dyck surface of 12 octagons. Septupling the cube by the so-called capra transformation followed by identifying the square faces to give a genus 3 surface gives the Klein surface of 24 heptagons. These geometrical transformations relate to the construction of low-density zeolite-like structures for carbon and boron nitride allotropes based on a cubic lattice     

About second kind continuous chirality measures. 1. Planar sets

The chirality index of a d-dimensional set of n points is defined as the sum of the n squared distances between the vertices of the set and those of its inverted image, normalized to 4T/d,T being the inertia of the set. The index is computed after minimization of the sum of the squared distances with respect to all rotations and translations and all permutations between equivalent vertices. The properties of the chiral index are examined for planar sets. The most achiral triangles are obtained analytically for all equivalence situations: one, two, and three equivalent vertices. These triangles are different from those obtained by Weinberg and Mislow with distance functions. This revised version was published online in July 2006 with corrections to the Cover Date.     

Morphometric approach to thermodynamic quantities of solvation of complex molecules: extension to multicomponent solvent

The morphometric approach (MA) is a powerful tool for calculating a solvation free energy (SFE) and related quantities of solvation thermodynamics of complex molecules. Here, we extend it to a solvent consisting of m components. In the integral equation theories, the SFE is expressed as the sum of m terms each of which comprises solute-component j correlation functions (j = 1,..., m). The MA is applied to each term in a formally separate manner: The term is expressed as a linear combination of the four geometric measures, excluded volume, solvent-accessible surface area, and integrated mean and Gaussian curvatures of the accessible surface, which are calculated for component j. The total number of the geometric measures or the coefficients in the linear combinations is 4m. The coefficients are determined in simple geometries, i.e., for spherical solutes with various diameters in the same multicomponent solvent. The SFE of the spherical solutes are calculated using the radial-symmetric integral equation theory. The extended version of the MA is illustrated for a protein modeled as a set of fused hard spheres immersed in a binary mixture of hard spheres. Several mixtures of different molecular-diameter ratios and compositions and 30 structures of the protein with a variety of radii of gyration are considered for the illustration purpose. The SFE calculated by the MA is compared with that by the direct application of the three-dimensional integral equation theory (3D-IET) to the protein. The deviations of the MA values from the 3D-IET values are less than 1.5%. The computation time required is over four orders of magnitude shorter than that in the 3D-IET. The MA thus developed is expected to be best suited to analyses concerning the effects of cosolvents such as urea on the structural stability of a protein.     

Hamiltonian circuits,Hamiltonian paths and branching graphs of benzenoid systems

A benzenoid systemH is a finite connected subgraph of the infinite hexagonal lattice with out cut bonds and non-hexagonal interior faces. The branching graphG ofH consists of all vertices ofH of degree 3 and bonds among them. In this paper, the following results are obtained:

1. A necessary condition for a benzenoid system to have a Hamiltonian circuit.
2. A necessary and sufficient condition for a benzenoid system to have a Hamiltonian path.
3. A characterization of connected subgraphs of the infinite hexagonal lattice which are branching graphs of benzenoid systems.
4. A proof that if a disconnected subgraph G of the infinite hexagonal lattice given along with the positions of its vertices is the branching graph of a benzenoid system H, then H is unique.

     

Grid inhomogeneous solvation theory: Hydration structure and thermodynamics of the miniature receptor cucurbit[7]uril

The displacement of perturbed water upon binding is believed to play a critical role in the thermodynamics of biomolecular recognition, but it is nontrivial to unambiguously define and answer questions about this process. We address this issue by introducing grid inhomogeneous solvation theory (GIST), which discretizes the equations of inhomogeneous solvation theory (IST) onto a three-dimensional grid situated in the region of interest around a solute molecule or complex. Snapshots from explicit solvent simulations are used to estimate localized solvation entropies, energies, and free energies associated with the grid boxes, or voxels, and properly summing these thermodynamic quantities over voxels yields information about hydration thermodynamics. GIST thus provides a smoothly varying representation of water properties as a function of position, rather than focusing on hydration sites where solvent is present at high density. It therefore accounts for full or partial displacement of water from sites that are highly occupied by water, as well as for partly occupied and water-depleted regions around the solute. GIST can also provide a well-defined estimate of the solvation free energy and therefore enables a rigorous end-states analysis of binding. For example, one may not only use a first GIST calculation to project the thermodynamic consequences of displacing water from the surface of a receptor by a ligand, but also account, in a second GIST calculation, for the thermodynamics of subsequent solvent reorganization around the bound complex. In the present study, a first GIST analysis of the molecular host cucurbit[7]uril is found to yield a rich picture of hydration structure and thermodynamics in and around this miniature receptor. One of the most striking results is the observation of a toroidal region of high water density at the center of the host's nonpolar cavity. Despite its high density, the water in this toroidal region is disfavored energetically and entropically, and hence may contribute to the known ability of this small receptor to bind guest molecules with unusually high affinities. Interestingly, the toroidal region of high water density persists even when all partial charges of the receptor are set to zero. Thus, localized regions of high solvent density can be generated in a binding site without strong, attractive solute-solvent interactions.     

Surface-integral QSPR models: local energy properties

Surface-integral models based on AM1 semiempirical molecular orbital calculations are presented for the free energies of solvation in water, n-octanol, and chloroform and for the enthalpy of solvation in water. A parametrized function of four local properties calculated at the isodensity surface (the molecular electrostatic potential, local ionization energy, electron affinity, and polarizability) is integrated over the triangulated surface area to obtain the target quantity. The resulting models give results only slightly less accurate than those reported for parametrized generalized Born/polar surface area models despite relying only on gas-phase calculations. The water and octanol free-energy models were validated by calculating the water-octanol partition coefficient for a test set of organic compounds with moderate success. The models lead to a local solvation energy, which can be projected onto the molecular isodensity surface and provides insight into "hot" areas for solvation in water or the other solvents.     

Solvent interaction energy calculations on molecular dynamics trajectories: increasing the efficiency using systematic frame selection

End-point methods such as linear interaction energy (LIE) analysis, molecular mechanics generalized Born solvent-accessible surface (MM/GBSA), and solvent interaction energy (SIE) analysis have become popular techniques to calculate the free energy associated with protein-ligand binding. Such methods typically use molecular dynamics (MD) simulations to generate an ensemble of protein structures that encompasses the bound and unbound states. The energy evaluation method (LIE, MM/GBSA, or SIE) is subsequently used to calculate the energy of each member of the ensemble, thus providing an estimate of the average free energy difference between the bound and unbound states. The workflow requiring both MD simulation and energy calculation for each frame and each trajectory proves to be computationally expensive. In an attempt to reduce the high computational cost associated with end-point methods, we study several methods by which frames may be intelligently selected from the MD simulation including clustering and address the question of how the number of selected frames influences the accuracy of the SIE calculations.     

Crucial importance of translational entropy of water in pressure denaturation of proteins

We present statistical thermodynamics of pressure denaturation of proteins, in which the three-dimensional integral equation theory is employed. It is applied to a simple model system focusing on the translational entropy of the solvent. The partial molar volume governing the pressure dependence of the structural stability of a protein is expressed for each structure in terms of the excluded volume for the solvent molecules, the solvent-accessible surface area (ASA), and a parameter related to the solvent-density profile formed near the protein surface. It is argued that the entropic effect originating from the translational movement of water molecules plays critical roles in the pressure-induced denaturation. We also show that the exceptionally small size of water molecules among dense liquids in nature is crucial for pressure denaturation. An unfolded structure, which is only moderately less compact than the native structure but has much larger ASA, is shown to turn more stable than the native one at an elevated pressure. The water entropy for the native structure is higher than that for the unfolded structure in the low-pressure region, whereas the opposite is true in the high-pressure region. Such a structure is characterized by the cleft and/or swelling and the water penetration into the interior. In another solvent whose molecular size is 1.5 times larger than that of water, however, the inversion of the stability does not occur any longer. The random coil becomes relatively more destabilized with rising pressure, irrespective of the molecular size of the solvent. These theoretical predictions are in qualitatively good agreement with the experimental observations.     

Lattice modeling: Accuracy of energy calculations

Energy calculations based on lattice models of protein chains are always approximate, because any such a model distorts distances between chain links and, consequently, the energies of interaction between them. The energetic errors of lattice models are examined here for 15 proteins of different sizes and types of secondary structure, for lattice spacings ranging from 0.25 to 2.5 Å. The lattice models are derived using previously described algorithms which insure a minimal root mean square (rms) deviation from the off-lattice structure for any given lattice-protein orientation. For each protein structure we computed a set of different lattice models with virtually equal rms deviations, and then compared their energies. Energy calculations were based on the pairwise potentials. We found that the energies of lattice models follows a normal distribution with a nonnegligible dispersion, even at a fine lattice spacing of 0.25 Å. For any lattice model of a protein, the lattice spacing must be 1.0 Å or less in order to be able to distinguish energetically between the folded and extended states. However, when an ensemble of lattice models is considered, this distinction can be made for lattice spacing up to 2.0 Å. We conclude that to attain a better approximation of the protein lattice model energies, one must adjust potentials to the lattice spacing. © 1996 by John Wiley & Sons, Inc.     

A precise boundary element method for macromolecular transport properties

A very precise boundary element numerical solution of the exact formulation of the hydrodynamic resistance problem with stick boundary conditions is presented. BEST, the Fortran 77 program developed for this purpose, computes the full transport tensors in the center of resistance or the center of diffusion for an arbitrarily shaped rigid body, including rotation-translation coupling. The input for this program is a triangulation of the solvent-defined surface of the molecule of interest, given by Connolly's MSROLL or other suitable triangulator. The triangulation is prepared for BEST by COALESCE, a program that allows user control over the quality and number of triangles to describe the surface. High numerical precision is assured by effectively exact integration of the Oseen tensor over triangular surface elements, and by scaling the hydrodynamic computation to the precise surface area of the molecule. Efficiency of computation is achieved by the use of public domain LAPACK routines that call BLAS Level 3 hardware-optimized subroutines available for most processors. A protein computation can be done in less than 10 min of CPU time in a modern Pentium IV processor. The present work includes a complete analysis of the sources of error in the numerical work and techniques to eliminate these errors. The operation of BEST is illustrated with applications to ellipsoids of revolution, and Lysozyme, a small protein. The typical numerical accuracy achieved is 0.05% compared to analytical theory. The numerical precision for a protein is better than 1%, much better than experimental errors in these quantities, and more than 10 times better than traditional bead-based methods.     

Metal triangles as building blocks in metal cluster chemistry

Many trinuclear metal clusters have structures based on isolated metal triangles with either single bonds (e.g.,M ₃(CO)₁₂ whereM = Fe, Ru, Os) or double bonds (e.g., Re₃ Cl ₁₂ ^3– ) along each edge of the triangle. Individual metal triangles can be joined in the following ways to form more complicated triangulated networks: (1) Bridging an edge of a triangle with a new vertex to give rafts in which adjacent triangles share edges; (2) Bridging a vertex of a triangle with a new edge to give bowties in which adjacent triangles share vertices; (3) Capping a triangular face with a new vertex to give a chain of tetrahedra in which adjacent tetrahedra share faces. Such triangulated metal networks are particularly common in osmium carbonyl chemistry and in mixed osmium/platinum carbonyl derivatives. Platinum triangles of the type Pt₃L₆ are analogous to cyclopropenyl rings and can form sandwiches with one or two mercury atoms in the center such as the mercuric derivative Hg[Pt.₃(µ₂-2,6-Me₂C₆H₃NC)₃] (2,6-Me₂C₆H₃NC)₃]₂ and the mercurous derivative Hg₂[Pt₃(µ₂-CO)₃L₃]₂. Platinum triangles can also be stacked in the absence of filling to give [Pt₃(µ₂-CO)₃(CO)₃] _n ^2– (n=2, 3, 4, 5, 6, 10). Metal triangles also form the faces of metal deltahera of which the octahedron, bicapped square antiprism, and icosahedron are found in globally delocalized transition metal clusters.This article is dedicated to Prof. L. F. Dahl in recognition of his many seminal contributions to metal cluster chemistry.