Total free energy analysis of fully hydrated proteins.

The total free energy of a hydrated biomolecule and its corresponding decomposition of energy and entropy provides detailed information about regions of thermodynamic stability or instability. The free energies of four hydrated globular proteins with different net charges are calculated from a molecular dynamics simulation, with the energy coming from the system Hamiltonian and entropy using multiscale cell correlation. Water is found to be most stable around anionic residues, intermediate around cationic and polar residues, and least stable near hydrophobic residues, especially when more buried, with stability displaying moderate entropy-enthalpy compensation. Conversely, anionic residues in the proteins are energetically destabilized relative to singly solvated amino acids, while trends for other residues are less clear-cut. Almost all residues lose intraresidue entropy when in the protein, enthalpy changes are negative on average but may be positive or negative, and the resulting overall stability is moderate for some proteins and negligible for others. The free energy of water around single amino acids is found to closely match existing hydrophobicity scales. Regarding the effect of secondary structure, water is slightly more stable around loops, of intermediate stability around ß strands and turns, and least stable around helices. An interesting asymmetry observed is that cationic residues stabilize a residue when bonded to its N-terminal side but destabilize it when on the C-terminal side, with a weaker reversed trend for anionic residues.

Energy-entropy multiscale cell correlation method to predict toluene-water log P in the SAMPL9 challenge.

The energy-entropy multiscale cell correlation (EE-MCC) method is used to calculate toluene-water log P values of 16 drug molecules in the SAMPL9 physical properties challenge. EE-MCC calculates the free energy, energy and entropy from molecular dynamics (MD) simulations of the water and toluene solutions. Specifically, MCC evaluates entropy by partitioning the system into cells of correlated atoms at multiple length scales and further partitioning the local coordinates into energy wells, yielding vibrational and topographical terms from the energy-well sizes and probabilities. The log P values calculated by EE-MCC using three 200 ns MD simulations have a mean average error of 0.82 and standard error of the mean of 0.97 versus experiment, which is comparable with the best methods entered in SAMPL9. The main contribution to log P is from energy. Less polar drugs have more favourable energies of transfer. The entropy of transfer consists of increased solute vibrational and conformational terms in toluene due to weaker interactions, fewer solute positions in the larger-molecule solvent, reduced water vibrational entropy, negligible change in toluene vibrational entropy, and gains in solvent orientational entropy. The solvent entropy contributions here may be slightly underestimated because software limitations and statistical fluctuations meant that only the first shell could be included while averaged over the whole solution. Nonetheless, such issues will be addressed in future software to offer a general method to calculate entropy directly from MD simulation and to provide molecular understanding or guide system design.

What Determines the Selectivity of Arginine Dihydroxylation by the Nonheme Iron Enzyme OrfP?

The nonheme iron enzyme OrfP reacts with l-Arg selectively to form the 3R,4R-dihydroxyarginine product, which in mammals can inhibit the nitric oxide synthase enzymes involved in blood pressure control. To understand the mechanisms of dioxygen activation of l-Arg by OrfP and how it enables two sequential oxidation cycles on the same substrate, we performed a density functional theory study on a large active site cluster model. We show that substrate binding and positioning in the active site guides a highly selective reaction through C3 -H hydrogen atom abstraction. This happens despite the fact that the C3 -H and C4 -H bond strengths of l-Arg are very similar. Electronic differences in the two hydrogen atom abstraction pathways drive the reaction with an initial C3 -H activation to a low-energy 5 σ-pathway, while substrate positioning destabilizes the C4 -H abstraction and sends it over the higher-lying 5 π-pathway. We show that substrate and monohydroxylated products are strongly bound in the substrate binding pocket and hence product release is difficult and consequently its lifetime will be long enough to trigger a second oxygenation cycle.

Energy-entropy prediction of octanol-water logP of SAMPL7 N-acyl sulfonamide bioisosters.

