A reversible minimum-to-minimum mapping method for the calculation of free-energy differences

A general method is introduced for the calculation of the free-energy difference between two systems, 0 and 1, with configuration spaces omega(0), omega(1) of the same dimensionality. The method relies upon establishing a objective mapping between disjoint subsets gamma(i)(0) of omega(0) and corresponding disjoint subsets gamma(i)(1) of omega(1), and averaging a function of the ratio of configurational integrals over gamma(i)(0) and gamma(i)(1) with respect to the probability densities of the two systems. The mapped subsets gamma(i)(0) and gamma(i)(1) need not span the entire configuration spaces omega(0) and omega(1). The method is applied for the calculation of the excess chemical potential mu(ex) in a Lennard-Jones (LJ) fluid. In this case, omega(0) is the configuration space of a (N-1) real molecule plus one ideal-gas molecule system, while omega(1) is the configuration space of a N real molecule system occupying the same volume. Gamma(i)(0) and gamma(i)(1) are constructed from hyperspheres of the same radius centered at minimum-energy configurations of a set of "active" molecules lying within distance a from the ideal-gas molecule and the last real molecule, respectively. An algorithm is described for sampling gamma(i)(0) and gamma(i)(1) given a point in omega(0) or in omega(1). The algorithm encompasses three steps: "quenching" (minimization with respect to the active-molecule degrees of freedom), "mutation" (gradual conversion of the ideal-gas molecule into a real molecule, with simultaneous minimization of the energy with respect to the active-molecule degrees of freedom), and "excitation" (generation of points on a hypersphere centered at the active-molecule energy minimum). These steps are also carried out in reverse, as required by the bijective nature of the mapping. The mutation step, which establishes a reversible mapping between energy minima with respect to the active degrees of freedom of systems 0 and 1, ensures that excluded volume interactions emerging in the process of converting the ideal-gas molecule into a real molecule are relieved through appropriate rearrangement of the surrounding active molecules. Thus, the insertion problem plaguing traditional methods for the calculation of chemical potential at high densities is alleviated. Results are presented at two state points of the LJ system for a variety of radii a of the active domain. It is shown that the estimated values of mu(ex) are correct in all cases and subject to an order of magnitude lower statistical uncertainty than values based on the same number of Widom [J. Chem. Phys. 39, 2808 (1963)] insertions at high fluid densities. Optimal settings for the new algorithm are identified and distributions of the quantities involved in it [number of active molecules, energy at the sampled minima of systems 0 and 1, and free-energy differences between subsets gamma(i)(0) and gamma(i)(1) that are mapped onto each other] are explored.     

Stereogenicity/Astereogenicity as Global/Local Permutation-Group Symmetry and Relevant Concepts for Restructuring Stereochemistry

Molecules of ligancy 4 that have been derived from an allene, an ethylene, a tetrahedral, and a square-planar skeleton have been investigated to show that their symmetries are dually and distinctly controlled by point groups and permutation groups. Insomuch as the point-group symmetry was exhibited to control the chirality/achirality of a molecule, sphericity in a molecule, and enantiomeric relationship between molecules [S. Fujita, J. Am. Chem. Soc. 112 (1990) 3390], the permutation-group symmetry has been now clarified to control the stereogenicity of a molecule, tropicity in a molecule, and diastereomeric relationship between molecules. To characterize permutation groups, proper and improper permutations have been defined by comparing proper and improper rotations. Thereby, such permutation groups are classified into stereogenic and astereogenic ones. After a coset representation (CR) of a permutation group has been ascribed to an orbit (equivalence class), the tropicity of the orbit has been defined in term of the global stereogenicity and the local stereogenicity of the CR. As a result, the conventional stereogenicity has now been replaced by the concept local stereogenicity of the present investigation. The terms homotropic, enantiotropic, and hemitropic are coined and used to characterize prostereogenicity. Thus, a molecule is defined as being prostereogenic if it has at least one enantiotropic orbit. Since this definition has been found to be parallel with the definition of prochirality, relevant concepts have been discussed with respect to the parallelism between stereogenicity and chirality in order to restructure the theoretical foundation of stereochemistry and stereoisomerism. The derivation of the skeletons has been characterized by desymmetrization due to the subduction of CRs. The Cahn–Ingold–Prelog (CIP) system has been discussed from the permutational point of view to show that it specifies diastereomeric relationships only. The apparent specification of enantiomeric relationships by the CIP system has been shown to stem from the fact that diastereomeric relationships and enantiomeric ones overlap occasionally in case of tetrahedral molecules.     

