Simulation of conformational transitions by the restricted perturbation-targeted molecular dynamics method

Molecular Simulations of Conformational Transitions within the Insulin Receptor Kinase Reveal Consensus Features in a Multistep Activation Pathway

Modulating the transitions between active and inactive conformations of protein kinases is the primary means of regulating their catalytic activity, achieved by phosphorylation of the activation loop (A-loop). To elucidate the mechanism of this conformational activation, we applied the string method to determine the conformational transition path of insulin receptor kinase between the active and inactive conformations and the corresponding free-energy profiles with and without A-loop phosphorylation. The conformational change was found to proceed in three sequential steps: first, the flipping of the DFG motif of the active site; second, rotation of the A-loop; finally, the inward movement of the αC helix. The main energetic bottleneck corresponds to the conformational change in the A-loop, while changes in the DFG motif and αC helix occur before and after A-loop conformational change, respectively. In accordance with this, two intermediate states are identified, the first state just after the DFG flipping and the second state after the A-loop rotation. These intermediates exhibit structural features characteristic of the corresponding inactive and active conformations of other protein kinases. To understand the impact of A-loop phosphorylation on kinase conformation, the free energies of A-loop phosphorylation were determined at several states along the conformational transition path using the free-energy perturbation simulations. The calculated free energies reveal that while the unphosphorylated kinase interconverts between the inactive and active conformations, A-loop phosphorylation restricts access to the inactive conformation, thereby increasing the active conformation population. Overall, this study suggests a consensus mechanism of conformational activation between different protein kinases.

Qualitative Prediction of Ligand Dissociation Kinetics from Focal Adhesion Kinase Using Steered Molecular Dynamics

Most early-stage drug discovery projects focus on equilibrium binding affinity to the target alongside selectivity and other pharmaceutical properties. Since many approved drugs have nonequilibrium binding characteristics, there has been increasing interest in optimizing binding kinetics early in the drug discovery process. As focal adhesion kinase (FAK) is an important drug target, we examine whether steered molecular dynamics (SMD) can be useful for identifying drug candidates with the desired drug-binding kinetics. In simulating the dissociation of 14 ligands from FAK, we find an empirical power–law relationship between the simulated time needed for ligand unbinding and the experimental rate constant for dissociation, with a strong correlation depending on the SMD force used. To improve predictions, we further develop regression models connecting experimental dissociation rate with various structural and energetic quantities derived from the simulations. These models can be used to predict dissociation rates from FAK for related compounds.

Exploring Conformational Space with Thermal Fluctuations Obtained by Normal-Mode Analysis

Proteins in their native states can be represented as ensembles of conformers in dynamical equilibrium. Thermal fluctuations are responsible for transitions between these conformers. Normal-modes analysis (NMA) using elastic network models (ENMs) provides an efficient procedure to explore global dynamics of proteins commonly associated with conformational transitions. In the present work, we present an iterative approach to explore protein conformational spaces by introducing structural distortions according to their equilibrium dynamics at room temperature. The approach can be used either to perform unbiased explorations of conformational space or to explore guided pathways connecting two different conformations, e.g., apo and holo forms. In order to test its performance, four proteins with different magnitudes of structural distortions upon ligand binding have been tested. In all cases, the conformational selection model has been confirmed and the conformational space between apo and holo forms has been encompassed. Different strategies have been tested that impact on the efficiency either to achieve a desired conformational change or to achieve a balanced exploration of the protein conformational multiplicity.

In silico Description of LAT1 Transport Mechanism at an Atomistic Level

The molecular mechanism of transport mediated by LAT1, a sodium-independent antiporter of large neutral amino acids, was investigated through in silico procedures, specifically making reference to two transported substrates, tyrosine (Tyr) and leucine methyl ester (LME), and to 3,5-diiodo-L-tyrosine (DIT), a well-known LAT1 inhibitor. Two models of the transporter were built by comparative modeling, with LAT1 either in an outward-facing (OF) or in an inward-facing (IF) conformation, based, respectively, on the crystal structure of AdiC and of GadC. As frequently classic Molecular Dynamics (MD) fails to monitor large-scale conformational transitions within a reasonable simulated time, the OF structure was equilibrated for 150 ns then processed through targeted MD (tMD). During this procedure, an elastic force pulled the OF structure to the IF structure and induced, at the same time, substrates/inhibitor to move through the transport channel. This elastic force was modulated by a spring constant (k) value; by decreasing its value from 100 to 70, it was possible to comparatively account for the propensity for transport of the three tested molecules. In line with our expectations, during the tMD simulations, Tyr and LME behaved as substrates, moving down the transport channel, or most of it, for all k values. On the contrary, DIT behaved as an inhibitor, being (almost) transported across the channel only at the highest k value (100). During their transit through the channel, Tyr and LME interacted with specific amino acids (first with Phe252 then with Thr345, Arg348, Tyr259, and Phe262); this suggests that a primary as well as a putative secondary gate may contribute to the transport of substrates. Quite on the opposite, DIT appeared to establish only transient interactions with side chains lining the external part of the transport channel. Our tMD simulations could thus efficiently discriminate between two transported substrates and one inhibitor, and therefore can be proposed as a benchmark for developing novel LAT1 inhibitors of pharmacological interest.