Partition coefficients quantify a molecule's distribution between two immiscible liquid phases. While there are many methods to compute them, there is not yet a method based on the free energy of each system in terms of energy and entropy, where entropy depends on the probability distribution of all quantum states of the system. Here we test a method in this class called Energy Entropy Multiscale Cell Correlation (EE-MCC) for the calculation of octanol-water logP values for 22 N-acyl sulfonamides in the SAMPL7 Physical Properties Challenge (Statistical Assessment of the Modelling of Proteins and Ligands). EE-MCC logP values have a mean error of 1.8 logP units versus experiment and a standard error of the mean of 1.0 logP units for three separate calculations. These errors are primarily due to getting sufficiently converged energies to give accurate differences of large numbers, particularly for the large-molecule solvent octanol. However, this is also an issue for entropy, and approximations in the force field and MCC theory also contribute to the error. Unique to MCC is that it explains the entropy contributions over all the degrees of freedom of all molecules in the system. A gain in orientational entropy of water is the main favourable entropic contribution, supported by small gains in solute vibrational and orientational entropy but offset by unfavourable changes in the orientational entropy of octanol, the vibrational entropy of both solvents, and the positional and conformational entropy of the solute.

Energy-entropy method using multiscale cell correlation to calculate binding free energies in the SAMPL8 host-guest challenge.

Free energy drives a wide range of molecular processes such as solvation, binding, chemical reactions and conformational change. Given the central importance of binding, a wide range of methods exist to calculate it, whether based on scoring functions, machine-learning, classical or electronic structure methods, alchemy, or explicit evaluation of energy and entropy. Here we present a new energy-entropy (EE) method to calculate the host-guest binding free energy directly from molecular dynamics (MD) simulation. Entropy is evaluated using Multiscale Cell Correlation (MCC) which uses force and torque covariance and contacts at two different length scales. The method is tested on a series of seven host-guest complexes in the SAMPL8 (Statistical Assessment of the Modeling of Proteins and Ligands) "Drugs of Abuse" Blind Challenge. The EE-MCC binding free energies are found to agree with experiment with an average error of 0.9 kcal mol-1. MCC makes clear the origin of the entropy changes, showing that the large loss of positional, orientational, and to a lesser extent conformational entropy of each binding guest is compensated for by a gain in orientational entropy of water released to bulk, combined with smaller decreases in vibrational entropy of the host, guest and contacting water.

Convergence behaviour of solvation shells in simulated liquids.

A convenient way to analyse solvent structure around a solute is to use solvation shells, whereby solvent position around the solute is discretised by the size of a solvent molecule, leading to multiple shells around the solute. The two main ways to define multiple shells around a solute are either directly with respect to the solute, called solute-centric, or locally for both solute and solvent molecules alike. It might be assumed that both methods lead to solvation shells with similar properties. However, our analysis suggests otherwise. Solvation shells are analysed in a series of simulations of five pure liquids of differing polarity. Shells are defined locally working outwards from each molecule treated as a reference molecule using two methods: the cutoff at the first minimum in the radial distribution function and the parameter-free Relative Angular Distance method (RAD). The molecular properties studied are potential energy, coordination number and coordination radius. Rather than converging to bulk values, as might be expected for pure solvents, properties are found to deviate as a function of shell index. This behaviour occurs because molecules with larger coordination numbers and radius have more neighbours, which make them more likely to be connected to the reference molecule via fewer shells. The effect is amplified for RAD because of its more variable coordination radii and for water with its more open structure and stronger interactions. These findings indicate that locally defined shells should not be thought of as directly comparable to solute-centric shells or to distance. As well as showing how box size and cutoff affect the non-convergence, to restore convergence we propose a hybrid method by defining a new set of shells with boundaries at the uppermost distance of each locally derived shell.

How Do Electrostatic Perturbations of the Protein Affect the Bifurcation Pathways of Substrate Hydroxylation versus Desaturation in the Nonheme Iron-Dependent Viomycin Biosynthesis Enzyme?