A new approach to analysis and display of local lipophilicity/hydrophilicity mapped on molecular surfaces

Summary A new method for display and analysis of lipophilic/hydrophilic properties on molecular surfaces is presented. The present approach is based on the concept of Crippen and coworkers that the overall hydrophobicity of a molecule (measured as the logarithm of the partition coefficient in an octanol/water system) can be obtained as a superposition of single atom contributions. It is also based on the concept of molecular lipophilicity potentials (MLP) first introduced by Audry and coworkers in order to establish a 3D lipophilicity potential profile in the molecular environment. Instead of using a l/r- or an exponential distance law between the atomic coordinates and a point on the molecular surface, a new distance dependency is introduced for the calculation of an MLP-value on the solvent-accessible surface of the molecule. In the present formalism the Crippen values (introduced for atoms in their characteristic structural environment) are projected onto the van der Waals surface of the molecule by a special weighting procedure. This guarantees that only those atomic fragments contribute significantly to the surface values that are in the close neighbourhood of the surface point. This procedure not only works for small molecules but also allows the characterization of the surfaces of biological macromolecules by means of local lipophilicity. Lipophilic and hydrophilic domains can be recognized by visual inspection of computer-generated images or by computational procedures using fuzzy logic strategies. Local hydrophobicities on different molecular surfaces can be quantitatively compared on the basis of the present approach.     

Rubber elasticity and network structure

By defining stretching-invariant subunits of deformation in terms of the van der Waals model it is shown that the density of these subsystems can be smaller than the density of the cross-links. This is observed in short-chain networks with rigid cross-linkages, where due to the broad chain length distribution, very short and therefore stiff network chains appear. It is in any case possible to relate the modulus to the weight fraction of the cross-linking agent if the production of chemical cross-limits per agent molecule is properly defined. The stress-strain behavior of the rubbers can then be fully described within the total range of strains embracing also the upturn at largest elongation.Moreover, on the use of a thermodynamics of swelling of van der Waals networks, a set of swelling experiments drawn from literature can be described: the key assumption is that configurational exclusion effects within the swollen short-chain network must be accounted for. In the case of swelling with flexible molecules there is a need for further empirical modification of the thermodynamical formalism indicating that the configurational abilities within these swollen networks seem also to depend on the relative ratio of the network chain length related to the length of the linear solvent molecules.Dedicated to Prof. Dr. R. Kosfeld on the occasion of his 60th birthday.     

Probing the transition from hydrophilic to hydrophobic solvation with atomic scale resolution

Picosecond and femtosecond X-ray absorption spectroscopy is used to probe the changes of the solvent shell structure upon electron abstraction of aqueous iodide using an ultrashort laser pulse. The transient L(1,3) edge EXAFS at 50 ps time delay points to the formation of an expanded water cavity around the iodine atom, in good agreement with classical and quantum mechanical/molecular mechanics (QM/MM) molecular dynamics (MD) simulations. These also show that while the hydrogen atoms pointed toward iodide, they predominantly point toward the bulk solvent in the case of iodine, suggesting a hydrophobic behavior. This is further confirmed by quantum chemical (QC) calculations of I(-)/I(0)(H(2)O)(n=1-4) clusters. The L(1) edge sub-picosecond spectra point to the existence of a transient species that is not present at 50 ps. The QC calculations and the QM/MM MD simulations identify this transient species as an I(0)(OH(2)) complex inside the cavity. The simulations show that upon electron abstraction most of the water molecules move away from iodine, while one comes closer to form the complex that lives for 3-4 ps. This time is governed by the reorganization of the main solvation shell, basically the time it takes for the water molecules to reform an H-bond network. Only then is the interaction with the solvation shell strong enough to pull the water molecule of the complex toward the bulk solvent. Overall, much of the behavior at early times is determined by the reorientational dynamics of water molecules and the formation of a complete network of hydrogen bonded molecules in the first solvation shell.     

