Escape of a Small Molecule from Inside T4 Lysozyme by Multiple Pathways

The T4 lysozyme L99A mutant is often used as a model system to study small-molecule binding to proteins, but pathways for ligand entry and exit from the buried binding site and the associated protein conformational changes have not been fully resolved. Here, molecular dynamics simulations were employed to model benzene exit from its binding cavity using the weighted ensemble (WE) approach to enhance sampling of low-probability unbinding trajectories. Independent WE simulations revealed four pathways for benzene exit, which correspond to transient tunnels spontaneously formed in previous simulations of apo T4 lysozyme. Thus, benzene unbinding occurs through multiple pathways partially created by intrinsic protein structural fluctuations. Motions of several α-helices and side chains were involved in ligand escape from metastable microstates. WE simulations also provided preliminary estimates of rate constants for each exit pathway. These results complement previous works and provide a semiquantitative characterization of pathway heterogeneity for binding of small molecules to proteins.     

Transmembrane helix orientation and dynamics: insights from ensemble dynamics with solid-state NMR observables

As the major component of membrane proteins, transmembrane helices embedded in anisotropic bilayer environments adopt preferential orientations that are characteristic or related to their functional states. Recent developments in solid-state nuclear magnetic resonance (SSNMR) spectroscopy have made it possible to measure NMR observables that can be used to determine such orientations in a native bilayer environment. A quasistatic single conformer model is frequently used to interpret the SSNMR observables, but important motional information can be missing or misinterpreted in the model. In this work, we have investigated the orientation of the single-pass transmembrane domain of viral protein ”u“ (VpuTM) from HIV-1 by determining an ensemble of structures using multiple conformer models based on the SSNMR ensemble dynamics technique. The resulting structure ensemble shows significantly larger orientational fluctuations while the ensemble-averaged orientation is compatible with the orientation based on the quasistatic model. This observation is further corroborated by comparison with the VpuTM orientation from comparative molecular dynamics simulations in explicit bilayer membranes. SSNMR ensemble dynamics not only reveals the importance of transmembrane helix dynamics in interpretation of SSNMR observables, but also provides a means to simultaneously extract both transmembrane helix orientation and dynamics information from the SSNMR measurements.     

Approximate quantum trajectory dynamics for reactive processes in condensed phase

A method of molecular dynamics with quantum corrections, practical for studies of large molecular systems, is reviewed. The approach is based on the Bohmian formulation of the time-dependent Schrödinger equation in which a wavefunction is represented by an ensemble of interdependent trajectories. The quantum effects come from the quantum potential acting on trajectories on par with the usual classical potential. The quantum potential is determined from the evolving nuclear wavefunction, i.e. from the quantum trajectory (QT) ensemble itself. For practical and conceptual reasons the quantum potential and corresponding quantum nuclear effect are computed only for the selected light nuclei. For studies of reactive chemical processes, the classical potential is computed on-the-fly using the density functional tight binding method of electronic structure. A massively parallel implementation, based on the message passing interface allows for efficient simulations of ensembles of thousands of trajectories describing systems of up to 200 atoms. As a biochemical application, the approximate QT approach is used to model the tunnelling-dominated proton transfer in soybean-lipoxygenase-1. A materials science application is represented by a study of the nuclear quantum effect on adsorption of hydrogen and deuterium on a C₃₇H₁₅ molecule, which is a model ‘flake’ of graphene.     

Classical molecular interaction potentials: improved setup procedure in molecular dynamics simulations of proteins.

The latest version of the classical molecular interaction potential (CMIP) has the ability to predict the position of crystallographic waters in several proteins with great accuracy. This article analyzes the ability of the CMIP functional to improve the setup procedure of the molecular system in molecular dynamics (MD) simulations of proteins. To this end, the CMIP strategy is used to include both water molecules and counterions in different protein systems. The structural details of the configurations sampled from trajectories obtained using the CMIP setup procedure are compared with those obtained from trajectories derived from a standard equilibration process. The results show that standard MD simulations can lead to artifactual results, which are avoided using the CMIP setup procedure. Because the CMIP is easy to implement at a low computational cost, it can be very useful in obtaining reliable MD trajectories.     