The viomycin biosynthesis enzyme VioC is a nonheme iron and α-ketoglutarate-dependent dioxygenase involved in the selective hydroxylation of l-arginine at the C3-position for antibiotics biosynthesis. Interestingly, experimental studies showed that using the substrate analogue, namely, l-homo-arginine, a mixture of products was obtained originating from C3-hydroxylation, C4-hydroxylation, and C3-C4-desaturation. To understand how the addition of one CH2 group to a substrate can lead to such a dramatic change in selectivity and activity, we decided to perform a computational study using quantum mechanical (QM) cluster models. We set up a large active-site cluster model of 245 atoms that includes the oxidant with its first- and second-coordination sphere influences as well as the substrate binding pocket. The model was validated against experimental work from the literature on related enzymes and previous computational studies. Thereafter, possible pathways leading to products and byproducts were investigated for a model containing l-Arg and one for l-homo-Arg as substrate. The calculated free energies of activation predict product distributions that match the experimental observation and give a low-energy C3-hydroxylation pathway for l-Arg, while for l-homo-Arg, several barriers are found to be close in energy leading to a mixture of products. We then analyzed the origins of the differences in product distributions using thermochemical, valence bond, and electrostatic models. Our studies show that the C3-H and C4-H bond strengths of l-Arg and l-homo-Arg are similar; however, external perturbations from an induced electric field of the protein affect the relative C-H bond strengths of l-Arg dramatically and make the C3-H bond the weakest and guide the reaction to a selective C3-hydroxylation channel. Therefore, the charge distribution in the protein and the induced electric dipole field of the active site of VioC guides the l-Arg substrate activation to C3-hydroxylation and disfavors the C4-hydroxylation pathway, while this does not occur for l-homo-Arg. Tight substrate positioning and electrostatic perturbations from the second-coordination sphere residues in VioC also result in a slower overall reaction for l-Arg; however, they enable a high substrate selectivity. Our studies highlight the importance of the second-coordination sphere in proteins that position the substrate and oxidant, perturb charge distributions, and enable substrate selectivity.

Lignin Biodegradation by a Cytochrome P450 Enzyme: A Computational Study into Syringol Activation by GcoA.

A recently characterized cytochrome P450 isozyme GcoA activates lignin components through a selective O-demethylation or alternatively an acetal formation reaction. These are important reactions in biotechnology and, because lignin is readily available; it being the main component in plant cell walls. In this work we present a density functional theory study on a large active site model of GcoA to investigate syringol activation by an iron(IV)-oxo heme cation radical oxidant (Compound I) leading to hemiacetal and acetal products. Several substrate-binding positions were tested and full energy landscapes calculated. The study shows that substrate positioning determines the product distributions. Thus, with the phenol group pointing away from the heme, an O-demethylation is predicted, whereas an initial hydrogen-atom abstraction of the weak phenolic O-H group would trigger a pathway leading to ring-closure to form acetal products. Predictions on how to engineer P450 GcoA to get more selective product distributions are given.

Model for Counterion Binding and Charge Reversal on Protein Surfaces.

The structural stability and solubility of proteins in liquid therapeutic formulations is important, especially since new generations of therapeutics are designed for efficacy before consideration of stability. We introduce an electrostatic binding model to measure the net charge of proteins with bound ions in solution. The electrostatic potential on a protein surface is used to separately group together acidic and basic amino acids into patches, which are then iteratively bound with oppositely charged counterions. This model is aimed toward formulation chemists for initial screening of a range of conditions prior to lab-work. Computed results compare well with experimental zeta potential measurements from the literature covering a range of solution conditions. Importantly, the binding model reproduces the charge reversal phenomenon that is observed with polyvalent ion binding to proteins and its dependence on ion charge and concentration. Intriguingly, protein sequence can be used to give similarly good agreement with experiment as protein structure, interpreted as resulting from the close proximity of charged side chains on a protein surface. Further, application of the model to human proteins suggests that polyanion binding and overcharging, including charge reversal for cationic proteins, is a general feature. These results add to evidence that addition of polyanions to protein formulations could be a general mechanism for modulating solution stability.