Chemical stereographs: An extension of chemical graphs for representing stereochemical and conformational structures

Chemical stereographs are presented as vehicles for representing qualitative three-dimensional features of molecules that put stereochemical and conformational distinctions in a common graph-theoretic formalism. They extend the concept of a chemical graph by adding tetrads, each qualitatively characterizing the three-dimensional arrangement of four atoms with respect to its clinicity and handedness components. The characterization is sufficiently precise to distinguish synperiplanar, synclinical (gauche), anticlinal and antiperiplanar relationships between vicinal atoms of various conformers. Collectively, the tetrads constitute the embedding graph which presents new possibilities in displaying the stereochemical and conformational features of a molecule. A chemical graph and one of its possible embedding graphs constitute a chemical stereograph. Potential applications of chemical stereographs in the areas of structural representations, molecular symmetry analysis, and stereo-specific substructure searching are discussed.Part of this work was presented by M.J. at the 1989 PaciChem Mathematical Chemistry Minisymposium. The work was funded by The Upjohn Company.     

Models of hydration and isomeric transitions of glucose molecules in aqueous solutions

The kinematic viscosity of aqueous glucose solutions is studied. It is found that the hydrodynamic radius of monosaccharide molecule in an aqueous solution depends on temperature in the range of 290–355 K. Using a bimodal model of the energy states of the volume in which the glucose molecule is located and local equilibrium is established, it is shown that the above-mentioned dependence can be attributed to disturbances in the equilibrium of isomeric transitions, induced by variations in temperature. The parameters of isomeric transitions for a glucose molecule in an aqueous solvent, the probability of “chair” and “boat” configurations occurring for glucose molecules, and the number of water molecules in the hydration shells of these configurations are calculated; the strain of the chemical bonds in the chair configuration of a glucose molecule is estimated.     

Kinetic Energy and the Covalent Bond in H2 +

A modified version of Ruedenbergs innovative analysis of the chemical bond in the hydrogen molecule ion is presented that factors the bond energy into bonding and nonbonding contributions. This simplified approach clearly illustrates Ruedenbergs main thesis: chemical bond formation is driven by a decrease in electron kinetic energy.     

Molecular bonding profiles

We present a refinement of recently proposed characterization of molecules based on a sequence of powers of interatomic separations referred to as molecular profiles. The molecular profiles and closely related shape profiles were based on the averaging contributions arising from different powers of interatomic distances for atoms in a molecule or atoms at the molecular periphery, respectively. Consequently, molecular models in which atoms have the same set of coordinates but different bonding patterns will result in identical molecular profiles. In this article we outline a refinement of molecular profiles in which the bonding pattern in a molecule is fully acknowledged. This is accomplished by adding ghost sites along chemical bonds. The distance-based invariants of the augmented matrix reflect the bonding pattern of a structure giving different molecular profiles for molecules having the same atomic coordinates but different bondings. The procedure is general and applies to two-dimensional and three-dimensional molecular skeletons. Equally, the approach can be applied to van der Waals-type molecular surfaces and molecular contours of equal electron densities in order to obtain characterization of more realistic molecular models.Dedicated to a leading graph theorist, Professor R. C. Read, for sharing his mathematical talents on chemical structures.     

Modeling and conformation analysis of β-cyclodextrin complexes

Summary A series of -cyclodextrin complexes containing various guest molecules was studied using computer-aided molecular modeling and conformation analysis techniques. The geometry of each complex was studied using crystallographic data. The positions of the glycosidic O4 atoms indicate that the -cyclodextrin molecules are elliptically distorted. This distortion can be related to the van der Waals volume of the guest molecules. This correlation is different for aromatic and non-aromatic guest compounds. Rigid body docking experiments demonstrated that in crystal structures the guest molecule occupies a position in the cavity of nearly minimum interaction energy when there are no other molecules having interactions with the guest molecule. From the crystallographic data several rules could be deduced which seem to determine the conformation of -cyclodextrin molecules in complexes. A procedure was developed to construct -cyclodextrin molecules that are able to encompass guest molecules having a given van der Waals volume.     

