Calculation of free-energy differences by confinement simulations. Application to peptide conformers

Optimal Transport for Free Energy Estimation

Optimal transport theory is a growing field of mathematics, which has recently found many applications. Here we take advantage of optimal transport for computational free energy estimation. We show analytically, and then via simulation, that this approach is effective in terms of optimizing the barriers of an alchemical transformation.

Dimer dissociation is a key energetic event in the fold-switch pathway of KaiB

Cyanobacteria possesses the simplest circadian clock, composed of three proteins that act as a phosphorylation oscillator: KaiA, KaiB, and KaiC. The timing of this oscillator is determined by the fold-switch of KaiB, a structural rearrangement of its C-terminal half that is accompanied by a change in the oligomerization state. During the day, KaiB forms a stable tetramer (gsKaiB), whereas it adopts a monomeric thioredoxin-like fold during the night (fsKaiB). Although the structures and functions of both native states are well studied, little is known about the sequence and structure determinants that control their structural interconversion. Here, we used confinement molecular dynamics (CCR-MD) and folding simulations using structure-based models to show that the dissociation of the gsKaiB dimer is a key energetic event for the fold-switch. Hydrogen-deuterium exchange mass spectrometry (HDXMS) recapitulates the local stability of protein regions reported by CCR-MD, with both approaches consistently indicating that the energy and backbone flexibility changes are solely associated with the region that fold-switches between gsKaiB and fsKaiB and that the localized regions that differentially stabilize gsKaiB also involve regions outside the dimer interface. Moreover, two mutants (R23C and R75C) previously reported to be relevant for altering the rhythmicity of the Kai clock were also studied by HDXMS. Particularly, R75C populates dimeric and monomeric states with a deuterium incorporation profile comparable to the one observed for fsKaiB, emphasizing the importance of the oligomerization state of KaiB for the fold-switch. These findings suggest that the information necessary to control the rhythmicity of the cyanobacterial biological clock is, to a great extent, encoded within the KaiB sequence.

Thermodynamics and Kinetics of Drug-Target Binding by Molecular Simulation

Computational studies play an increasingly important role in chemistry and biophysics, mainly thanks to improvements in hardware and algorithms. In drug discovery and development, computational studies can reduce the costs and risks of bringing a new medicine to market. Computational simulations are mainly used to optimize promising new compounds by estimating their binding affinity to proteins. This is challenging due to the complexity of the simulated system. To assess the present and future value of simulation for drug discovery, we review key applications of advanced methods for sampling complex free-energy landscapes at near nonergodicity conditions and for estimating the rate coefficients of very slow processes of pharmacological interest. We outline the statistical mechanics and computational background behind this research, including methods such as steered molecular dynamics and metadynamics. We review recent applications to pharmacology and drug discovery and discuss possible guidelines for the practitioner. Recent trends in machine learning are also briefly discussed. Thanks to the rapid development of methods for characterizing and quantifying rare events, simulation's role in drug discovery is likely to expand, making it a valuable complement to experimental and clinical approaches.

About the Protein Space Vastness

An accurate estimation of the Protein Space size, in light of the factors that govern it, is a long-standing problem and of paramount importance in evolutionary biology, since it determines the nature of protein evolvability. A simple analysis will enable us to, firstly, reduce an unrealistic Protein Space size of ~ 10¹³⁰ sequences, for a 100-residues polypeptide chain, to ~ 10⁹ functional proteins and, secondly, estimate a robust average-mutation rate per amino acid (ξ ~ 1.23) and infer from it, in light of the protein marginal stability, that only a fraction of the sequence will be available at any one time for a functional protein to evolve. Although this result does not solve the Protein Space vastness problem frames it in a more rational one and illustrates the impact of the marginal stability on protein evolvability.

Smooth orientation-dependent scoring function for coarse-grained protein quality assessment

MOTIVATION

Protein quality assessment (QA) is a crucial element of protein structure prediction, a fundamental and yet open problem in structural bioinformatics. QA aims at ranking predicted protein models to select the best candidates. The assessment can be performed based either on a single model or on a consensus derived from an ensemble of models. The latter strategy can yield very high performance but substantially depends on the pool of available candidate models, which limits its applicability. Hence, single-model QA methods remain an important research target, also because they can assist the sampling of candidate models.

RESULTS