The MUMO (minimal under-restraining minimal over-restraining) method for the determination of native state ensembles of proteins

While reliable procedures for determining the conformations of proteins are available, methods for generating ensembles of structures that also reflect their flexibility are much less well established. Here we present a systematic assessment of the ability of ensemble-averaged molecular dynamics simulations with ensemble-averaged NMR restraints to simultaneously reproduce the average structure of proteins and their associated dynamics. We discuss the effects that under-restraining (overfitting) and over-restraining (underfitting) have on the structures generated in ensemble-averaged molecular simulations. We then introduce the MUMO (minimal under-restraining minimal over-restraining) method, a procedure in which different observables are averaged over a different number of molecules. As both over-restraining and under-restraining are significantly reduced in the MUMO method, it is possible to generate ensembles of conformations that accurately characterize both the structure and the dynamics of native states of proteins. The application of the MUMO method to the protein ubiquitin yields a high-resolution structural ensemble with an RDC Q-factor of 0.19.     

Solid-State NMR-Restrained Ensemble Dynamics of a Membrane Protein in Explicit Membranes

Solid-state NMR has been used to determine the structures of membrane proteins in native-like lipid bilayer environments. Most structure calculations based on solid-state NMR observables are performed using simulated annealing with restrained molecular dynamics and an energy function, where all nonbonded interactions are represented by a single, purely repulsive term with no contributions from van der Waals attractive, electrostatic, or solvation energy. To our knowledge, this is the first application of an ensemble dynamics technique performed in explicit membranes that uses experimental solid-state NMR observables to obtain the refined structure of a membrane protein together with information about its dynamics and its interactions with lipids. Using the membrane-bound form of the fd coat protein as a model membrane protein and its experimental solid-state NMR data, we performed restrained ensemble dynamics simulations with different ensemble sizes in explicit membranes. For comparison, a molecular dynamics simulation of fd coat protein was also performed without any restraints. The average orientation of each protein helix is similar to a structure determined by traditional single-conformer approaches. However, their variations are limited in the resulting ensemble of structures with one or two replicas, as they are under the strong influence of solid-state NMR restraints. Although highly consistent with all solid-state NMR observables, the ensembles of more than two replicas show larger orientational variations similar to those observed in the molecular dynamics simulation without restraints. In particular, in these explicit membrane simulations, Lys⁴⁰, residing at the C-terminal side of the transmembrane helix, is observed to cause local membrane curvature. Therefore, compared to traditional single-conformer approaches in implicit environments, solid-state NMR restrained ensemble simulations in explicit membranes readily characterize not only protein dynamics but also protein-lipid interactions in detail.     

Ensemble MD simulations restrained via crystallographic data: Accurate structure leads to accurate dynamics

Currently, the best existing molecular dynamics (MD) force fields cannot accurately reproduce the global free‐energy minimum which realizes the experimental protein structure. As a result, long MD trajectories tend to drift away from the starting coordinates (e.g., crystallographic structures). To address this problem, we have devised a new simulation strategy aimed at protein crystals. An MD simulation of protein crystal is essentially an ensemble simulation involving multiple protein molecules in a crystal unit cell (or a block of unit cells). To ensure that average protein coordinates remain correct during the simulation, we introduced crystallography‐based restraints into the MD protocol. Because these restraints are aimed at the ensemble‐average structure, they have only minimal impact on conformational dynamics of the individual protein molecules. So long as the average structure remains reasonable, the proteins move in a native‐like fashion as dictated by the original force field. To validate this approach, we have used the data from solid‐state NMR spectroscopy, which is the orthogonal experimental technique uniquely sensitive to protein local dynamics. The new method has been tested on the well‐established model protein, ubiquitin. The ensemble‐restrained MD simulations produced lower crystallographic R factors than conventional simulations; they also led to more accurate predictions for crystallographic temperature factors, solid‐state chemical shifts, and backbone order parameters. The predictions for ¹⁵N relaxation rates are at least as accurate as those obtained from conventional simulations. Taken together, these results suggest that the presented trajectories may be among the most realistic protein MD simulations ever reported. In this context, the ensemble restraints based on high‐resolution crystallographic data can be viewed as protein‐specific empirical corrections to the standard force fields.     

Learning chemistry with multiple first-principles simulations