Advanced molecular dynamics simulation methods for kinase drug discovery

Interest in the application of molecular dynamics (MD) simulations has increased in the field of protein kinase (PK) drug discovery. PKs belong to an important drug target class because they are directly involved in a number of diseases, including cancer. MD methods simulate dynamic biological and chemical events at an atomic level. This information can be combined with other in silico and experimental methods to efficiently target selected receptors. In this review, we present common and advanced methods of MD simulations and we focus on the recent applications of MD-based methodologies that provided significant insights into the elucidation of biological mechanisms involving PKs and into the discovery of novel kinase inhibitors.

Computing ensembles of transitions with molecular dynamics simulations

A molecular understanding of conformational change is important for connecting structure and function. Without the ability to sample on the meaningful large-scale conformational changes, the ability to infer biological function and to understand the effect of mutations and changes in environment is not possible. Our Dynamic Importance Sampling method (DIMS), part of the CHARMM simulation package, is a method that enables sampling over ensembles of transition intermediates. This chapter outlines the context for the method and the usage within the program.

Molecular insight into conformational transmission of human P-glycoprotein

P-glycoprotein (P-gp), a kind of ATP-binding cassette transporter, can export candidates through a channel at the two transmembrane domains (TMDs) across the cell membranes using the energy released from ATP hydrolysis at the two nucleotide-binding domains (NBDs). Considerable evidence has indicated that human P-gp undergoes large-scale conformational changes to export a wide variety of anti-cancer drugs out of the cancer cells. However, molecular mechanism of the conformational transmission of human P-gp from the NBDs to the TMDs is still unclear. Herein, targeted molecular dynamics simulations were performed to explore the atomic detail of the conformational transmission of human P-gp. It is confirmed that the conformational transition from the inward- to outward-facing is initiated by the movement of the NBDs. It is found that the two NBDs move both on the two directions (x and y). The movement on the x direction leads to the closure of the NBDs, while the movement on the y direction adjusts the conformations of the NBDs to form the correct ATP binding pockets. Six key segments (KSs) protruding from the TMDs to interact with the NBDs are identified. The relative movement of the KSs along the y axis driven by the NBDs can be transmitted through α-helices to the rest of the TMDs, rendering the TMDs to open towards periplasm in the outward-facing conformation. Twenty eight key residue pairs are identified to participate in the interaction network that contributes to the conformational transmission from the NBDs to the TMDs of human P-gp. In addition, 9 key residues in each NBD are also identified. The studies have thus provided clear insight into the conformational transmission from the NBDs to the TMDs in human P-gp.

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.

The study of interactions between DNA and PcrA DNA helicase by using targeted molecular dynamic simulations

DNA helicases are important enzymes involved in all aspects of nucleic acid metabolism, ranging from DNA replication and repair to recombination, rescue of stalled replication and translation. DNA helicases are molecular motors. Through conformational changes caused by ATP hydrolysis and binding, they move along the template double helix, break the hydrogen bonds between the two strands and separate the template chains, so that the genetic information can be accessed. In this paper, targeted molecular dynamic simulations were performed to study the important interactions between DNA and PcrA DNA helicase, which can not be observed from the crystal structures. The key residues on PcrA DNA helicase that have strong interactions with both double stranded DNA (ds-DNA) and single stranded DNA (ss-DNA) have been identified, and it was found that such interactions mostly exist between the protein and DNA backbone, which indicates that the translocation of PcrA is independent of the DNA sequence. The simulations indicate that the ds-DNA is separated upon ATP rebinding, rather than ATP hydrolysis, which suggests that the two strokes in the mechanism have two different major roles. Firstly, in the power stroke (ATP hydrolysis), most of the translocations of the bases from one pocket to the next occur. In the relaxation stroke (ATP binding), most of the 'work' is being done to 'melt' the DNA at the separation fork. Therefore, we propose a mechanism whereby the translocation of the ss-DNA is powered by ATP hydrolysis and the separation of the ds-DNA is powered by ATP binding.

Frequency response of a protein to local conformational perturbations

Signals created by local perturbations are known to propagate long distances through proteins via backbone connectivity and nonbonded interactions. In the current study, signal propagation from the flexible ligand binding loop to the rest of Protein Tyrosine Phosphatase 1B (PTP1B) was investigated using frequency response techniques. Using restrained Targeted Molecular Dynamics (TMD) potential on WPD and R loops, PTP1B was driven between its crystal structure conformations at different frequencies. Propagation of the local perturbation signal was manifested via peaks at the fundamental frequency and upper harmonics of 1/f distributed spectral density of atomic variables, such as Cα atoms, dihedral angles, or polar interaction distances. Frequency of perturbation was adjusted high enough (simulation length >∼10×period of a perturbation cycle) not to be clouded by random diffusional fluctuations, and low enough (<∼0.8 ns(-1)) not to attenuate the propagating signal and enhance the contribution of the side-chains to the dissipation of the signals. Employing Discrete Fourier Transform (DFT) to TMD simulation trajectories of 16 cycles of conformational transitions at periods of 1.2 to 5 ns yielded Cα displacements consistent with those obtained from crystal structures. Identification of the perturbed atomic variables by statistical t-tests on log-log scale spectral densities revealed the extent of signal propagation in PTP1B, while phase angles of the filtered trajectories at the fundamental frequency were used to cluster collectively fluctuating elements. Hydrophobic interactions were found to have a higher contribution to signal transduction between side-chains compared to the role of polar interactions. Most of in-phase fluctuating residues on the signaling pathway were found to have high identity among PTP domains, and located over a wide region of PTP1B including the allosteric site. Due to its simplicity and efficiency, the suggested technique may find wide applications in identification of signaling pathways of different proteins.