Cross-linking of aromatic phenolate groups by cytochrome P450 enzymes: a model for the biosynthesis of vancomycin by OxyB.

The cytochromes P450 are a versatile class of enzymes involved in many chemical reactions in biosystems and as such they take part in biodegradation as well as biosynthesis pathways in many organisms. These enzymes use molecular oxygen on a heme centre and often react as mono-oxygenases. Lesser known reactions catalyzed by the P450s include desaturation pathways and ring-closure reactions. In this work we study the aromatic cross-linking of glycopeptide units as, for instance, performed by the P450 isozyme OxyB as part of vancomycin biosynthesis. A series of density functional theory studies are reported on a large active site cluster model of 258 atoms containing the heme with its coordinated ligands, a representative substrate and its interacting protein residues. We show that the catalytic cycle intermediates Compound I and Compound II of P450 can rapidly and successively abstract a phenolic hydrogen atom from adjacent peptide groups to give a biradical intermediate with small reaction barriers. The latter can form the ether cross-link between the two aromatic residues, which is the rate-determining step in the reaction mechanism and involves a simultaneous proton transfer from the ipso-position to the ketone. A thermochemical analysis reveals that weak phenolic O-H bonds lead to hydrogen atom abstraction easily by Compound I and Compound II, enabling a selective aromatic cross-linking reaction.

Entropy of Proteins Using Multiscale Cell Correlation.

A new multiscale method is presented to calculate the entropy of proteins from molecular dynamics simulations. Termed Multiscale Cell Correlation (MCC), the method decomposes the protein into sets of rigid-body units based on their covalent-bond connectivity at three levels of hierarchy: molecule, residue, and united atom. It evaluates the vibrational and topographical entropy from forces, torques, and dihedrals at each level, taking into account correlations between sets of constituent units that together make up a larger unit at the coarser length scale. MCC gives entropies in close agreement with normal-mode analysis and smaller than those using quasiharmonic analysis as well as providing much faster convergence. Moreover, MCC provides an insightful decomposition of entropy at each length scale and for each type of amino acid according to their solvent exposure and whether they are terminal residues. While the residue entropy depends weakly on solvent exposure, there is greater variation in entropy components for larger, more polar amino acids, which have increased conformational entropy but reduced vibrational entropy with greater solvent exposure.

Extending the Martini Coarse-Grained Force Field to N-Glycans.

Glycans play a vital role in a large number of cellular processes. Their complex and flexible nature hampers structure-function studies using experimental techniques. Molecular dynamics (MD) simulations can help in understanding dynamic aspects of glycans if the force field parameters used can reproduce key experimentally observed properties. Here, we present optimized coarse-grained (CG) Martini force field parameters for N-glycans, calibrated against experimentally derived binding affinities for lectins. The CG bonded parameters were obtained from atomistic (ATM) simulations for different glycan topologies including high mannose and complex glycans with various branching patterns. In the CG model, additional elastic networks are shown to improve maintenance of the overall conformational distribution. Solvation free energies and octanol-water partition coefficients were also calculated for various N-glycan disaccharide combinations. When using standard Martini nonbonded parameters, we observed that glycans spontaneously aggregated in the solution and required down-scaling of their interactions for reproduction of ATM model radial distribution functions. We also optimized the nonbonded interactions for glycans interacting with seven lectin candidates and show that a relatively modest scaling down of the glycan-protein interactions can reproduce free energies obtained from experimental studies. These parameters should be of use in studying the role of glycans in various glycoproteins and carbohydrate binding proteins as well as their complexes, while benefiting from the efficiency of CG sampling.

Entropy of Simulated Liquids Using Multiscale Cell Correlation.