Conformational energy downward driver (CEDD): Characterization and calibration of the method

Summary A method has been developed that allows one to drive a molecule to conformations of lowest energy given the starting conformation, the identity of the rotatable bonds and the step size. This method has proved useful in our hands in the drug design arena where it is frequently more important to get low-energy conformers of a molecule that match some other (e.g. pharmacophoric or enzyme pocket) requirements than to exhaustively enumerate all possible low-energy conformations for each of the molecules to be studied. The method has been shown to work in the test cases studied to date. Furthermore, so far it has been shown to be sufficiently fast to be used for molecules containing up to 70 rotatable bonds.     

Maximum bond order hybrid orbitals

Summary Based on the simplified calculation scheme of the maximum bond order principle and the basic idea of the maximum overlap symmetry orbital method, a simple procedure is suggested for constructing systematically the bonding hybrid orbitals, called maximum bond order hybrid orbitals, for a given molecule from the first-order density matrix obtained from a molecular orbital calculation. As an example, the proposed procedure is performed for some typical small molecules by use of the density matrix obtained from CNDO/2 calculation. It is shown that the bonding hybrid orbitals constructed by using the procedure are extremely close to those by using the natural hybrid orbital procedure and in good agreement with chemical intuition, and that the proposed procedure can be performed more easily than the natural hybrid orbital procedure and can given simultaneously the values of the maximum bond order for all bonds in molecules.The project was supported by National Natural Science Foundation of China and also supported partly by Foundation of Hubei Education Commission     

Probing intrazeolite space

Molecular sieves are inorganic framework structures generally composed of crystalline aluminosilicate tetrahedra which are arranged to form channels of 2–10 Å, diameters and cages with dimensions from 6–15 Å. Absorption of probe molecules of varying geometries and sizes is used to characterize the framework dimensions and topography in concert with X-ray diffraction identification of the specific structure. From unit cell dimensions and assumptions about the size of the framework forming species, a pore volume can be calculated. The volumes of the absorbed probe molecules, using their liquid densities, are then compared to the calculated pore volume. The constraint on the packing of the absorbed molecules is quantified by comparing their packing density to their density in the liquid state. Further, the packing of different probe molecules into the same pore volume is compared via a ratio technique called the packing ratio. The effect of the lattice geometry and framework dimensions on the packing ratios is to provide a set of characteristic values for a given molecular sieve. The packing ratios for zeolites rho and ZSM-5 are presented as expectation values for other scientists to use as bases of comparison.     

Classification of water molecules in protein binding sites

Water molecules play a crucial role in mediating the interaction between a ligand and a macromolecular receptor. An understanding of the nature and role of each water molecule in the active site of a protein could greatly increase the efficiency of rational drug design approaches: if the propensity of a water molecule for displacement can be determined, then synthetic effort may be most profitably applied to the design of specific ligands with the displacement of this water molecule in mind. In this paper, a thermodynamic analysis of water molecules in the binding sites of six proteins, each complexed with a number of inhibitors, is presented. Two classes of water molecules were identified: those conserved and not displaced by any of the ligands, and those that are displaced by some ligands. The absolute binding free energies of 54 water molecules were calculated using the double decoupling method, with replica exchange thermodynamic integration in Monte Carlo simulations. It was found that conserved water molecules are on average more tightly bound than displaced water molecules. In addition, Bayesian statistics is used to calculate the probability that a particular water molecule may be displaced by an appropriately designed ligand, given the calculated binding free energy of the water molecule. This approach therefore allows the numerical assessment of whether or not a given water molecule should be targeted for displacement as part of a rational drug design strategy.     

Structure of hydrated Na-Nafion polymer membranes