Allosteric activation transitions in enzymes and biomolecular motors: insights from atomistic and coarse-grained simulations

The chemical step in enzymes is usually preceded by a kinetically distinct activation step that involves large-scale conformational transitions. In "simple" enzymes this step corresponds to the closure of the active site; in more complex enzymes, such as biomolecular motors, the activation step is more complex and may involve interactions with other biomolecules. These activation transitions are essential to the function of enzymes and perturbations in the scale and/or rate of these transitions are implicated in various serious human diseases; incorporating key flexibilities into engineered enzymes is also considered a major remaining challenge in rational enzyme design. Therefore it is important to understand the underlying mechanism of these transitions. This is a significant challenge to both experimental and computational studies because of the allosteric and multi-scale nature of such transitions. Using our recent studies of two enzyme systems, myosin and adenylate kinase (AK), we discuss how atomistic and coarse-grained simulations can be used to provide insights into the mechanism of activation transitions in realistic systems. Collectively, the results suggest that although many allosteric transitions can be viewed as domain displacements mediated by flexible hinges, there are additional complexities and various deviations. For example, although our studies do not find any evidence for "cracking" in AK, our results do underline the contribution of intra-domain properties (e.g., dihedral flexibility) to the rate of the transition. The study of mechanochemical coupling in myosin highlights that local changes important to chemistry require stabilization from more extensive structural changes; in this sense, more global structural transitions are needed to activate the chemistry in the active site. These discussions further emphasize the importance of better understanding factors that control the degree of co-operativity for allosteric transitions, again hinting at the intimate connection between protein stability and functional flexibility. Finally, a number of topics of considerable future interest are briefly discussed.

Analysis and elimination of a bias in targeted molecular dynamics simulations of conformational transitions: application to calmodulin

The popular targeted molecular dynamics (TMD) method for generating transition paths in complex biomolecular systems is revisited. In a typical TMD transition path, the large-scale changes occur early and the small-scale changes tend to occur later. As a result, the order of events in the computed paths depends on the direction in which the simulations are performed. To identify the origin of this bias, and to propose a method in which the bias is absent, variants of TMD in the restraint formulation are introduced and applied to the complex open ↔ closed transition in the protein calmodulin. Due to the global best-fit rotation that is typically part of the TMD method, the simulated system is guided implicitly along the lowest-frequency normal modes, until the large spatial scales associated with these modes are near the target conformation. The remaining portion of the transition is described progressively by higher-frequency modes, which correspond to smaller-scale rearrangements. A straightforward modification of TMD that avoids the global best-fit rotation is the locally restrained TMD (LRTMD) method, in which the biasing potential is constructed from a number of TMD potentials, each acting on a small connected portion of the protein sequence. With a uniform distribution of these elements, transition paths that lack the length-scale bias are obtained. Trajectories generated by steered MD in dihedral angle space (DSMD), a method that avoids best-fit rotations altogether, also lack the length-scale bias. To examine the importance of the paths generated by TMD, LRTMD, and DSMD in the actual transition, we use the finite-temperature string method to compute the free energy profile associated with a transition tube around a path generated by each algorithm. The free energy barriers associated with the paths are comparable, suggesting that transitions can occur along each route with similar probabilities. This result indicates that a broad ensemble of paths needs to be calculated to obtain a full description of conformational changes in biomolecules. The breadth of the contributing ensemble suggests that energetic barriers for conformational transitions in proteins are offset by entropic contributions that arise from a large number of possible paths.

Towards the prediction of order parameters from molecular dynamics simulations in proteins

A molecular understanding of how protein function is related to protein structure requires an ability to understand large conformational changes between multiple states. Unfortunately these states are often separated by high free energy barriers and within a complex energy landscape. This makes it very difficult to reliably connect, for example by all-atom molecular dynamics calculations, the states, their energies, and the pathways between them. A major issue needed to improve sampling on the intermediate states is an order parameter--a reduced descriptor for the major subset of degrees of freedom--that can be used to aid sampling for the large conformational change. We present a method to combine information from molecular dynamics using non-linear time series and dimensionality reduction, in order to quantitatively determine an order parameter connecting two large-scale conformationally distinct protein states. This new method suggests an implementation for molecular dynamics calculations that may be used to enhance sampling of intermediate states.

Strandwise translocation of a DNA glycosylase on undamaged DNA