Accurately calculating the entropy of liquids is an important goal, given that many processes take place in the liquid phase. Of almost equal importance is understanding the values obtained. However, there are few methods that can calculate the entropy of such systems, and fewer still to make sense of the values obtained. We present our multiscale cell correlation (MCC) method to calculate the entropy of liquids from molecular dynamics simulations. The method uses forces and torques at the molecule and united-atom levels and probability distributions of molecular coordinations and conformations. The main differences with previous work are the consistent treatment of the mean-field cell approximation to the approriate degrees of freedom, the separation of the force and torque covariance matrices, and the inclusion of conformation correlation for molecules with multiple dihedrals. MCC is applied to a broader set of 56 important industrial liquids modeled using the Generalized AMBER Force Field (GAFF) and Optimized Potentials for Liquid Simulations (OPLS) force fields with 1.14*CM1A charges. Unsigned errors versus experimental entropies are 8.7 J K - 1 mol - 1 for GAFF and 9.8 J K - 1 mol - 1 for OPLS. This is significantly better than the 2-Phase Thermodynamics method for the subset of molecules in common, which is the only other method that has been applied to such systems. MCC makes clear why the entropy has the value it does by providing a decomposition in terms of translational and rotational vibrational entropy and topographical entropy at the molecular and united-atom levels.

Overcoming the limitations of cutoffs for defining atomic coordination in multicomponent systems.

A common way to understand structure in multimolecular systems is the coordination shell which comprises all the neighbors of an atom. Coordination, however, is nontrivial to determine because there is no obvious way to determine when atoms are neighbors. A common solution is to take all atoms within a cutoff at the first minimum of the radial distribution function, g(r). We show that such an approach cannot be consistently applied to model multicomponent systems, namely mixtures of atoms differing in size or charge. Coordination shells using the total g(r) are found to be too restrictive for atoms of different size while those using pairwise g(r)s are excessive for charged mixtures. The recently introduced relative angular distance algorithm, however, which defines coordination instantaneously from atomic positions, is consistently able to define coordination shells containing the expected neighboring atoms for all these systems. This more robust way to determine coordination should in turn make coordination a more robust way to understand structure. © 2017 Wiley Periodicals, Inc.

Protons and Hydroxide Ions in Aqueous Systems.

Understanding the structure and dynamics of water's constituent ions, proton and hydroxide, has been a subject of numerous experimental and theoretical studies over the last century. Besides their obvious importance in acid-base chemistry, these ions play an important role in numerous applications ranging from enzyme catalysis to environmental chemistry. Despite a long history of research, many fundamental issues regarding their properties continue to be an active area of research. Here, we provide a review of the experimental and theoretical advances made in the last several decades in understanding the structure, dynamics, and transport of the proton and hydroxide ions in different aqueous environments, ranging from water clusters to the bulk liquid and its interfaces with hydrophobic surfaces. The propensity of these ions to accumulate at hydrophobic surfaces has been a subject of intense debate, and we highlight the open issues and challenges in this area. Biological applications reviewed include proton transport along the hydration layer of various membranes and through channel proteins, problems that are at the core of cellular bioenergetics.

Locally adaptive method to define coordination shell.

An algorithm is presented to define a particle's coordination shell for any collection of particles. It requires only the particles' positions and no pre-existing knowledge or parameters beyond those already in the force field. A particle's shell is taken to be all particles that are not blocked by any other particle and not further away than a blocked particle. Because blocking is based on two distances and an angle for triplets of particles, it is called the relative angular distance (RAD) algorithm. RAD is applied to Lennard-Jones particles in molecular dynamics simulations of crystalline, liquid, and gaseous phases at various temperatures and densities. RAD coordination shells agree well with those from a cut-off in the radial distribution function for the crystals and liquids and are slightly higher for the gas.

Instantaneous, parameter-free methods to define a solute's hydration shell.