We present a novel single-model QA method called SBROD. The SBROD (Smooth Backbone-Reliant Orientation-Dependent) method uses only the backbone protein conformation, and hence it can be applied to scoring coarse-grained protein models. The proposed method deduces its scoring function from a training set of protein models. The SBROD scoring function is composed of four terms related to different structural features: residue-residue orientations, contacts between backbone atoms, hydrogen bonding and solvent-solute interactions. It is smooth with respect to atomic coordinates and thus is potentially applicable to continuous gradient-based optimization of protein conformations. Furthermore, it can also be used for coarse-grained protein modeling and computational protein design. SBROD proved to achieve similar performance to state-of-the-art single-model QA methods on diverse datasets (CASP11, CASP12 and MOULDER).

AVAILABILITY AND IMPLEMENTATION

The standalone application implemented in C++ and Python is freely available at https://gitlab.inria.fr/grudinin/sbrod and supported on Linux, MacOS and Windows.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Solubility Advantage of Amorphous Ketoprofen. Thermodynamic and Kinetic Aspects by Molecular Dynamics and Free Energy Approaches

Thermodynamic and kinetic aspects of crystalline (c-KTP) and amorphous (a-KTP) ketoprofen dissolution in water have been investigated by molecular dynamics simulation focusing on free energy properties. Absolute free energies of all relevant species and phases have been determined by thermodynamic integration on a novel path, first connecting the harmonic to the anharmonic system Hamiltonian at low T and then extending the result to the temperature of interest. The free energy required to transfer one ketoprofen molecule from the crystal to the solution is in fair agreement with the experimental value. The absolute free energy of the amorphous form is 19.58 kJ/mol higher than for the crystal, greatly enhancing the ketoprofen concentration in water, although as a metastable species in supersaturated solution. The kinetics of the dissolution process has been analyzed by computing the free energy profile along a reaction coordinate bringing one ketoprofen molecule from the crystal or amorphous phase to the solvated state. This computation confirms that, compared to the crystal form, the dissolution rate is nearly 7 orders of magnitude faster for the amorphous form, providing one further advantage to the latter in terms of bioavailability. The problem of drug solubility, of great practical importance, is used here as a test bed for a refined method to compute absolute free energies, which could be of great interest in biophysics and drug discovery in particular.

Differential Local Stability Governs the Metamorphic Fold Switch of Bacterial Virulence Factor RfaH

RfaH, a two-domain protein from a universally conserved NusG/Spt5 family of regulators, is required for the transcription and translation of long virulence and conjugation operons in many Gram-negative bacterial pathogens. Escherichia coli RfaH action is controlled by a unique large-scale structural rearrangement triggered by recruitment to transcription elongation complexes through a specific DNA element. Upon recruitment, the C-terminal domain of RfaH refolds from an α-hairpin, which is bound to RNA polymerase binding site within the N-terminal domain, into an unbound β-barrel that interacts with the ribosome. Although structures of the autoinhibited (α-hairpin) and active (β-barrel) states and plausible refolding pathways have been reported, how this reversible switch is encoded within RfaH sequence and structure is poorly understood. Here, we combined hydrogen-deuterium exchange measurements by mass spectrometry and nuclear magnetic resonance with molecular dynamics to evaluate the differential local stability between both RfaH folds. Deuteron incorporation reveals that the tip of the C-terminal hairpin (residues 125-145) is stably folded in the autoinhibited state (∼20% deuteron incorporation), whereas the rest of this domain is highly flexible (>40% deuteron incorporation), and its flexibility only decreases in the β-folded state. Computationally predicted ΔG agree with these results by displaying similar anisotropic stability within the tip of the α-hairpin and on neighboring N-terminal domain residues. Remarkably, the β-folded state shows comparable structural flexibility than nonmetamorphic homologs. Our findings provide information critical for understanding the metamorphic behavior of RfaH and other chameleon proteins and for devising targeted strategies to combat bacterial infections.

Free-Energy Simulations Resolve the Low-Affinity Na+-High-Affinity Asp Binding Paradox in GltPh