Base excision repair of genotoxic nucleobase lesions in the genome is critically dependent upon the ability of DNA glycosylases to locate rare sites of damage embedded in a vast excess of undamaged DNA, using only thermal energy to fuel the search process. Considerable interest surrounds the question of how DNA glycosylases translocate efficiently along DNA while maintaining their vigilance for target damaged sites. Here, we report the observation of strandwise translocation of 8-oxoguanine DNA glycosylase, MutM, along undamaged DNA. In these complexes, the protein is observed to translocate by one nucleotide on one strand while remaining untranslocated on the complementary strand. We further report that alterations of single base-pairs or a single amino acid substitution (R112A) can induce strandwise translocation. Molecular dynamics simulations confirm that MutM can translocate along DNA in a strandwise fashion. These observations reveal a previously unobserved mode of movement for a DNA-binding protein along the surface of DNA.

Hysteresis-based mechanism for the directed motility of the Ncd motor

Ncd is a Kinesin-14 family protein that walks to the microtubule's minus end. Although available structures show its α-helical neck in either pre- or post-stroke orientations, little is known about the transition between these two states. Using a combination of molecular dynamics simulations and structural analyses, we find that the neck sequentially makes intermediate contacts with the motor head along its mostly longitudinal path, and it develops a 24° twist in the post-stroke orientation. The forward (pre-stroke to post-stroke) motion has an ∼4.5 k(B)T (where k(B) is the Boltzmann constant, and T=300 K) free-energy barrier and is a diffusion guided by the intermediate contacts. The post-stroke free-energy minimum is higher and is formed ∼10° before reaching the orientation in the post-stroke crystal structure, consistent with previous structural data. The importance of intermediate contacts correlates with the existing motility data, including those for mutant Ncds. Unlike the forward motion, the recovery stroke goes nearly downhill in free energy, powered in part by torsional relaxation of the neck. The hysteresis in the energetics of the neck motion arises from the mechanical compliance of the protein, and together with guided diffusion, it may be key to the directed motility of Ncd.

Importance of domain closure for the autoactivation of ERK2

Extracellular signal-regulated kinases 1 and 2 (ERK1 and -2, respectively) play a critical role in regulating cell division and have been implicated in cancer. In addition to activation by MAPK/ERK kinases 1 and 2 (MEK1 and -2, respectively), certain mutants of ERK2 can be activated by autophosphorylation. To identify the mechanism of autoactivation, we have performed a series of molecular dynamics simulations of ERK1 and -2 in various stages of activation as well as the constitutively active Q103A, I84A, L73P, and R65S ERK2 mutants. Our simulations indicate the importance of domain closure for autoactivation and activity regulation, with that event occurring prior to folding of the activation lip and of loop L16. Results indicate that the second phosphorylation event, that of T183, disrupts hydrogen bonding involving D334, thereby allowing the kinase to lock into the active conformation. On the basis of the simulations, three predictions were made. G83A was suggested to impede activation; K162M was suggested to perturb the interface between the N- and C-domains leading to activation, and Q64C was hypothesized to stop folding of loop L16, thereby perturbing the homodimerization interface. Functional analysis of the mutants validated the predictions concerning the G83A and Q64C mutants. The K162M mutant did not autoactivate as predicted, however, which may be due to the location of the residue on the protein surface near the ED substrate docking domain.

A normal mode-based geometric simulation approach for exploring biologically relevant conformational transitions in proteins

A three-step approach for multiscale modeling of protein conformational changes is presented that incorporates information about preferred directions of protein motions into a geometric simulation algorithm. The first two steps are based on a rigid cluster normal-mode analysis (RCNMA). Low-frequency normal modes are used in the third step (NMSim) to extend the recently introduced idea of constrained geometric simulations of diffusive motions in proteins by biasing backbone motions of the protein, whereas side-chain motions are biased toward favorable rotamer states. The generated structures are iteratively corrected regarding steric clashes and stereochemical constraint violations. The approach allows performing three simulation types: unbiased exploration of conformational space; pathway generation by a targeted simulation; and radius of gyration-guided simulation. When applied to a data set of proteins with experimentally observed conformational changes, conformational variabilities are reproduced very well for 4 out of 5 proteins that show domain motions, with correlation coefficients r > 0.70 and as high as r = 0.92 in the case of adenylate kinase. In 7 out of 8 cases, NMSim simulations starting from unbound structures are able to sample conformations that are similar (root-mean-square deviation = 1.0-3.1 Å) to ligand bound conformations. An NMSim generated pathway of conformational change of adenylate kinase correctly describes the sequence of domain closing. The NMSim approach is a computationally efficient alternative to molecular dynamics simulations for conformational sampling of proteins. The generated conformations and pathways of conformational transitions can serve as input to docking approaches or as starting points for more sophisticated sampling techniques.

Cy3-DNA stacking interactions strongly depend on the identity of the terminal basepair

We characterized the effect of the first basepair on the conformational dynamics of the fluorescent dye Cy3 attached to the 5' end of double-stranded DNA using gaussian-mixture adaptive umbrella sampling simulations. In the simulations, the sampling of all five dihedral angles along the linker was enhanced, so that both stacked and unstacked states were sampled. The affinity of Cy3 for a T·A basepair (with the dye attached to T) was found to be significantly less than for the other basepairs. This was verified experimentally by measuring the activation energies for cis-trans isomerization of the dye. The simulation and experimental results indicate the existence of partially unstacked conformations amenable to photoisomerization. The simulations also showed that stacking of Cy3 straightens the DNA while stabilizing the first basepair. Our findings indicate that fluorescence is modulated by Cy3-DNA interactions in a sequence-dependent manner.