Huge parallel high-performance computing (HPC) architectures are today available laboratories for modelling atomic forces with high accuracy and for large samples of atoms. Modern statistical tools allow to simulate the statistics of these samples, while first-principles molecular dynamics (MD) can probe the interactions within large atomic samples, including possible chemical reactions. But a proper statistical convergence for the ensemble, represented in terms of a bundle of trajectories, is still unsatisfactory in terms of comparisons with experiments. Can we learn something by these HPC experiments? In this contribution, we show one example, where the occurrence of a chemical reaction in a disordered system is probed. The complex of the copper ion and a segment of the amyloid-β peptide, of wide interest in understanding the progress of Alzheimer's disease, was modelled combining constructions based on empirical force fields with first-principles MD simulations. We simulate a bundle of 16 different structures, biasing different Cu coordination numbers and changing the charge (oxidation state) of the assembly. Even within the given approximations for forces and the poor sampling, we could identify the structures of the complex that are able to react with hydrogen peroxide. The observation explains, at a molecular level, one important linkage between Alzheimer's disease and oxidative stress. This is an example of a general strategy for exploiting reactive configurations within a large set of possible reasonable candidates.     

Lipid tempering simulation of model biological membranes on parallel platforms

In this report we have tested a parallel implementation for the simulation of lipid bilayers at the atomistic level, based on a generalized ensemble protocol where only the torsional degrees of freedom of the alkyl chains of the lipids are heated. The results in terms of configurational sampling enhancement have been compared with a conventional simulation produced with a widespread molecular dynamics code. Results show that the proposed thermodynamic-based multiple trajectories parallel protocol for membrane simulations allows for an efficient use of CPU resources with respect to the conventional single trajectory, providing accurate results for area and volume per lipid, membrane thickness, undulation spectra and boosting significantly diffusion and mixing in lipid bilayers due to the sampling enhancement of gauche/trans ratios of the alkyl chain dihedral angles.     

Blue Moon Approach to Rare Events

The "Blue Moon" ensemble is a computationally efficient molecular dynamics method to estimate the rate constants of rare activated events when the process can be described by a reaction coordinate ξ(r), a well-defined function in configuration space. By means of holonomic constraints a number of values of ξ(r) can be prescribed along the relevant path to identify the "bottleneck" region first and to sample an ensemble of starting conditions to generate activated trajectories. These MD trajectories sample phase space according to a biased configurational distribution. With a suitable re-weighting of averages from such ensemble of trajectories one can characterize completely rare events.     

The Folding Transition-State Ensemble of a Four-Helix Bundle Protein: Helix Propensity as a Determinant and Macromolecular Crowding as a Probe

The four-helix bundle protein Rd-apocyt b₅₆₂, a redesigned stable variant of apocytochrome b₅₆₂, exhibits two-state folding kinetics. Its transition-state ensemble has been characterized by Φ-value analysis. To elucidate the molecular basis of the transition-state ensemble, we have carried out high-temperature molecular dynamics simulations of the unfolding process. In six parallel simulations, unfolding started with the melting of helix I and the C-terminal half of helix IV, and followed by helix III, the N-terminal half of helix IV and helix II. This ordered melting of the helices is consistent with the conclusion from native-state hydrogen exchange, and can be rationalized by differences in intrinsic helix propensity. Guided by experimental Φ-values, a putative transition-state ensemble was extracted from the simulations. The residue helical probabilities of this transition-state ensemble show good correlation with the Φ-values. To further validate the putative transition-state ensemble, the effect of macromolecular crowding on the relative stability between the unfolded ensemble and the transition-state ensemble was calculated. The resulting effect of crowding on the folding kinetics agrees well with experimental observations. This study shows that molecular dynamics simulations combined with calculation of crowding effects provide an avenue for characterize the transition-state ensemble in atomic details.     

GRID-MD-A tool for massive simulation of protein channels

We present here a fast method for the exploration of channels in proteins based on molecular dynamics simulations of probe particles in a discrete grid space defined by an ensemble of protein conformations obtained either experimentally or by out-of-the-grid atomistic molecular dynamics simulations. The method is able to provide millisecond-long trajectories with a small computational effort, requires no human intervention in defining possible exit pathways and can detect both major and minor channels, giving a correct balance to the relative flux between them. The Grid-Molecular-Dynamics approach is then a suitable method for massive exploration of channels in proteins, even of those with unknown functional annotation.     