Glutamate transporters clear up excess extracellular glutamate by cotransporting three Na⁺ and one H⁺ with the countertransport of one K⁺. The archaeal homologs are selective to aspartate and only cotransport three Na⁺. The crystal structures of Glt_Ph from archaea have been used in computational studies to understand the transport mechanism. Although some progress has been made with regard to the ligand-binding sites, a consistent picture of transport still eludes us. A major concern is the discrepancy between the computed binding free energies, which predict high-affinity Na⁺-low-affinity aspartate binding, and the experimental results in which the opposite is observed. Here, we show that the binding of the first two Na⁺ ions involves an intermediate state near the Na1 site, where two Na⁺ ions coexist and couple to aspartate with similar strengths, boosting its affinity. Binding free energies for Na⁺ and aspartate obtained using this intermediate state are in good agreement with the experimental values. Thus, the paradox in binding affinities arises from the assumption that the ligands bind to the sites observed in the crystal structure following the order dictated by their binding free energies with no intermediate states. In fact, the presence of an intermediate state eliminates such a correlation between the binding free energies and the binding order. The intermediate state also facilitates transition of the first Na⁺ ion to its final binding site via a knock-on mechanism, which induces substantial conformational changes in the protein consistent with experimental observations.

Allosterism and signal transfer in DNA

We analysed the basic mechanisms of signal transmission in DNA and the origins of the allostery exhibited by systems such as the ternary complex BAMHI-DNA-GRDBD. We found that perturbation information generated by a primary protein binding event travels as a wave to distant regions of DNA following a hopping mechanism. However, such a structural perturbation is transient and does not lead to permanent changes in the DNA geometry and interaction properties at the secondary binding site. The BAMHI-DNA-GRDBD allosteric mechanism does not occur through any traditional models: direct (protein-protein), indirect (reorganization of the secondary site) readout or solvent-release. On the contrary, it is generated by a subtle and less common entropy-mediated mechanism, which might have an important role to explain other DNA-mediated cooperative effects.

Harnessing complexity in molecular self-assembly using computer simulations

In molecular self-assembly, hundreds of thousands of freely-diffusing molecules associate to form ordered and functional architectures in the absence of an actuator. This intriguing phenomenon plays a critical role in biology and has become a powerful tool for the fabrication of advanced nanomaterials. Due to the limited spatial and temporal resolutions of current experimental techniques, computer simulations offer a complementary strategy to explore self-assembly with atomic resolution. Here, we review recent computational studies focusing on both thermodynamic and kinetic aspects. As we shall see, thermodynamic approaches based on modeling and statistical mechanics offer initial guidelines to design nanostructures with modest computational effort. Computationally more intensive analyses based on molecular dynamics simulations and kinetic network models (KNMs) reach beyond it, opening the door to the rational design of self-assembly pathways. Current limitations of these methodologies are discussed. We anticipate that the synergistic use of thermodynamic and kinetic analyses based on computer simulations will provide an important contribution to the de novo design of self-assembly.

Modeling the adsorption equilibrium of small-molecule gases on graphene: effect of the volume to surface ratio

The adsorption probability of molecules at surfaces depends on the available volume-to-surface ratio with important implications for 2D self-assembly.

Computation of Hemagglutinin Free Energy Difference by the Confinement Method

Hemagglutinin (HA) mediates membrane fusion, a crucial step during influenza virus cell entry. How many HAs are needed for this process is still subject to debate. To aid in this discussion, the confinement free energy method was used to calculate the conformational free energy difference between the extended intermediate and postfusion state of HA. Special care was taken to comply with the general guidelines for free energy calculations, thereby obtaining convergence and demonstrating reliability of the results. The energy that one HA trimer contributes to fusion was found to be 34.2 ± 3.4kBT, similar to the known contributions from other fusion proteins. Although computationally expensive, the technique used is a promising tool for the further energetic characterization of fusion protein mechanisms. Knowledge of the energetic contributions per protein, and of conserved residues that are crucial for fusion, aids in the development of fusion inhibitors for antiviral drugs.

Binding free energies in the SAMPL5 octa-acid host-guest challenge calculated with DFT-D3 and CCSD(T)