Computing ensembles of transitions from stable states: Dynamic importance sampling

There is an increasing dataset of solved biomolecular structures in more than one conformation and increasing evidence that large-scale conformational change is critical for biomolecular function. In this article, we present our implementation of a dynamic importance sampling (DIMS) algorithm that is directed toward improving our understanding of important intermediate states between experimentally defined starting and ending points. This complements traditional molecular dynamics methods where most of the sampling time is spent in the stable free energy wells defined by these initial and final points. As such, the algorithm creates a candidate set of transitions that provide insights for the much slower and probably most important, functionally relevant degrees of freedom. The method is implemented in the program CHARMM and is tested on six systems of growing size and complexity. These systems, the folding of Protein A and of Protein G, the conformational changes in the calcium sensor S100A6, the glucose-galactose-binding protein, maltodextrin, and lactoferrin, are also compared against other approaches that have been suggested in the literature. The results suggest good sampling on a diverse set of intermediates for all six systems with an ability to control the bias and thus to sample distributions of trajectories for the analysis of intermediate states.

Putative active states of a prototypic g-protein-coupled receptor from biased molecular dynamics

A major current focus of structural work on G-protein-coupled receptors (GPCRs) pertains to the investigation of their active states. However, for virtually all GPCRs, active agonist-bound intermediate states have been difficult to characterize experimentally owing to their higher conformational flexibility, and thus intrinsic instability, as compared to inactive inverse agonist-bound states. In this work, we explored possible activation pathways of the prototypic GPCR bovine rhodopsin by means of biased molecular dynamics simulations. Specifically, we used an explicit atomistic representation of the receptor and its environment, and sampled the conformational transition from the crystal structure of a photoactivated deprotonated state of rhodopsin to the low pH crystal structure of opsin in the presence of 11-trans-retinal, using adiabatic biased molecular dynamics simulations. We then reconstructed the system free-energy landscape along the predetermined transition trajectories using a path collective variable approach based on metadynamics. Our results suggest that the two experimental endpoints of rhodopsin/opsin are connected by at least two different pathways, and that the conformational transition is populated by at least four metastable states of the receptor, characterized by a different amplitude of the outward movement of transmembrane helix 6.

Equilibrium sampling for biomolecules under mechanical tension

In the studies of force-induced conformational transitions of biomolecules, the large timescale difference from experiments presents the challenge of obtaining convergent sampling for molecular dynamics simulations. To circumvent this fundamental problem, an approach combining the replica-exchange method and umbrella sampling (REM-US) was developed to simulate mechanical stretching of biomolecules under equilibrium conditions. Equilibrium properties of conformational transitions can be obtained directly from simulations without further assumptions. To test the performance, we carried out REM-US simulations of atomic force microscope (AFM) stretching and relaxing measurements on the polysaccharide pustulan, a (1-->6)-beta-D-glucan, which undergoes well-characterized rotameric transitions in the backbone bonds. With significantly enhanced sampling convergence and efficiency, the REM-US approach closely reproduced the equilibrium force-extension curves measured in AFM experiments. Consistent with the reversibility in the AFM measurements, the new approach generated identical force-extension curves in both stretching and relaxing simulations-an outcome not reported in previous studies, proving that equilibrium conditions were achieved in the simulations. REM-US may provide a robust approach to modeling of mechanical stretching on polysaccharides and even nucleic acids.

Mechanical coupling in myosin V: a simulation study

Myosin motor function depends on the interaction between different domains that transmit information from one part of the molecule to another. The interdomain coupling in myosin V is studied with restrained targeted molecular dynamics using an all-atom representation in explicit solvent. To elucidate the origin of the conformational change due to the binding of ATP, targeting forces are applied to small sets of atoms (the forcing sets, FSs) in the direction of their displacement from the rigor conformation, which has a closed actin-binding cleft, to the post-rigor conformation, in which the cleft is open. The "minimal" FS that results in extensive structural changes in the overall myosin conformation is composed of ATP, switch 1, and the nearby HF, HG, and HH helices. Addition of switch 2 to the FS is required to achieve a complete opening of the actin-binding cleft. The restrained targeted molecular dynamics simulations reveal the mechanical coupling pathways between (i) the nucleotide-binding pocket (NBP) and the actin-binding cleft, (ii) the NBP and the converter, and (iii) the actin-binding cleft and the converter. Closing of the NBP due to ATP binding is tightly coupled to the opening of the cleft and leads to the rupture of a key hydrogen bond (F441N/A684O) between switch 2 and the SH1 helix. The actin-binding cleft may mediate the rupture of this bond via a connection between the HW helix, the relay helix, and switch 2. The findings are consistent with experimental studies and a recent normal mode analysis. The present method is expected to be useful more generally in studies of interdomain coupling in proteins.

Encounter and extrusion of an intrahelical lesion by a DNA repair enzyme