A range of methods are presented to calculate a solute's hydration shell from computer simulations of dilute solutions of monatomic ions and noble gas atoms. The methods are designed to be parameter-free and instantaneous so as to make them more general, accurate, and consequently applicable to disordered systems. One method is a modified nearest-neighbor method, another considers solute-water Lennard-Jones overlap followed by hydrogen-bond rearrangement, while three methods compare various combinations of water-solute and water-water forces. The methods are tested on a series of monatomic ions and solutes and compared with the values from cutoffs in the radial distribution function, the nearest-neighbor distribution functions, and the strongest-acceptor hydrogen bond definition for anions. The Lennard-Jones overlap method and one of the force-comparison methods are found to give a hydration shell for cations which is in reasonable agreement with that using a cutoff in the radial distribution function. Further modifications would be required, though, to make them capture the neighboring water molecules of noble-gas solutes if these weakly interacting molecules are considered to constitute the hydration shell.

Long-range hydrogen-bond structure in aqueous solutions and the vapor-water interface.

There is a considerable disagreement about the extent to which solutes perturb water structure. On the one hand, studies that analyse structure directly only show local structuring in a solute's first and possibly second hydration shells. On the other hand, thermodynamic and kinetic data imply indirectly that structuring occurs much further away. Here, the hydrogen-bond structure of water around halide anions, alkali cations, noble-gas solutes, and at the vapor-water interface is examined using molecular dynamics simulations. In addition to the expected perturbation in the first hydration shell, deviations from bulk behavior are observed at longer range in the rest of the simulation box. In particular, at the longer range, there is an excess of acceptors around halide anions, an excess of donors around alkali cations, weakly enhanced tetrahedrality and an oscillating excess and deficiency of donors and acceptors around noble-gas solutes, and enhanced tetrahedrality at the vapor-water interface. The structuring compensates for the short-range perturbation in water-water hydrogen bonds induced by the solute. Rather than being confined close to the solute, it is spread over as many water molecules as possible, presumably to minimize the perturbation to each water molecule.

Thermodynamic Origin of Differential Excipient-Lysozyme Interactions.

Understanding the intricate interplay of interactions between proteins, excipients, ions and water is important to achieve the effective purification and stable formulation of protein therapeutics. The free energy of lysozyme interacting with two kinds of polyanionic excipients, citrate and tripolyphosphate, together with sodium chloride and TRIS-buffer, are analysed in multiple-walker metadynamics simulations to understand why tripolyphosphate causes lysozyme to precipitate but citrate does not. The resulting multiscale decomposition of energy and entropy components for water, sodium chloride, excipients and lysozyme reveals that lysozyme is more stabilised by the interaction of tripolyphosphate with basic residues. This is accompanied by more sodium ions being released into solution from tripolyphosphate than for citrate, whilst the latter instead has more water molecules released into solution. Even though lysozyme aggregation is not directly probed in this study, these different mechanisms are suspected to drive the cross-linking between lysozyme molecules with vacant basic residues, ultimately leading to precipitation.

Inhibitors of PIM-1 kinase: a computational analysis of the binding free energies of a range of imidazo [1,2-b] pyridazines.

The binding of a selection of competitive imidazo [1,2-b] pyridazine inhibitors of PIM-1 kinase with nanomolar activity has been analyzed using computational methods. Molecular dynamics simulations using umbrella sampling to determine a potential of mean force have been used to accurately predict the relative free energies of binding of these inhibitors, from -4.3 to -9.5 kcal mol(-1), in excellent agreement with the trends observed in previous experimental assays. The relative activity of the inhibitors could not be accounted for by any single effect or interaction within the active site and could only be fully reproduced when the overall free energies were considered, including important contributions from interactions outside the hinge region and using explicit solvent in the active site. The potential of mean force for the displacement of the glycine-rich phosphate binding loop (P-loop) has also been estimated and shown to be an important feature in the binding of these ligands.