We have tried to calculate the free energy for the binding of six small ligands to two variants of the octa-acid deep cavitand host in the SAMPL5 blind challenge. We employed structures minimised with dispersion-corrected density-functional theory with small basis sets and energies were calculated using large basis sets. Solvation energies were calculated with continuum methods and thermostatistical corrections were obtained from frequencies calculated at the HF-3c level. Care was taken to minimise the effects of the flexibility of the host by keeping the complexes as symmetric and similar as possible. In some calculations, the large net charge of the host was reduced by removing the propionate and benzoate groups. In addition, the effect of a restricted molecular dynamics sampling of structures was tested. Finally, we tried to improve the energies by using the DLPNO-CCSD(T) approach. Unfortunately, results of quite poor quality were obtained, with no correlation to the experimental data, systematically too positive affinities (by ~50 kJ/mol) and a mean absolute error (after removal of the systematic error) of 11-16 kJ/mol. DLPNO-CCSD(T) did not improve the results, so the accuracy is not limited by the energy function. Instead, four likely sources of errors were identified: first, the minimised structures were often incorrect, owing to the omission of explicit solvent. They could be partly improved by performing the minimisations in a continuum solvent with four water molecules around the charged groups of the ligands. Second, some ligands could bind in several different conformations, requiring sampling of reasonable structures. Third, there is an indication the continuum-solvation model has problems to accurately describe the binding of both the negatively and positively charged guest molecules. Fourth, different methods to calculate the thermostatistical corrections gave results that differed by up to 30 kJ/mol and there is an indication that HF-3c overestimates the entropy term. In conclusion, it is a challenge to calculate binding affinities for this octa-acid system with quantum-mechanical methods.

Balancing Accuracy and Cost of Confinement Simulations by Interpolation and Extrapolation of Confinement Energies

Improvements to the confinement method for the calculation of conformational free energy differences are presented. By taking advantage of phase space overlap between simulations at different frequencies, significant gains in accuracy and speed are reached. The optimal frequency spacing for the simulations is obtained from extrapolations of the confinement energy, and relaxation time analysis is used to determine time steps, simulation lengths, and friction coefficients. At postprocessing, interpolation of confinement energies is used to significantly reduce discretization errors in the calculation of conformational free energies. The efficiency of this protocol is illustrated by applications to alanine n-peptides and lactoferricin. For the alanine-n-peptide, errors were reduced between 2- and 10-fold and sampling times between 8- and 67-fold, while for lactoferricin the long sampling times at low frequencies were reduced 10-100-fold.

A Simple and Accurate Method To Calculate Free Energy Profiles and Reaction Rates from Restrained Molecular Simulations of Diffusive Processes

A method is developed to obtain simultaneously free energy profiles and diffusion constants from restrained molecular simulations in diffusive systems. The method is based on low-order expansions of the free energy and diffusivity as functions of the reaction coordinate. These expansions lead to simple analytical relationships between simulation statistics and model parameters. The method is tested on 1D and 2D model systems; its accuracy is found to be comparable to or better than that of the existing alternatives, which are briefly discussed. An important aspect of the method is that the free energy is constructed by integrating its derivatives, which can be computed without need for overlapping sampling windows. The implementation of the method in any molecular simulation program that supports external umbrella potentials (e.g., CHARMM) requires modification of only a few lines of code. As a demonstration of its applicability to realistic biomolecular systems, the method is applied to model the α-helix ↔ β-sheet transition in a 16-residue peptide in implicit solvent, with the reaction coordinate provided by the string method. Possible modifications of the method are briefly discussed; they include generalization to multidimensional reaction coordinates [in the spirit of the model of Ermak and McCammon (Ermak, D. L.; McCammon, J. A. J. Chem. Phys. 1978, 69, 1352-1360)], a higher-order expansion of the free energy surface, applicability in nonequilibrium systems, and a simple test for Markovianity. In view of the small overhead of the method relative to standard umbrella sampling, we suggest its routine application in the cases where umbrella potential simulations are appropriate.

Calculation of the Local Free Energy Landscape in the Restricted Region by the Modified Tomographic Method

The free energy landscape is the most important information in the study of the reaction mechanisms of the molecules. However, it is difficult to calculate. In a large collective variable space, a molecule must take a long time to obtain the sufficient sampling during the simulation. To save the calculation quantity, decreasing the sampling region and constructing the local free energy landscape is required in practice. However, the restricted region in the collective variable space may have an irregular shape. Simply restricting one or more collective variables of the molecule cannot satisfy the requirement. In this paper, we propose a modified tomographic method to perform the simulation. First, it divides the restricted region by some hyperplanes and connects the centers of hyperplanes together by a curve. Second, it forces the molecule to sample on the curve and the hyperplanes in the simulation and calculates the free energy data on them. Finally, all the free energy data are combined together to form the local free energy landscape. Without consideration of the area outside the restricted region, this free energy calculation can be more efficient. By this method, one can further optimize the path quickly in the collective variable space.

Assessment of Mutational Effects on Peptide Stability through Confinement Simulations