How living systems detect the presence of genotoxic damage embedded in a million-fold excess of undamaged DNA is an unresolved question in biology. Here we have captured and structurally elucidated a base-excision DNA repair enzyme, MutM, at the stage of initial encounter with a damaged nucleobase, 8-oxoguanine (oxoG), nested within a DNA duplex. Three structures of intrahelical oxoG-encounter complexes are compared with sequence-matched structures containing a normal G base in place of an oxoG lesion. Although the protein-DNA interfaces in the matched complexes differ by only two atoms-those that distinguish oxoG from G-their pronounced structural differences indicate that MutM can detect a lesion in DNA even at the earliest stages of encounter. All-atom computer simulations show the pathway by which encounter of the enzyme with the lesion causes extrusion from the DNA duplex, and they elucidate the critical free energy difference between oxoG and G along the extrusion pathway.

The mechanism of the translocation step in DNA replication by DNA polymerase I: a computer simulation analysis

High-fidelity DNA polymerases copy DNA rapidly and accurately by adding correct deoxynucleotide triphosphates to a growing primer strand of DNA. Following nucleotide incorporation, a series of conformational changes translocate the DNA substrate by one base pair step, readying the polymerase for the next round of incorporation. Molecular dynamics simulations indicate that the translocation consists globally of a polymerase fingers-opening transition, followed by the DNA displacement and the insertion of the template base into the preinsertion site. They also show that the pyrophosphate release facilitates the opening transition and that the universally conserved Y714 plays a key role in coupling polymerase opening to DNA translocation. The transition involves several metastable intermediates in one of which the O helix is bent in the vicinity of G711. Completion of the translocation appears to require a gating motion of the O1 helix, perhaps facilitated by the presence of G715. These roles are consistent with the high level of conservation of Y714 and the two glycine residues at these positions. It is likely that a corresponding mechanism is applicable to other polymerases.

Transition path sampling study of the conformational fluctuation of His-64 in human carbonic anhydrase II

We report here a transition path sampling study of the conformational fluctuation of His-64 that is known to be important in the enzymatic catalysis of human carbonic anhydrase II. The dynamical transition between experimentally detected conformations of His-64 could not be observed using classical molecular dynamics trajectories extended to 3.5 ns, indicating the transition to be rare on the time scale of molecular dynamics. Using the transition path sampling method, an ensemble of transition paths between these two conformers has been generated and analyzed in detail to identify the mechanism of coupling of His-64 to its neighboring residues during the conformational transition. It is found that both Asn-62 and Tyr-7 may contribute toward retaining the His-64 residue in its outward conformation. Trp-5, on the other hand, shows marked motions at the transition state. The number of water molecules inside a part of the active site cavity and the corresponding cavity volume are also found to vary coupled to the His-64 conformational dynamics.

Targeted molecular dynamics reveals overall common conformational changes upon hybrid domain swing-out in beta3 integrins

The beta3 integrin family members alphaIIbeta3 and alphaVbeta3 signal bidirectionally through long-range allosteric changes, including a transition from a bent unliganded-closed low-affinity state to an extended liganded-open high-affinity state. To obtain an atomic-level description of this transition in an explicit solvent, we carried out targeted molecular dynamics simulations of the headpieces of alphaIIbeta3 and alphaVbeta3 integrins. Although minor differences were observed between these receptors, our results suggest a common transition pathway in which the hybrid domain swing-out is accompanied by conformational changes within the beta3 betaA (I-like) domain that propagate through the alpha7 helix C-terminus, and are followed by the alpha7 helix downward motion and the opening of the beta6-alpha7 loop. Breaking of contact interactions between the beta6-alpha7 loop and the alpha1 helix N-terminus results in helix straightening, internal rearrangements of the specificity determining loop (SDL), movement of the beta1-alpha1 loop toward the metal ion dependent adhesion site (MIDAS), and final changes at the interfaces between the beta3 betaA (I-like) domain and either the hybrid or the alpha beta-propeller domains. Taken together, our results suggest novel testable hypotheses of intradomain and interdomain interactions responsible for beta3 integrin activation.

A computational analysis of ATP binding of SV40 large tumor antigen helicase motor

Simian Virus 40 Large Tumor Antigen (LTag) is an efficient helicase motor that unwinds and translocates DNA. The DNA unwinding and translocation of LTag is powered by ATP binding and hydrolysis at the nucleotide pocket between two adjacent subunits of an LTag hexamer. Based on the set of high-resolution hexameric structures of LTag helicase in different nucleotide binding states, we simulated a conformational transition pathway of the ATP binding process using the targeted molecular dynamics method and calculated the corresponding energy profile using the linear response approximation (LRA) version of the semi-macroscopic Protein Dipoles Langevin Dipoles method (PDLD/S). The simulation results suggest a three-step process for the ATP binding from the initial interaction to the final tight binding at the nucleotide pocket, in which ATP is eventually “locked” by three pairs of charge-charge interactions across the pocket. Such a “cross-locking” ATP binding process is similar to the binding zipper model reported for the F1-ATPase hexameric motor. The simulation also shows a transition mechanism of Mg²⁺ coordination to form the Mg-ATP complex during ATP binding, which is accompanied by the large conformational changes of LTag. This simulation study of the ATP binding process to an LTag and the accompanying conformational changes in the context of a hexamer leads to a refined cooperative iris model that has been proposed previously.