We use molecular dynamics simulations to investigate the structure of the hydrated Na-Nafion membranes. The membrane is "prepared" by starting with the Nafion chains placed on a cylinder having the water inside it. Minimizing the energy of the system leads to a filamentary hydrophilic domain whose structure depends on the degree of hydration. At 5 wt % water the system does not have enough water molecules to solvate all the ions that could be formed by the dissociation of the -SO3Na groups. As a result, the -SO3Na groups aggregate with the water to form very small droplets that do not join into a continuous phase. The size of the droplets is between 5 and 8 A. As the amount of water present in the membrane is increased, the membrane swells, and SO3Na has an increasing tendency to dissociate into ions. Furthermore, a transition to a percolating hydrophilic network is observed. In the percolating structure, the water forms irregular curvilinear channels branching in all directions. The typical dimension of the cross section of these channels is about 10-20 A. Calculated neutron scattering from the simulated system is in qualitative agreement with experiment. In all simulations, the pendant sulfonated perfluorovinyl side chains of the Nafion hug the walls of the hydrophilic channel, while the sulfonate groups point toward the center of the hydrophilic phase. The expulsion of the side chains from the hydrophilic domain is favored because it allows better interaction between the water molecules. We have also examined the probability of finding water molecules around the Na+ and the -SO3(-) ions as well as the probability of finding other water molecules next to a given water molecule. These probabilities are much broader than those found in bulk water or for one ion in bulk water (calculated with the potentials used in the present simulation). This is due to the highly inhomogeneous nature of the material contained in the small hydrophilic pores.     

An approach to computer-aided inhibitor design: Application to cathepsin L

Summary We have developed an approach to search for molecules that can be used as lead compounds in designing an inhibitor for a given proteolytic enzyme when the 3D structure of a homologous protein is known. This approach is based on taking the cast of the binding pocket of the protease and comparing its dimensions with that of the dimensions of small molecules. Herein the 3D structure of papain is used to model cathepsin L using the comparative modeling technique. The cast of the binding pocket is computed using the crystal structure of papain because the structures of papain and the model of cathepsin L are found to be similar at the binding site. The dimensions of the cast of the binding site of papain are used to screen for molecules from the Cambridge Structural Database (CSD) of small molecules. Twenty molecules out of the 80 000 small molecules in the CSD are found to have dimensions that are accommodated by the papain binding pocket. Visual comparison of the shapes of the cast and the 20 screened molecules resulted in identifying brevotoxin b, a toxin isolated from the red tide dinoflagellate Ptycho brevis (previously classified as Gymonodium breve), as the structure that best fits the binding pocket of papain. We tested the proteolytic activity of papain and cathepsin L in the presence of brevotoxin b and found inhibition of papain and cathepsin L with K_is of 25 M and 0.6 M, respectively. We also compare our method with a more elaborate method in the literature, by presenting our results on the computer search for inhibitors of the HIV-1 protease.     

Effect of Cyclodextrins on the Chemical Stability of ST1435, a Contraceptive Steroid Progestin, in Aqueous Solution

The influence of cyclodextrins (CDs) on the chemical stability of the contraceptive steroid progestin, ST1435, in aqueous solution has been studied using reversed phase high performance liquid chromatography. The effects of CD structure, temperature, and CD concentration on the rate of degradation were investigated. It was found that the drug degraded to different extents following a pseudo-first order reaction mechanism. The presence of the host molecules affected the degradation rate as a result of complexation which might result in protection of the labile moiety of the drug molecule against degradation. Hydroxypropyl--cyclodextrin (HP--CD) and hydroxyethyl--cyclodextrin (HE--CD) retarded the degradation in contrast to -cyclodextrin (-CD) which accelerated the steroid degradation. The stabilizing action of HP--CD is larger than that of HE--CD. The degradation rate increased upon increasing temperature and the Arrhenius equation is valid. Lineweaver-Burk equation analysis indicated that the steroid included inside the CD cavity degraded three times more slowly than did the free ST1435 in solution. This equation further supported the formation of a 1 : 1 inclusion complex between ST1435 and HP--CD with a stability constant of 934.5 M-1 at 65°C.     

The Energy of Interaction between Water Monolayer and Carbon Surface and the Determination of H-Bond Energy from the Isotherms of Water Desorption by Organic Molecules into an Aqueous Solution