The evaluation of free energy differences between specific states of a system is of fundamental interest in the study of (bio)chemical systems. Herein, we examine the use of the recently introduced confinement method (CM) to evaluate relative free energy changes upon protein/peptide mutations. CM is a path-independent technique that involves the transformation of a configurational state of the system into an ideal crystal permitting the direct computation of free energy differences. We illustrate the method by evaluating the differential stabilities between native and mutant sequences of a model peptide that has been extensively characterized by experimental approaches, the GB1 hairpin. We show a good correlation between calculated and experimental relative stabilities and discuss other possible applications of this method in the context of complex molecular conversions.

Predicting molecular self-assembly at surfaces: a statistical thermodynamics and modeling approach

A self-consistent framework based on modeling and statistical mechanics for the theoretical interpretation of self-assembly at surfaces and interfaces is presented.

Constraint methods that accelerate free-energy simulations of biomolecules

Atomistic molecular dynamics simulations of biomolecules are critical for generating narratives about biological mechanisms. The power of atomistic simulations is that these are physics-based methods that satisfy Boltzmann's law, so they can be used to compute populations, dynamics, and mechanisms. But physical simulations are computationally intensive and do not scale well to the sizes of many important biomolecules. One way to speed up physical simulations is by coarse-graining the potential function. Another way is to harness structural knowledge, often by imposing spring-like restraints. But harnessing external knowledge in physical simulations is problematic because knowledge, data, or hunches have errors, noise, and combinatoric uncertainties. Here, we review recent principled methods for imposing restraints to speed up physics-based molecular simulations that promise to scale to larger biomolecules and motions.

Investigations of α-helix↔β-sheet transition pathways in a miniprotein using the finite-temperature string method

A parallel implementation of the finite-temperature string method is described, which takes into account the invariance of coordinates with respect to rigid-body motions. The method is applied to the complex α-helix↔β-sheet transition in a β-hairpin miniprotein in implicit solvent, which exhibits much of the complexity of conformational changes in proteins. Two transition paths are considered, one derived from a linear interpolant between the endpoint structures and the other derived from a targeted dynamics simulation. Two methods for computing the conformational free energy (FE) along the string are compared, a restrained method, and a tessellation method introduced by E. Vanden-Eijnden and M. Venturoli [J. Chem. Phys. 130, 194103 (2009)]. It is found that obtaining meaningful free energy profiles using the present atom-based coordinates requires restricting sampling to a vicinity of the converged path, where the hyperplanar approximation to the isocommittor surface is sufficiently accurate. This sampling restriction can be easily achieved using restraints or constraints. The endpoint FE differences computed from the FE profiles are validated by comparison with previous calculations using a path-independent confinement method. The FE profiles are decomposed into the enthalpic and entropic contributions, and it is shown that the entropy difference contribution can be as large as 10 kcal/mol for intermediate regions along the path, compared to 15-20 kcal/mol for the enthalpy contribution. This result demonstrates that enthalpic barriers for transitions are offset by entropic contributions arising from the existence of different paths across a barrier. The possibility of using systematically coarse-grained representations of amino acids, in the spirit of multiple interaction site residue models, is proposed as a means to avoid ad hoc sampling restrictions to narrow transition tubes.

Computing the relative stabilities and the per-residue components in protein conformational changes

Protein molecules often undergo conformational changes. In order to gain insights into the forces that drive such changes, it would be useful to have a method that computes the per-residue contributions to the conversion free energy. Here, we describe the "confine-convert-release" (CCR) method, which is applicable to large conformational changes. We show that CCR correctly predicts the stable states of several "chameleon" sequences that have previously been challenging for molecular simulations. CCR can often discriminate better from worse predictions of native protein models in critical assessment of protein structure prediction (CASP). We show how the total conversion free energies can be parsed into per-residue free-energy components. Such parsing gives insights into which amino acids are most responsible for given transformations. For example, here we are able to "reverse-engineer" the known design principles of the chameleon proteins. This opens up the possibility for systematic improvements in structure-prediction scoring functions, in the design of protein conformational switches, and in interpreting protein mechanisms at the amino-acid level.

String method for calculation of minimum free-energy paths in Cartesian space in freely-tumbling systems