The Large Tumor antigen (LTag) encoded by Simian Virus 40 (SV40) is a marvelous molecule that is not only a viral oncogene, but also an efficient molecular machine as a helicase that unwinds double helix DNA for genome replication, an essential process in all living organisms. LTag hexameric helicase uses the energy of ATP to power its conformational switch for DNA unwinding. Understanding how the LTag conformational switch is coupled to the energy from ATP usage by LTag to do the mechanical work of unwinding DNA is of great interest to biologists, and yet remains to be established. Based on our previous high-resolution structures of LTag helicase in different conformational states, we simulated an LTag conformational transition pathway in the ATP binding process using the targeted molecular dynamics method. Our simulation results suggest a three-step process for the ATP binding to the nucleotide pocket, in which ATP is eventually “locked” into the pocket by three pairs of “locker” interactions. We have also quantitatively evaluated the energy profile of ATP binding using a special computational simulation technique. Additionally, our simulation study of ATP binding by LTag and the accompanying conformational switches in the context of a hexamer leads to a refined cooperative iris model that may be used for DNA unwinding.

CHARMM: the biomolecular simulation program

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983.

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

Conformational free-energy differences are key quantities for understanding important phenomena in molecular biology that involve large structural changes of macromolecules. In this paper, an improved version of the confinement approach, which is based on earlier developments, is used to determine the free energy of individual molecular states by progressively restraining the corresponding molecular structures to pure harmonic basins, whose absolute free energy can be computed by normal-mode analysis. The method is used to calculate the free-energy difference between two conformational states of the alanine dipeptide in vacuo, and of the beta-hairpin from protein G with an implicit solvation model. In all cases, the confinement results are in excellent agreement with the ones obtained from converged equilibrium molecular dynamics simulations, which have a much larger computational cost. The systematic and statistical errors of the results are evaluated and the origin of the errors is identified. The sensitivity of the calculated free-energy differences to structure-based definitions of the molecular states is discussed. A variant of the method, which closes the thermodynamic cycle by a quasi-harmonic rather than harmonic analysis, is introduced. The latter is proposed for possible use with explicit solvent simulations.

An optimized replica exchange molecular dynamics method

We introduce a new way to perform swaps between replicas in replica exchange molecular dynamics simulations. The method is based on a generalized canonical probability distribution function and flattens the potential of mean force along the temperature coordinate, such that a random walk in temperature space is achieved. Application to a Go model of protein A showed that the method is more efficient than conventional replica exchange. The method results in a constant probability distribution of the replicas over the thermostats, yields a minimum round-trip time between extremum temperatures, and leads to faster ergodic convergence.

Computational modeling of structurally conserved cancer mutations in the RET and MET kinases: the impact on protein structure, dynamics, and stability

Structural and biochemical characterization of protein kinases that confer oncogene addiction and harbor a large number of disease-associated mutations, including RET and MET kinases, have provided insights into molecular mechanisms associated with the protein kinase activation in human cancer. In this article, structural modeling, molecular dynamics, and free energy simulations of a structurally conserved mutational hotspot, shared by M918T in RET and M1250T in MET kinases, are undertaken to quantify the molecular mechanism of activation and the functional role of cancer mutations in altering protein kinase structure, dynamics, and stability. The mechanistic basis of the activating RET and MET cancer mutations may be driven by an appreciable free energy destabilization of the inactive kinase state in the mutational forms. According to our results, the locally enhanced mobility of the cancer mutants and a higher conformational entropy are counterbalanced by a larger enthalpy loss and result in the decreased thermodynamic stability. The computed protein stability differences between the wild-type and cancer kinase mutants are consistent with circular dichroism spectroscopy and differential scanning calorimetry experiments. These results support the molecular mechanism of activation, which causes a detrimental imbalance in the dynamic equilibrium shifted toward the active form of the enzyme. Furthermore, computer simulations of the inhibitor binding with the oncogenic and drug-resistant RET mutations have also provided a plausible molecular rationale for the observed differences in the inhibition profiles, which is consistent with the experimental data. Finally, structural mapping of RET and MET cancer mutations and the computed protein stability changes suggest a similar mechanism of activation, whereby the cancer mutations which display the higher oncogenic activity tend to have the greatest destabilization effect on the inactive kinase structure.

Mapping the conformational transition in Src activation by cumulating the information from multiple molecular dynamics trajectories

The Src-family kinases are allosteric enzymes that play a key role in the regulation of cell growth and proliferation. In response to cellular signals, they undergo large conformational changes to switch between distinct inactive and active states. A computational strategy for characterizing the conformational transition pathway is presented to bridge the inactive and active states of the catalytic domain of Hck. The information from a large number (78) of independent all-atom molecular dynamics trajectories with explicit solvent is combined together to assemble a connectivity map of the conformational transition. Two intermediate states along the activation pathways are identified, and their structural features are characterized. A coarse free-energy landscape is built in terms of the collective motions corresponding to the opening of the activation loop (A-loop) and the rotation of the alphaC helix. This landscape shows that the protein can adopt a multitude of conformations in which the A-loop is partially open, while the alphaC helix remains in the orientation characteristic of the inactive conformation. The complete transition leading to the active conformation requires a concerted movement involving further opening of the A-loop, the relative alignment of N-lobe and C-lobe, and the rotation of the alphaC helix needed to recruit the residues necessary for catalysis in the active site. The analysis leads to a dynamic view of the full-length kinase activation, whereby transitions of the catalytic domain to intermediate configurations with a partially open A-loop are permitted, even while the SH2-SH3 clamp remains fully engaged. These transitions would render Y416 available for the transphosphorylation event that ultimately locks down the active state. The results provide a broad framework for picturing the conformational transitions leading to kinase activation.