Efficient sampling in collective coordinate space

Collective motions in biological macromolecules have been shown to be important for function. The most important collective motions occur on slow time scales, which poses a sampling problem in dynamic simulation of biomolecules. We present a novel method for efficient conformational sampling. The method combines the simulation of an ensemble of concurrent trajectories with restraints acting on the ensemble of structures as a whole. Two properties of the ensemble may be restrained: (i) the variance of the ensemble and (ii) the average position of the ensemble. Both properties are defined in a subspace of collective coordinate space spanned by an arbitrary number of modes. We show that weak restraints on the ensemble variance suffice for an increase in sampling efficiency along soft modes by two orders of magnitudes. The resulting trajectories exhibit virtually the same structural quality as trajectories generated by restraint-free-molecular dynamics simulation, as judged by standard structure validation tools. The method is used to probe the resistance of a structure against conformational changes along collective modes and clearly distinguishes soft from stiff modes. Further applications are discussed. Proteins 2000;39:82-88.     

Highly sampled tetranucleotide and tetraloop motifs enable evaluation of common RNA force fields

Recent modifications and improvements to standard nucleic acid force fields have attempted to fix problems and issues that have been observed as longer timescale simulations have become routine. Although previous work has shown the ability to fold the UUCG stem–loop structure, until now no group has attempted to quantify the performance of current force fields using highly converged structural populations of the tetraloop conformational ensemble. In this study, we report the use of multiple independent sets of multidimensional replica exchange molecular dynamics (M-REMD) simulations with different initial conditions to generate well-converged conformational ensembles for the tetranucleotides r(GACC) and r(CCCC), as well as the larger UUCG tetraloop motif. By generating what is to our knowledge the most complete RNA structure ensembles reported to date for these systems, we remove the coupling between force field errors and errors due to incomplete sampling, providing a comprehensive comparison between current top-performing MD force fields for RNA. Of the RNA force fields tested in this study, none demonstrate the ability to correctly identify the most thermodynamically stable structure for all three systems. We discuss the deficiencies present in each potential function and suggest areas where improvements can be made. The results imply that although “short” (nsec-μsec timescale) simulations may stay close to their respective experimental structures and may well reproduce experimental observables, inevitably the current force fields will populate alternative incorrect structures that are more stable than those observed via experiment.     

Solid-state NMR ensemble dynamics as a mediator between experiment and simulation

Solid-state NMR (SSNMR) is a powerful technique to describe the orientations of membrane proteins and peptides in their native membrane bilayer environments. The deuterium (²H) quadrupolar splitting (DQS), one of the SSNMR observables, has been used to characterize the orientations of various single-pass transmembrane (TM) helices using a semistatic rigid-body model such as the geometric analysis of labeled alanine (GALA) method. However, dynamic information of these TM helices, which could be related to important biological function, can be missing or misinterpreted with the semistatic model. We have investigated the orientation of WALP23 in an implicit membrane of dimyristoylglycerophosphocholine by determining an ensemble of structures using multiple conformer models with a DQS restraint potential. When a single conformer is used, the resulting helix orientation (tilt angle (τ) of 5.6 ± 3.2° and rotation angle (ρ) of 141.8 ± 40.6°) is similar to that determined by the GALA method. However, as the number of conformers is increased, the tilt angles of WALP23 ensemble structures become larger (26.9 ± 6.7°), which agrees well with previous molecular dynamics simulation results. In addition, the ensemble structure distribution shows excellent agreement with the two-dimensional free energy surface as a function of WALP23's τ and ρ. These results demonstrate that SSNMR ensemble dynamics provides a means to extract orientational and dynamic information of TM helices from their SSNMR observables and to explain the discrepancy between molecular dynamics simulation and GALA-based interpretation of DQS data.     

Kinetic studies of folding of the B-domain of staphylococcal protein A with molecular dynamics and a united-residue (UNRES) model of polypeptide chains