The string method is a molecular-simulation technique that aims to calculate the minimum free-energy path of a chemical reaction or conformational transition, in the space of a pre-defined set of reaction coordinates that is typically highly dimensional. Any descriptor may be used as a reaction coordinate, but arguably the Cartesian coordinates of the atoms involved are the most unprejudiced and intuitive choice. Cartesian coordinates, however, present a non-trivial problem, in that they are not invariant to rigid-body molecular rotations and translations, which ideally ought to be unrestricted in the simulations. To overcome this difficulty, we reformulate the framework of the string method to integrate an on-the-fly structural-alignment algorithm. This approach, referred to as SOMA (String method with Optimal Molecular Alignment), enables the use of Cartesian reaction coordinates in freely tumbling molecular systems. In addition, this scheme permits the dissection of the free-energy change along the most probable path into individual atomic contributions, thus revealing the dominant mechanism of the simulated process. This detailed analysis also provides a physically-meaningful criterion to coarse-grain the representation of the path. To demonstrate the accuracy of the method we analyze the isomerization of the alanine dipeptide in vacuum and the chair-to-inverted-chair transition of β-D mannose in explicit water. Notwithstanding the simplicity of these systems, the SOMA approach reveals novel insights into the atomic mechanism of these isomerizations. In both cases, we find that the dynamics and the energetics of these processes are controlled by interactions involving only a handful of atoms in each molecule. Consistent with this result, we show that a coarse-grained SOMA calculation defined in terms of these subsets of atoms yields nearidentical minimum free-energy paths and committor distributions to those obtained via a highly-dimensional string.

Conformational contribution to thermodynamics of binding in protein-peptide complexes through microscopic simulation

We extract the thermodynamics of conformational changes in biomacromolecular complexes from the distributions of the dihedral angles of the macromolecules. These distributions are obtained from the equilibrium configurations generated via all-atom molecular dynamics simulations. The conformational thermodynamics data we obtained for calmodulin-peptide complexes using our methodology corroborate well with the experimentally observed conformational and binding entropies. The conformational free-energy changes and their contributions for different peptide-binding regions of calmodulin are evaluated microscopically.

A simplified confinement method for calculating absolute free energies and free energy and entropy differences

A simple and robust formulation of the path-independent confinement method for the calculation of free energies is presented. The simplified confinement method (SCM) does not require matrix diagonalization or switching off the molecular force field, and has a simple convergence criterion. The method can be readily implemented in molecular dynamics programs with minimal or no code modifications. Because the confinement method is a special case of thermodynamic integration, it is trivially parallel over the integration variable. The accuracy of the method is demonstrated using a model diatomic molecule, for which exact results can be computed analytically. The method is then applied to the alanine dipeptide in vacuum, and to the α-helix ↔ β-sheet transition in a 16-residue peptide modeled in implicit solvent. The SCM requires less effort for the calculation of free energy differences than previous formulations because it does not require computing normal modes. The SCM has a diminished advantage for determining absolute free energy values, because it requires decreasing the MD integration step to obtain accurate results. An approximate confinement procedure is introduced, which can be used to estimate directly the configurational entropy difference between two macrostates, without the need for additional computation of the difference in the free energy or enthalpy. The approximation has convergence properties similar to those of the standard confinement method for the calculation of free energies. The use of the approximation requires about 5 times less wall-clock simulation time than that needed to compute enthalpy differences to similar precision from an MD trajectory. For the biomolecular systems considered in this study, the errors in the entropy approximation are under 10%. Practical applications of the methods to proteins are currently limited to implicit solvent simulations.

Free energy simulations of ligand binding to the aspartate transporter Glt(Ph)

Glutamate/Aspartate transporters cotransport three Na(+) and one H(+) ions with the substrate and countertransport one K(+) ion. The binding sites for the substrate and two Na(+) ions have been observed in the crystal structure of the archeal homolog Glt(Ph), while the binding site for the third Na(+) ion has been proposed from computational studies and confirmed by experiments. Here we perform detailed free energy simulations of Glt(Ph), giving a comprehensive characterization of the substrate and ion binding sites, and calculating their binding free energies in various configurations. Our results show unequivocally that the substrate binds after the binding of two Na(+) ions. They also shed light into Asp/Glu selectivity of Glt(Ph), which is not observed in eukaryotic glutamate transporters.

A conformational transition in the myosin VI converter contributes to the variable step size