It is proven that, when passing from a liquid into an adsorption phase on a carbon surface, the maximal number of H-bonds between water molecules decreases from four to three because of the screening of one unpaired electron of the oxygen atom of an adsorbed water molecule by aromatic rings of the carbon surface. An energy gain equal to the energy of one H-bond arises upon water desorption by organic molecules adsorbed from an aqueous solution. The ratio between the number of H-bonds of a group of water molecules, which is displaced into a solution by one organic molecule, in the solution and in the adsorption phase is independent of the number of molecules in this group and is, on the average, equal to 2.038 for all possible structures of H-bonds in both phases. The allowance for this ratio in the isotherm of water desorption into a solution and the introduction of a coefficient, which depends on the relative water content ( ) in the adsorption phase, in the form of into the equation of the desorption isotherm make it possible to determine the balance of the change in the Gibbs energy at the desorption equilibrium and the standard Gibbs energy = –1.76 kJ/mol of water desorption into a solution from a carbon surface. This value determined by an independent method is = –1.79 kJ/mol; i.e., both values are close to each other. The RTlnf energy of the additional H-bond, which is formed between water molecules upon passing from the adsorption phase into the solution, is found by the extrapolation of the isotherms of water desorption by molecules of several benzene derivatives. This energy ranges from 9.13 to 9.24 kJ/mol, thus corresponding to the energy of one H-bond, as measured by IR spectroscopy and NMR.     

Interaction of low-energy ions and atoms of light elements with a fluorinated carbon molecular lattice

The mechanism of interaction of low-energy atoms and ions of light elements (H, H+, He, Li, the kinetic energy of the particles 2-40 eV) with C6H6, C6F12, C60, and C60F48 molecules was studied by ab initio MD simulations and quantum-chemical calculations. It was shown that starting from 6 A from the carbon skeleton for the "C6H6 + proton" and "C60 + proton" systems, the electronic charge transfer from the aromatic molecule to H+ occurs with a probability close to 1. The process transforms the H+ to a hydrogen atom and the neutral C6H6 and C60 molecules to cation radicals. The mechanism of interaction of low-energy protons with C6F12 and C60F48 molecules has a substantially different character and can be considered qualitatively as the interaction between a neutral molecule and a point charge. The Coulomb perturbation of the system arising from the interaction of the uncompensated proton charge with the Mulliken charges of fluorine atoms results in an inversion of the energies of the electronic states localized on the proton and on the C6F12 and C60F48 molecules and makes the electronic charge transfer energetically unfavorable. On the different levels of theory, the barriers of the proton penetration for the C6F12 and C60F48 molecules are from two to four times lower than those for the corresponding parent systems (C6H6 and C60). The penetration barriers of the He atom and Li+ ion depend mainly on the effective radii of the bombarding particles. The theoretical penetration and escaped barriers for the "Li+ + C60" process qualitatively explain the experimental conditions of synthesis of the Li@C60 complex.     

A Monte Carlo simulation of energy transfer processes in thermal unimolecular reactions. II. An applications to polyatomic dissociation

The Monte Carlo method has been used to provide a numerical solution to the ro-vibrational master equation for the low pressure unimolecular decomposition of a polyatomic molecule. This type of solution is made possible through the use of a simple exponential transition probability function, that represents the efficiency with which energy transfer takes place between the reactant molecule and an unspecified heat bath gas. The Monte Carlo technique is used to generate random variables that are distributed in a manner prescribed by the transition probability function. In the case of the present simulation, these variables correspond to random energy jumps induced in the molecule through single collision events. In order to account for the energy dependence of the vibrational state densities, we have proposed that vibrational relaxation in the polyatomic takes place from a single vibrational mode. Under equilibrium conditions we are able to show that with this assumption, the Monte Carlo model is capable of reproducing molecular quantities, such as the average vibrational energy per molecule and the vibrational specific heat, that compare favourable with the corresponding values calculated from equilibrium statistical mechanics. The model has been applied to a study of the low pressure unimolecular decomposition of a series of polyatomics. For three of the molecules, CH₄, CD₄, and C₂H₆ the agreement between the calculated and the high temperature experimental rate constants is very good. The calculations indicate that a significant proportion of the molecules that dissociate are rotationally as well as vibrationally excited. Very few of the reactive molecules have a vibrational energy content equal to or greater than E₀, the dissociation energy. The extent of rotational excitation is found to be temperature dependent.