Langevin dynamics is used with our physics-based united-residue (UNRES) force field to study the folding pathways of the B-domain of staphylococcal protein A (1BDD (alpha; 46 residues)). With 400 trajectories of protein A started from the extended state (to gather meaningful statistics), and simulated for more than 35 ns each, 380 of them folded to the native structure. The simulations were carried out at the optimal folding temperature of protein A with this force field. To the best of our knowledge, this is the first simulation study of protein-folding kinetics with a physics-based force field in which reliable statistics can be gathered. In all the simulations, the C-terminal alpha-helix forms first. The ensemble of the native basin has an average RMSD value of 4 A from the native structure. There is a stable intermediate along the folding pathway, in which the N-terminal alpha-helix is unfolded; this intermediate appears on the way to the native structure in less than one-fourth of the folding pathways, while the remaining ones proceed directly to the native state. Non-native structures persist until the end of the simulations, but the native-like structures dominate. To express the kinetics of protein A folding quantitatively, two observables were used: (i) the average alpha-helix content (averaged over all trajectories within a given time window); and (ii) the fraction of conformations (averaged over all trajectories within a given time window) with Calpha RMSD values from the native structure less than 5 A (fraction of completely folded structures). The alpha-helix content grows quickly with time, and its variation fits well to a single-exponential term, suggesting fast two-state kinetics. On the other hand, the fraction of folded structures changes more slowly with time and fits to a sum of two exponentials, in agreement with the appearance of the intermediate, found when analyzing the folding pathways. This observation demonstrates that different qualitative and quantitative conclusions about folding kinetics can be drawn depending on which observable is monitored.     

Formation of the folding nucleus of an SH3 domain investigated by loosely coupled molecular dynamics simulations

The experimentally well-established folding mechanism of the src-SH3 domain, and in particular the phi-value analysis of its transition state, represents a sort of testing table for computational investigations of protein folding. Here, parallel molecular dynamics simulations of the src-SH3 domain have been performed starting from denatured conformations. By rescuing and restarting only trajectories approaching the folding transition state, an ensemble of conformations was obtained with a completely structured central beta-sheet and a native-like packing of residues Ile-110, Ala-121, and Ile-132. An analysis of the trajectories shows that there are several pathways leading to the formation of the central beta-sheet whereas its two hairpins form in a different but consistent way.     

Thermodynamic and structural effect of urea and guanidine chloride on the helical and on a hairpin fragment of GB1 from molecular simulations

With the help of molecular‐dynamics simulations, we studied the effect of urea and guanidine chloride on the thermodynamic and structural properties of the helical fragment of protein GB1, comparing them with those of its second beta hairpin. We showed that the helical fragment in different solvents populates an ensemble of states that is more complex than that of the hairpin, and thus the associated experimental observables (circular‐dichroism spectra, secondary chemical shifts, m values), that we back‐calculated from the simulations and compared with the actual data, are more difficult to interpret. We observed that in the case of both peptides, urea binds tightly to their backbone, while guanidine exerts its denaturing effect in a more subtle way, strongly affecting the electrostatic properties of the solution. This difference can have consequences in the way denaturation experiments are interpreted. Proteins 2017; 85:753–763. © 2016 Wiley Periodicals, Inc.     

The folding transition state of the cold shock protein is strongly polarized

It has been shown recently that an 11-residue peptide fragment of transthyretin, TTR(105-115), can form amyloid fibrils in vitro by adopting an extended beta-strand conformation. We used molecular dynamics simulations on systems of TTR(105-115) peptides, for a total length of about 5 micros, to explore the process of self-assembly and the structures of the resulting aggregates. Our results suggest that an antiparallel association of the beta-strands is more probable than a parallel one and that the central residues (T106-L111) in a beta-strand have a high propensity to form inter-peptide hydrogen bonds. The study of the dynamics of self-association indicated that, for this peptide, trajectories leading to conformations with high alpha-helical content are off-pathway from those leading to aggregates with high beta-structure content. We also show that the diverse oligomeric structures that form spontaneously in the molecular dynamics simulations are, to a large extent, compatible with solid-state NMR experimental measurements, including chemical shifts, on fully formed fibrils. The strategy that we present may therefore be used in the design of new experiments to determine the structure of amyloid fibrils, such as those involving site-specific isotope labelling schemes to measure key inter-atomic distances.