Myosin VI (MVI) is a dimeric molecular motor that translocates backwards on actin filaments with a surprisingly large and variable step size, given its short lever arm. A recent x-ray structure of MVI indicates that the large step size can be explained in part by a novel conformation of the converter subdomain in the prepowerstroke state, in which a 53-residue insert, unique to MVI, reorients the lever arm nearly parallel to the actin filament. To determine whether the existence of the novel converter conformation could contribute to the step-size variability, we used a path-based free-energy simulation tool, the string method, to show that there is a small free-energy difference between the novel converter conformation and the conventional conformation found in other myosins. This result suggests that MVI can bind to actin with the converter in either conformation. Models of MVI/MV chimeric dimers show that the variability in the tilting angle of the lever arm that results from the two converter conformations can lead to step-size variations of ∼12 nm. These variations, in combination with other proposed mechanisms, could explain the experimentally determined step-size variability of ∼25 nm for wild-type MVI. Mutations to test the findings by experiment are suggested.

Accurate determination of the binding free energy for KcsA-charybdotoxin complex from the potential of mean force calculations with restraints

Free energy calculations for protein-ligand dissociation have been tested and validated for small ligands (50 atoms or less), but there has been a paucity of studies for larger, peptide-size ligands due to computational limitations. Previously we have studied the energetics of dissociation in a potassium channel-charybdotoxin complex by using umbrella sampling molecular-dynamics simulations, and established the need for carefully chosen coordinates and restraints to maintain the physiological ligand conformation. Here we address the ligand integrity problem further by constructing additional potential of mean forces for dissociation of charybdotoxin using restraints. We show that the large discrepancies in binding free energy arising from simulation artifacts can be avoided by using appropriate restraints on the ligand, which enables determination of the binding free energy within the chemical accuracy. We make several suggestions for optimal choices of harmonic potential parameters and restraints to be used in binding studies of large ligands.

Modeling intrinsically disordered proteins with bayesian statistics

The characterization of intrinsically disordered proteins is challenging because accurate models of these systems require a description of both their thermally accessible conformers and the associated relative stabilities or weights. These structures and weights are typically chosen such that calculated ensemble averages agree with some set of prespecified experimental measurements; however, the large number of degrees of freedom in these systems typically leads to multiple conformational ensembles that are degenerate with respect to any given set of experimental observables. In this work we demonstrate that estimates of the relative stabilities of conformers within an ensemble are often incorrect when one does not account for the underlying uncertainty in the estimates themselves. Therefore, we present a method for modeling the conformational properties of disordered proteins that estimates the uncertainty in the weights of each conformer. The Bayesian weighting (BW) formalism incorporates information from both experimental data and theoretical predictions to calculate a probability density over all possible ways of weighting the conformers in the ensemble. This probability density is then used to estimate the values of the weights. A unique and powerful feature of the approach is that it provides a built-in error measure that allows one to assess the accuracy of the ensemble. We validate the approach using reference ensembles constructed from the five-residue peptide met-enkephalin and then apply the BW method to construct an ensemble of the K18 isoform of the tau protein. Using this ensemble, we indentify a specific pattern of long-range contacts in K18 that correlates with the known aggregation properties of the sequence.

Accounting for ligand conformational restriction in calculations of protein-ligand binding affinities

The conformation adopted by a ligand on binding to a receptor may differ from its lowest-energy conformation in solution. In addition, the bound ligand is more conformationally restricted, which is associated with a configurational entropy loss. The free energy change due to these effects is often neglected or treated crudely in current models for predicting binding affinity. We present a method for estimating this contribution, based on perturbation theory using the quasi-harmonic model of Karplus and Kushick as a reference system. The consistency of the method is checked for small model systems. Subsequently we use the method, along with an estimate for the enthalpic contribution due to ligand-receptor interactions, to calculate relative binding affinities. The AMBER force field and generalized Born implicit solvent model is used. Binding affinities were estimated for a test set of 233 protein-ligand complexes for which crystal structures and measured binding affinities are available. In most cases, the ligand conformation in the bound state was significantly different from the most favorable conformation in solution. In general, the correlation between measured and calculated ligand binding affinities including the free energy change due to ligand conformational change is comparable to or slightly better than that obtained by using an empirically-trained docking score. Both entropic and enthalpic contributions to this free energy change are significant.

Energy landscape of a hydrogen-bonded non-degenerate molecular shuttle

In rotaxane 1, two co-conformations are populated in CDCl(3) at temperatures between 250 and 330 K. The thermodynamic parameters show strong enthalpy-entropy compensation, and a non-negligible heat capacity difference between the two forms.