Diffusive reaction dynamics on invariant free energy profiles

A fundamental problem in the analysis of protein folding and other complex reactions in which the entropy plays an important role is the determination of the activation free energy from experimental measurements or computer simulations. This article shows how to combine minimum-cut-based free energy profiles (F(C)), obtained from equilibrium molecular dynamics simulations, with conventional histogram-based free energy profiles (F(H)) to extract the coordinate-dependent diffusion coefficient on the F(C) (i.e., the method determines free energies and a diffusive preexponential factor along an appropriate reaction coordinate). The F(C), in contrast to the F(H), is shown to be invariant with respect to arbitrary transformations of the reaction coordinate, which makes possible partition of configuration space into basins in an invariant way. A "natural coordinate," for which F(H) and F(C) differ by a multiplicative constant (constant diffusion coefficient), is introduced. The approach is illustrated by a model one-dimensional system, the alanine dipeptide, and the folding reaction of a double beta-hairpin miniprotein. It is shown how the results can be used to test whether the putative reaction coordinate is a good reaction coordinate.

Conformational transition pathways explored by Monte Carlo simulation integrated with collective modes

Conformational transitions between open/closed or free/bound states in proteins possess functional importance. We propose a technique in which the collective modes obtained from an anisotropic network model (ANM) are used in conjunction with a Monte Carlo (MC) simulation approach, to investigate conformational transition pathways and pathway intermediates. The ANM-MC technique is applied to adenylate kinase (AK) and hemoglobin. The iterative method, in which normal modes are continuously updated during the simulation, proves successful in accomplishing the transition between open-closed conformations of AK and tense-relaxed forms of hemoglobin (C(alpha)-root mean square deviations between two end structures of 7.13 A and 3.55 A, respectively). Target conformations are reached by root mean-square deviations of 2.27 A and 1.90 A for AK and hemoglobin, respectively. The intermediate conformations overlap with crystal structures from the AK family within a 3.0-A root mean-square deviation. In the case of hemoglobin, the transition of tense-to-relaxed passes through the relaxed state. In both cases, the lowest-frequency modes are effective during transitions. The targeted Monte Carlo approach is used without the application of collective modes. Both the ANM-MC and targeted Monte Carlo techniques can explore sequences of events in transition pathways with an efficient yet realistic conformational search.

pFlexAna: detecting conformational changes in remotely related proteins

The pFlexAna (protein flexibility analyzer) web server detects and displays conformational changes in remotely related proteins, without relying on sequence homology. To do so, it first applies a reliable statistical test to align core protein fragments that are structurally similar and then clusters these aligned fragment pairs into 'super-alignments', according to the similarity of geometric transformations that align them. The result is that the dominant conformational changes occur between the clusters, while the smaller conformational changes occur within a cluster. pFlexAna is available at http://bigbird.comp.nus.edu.sg/pfa2/.

Geometry-based sampling of conformational transitions in proteins

The fast and accurate prediction of protein flexibility is one of the major challenges in protein science. Enzyme activity, signal transduction, and ligand binding are dynamic processes involving essential conformational changes ranging from small side chain fluctuations to reorientations of entire domains. In the present work, we describe a reimplementation of the CONCOORD approach, termed tCONCOORD, which allows a computationally efficient sampling of conformational transitions of a protein based on geometrical considerations. Moreover, it allows for the extraction of the essential degrees of freedom, which, in general, are the biologically relevant ones. The method rests on a reliable estimate of the stability of interactions observed in a starting structure, in particular those interactions that change during a conformational transition. Applications to adenylate kinase, calmodulin, aldose reductase, T4-lysozyme, staphylococcal nuclease, and ubiquitin show that experimentally known conformational transitions are faithfully predicted.

How subunit coupling produces the gamma-subunit rotary motion in F1-ATPase

F(o)F(1)-ATP synthase manufactures the energy "currency," ATP, of living cells. The soluble F(1) portion, called F(1)-ATPase, can act as a rotary motor, with ATP binding, hydrolysis, and product release, inducing a torque on the gamma-subunit. A coarse-grained plastic network model is used to show at a residue level of detail how the conformational changes of the catalytic beta-subunits act on the gamma-subunit through repulsive van der Waals interactions to generate a torque that drives unidirectional rotation, as observed experimentally. The simulations suggest that the calculated 85 degrees substep rotation is driven primarily by ATP binding and that the subsequent 35 degrees substep rotation is produced by product release from one beta-subunit and a concomitant binding pocket expansion of another beta-subunit. The results of the simulation agree with single-molecule experiments [see, for example, Adachi K, et al. (2007) Cell 130:309-321] and support a tri-site rotary mechanism for F(1)-ATPase under physiological condition.